From 57847266c05bc38218bf162efdb08e8dd40894d7 Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 3 May 2007 16:40:02 +0000 Subject: [PATCH] put repulsive part to 3bp/j2 function --- moldyn.c | 2 +- potentials/tersoff.c | 40 ++++++++++++++++++++++++---------------- sic.c | 2 +- 3 files changed, 26 insertions(+), 18 deletions(-) diff --git a/moldyn.c b/moldyn.c index 4d73617..705a28c 100644 --- a/moldyn.c +++ b/moldyn.c @@ -1711,7 +1711,7 @@ int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) { } /* - * periodic boundayr checking + * periodic boundary checking */ //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { diff --git a/potentials/tersoff.c b/potentials/tersoff.c index b6533e1..a4a87e3 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -94,6 +94,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { double s_r; double arg; + printf("WARNING! - tersoff_mult_2bp is obsolete.\n"); + printf("WARNING! - repulsive part handled in 3bp/j2 routine.\n"); + /* use newtons third law */ if(aivirial.xx-=force.x*dist_ij.x; - //ai->virial.yy-=force.y*dist_ij.y; - //ai->virial.zz-=force.z*dist_ij.z; - //ai->virial.xy-=force.x*dist_ij.y; - //ai->virial.xz-=force.x*dist_ij.z; - //ai->virial.yz-=force.y*dist_ij.z; + virial_calc(ai,&force,&dist_ij); /* energy 2bp contribution */ moldyn->energy+=f_r*f_c; @@ -336,7 +332,10 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_tersoff_exchange *exchange; t_3dvec force; double f_a,df_a,b,db,f_c,df_c; + double f_r,df_r; + double scale; double mu,B,chi; + double lambda,A; double d_ij; unsigned char brand; double ni,tmp; @@ -350,14 +349,18 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { S=params->S[brand]; R=params->R[brand]; B=params->B[brand]; + A=params->A[brand]; mu=params->mu[brand]; + lambda=params->lambda[brand]; chi=1.0; } else { S=params->Smixed; R=params->Rmixed; B=params->Bmixed; + A=params->Amixed; mu=params->mu_m; + lambda=params->lambda_m; chi=params->chi; } @@ -379,6 +382,10 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { f_a=-B*exp(-mu*d_ij); df_a=mu*f_a/d_ij; + /* f_r, df_r */ + f_r=A*exp(-lambda*d_ij); + df_r=lambda*f_r/d_ij; + /* b, db */ if(exchange->zeta_ij==0.0) { b=chi; @@ -394,8 +401,8 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } /* force contribution */ - v3_scale(&force,&(exchange->dist_ij),df_a*f_c+f_a*df_c); - v3_scale(&force,&force,-0.5*b); + scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_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 @@ -411,14 +418,13 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { #endif /* virial */ - //virial_calc(ai,&force,&(exchange->dist_ij)); - //virial_calc(aj,&force,&(exchange->dist_ij)); + 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; /* energy contribution */ - moldyn->energy+=0.5*f_c*b*f_a; + moldyn->energy+=0.5*f_c*(b*f_a+f_r); /* reset k counter for second k loop */ exchange->kcount=0; @@ -537,10 +543,11 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, #endif /* virial */ - //virial_calc(aj,&force,&dist_ij); + //v3_scale(&force,&force,-1.0); + virial_calc(ai,&force,&dist_ij); /* derivative wrt k */ - v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rij = - drj rij + v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik v3_scale(&tmp,&dcosdrk,fcdg); v3_add(&force,&force,&tmp); v3_scale(&force,&force,pre_dzeta); @@ -557,7 +564,8 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn, #endif /* virial */ - virial_calc(ak,&force,&dist_ik); + //v3_scale(&force,&force,-1.0); + virial_calc(ai,&force,&dist_ik); /* increase k counter */ exchange->kcount++; diff --git a/sic.c b/sic.c index 714710b..2221025 100644 --- a/sic.c +++ b/sic.c @@ -97,7 +97,7 @@ int main(int argc,char **argv) { /* choose potential */ set_potential1b(&md,tersoff_mult_1bp); - set_potential2b(&md,tersoff_mult_2bp); + //set_potential2b(&md,tersoff_mult_2bp); //set_potential3b_j1(&md,tersoff_mult_2bp); //set_potential3b_k1(&md,tersoff_mult_3bp); //set_potential3b_j3(&md,tersoff_mult_post_2bp); -- 2.20.1