optimized force eval in k loop!

diff --git a/potentials/albe.c b/potentials/albe.c
index bab0e27..677d6e3 100644 (file)
--- 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 */