X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe.c;h=bab0e27ce156f74c8d1a8771980ee69a0d00e7bb;hb=f828bb665c033d828c520c68a98abdb6fa63a83a;hp=28008b0541b1d20a68351a46669d77f314b00bf6;hpb=20409ee0c545235c9246edde2d3cda5de0ddabde;p=physik%2Fposic.git diff --git a/potentials/albe.c b/potentials/albe.c index 28008b0..bab0e27 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)); 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); 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; @@ -406,11 +400,15 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, v3_add(&(ai->f),&(ai->f),&force); #ifdef DEBUG - if(ai==&(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 i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); + 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 */ @@ -423,11 +421,15 @@ 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 */ @@ -444,11 +446,15 @@ 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 */