finished basic tersoff stuff, now start fixing ... !

diff --git a/moldyn.c b/moldyn.c
index 830b8f0..14ae320 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -925,6 +925,7 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
        exchange->d=&(params->d[num]);
        exchange->h=&(params->h[num]);
 
+       exchange->betan=pow(*(exchange->beta),*(exchange->n));
        exchange->c2=params->c[num]*params->c[num];
        exchange->d2=params->d[num]*params->d[num];
        exchange->c2d2=exchange->c2/exchange->d2;
@@ -1035,11 +1036,12 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        t_tersoff_mult_params *params;
        t_tersoff_exchange *exchange;
        t_3dvec dist_ij,dist_ik,dist_jk;
-       t_3dvec db_ij,temp,force;
+       t_3dvec temp,force;
        double R,S,s_r;
        double d_ij,d_ik,d_jk;
        double f_c,df_c,b_ij,f_a,df_a;
-       double B,mu;
+       double n,c,d,h,neta,betan,betan_1;
+       double theta,cos_theta,sin_theta;
        int num;
 
        params=moldyn->pot_params;
@@ -1058,9 +1060,6 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        d_ij=exchange->d_ij;
        d_ij2=exchange->d_ij2;
 
-       B=*(params->exchange.B);
-       mu=*(params->exchange.mu);
-
        f_a=params->exchange.f_a;
        df_a=params->exchange.df_a;
        
@@ -1111,17 +1110,68 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        if(bc) check_per_bound(moldyn,&dist_jk);
        d_jk=v3_norm(&dist_jk);
 
-
-       // GO ON HERE !!!
-
-       cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/(2*d_ij*d_ik);
+       beta=*(exchange->beta);
+       betan=exchange->betan;
+       n=*(exchange->n);
+       c=*(exchange->c);
+       d=*(exchange->d);
+       h=*(exchange->h);
+       c2=exchange->c2;
+       d2=exchange->d2;
+       c2d2=exchange->c2d2;
+
+       numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
+       denom=2*d_ij*d_ik;
+       cos_theta=numer/denom;
        sin_theta=sqrt(1.0-(cos_theta*cos_theta));
        theta=arccos(cos_theta);
+       d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom);
+       d_theta1=2*denom-numer*2*d_ik/d_ij;
+       d_theta2=2*denom-numer*2*d_ij/d_ik;
+       d_theta1*=d_theta;
+       d_theta2*=d_theta;
 
        h_cos=(h-cos_theta);
        h_cos2=h_cos*h_cos;
        d2_h_cos2=d2-h_cos2;
 
+       /* some usefull expressions */
+       frac1=c2/(d2-h_cos2);
+       bracket1=1+c2d2-frac1;
+       bracket2=f_c_ik*bracket1;
+       bracket2_n_1=pow(bracket2,n-1.0);
+       bracket2_n=bracket2_n_1*bracket2;
+       bracket3=1+betan*bracket2_n;
+       bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
+       bracket3_pow=bracket3_pow_1*bracket3;
+
+       /* now go on with calc of b_ij and derivation of b_ij */
+       b_ij=chi*bracket3_pow;
+
+       /* derivation of theta */
+       v3_scale(&force,&dist_ij,d1_theta);
+       v3_scale(&temp,&dist_ik,d_theta2);
+       v3_add(&force,&force,&temp);
+
+       /* part 1 of derivation of b_ij */
+       v3_scale(&force,sin_theta*2*h_cos*f_c_ik*frac1);
+
+       /* part 2 of derivation of b_ij */
+       v3_scale(&temp,&dist_ik,df_c_ik*bracket1);
+
+       /* sum up and scale ... */
+       v3_add(&temp,&temp,&force);
+       scale=bracket2_n_1*n*betan*(1+betan*bracket3_pow_1)*chi*(1.0/(2.0*n));
+       v3_scale(&temp,&temp,scale);
+
+       /* now construct an energy and a force out of that */
+       v3_scale(&temp,&temp,f_a);
+       v3_scale(&force,&dist_ij,df_a*b_ij);
+       v3_add(&temp,&temp,&force);
+       v3_scale(&temp,&temp,f_c);
+       v3_scale(&force,&dist_ij,df_c*b_ij*f_a);
+       v3_add(&force,&force,&temp);
+
        /* add forces */
        v3_add(&(ai->f),&(ai->f),&force);
        /* energy is 0.5 f_r f_c, but we will sum it up twice ... */

diff --git a/moldyn.h b/moldyn.h
index 9c63307..0a6b4f2 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -145,6 +145,7 @@ typedef struct s_tersoff_exchange {
        double c2;
        double d2;
        double c2d2;
+       double betan;
 
        u8 run3bp;
 } t_tersoff_exchange;