From da2d9866e05b1b7a408ecda2b1695e07c30b0533 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 3 Apr 2008 18:52:10 +0200
Subject: [PATCH] optimized force eval in k loop!

---
 potentials/albe.c | 44 ++++++++++++--------------------------------
 1 file changed, 12 insertions(+), 32 deletions(-)

diff --git a/potentials/albe.c b/potentials/albe.c
index bab0e27..677d6e3 100644
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -276,7 +276,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	}
 
 	/* force contribution */
-	scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a));
+	scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
 	v3_scale(&force,&(exchange->dist_ij),scale);
 	v3_add(&(ai->f),&(ai->f),&force);
 	v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij
@@ -304,7 +304,7 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 	exchange->pre_dzeta=0.5*f_a*f_c*db;
 
 	/* energy contribution */
-	energy=0.5*f_c*(f_r-b*f_a);
+	energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 	moldyn->energy+=energy;
 	ai->e+=energy;
 
@@ -329,7 +329,7 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 	double pre_dzeta;
 	double f_c_ik,df_c_ik;
 	double dijdik_inv,fcdg,dfcg;
-	t_3dvec dcosdri,dcosdrj,dcosdrk;
+	t_3dvec dcosdrj,dcosdrk;
 	t_3dvec force,tmp;
 
 	params=moldyn->pot_params;
@@ -375,7 +375,7 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 	dg=exchange->dg[kcount];
 	cos_theta=exchange->cos_theta[kcount];
 
-	/* cos_theta derivatives wrt i,j,k */
+	/* cos_theta derivatives wrt j,k */
 	dijdik_inv=1.0/(d_ij*d_ik);
 	v3_scale(&dcosdrj,&dist_ik,dijdik_inv);		// j
 	v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
@@ -383,37 +383,11 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 	v3_scale(&dcosdrk,&dist_ij,dijdik_inv);		// k
 	v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 	v3_add(&dcosdrk,&dcosdrk,&tmp);
-	v3_add(&dcosdri,&dcosdrj,&dcosdrk);		// i
-	v3_scale(&dcosdri,&dcosdri,-1.0);
 
 	/* f_c_ik * dg, df_c_ik * g */
 	fcdg=f_c_ik*dg;
 	dfcg=df_c_ik*g;
 
-	/* derivative wrt i */
-	v3_scale(&force,&dist_ik,dfcg);
-	v3_scale(&tmp,&dcosdri,fcdg);
-	v3_add(&force,&force,&tmp);
-	v3_scale(&force,&force,pre_dzeta);
-
-	/* force contribution */
-	v3_add(&(ai->f),&(ai->f),&force);
-	
-#ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
-	if(ai==&(moldyn->atom[DATOM])) {
-		printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-		printf("  adding %f %f %f\n",force.x,force.y,force.z);
-		printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
-		printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
-		printf("    d ij ik = %f %f\n",d_ij,d_ik);
-	}
-}
-#endif
-
-	/* virial */
-	//virial_calc(ai,&force,&dist_ij);
-
 	/* derivative wrt j */
 	v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 
@@ -432,8 +406,11 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 }
 #endif
 
-	/* virial */
+	/* force contribution to atom i */
 	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
+
+	/* virial */
 	virial_calc(ai,&force,&dist_ij);
 
 	/* derivative wrt k */
@@ -457,8 +434,11 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 }
 #endif
 
-	/* virial */
+	/* force contribution to atom i */
 	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
+
+	/* virial */
 	virial_calc(ai,&force,&dist_ik);
 	
 	/* increase k counter */
-- 
2.20.1