X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe.c;h=677d6e3b6315d48226a4b3da950378b1c0fa865a;hb=da2d9866e05b1b7a408ecda2b1695e07c30b0533;hp=28008b0541b1d20a68351a46669d77f314b00bf6;hpb=20409ee0c545235c9246edde2d3cda5de0ddabde;p=physik%2Fposic.git diff --git a/potentials/albe.c b/potentials/albe.c index 28008b0..677d6e3 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -98,7 +98,7 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* albe 3 body potential function (first k loop) */ int albe_mult_3bp_k1(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; @@ -191,14 +191,6 @@ int albe_mult_3bp_k1(t_moldyn *moldyn, exchange->zeta_ij+=f_c_ik*g; } -#ifdef DEBUG - if((ai==&(moldyn->atom[0]))| - (aj==&(moldyn->atom[864]))| - (ak==&(moldyn->atom[1003]))) { - printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos); - } -#endif - /* store even more data for second k loop */ exchange->g[kcount]=g; exchange->dg[kcount]=dg; @@ -284,33 +276,35 @@ 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 #ifdef DEBUG - if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) { +if(moldyn->time>DSTART&&moldyn->timeatom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) { printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag); - printf("adding %f %f %f\n",force.x,force.y,force.z); + printf(" adding %f %f %f\n",force.x,force.y,force.z); if(ai==&(moldyn->atom[0])) - printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); + printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); if(aj==&(moldyn->atom[0])) - printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); - printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r), + printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); + printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r), f_c,b,f_a,f_r); - printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); + printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); } +} #endif /* virial */ virial_calc(ai,&force,&(exchange->dist_ij)); - /* dzeta prefactor = - 0.5 f_c f_a db */ - exchange->pre_dzeta=-0.5*f_a*f_c*db; + /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ + 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; @@ -322,7 +316,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* albe 3 body potential function (second k loop) */ int albe_mult_3bp_k2(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; @@ -335,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; @@ -381,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); @@ -389,33 +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(ai==&(moldyn->atom[0])) { - 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); - } -#endif - - /* virial */ - //virial_calc(ai,&force,&dist_ij); - /* derivative wrt j */ v3_scale(&force,&dcosdrj,fcdg*pre_dzeta); @@ -423,15 +395,22 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, v3_add(&(aj->f),&(aj->f),&force); #ifdef DEBUG - if(aj==&(moldyn->atom[0])) { +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 j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + printf(" total j: %f %f %f\n",aj->f.x,aj->f.y,aj->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 */ + /* 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 */ @@ -444,15 +423,22 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, v3_add(&(ak->f),&(ak->f),&force); #ifdef DEBUG - if(ak==&(moldyn->atom[0])) { +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 k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + printf(" total k: %f %f %f\n",ak->f.x,ak->f.y,ak->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 */ + /* 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 */