From: hackbard <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/cgi-bin/gitweb.cgi?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);