fixed sign

[physik/posic.git] / potentials / tersoff.c

diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index 531d292..b6533e1 100644 (file)
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -153,7 +153,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* add forces */
        v3_add(&(ai->f),&(ai->f),&force);
-       v3_sub(&(aj->f),&(aj->f),&force);
+       v3_sub(&(aj->f),&(aj->f),&force); // reason: dri rij = - drj rij
 
 #ifdef DEBUG
        if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
@@ -388,7 +388,7 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                ni=*(exchange->n_i);
                tmp=exchange->betaini*pow(exchange->zeta_ij,ni-1.0);
                b=(1.0+exchange->zeta_ij*tmp);
-               db=chi*pow(b,-1.0/(2*ni)-1.0);
+               db=chi*pow(b,-1.0/(2.0*ni)-1.0);
                b=db*b;
                db*=-0.5*tmp;
        }
@@ -397,7 +397,7 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        v3_scale(&force,&(exchange->dist_ij),df_a*f_c+f_a*df_c);
        v3_scale(&force,&force,-0.5*b);
        v3_add(&(ai->f),&(ai->f),&force);
-       v3_sub(&(aj->f),&(aj->f),&force);
+       v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
 
 #ifdef DEBUG
        if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
@@ -522,7 +522,7 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
        /* virial */
        //virial_calc(ai,&force,&dist_ij);
 
-       /* derivatice wrt j */
+       /* derivative wrt j */
        v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 
        /* force contribution */
@@ -540,7 +540,7 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
        //virial_calc(aj,&force,&dist_ij);
 
        /* derivative wrt k */
-       v3_scale(&force,&dist_ik,dfcg);
+       v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rij = - drj rij
        v3_scale(&tmp,&dcosdrk,fcdg);
        v3_add(&force,&force,&tmp);
        v3_scale(&force,&force,pre_dzeta);