From: hackbard Date: Fri, 8 Dec 2006 00:21:15 +0000 (+0000) Subject: dist vectors changed, lj and ho forces adjusted, tersoff needs redesign X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=48c565cf62ba65a94b757b1c4118fe8938d736f8;p=physik%2Fposic.git dist vectors changed, lj and ho forces adjusted, tersoff needs redesign anyways ... *sigh* --- diff --git a/moldyn.c b/moldyn.c index f09c93a..880e56c 100644 --- 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);