dist vectors changed, lj and ho forces adjusted, tersoff needs redesign

anyways ... *sigh*

diff --git a/moldyn.c b/moldyn.c
index f09c93a27e10fb00dc7a4dc6087bb1cb3f921948..880e56c1cdfa5809741bfe3a9d07792a907b3616 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -983,14 +983,15 @@ int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        sc=params->spring_constant;
        equi_dist=params->equilibrium_distance;
 
-       v3_sub(&distance,&(ai->r),&(aj->r));
+       v3_sub(&distance,&(aj->r),&(ai->r));
        
        if(bc) check_per_bound(moldyn,&distance);
        d=v3_norm(&distance);
        if(d<=moldyn->cutoff) {
                /* energy is 1/2 (d-d0)^2, but we will add this twice ... */
                moldyn->energy+=(0.25*sc*(d-equi_dist)*(d-equi_dist));
-               v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d)));
+               /* f = -grad E; grad r_ij = -1 1/r_ij distance */
+               v3_scale(&force,&distance,sc*(1.0-(equi_dist/d)));
                v3_add(&(ai->f),&(ai->f),&force);
        }
 
@@ -1011,7 +1012,7 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        sig6=params->sigma6;
        sig12=params->sigma12;
 
-       v3_sub(&distance,&(ai->r),&(aj->r));
+       v3_sub(&distance,&(aj->r),&(ai->r));
        if(bc) check_per_bound(moldyn,&distance);
        d=v3_absolute_square(&distance);        /* 1/r^2 */
        if(d<=moldyn->cutoff_square) {
@@ -1027,7 +1028,7 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                h1*=12*sig12;
                d=+h1-h2;
                d*=eps;
-               v3_scale(&force,&distance,d);
+               v3_scale(&force,&distance,-1.0*d); /* f = - grad E */
                v3_add(&(ai->f),&(ai->f),&force);
        }
 
@@ -1126,7 +1127,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
         * therefore we need: R, S, A, lambda
         */
 
-       v3_sub(&dist_ij,&(ai->r),&(aj->r));
+       v3_sub(&dist_ij,&(aj->r),&(ai->r));
 
        if(bc) check_per_bound(moldyn,&dist_ij);
 
@@ -1319,7 +1320,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
         *
         */
 
-       v3_sub(&dist_ik,&(ai->r),&(ak->r));
+       v3_sub(&dist_ik,&(ak->r),&(ai->r));
        if(bc) check_per_bound(moldyn,&dist_ik);
        d_ik=v3_norm(&dist_ik);