From 8bb8358a7dca3d77f5c504f9da57002d2cf30303 Mon Sep 17 00:00:00 2001 From: hackbard Date: Sun, 10 Dec 2006 04:06:24 +0000 Subject: [PATCH] nearly finished 3bp tersoff routine (2bp post needs to be done after that!) --- moldyn.c | 357 ++++++++++++++++++++++++++++++++++++++----------------- moldyn.h | 37 ++++-- 2 files changed, 270 insertions(+), 124 deletions(-) diff --git a/moldyn.c b/moldyn.c index 880e56c..8d1cb23 100644 --- a/moldyn.c +++ b/moldyn.c @@ -1084,17 +1084,16 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { * their right values */ - exchange->beta=&(params->beta[num]); - exchange->n=&(params->n[num]); - exchange->c=&(params->c[num]); - exchange->d=&(params->d[num]); - exchange->h=&(params->h[num]); - - exchange->betan=pow(*(exchange->beta),*(exchange->n)); - exchange->n_betan=*(exchange->n)*exchange->betan; - exchange->c2=params->c[num]*params->c[num]; - exchange->d2=params->d[num]*params->d[num]; - exchange->c2d2=exchange->c2/exchange->d2; + exchange->beta_i=&(params->beta[num]); + exchange->n_i=&(params->n[num]); + exchange->c_i=&(params->c[num]); + exchange->d_i=&(params->d[num]); + exchange->h_i=&(params->h[num]); + + exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i)); + exchange->ci2=params->c[num]*params->c[num]; + exchange->di2=params->d[num]*params->d[num]; + exchange->ci2di2=exchange->ci2/exchange->di2; return 0; } @@ -1115,22 +1114,33 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { double scale; params=moldyn->pot2b_params; - num=ai->bnum; + num=aj->bnum; exchange=&(params->exchange); - exchange->run3bp=0; - exchange->run2bp_post=0; + /* clear 3bp and 2bp post run */ + exchange->run3bp_ij=0; + exchange->run3bp_ji=0; + exchange->run3bp_jk=0; + exchange->run2bp_post_ij=0; + exchange->run2bp_post_ji=0; + exchange->run2bp_post_jk=0; /* - * we need: f_c, df_c, f_r, df_r + * calc of 2bp contribution of V_ij and dV_ij/ji + * + * for Vij and dV_ij we need: + * - f_c_ij, df_c_ij + * - f_r_ij, df_r_ij + * + * for dV_ji we need: + * - f_c_ji = f_c_ij, df_c_ji = df_c_ij + * - f_r_ji = f_r_ij; df_r_ji = df_r_ij * - * therefore we need: R, S, A, lambda */ + /* dist_ij, d_ij */ v3_sub(&dist_ij,&(aj->r),&(ai->r)); - if(bc) check_per_bound(moldyn,&dist_ij); - d_ij=v3_norm(&dist_ij); /* save for use in 3bp */ @@ -1138,14 +1148,14 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->dist_ij=dist_ij; /* constants */ - if(num==aj->bnum) { + if(num==ai->bnum) { S=params->S[num]; R=params->R[num]; A=params->A[num]; B=params->B[num]; lambda=params->lambda[num]; mu=params->mu[num]; - params->exchange.chi=1.0; + exchange->chi=1.0; } else { S=params->Smixed; @@ -1156,34 +1166,60 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { mu=params->mu_m; params->exchange.chi=params->chi; } + + /* if d_ij > S => no force & potential energy contribution */ if(d_ij>S) return 0; + /* more constants */ + exchange->beta_j=&(params->beta[num]); + exchange->n_j=&(params->n[num]); + exchange->c_j=&(params->c[num]); + exchange->d_j=&(params->d[num]); + exchange->h_j=&(params->h[num]); + if(num==ai->bnum) { + exchange->betajnj=exchange->betaini; + exchange->cj2=exchange->ci2; + exchange->dj2=exchange->di2; + exchange->cj2dj2=exchange->ci2di2; + } + else { + exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j)); + exchange->cj2=params->c[num]*params->c[num]; + exchange->dj2=params->d[num]*params->d[num]; + exchange->cj2dj2=exchange->cj2/exchange->dj2; + } + + /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */ f_r=A*exp(-lambda*d_ij); df_r=-lambda*f_r/d_ij; - /* f_a, df_a calc + save for later use */ + /* f_a, df_a calc (again, same for ij and ji) | save for later use! */ exchange->f_a=-B*exp(-mu*d_ij); exchange->df_a=-mu*exchange->f_a/d_ij; + /* f_c, df_c calc (again, same for ij and ji) */ if(d_ijf),&(ai->f),&force); - /* energy is 0.5 f_r f_c ... */ + + /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */ moldyn->energy+=(0.5*f_r*f_c); /* save for use in 3bp */ @@ -1195,8 +1231,12 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->run2bp_post=1; /* reset 3bp sums */ - exchange->zeta=0.0; + exchange->zeta_ij=0.0; + exchange->zeta_ji=0.0; + exchange->zeta_kl=0.0; v3_zero(&(exchange->db_ij)); + v3_zero(&(exchange->db_ji)); + v3_zero(&(exchange->db_jk)); return 0; } @@ -1272,11 +1312,11 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_tersoff_mult_params *params; t_tersoff_exchange *exchange; - t_3dvec dist_ij,dist_ik; - t_3dvec temp,force; + t_3dvec dist_ij,dist_ik,dist_jk; + t_3dvec temp1,temp2; double R,S,s_r; - double d_ij,d_ik; - double rijrik,dijdik; + double d_ij,d_ik,d_jk; + double rxxryy,dxxdyy; double f_c,df_c,f_a,df_a; double f_c_ik,df_c_ik,arg; double n,c,d,h; @@ -1288,43 +1328,56 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { int num; params=moldyn->pot3b_params; - num=ai->bnum; exchange=&(params->exchange); if(!(exchange->run3bp)) return 0; /* - * we need: f_c, d_fc, b_ij, db_ij, f_a, df_a + * calc of 3bp contribution of V_ij and dV_ij/ji/jk & + * 2bp contribution of dV_jk + * + * for Vij and dV_ij we still need: + * - b_ij, db_ij (zeta_ij) + * - f_c_ik, df_c_ik, constants_i, cos_theta_ijk, d_costheta_ijk + * + * for dV_ji we still need: + * - b_ji, db_ji (zeta_ji) + * - f_c_jk, d_c_jk, constants_j, cos_theta_jik, d_costheta_jik + * + * for dV_jk we need: + * - f_c_jk + * - f_a_jk + * - db_jk (zeta_jk) + * - f_c_ji, df_c_ji, constants_j, cos_theta_jki, d_costheta_jki * - * we got f_c, df_c, f_a, df_a from 2bp calculation */ - d_ij=exchange->d_ij; - dist_ij=exchange->dist_ij; + /* + * get exchange data + */ - f_a=params->exchange.f_a; - df_a=params->exchange.df_a; + /* dist_ij, d_ij - this is < S_ij ! */ + dist_ij=exchange->dist_ij; + d_ij=exchange->d_ij; + /* f_c_ij, df_c_ij (same for ji) */ f_c=exchange->f_c; df_c=exchange->df_c; - - /* d_ij is <= S, as we didn't return so far! */ /* - * calc of b_ij (scalar) and db_ij (vector) - * - * - for b_ij: chi, beta, f_c_ik, w(=1), c, d, h, n, cos_theta - * - * - for db_ij: d_costheta, cos_theta, f_c_ik, df_c_ik, w_ik - * + * calculate unknown values now ... */ + /* V_ij and dV_ij stuff (in b_ij there is f_c_ik) */ + + /* 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); - /* constants */ + /* ik constants */ + num=ai->bnum; if(num==ak->bnum) { R=params->R[num]; S=params->S[num]; @@ -1334,84 +1387,166 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { S=params->Smixed; } - /* there is no contribution if f_c_ik = 0 */ - if(d_ik>S) - return 0; - - /* get exchange data */ - n=*(exchange->n); - c=*(exchange->c); - d=*(exchange->d); - h=*(exchange->h); - c2=exchange->c2; - d2=exchange->d2; - c2d2=exchange->c2d2; - - /* cosine of theta by scalaproduct */ - rijrik=v3_scalar_product(&dist_ij,&dist_ik); - dijdik=d_ij*d_ik; - cos_theta=rijrik/dijdik; - - /* hack - cos theta machine accuracy problems! */ - if(cos_theta>1.0||cos_theta<-1.0) { - printf("THETA CORRECTION\n"); - moldyn->debug++; - if(fabs(cos_theta)>1.0+ACCEPTABLE_ERROR) - printf("[moldyn] WARNING: cos theta failure!\n"); - if(cos_theta<0) { - cos_theta=-1.0; + /* zeta_ij/dzeta_ij contribution only for d_ik < S */ + if(d_ikn_i); + c=*(exchange->c_i); + d=*(exchange->d_i); + h=*(exchange->h_i); + c2=exchange->ci2; + d2=exchange->di2; + c2d2=exchange->ci2di2; + + /* cosine of theta_ijk by scalaproduct */ + rijrik=v3_scalar_product(&dist_ij,&dist_ik); + dijdik=d_ij*d_ik; + cos_theta=rijrik/dijdik; + + /* d_costheta */ + tmp=1.0/dijdik; + d_costheta1=cos_theta/(d_ij*d_ij)-tmp; + d_costheta2=cos_theta/(d_ik*d_ik)-tmp; + + /* some usefull values */ + h_cos=(h-cos_theta); + d2_h_cos2=d2+(h_cos*h_cos); + frac=c2/(d2_h_cos2); + + /* g(cos_theta) */ + g=1.0+c2d2-frac; + + /* d_costheta_ij and dg(cos_theta) - needed in any case! */ + v3_scale(&temp1,&dist_ij,d_costheta1); + v3_scale(&temp2,&dist_ik,d_costheta2); + v3_add(&temp1,&temp1,&temp2); + v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ + + /* f_c_ik & df_c_ik + {d,}zeta contribution */ + if(d_ik f_c_ik=1.0; + // => df_c_ik=0.0; of course we do not set this! + + /* zeta_ij */ + exchange->zeta_ij+=g; + + /* dzeta_ij */ + v3_add(dzeta_ij,dzeta_ij,&temp1); } else { - cos_theta=1.0; + /* {d,}f_c_ik */ + s_r=S-R; + arg=M_PI*(d_ik-R)/s_r; + f_c_ik=0.5+0.5*cos(arg); + df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); + + /* zeta_ij */ + exchange->zeta_ij+=f_c_ik*g; + + /* dzeta_ij */ + v3_scale(&temp1,&temp1,f_c_ik); + v3_scale(&temp2,&dist_ik,g*df_c_ik); + v3_add(dzeta_ij,&temp2,&temp1); } } - d_costheta1=dijdik-rijrik*d_ik/d_ij; - d_costheta2=dijdik-rijrik*d_ij/d_ik; + /* dV_ji stuff (in b_ji there is f_c_jk) + dV_jk stuff! */ + + /* 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); - h_cos=(h-cos_theta); - d2_h_cos2=d2+(h_cos*h_cos); + /* jk constants */ + num=aj->bnum; + if(num==ak->bnum) { + R=params->R[num]; + S=params->S[num]; + } + else { + R=params->Rmixed; + S=params->Smixed; + } - frac=c2/(d2_h_cos2); - g=1.0+c2d2-frac; + /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */ + if(d_jkn_j); + c=*(exchange->c_j); + d=*(exchange->d_j); + h=*(exchange->h_j); + c2=exchange->cj2; + d2=exchange->dj2; + c2d2=exchange->cj2dj2; + + /* cosine of theta_jik by scalaproduct */ + rxxryy=v3_scalar_product(&dist_ij,&dist_jk); /* times -1 */ + dxxdyy=d_ij*d_jk; + cos_theta=rxxryy/dxxdyy; + + /* d_costheta */ + d_costheta1=1.0/(d_jk*d_ij); + d_costheta2=cos_theta/(d_ij*d_ij); /* in fact -cos(), but ^ */ + + /* some usefull values */ + h_cos=(h-cos_theta); + d2_h_cos2=d2+(h_cos*h_cos); + frac=c2/(d2_h_cos2); + + /* g(cos_theta) */ + g=1.0+c2d2-frac; + + /* d_costheta_ij 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_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ + + /* dV_jk stuff | add force contribution on atom i immediately */ + if(exchange->d_ij_between_rs) { + tmp=pow(f_c_ij*g,(n_j-1.0)); /* zeta_jk ^ n_j-1 */ + v3_scale(&temp2,&temp1,f_c_ij) + v3_scale(&temp3,&dist_ij,df_c_ij); + v3_add(&temp3,&temp3,&temp2); /* dzeta_jk */ + } + else { + /* f_c_ij = 1, df_c_ij = 0 */ + tmp=pow(g,(n_j-1.0)); /* zeta_jk ^ n_j-1 */ + tmp + /* dzeta_jk in temp1 */ + /* HIER WEITER !!! */ + } - /* d_costheta contrib to db_ij (needed in all remaining cases) */ - v3_scale(&temp,&dist_ij,d_costheta1); - v3_scale(&force,&dist_ik,d_costheta2); - v3_add(&force,&force,&temp); - v3_scale(&force,&force,-2.0*frac*h_cos/d2_h_cos2); /* f_c_ik missing */ + /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */ + if(d_jk only d_costheta contrib to db_ij */ - // => do nothing ... + /* zeta_ji */ + exchange->zeta_ji+=g; - /* zeta, f_c_ik = 1 */ - exchange->zeta+=g; - } - else { - s_r=S-R; - arg=M_PI*(d_ik-R)/s_r; - f_c_ik=0.5+0.5*cos(arg); - df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); - /* scale d_costheta contrib with f_c_ik */ - v3_scale(&force,&force,f_c_ik); + /* dzeta_ji */ + v3_add(dzeta_ji,dzeta_ji,&temp1); + } + else { + /* f_c_jk */ + s_r=S-R; + arg=M_PI*(d_jk-R)/s_r; + f_c_jk=0.5+0.5*cos(arg); - /* df_c_ik contrib to db_ij */ - v3_scale(&temp,&dist_ik,df_c_ik*g); + /* zeta_ji */ + exchange->zeta_ji+=f_c_jk*g; - /* sum up both parts */ - v3_add(&force,&force,&temp); - - /* zeta */ - exchange->zeta+=f_c_ik*g; + /* dzeta_ij */ + v3_scale(&temp1,&temp1,f_c_jk); + v3_add(dzeta_ji,dzeta_ji,&temp1); + } } -printf("%.30f\n",exchange->zeta); - - /* add to db_ij */ - v3_add(&(exchange->db_ij),&(exchange->db_ij),&force); - + return 0; } diff --git a/moldyn.h b/moldyn.h index 32025fd..d62ec47 100644 --- a/moldyn.h +++ b/moldyn.h @@ -183,23 +183,34 @@ typedef struct s_tersoff_exchange { double chi; - double *beta; - double *n; - double *c; - double *d; - double *h; - - double c2; - double d2; - double c2d2; - double betan; - double n_betan; + double *beta_i; + double *beta_j; + double *n_i; + double *n_j; + double *c_i; + double *c_j; + double *d_i; + double *d_j; + double *h_i; + double *h_j; + + double ci2; + double cj2; + double di2; + double dj2; + double ci2di2; + double cj2dj2; + double betaini; + double betajnj; u8 run3bp; u8 run2bp_post; + u8 d_ij_between_rs; - t_3dvec db_ij; - double zeta; + double zeta_ij; + double zeta_ji; + t_3dvec dzeta_ij; + t_3dvec dzeta_ji; } t_tersoff_exchange; /* tersoff multi (2!) potential parameters */ -- 2.20.1