From 115ab83cedba54af2d165b8900781b98c5326b55 Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 28 Dec 2006 14:15:13 +0000 Subject: [PATCH] bugfix dg of dV_ji,jk --- moldyn.c | 70 +++++++++++++++++++++++++++++++++++++++++++++++--------- moldyn.h | 1 + 2 files changed, 60 insertions(+), 11 deletions(-) diff --git a/moldyn.c b/moldyn.c index faa6b7c..c21f59f 100644 --- a/moldyn.c +++ b/moldyn.c @@ -1129,6 +1129,9 @@ int potential_force_calc(t_moldyn *moldyn) { } #ifdef DEBUG printf("\n\n"); +#endif +#ifdef VDEBUG +printf("\n\n"); #endif return 0; @@ -1376,6 +1379,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* save for use in 3bp */ exchange->d_ij=d_ij; + exchange->d_ij2=d_ij2; exchange->dist_ij=dist_ij; /* more constants */ @@ -1443,6 +1447,14 @@ if(ai==&(moldyn->atom[0])) { printf("%f | %f\n",force.y,ai->f.y); printf("%f | %f\n",force.z,ai->f.z); } +#endif +#ifdef VDEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVij, dVji (2bp) contrib:\n"); + printf("%f | %f\n",force.x*dist_ij.x,ai->virial.xx); + printf("%f | %f\n",force.y*dist_ij.y,ai->virial.yy); + printf("%f | %f\n",force.z*dist_ij.z,ai->virial.zz); +} #endif /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */ @@ -1551,6 +1563,14 @@ if(ai==&(moldyn->atom[0])) { printf("%f | %f\n",force.y,ai->f.y); printf("%f | %f\n",force.z,ai->f.z); } +#endif +#ifdef VDEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVij (3bp) contrib:\n"); + printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx); + printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy); + printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz); +} #endif /* add energy of 3bp sum */ @@ -1596,6 +1616,14 @@ if(ai==&(moldyn->atom[0])) { printf("%f | %f\n",force.y,ai->f.y); printf("%f | %f\n",force.z,ai->f.z); } +#endif +#ifdef VDEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVji (3bp) contrib:\n"); + printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx); + printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy); + printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz); +} #endif return 0; @@ -1610,9 +1638,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_3dvec dist_ij,dist_ik,dist_jk; t_3dvec temp1,temp2; t_3dvec *dzeta; - double R,S,s_r; + double R,S,S2,s_r; double B,mu; - double d_ij,d_ik,d_jk; + double d_ij,d_ik,d_jk,d_ij2,d_ik2,d_jk2; double rr,dd; double f_c,df_c; double f_c_ik,df_c_ik,arg; @@ -1658,6 +1686,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* dist_ij, d_ij - this is < S_ij ! */ dist_ij=exchange->dist_ij; d_ij=exchange->d_ij; + d_ij2=exchange->d_ij2; /* f_c_ij, df_c_ij (same for ji) */ f_c=exchange->f_c; @@ -1672,21 +1701,26 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* dist_ik, d_ik */ v3_sub(&dist_ik,&(ak->r),&(ai->r)); if(bc) check_per_bound(moldyn,&dist_ik); - d_ik=v3_norm(&dist_ik); + d_ik2=v3_absolute_square(&dist_ik); /* ik constants */ brand=ai->brand; if(brand==ak->brand) { R=params->R[brand]; S=params->S[brand]; + S2=params->S2[brand]; } else { R=params->Rmixed; S=params->Smixed; + S2=params->S2mixed; } /* zeta_ij/dzeta_ij contribution only for d_ik < S */ - if(d_ikn_i); @@ -1704,8 +1738,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* d_costheta */ tmp=1.0/dd; - d_costheta1=cos_theta/(d_ij*d_ij)-tmp; - d_costheta2=cos_theta/(d_ik*d_ik)-tmp; + d_costheta1=cos_theta/d_ij2-tmp; + d_costheta2=cos_theta/d_ik2-tmp; /* some usefull values */ h_cos=(h-cos_theta); @@ -1757,13 +1791,14 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* dist_jk, d_jk */ v3_sub(&dist_jk,&(ak->r),&(aj->r)); if(bc) check_per_bound(moldyn,&dist_jk); - d_jk=v3_norm(&dist_jk); + d_jk2=v3_absolute_square(&dist_jk); /* jk constants */ brand=aj->brand; if(brand==ak->brand) { R=params->R[brand]; S=params->S[brand]; + S2=params->S2[brand]; B=params->B[brand]; mu=params->mu[brand]; chi=1.0; @@ -1771,13 +1806,17 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { else { R=params->Rmixed; S=params->Smixed; + S2=params->S2mixed; B=params->Bmixed; mu=params->mu_m; chi=params->chi; } /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */ - if(d_jkn_j); @@ -1795,7 +1834,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* d_costheta */ d_costheta1=1.0/dd; - d_costheta2=cos_theta/(d_ij*d_ij); + d_costheta2=cos_theta/d_ij2; /* some usefull values */ h_cos=(h-cos_theta); @@ -1805,10 +1844,11 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* g(cos_theta) */ g=1.0+c2d2-frac; - /* d_costheta_ij and dg(cos_theta) - needed in any case! */ + /* d_costheta_jik and dg(cos_theta) - needed in any case! */ v3_scale(&temp1,&dist_jk,d_costheta1); v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */ - v3_add(&temp1,&temp1,&temp2); + //v3_add(&temp1,&temp1,&temp2); + v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */ v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ /* store dg in temp2 and use it for dVjk later */ @@ -1873,6 +1913,14 @@ if(ai==&(moldyn->atom[0])) { printf("%f | %f\n",temp2.y,ai->f.y); printf("%f | %f\n",temp2.z,ai->f.z); } +#endif +#ifdef VDEBUG +if(ai==&(moldyn->atom[0])) { + printf("dVjk (3bp) contrib:\n"); + printf("%f | %f\n",temp2.x*dist_jk.x,ai->virial.xx); + printf("%f | %f\n",temp2.y*dist_jk.y,ai->virial.yy); + printf("%f | %f\n",temp2.z*dist_jk.z,ai->virial.zz); +} #endif } diff --git a/moldyn.h b/moldyn.h index dd94b0e..23b5100 100644 --- a/moldyn.h +++ b/moldyn.h @@ -195,6 +195,7 @@ typedef struct s_tersoff_exchange { double f_a,df_a; t_3dvec dist_ij; + double d_ij2; double d_ij; double chi; -- 2.20.1