From: hackbard Date: Thu, 3 Apr 2008 16:52:10 +0000 (+0200) Subject: optimized force eval in k loop! X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=da2d9866e05b1b7a408ecda2b1695e07c30b0533;p=physik%2Fposic.git optimized force eval in k loop! --- 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->timepre_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->timeatom[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->timef),&(ai->f),&force); + + /* virial */ virial_calc(ai,&force,&dist_ij); /* derivative wrt k */ @@ -457,8 +434,11 @@ if(moldyn->time>DSTART&&moldyn->timef),&(ai->f),&force); + + /* virial */ virial_calc(ai,&force,&dist_ik); /* increase k counter */