From b30d3d558354102cf5d0c06392961032b21bf9a9 Mon Sep 17 00:00:00 2001 From: hackbard Date: Sun, 10 Dec 2006 22:44:50 +0000 Subject: [PATCH] mods to corect 3bp derivative, still segfaulting ... --- moldyn.c | 182 +++++++++++++++++++++++++++++++------------------------ sic.c | 4 +- 2 files changed, 106 insertions(+), 80 deletions(-) diff --git a/moldyn.c b/moldyn.c index 8d1cb23..6736e45 100644 --- a/moldyn.c +++ b/moldyn.c @@ -697,7 +697,6 @@ int moldyn_integrate(t_moldyn *moldyn) { for(i=0;itime_steps;i++) { /* integration step */ -printf("MOVE\n"); moldyn->integrate(moldyn); /* p/t scaling */ @@ -921,14 +920,15 @@ int potential_force_calc(t_moldyn *moldyn) { /* 2bp post function */ if(moldyn->func2b_post) { +printf("DEBUG: vor 2bp post\n"); moldyn->func2b_post(moldyn, &(itom[i]), jtom,bc_ij); +printf("DEBUG: nach 2bp post\n"); } } } -printf("debug atom %d: %.15f %.15f %.15f\n",i,itom[i].r.x,itom[i].v.x,itom[i].f.x); } return 0; @@ -1111,19 +1111,14 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { int num; double s_r; double arg; - double scale; params=moldyn->pot2b_params; num=aj->bnum; exchange=&(params->exchange); /* 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; + exchange->run3bp=0; + exchange->run2bp_post=0; /* * calc of 2bp contribution of V_ij and dV_ij/ji @@ -1233,10 +1228,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* reset 3bp sums */ 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)); + v3_zero(&(exchange->dzeta_ij)); + v3_zero(&(exchange->dzeta_ji)); return 0; } @@ -1245,16 +1238,27 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { - /* here we have to allow for the 3bp sums */ + /* + * here we have to allow for the 3bp sums + * + * that is: + * - zeta_ij, dzeta_ij + * - zeta_ji, dzeta_ji + * + * to compute the 3bp contribution to: + * - Vij, dVij + * - dVji + * + */ t_tersoff_mult_params *params; t_tersoff_exchange *exchange; - t_3dvec force,temp,*db_ij,*dist_ij; - double db_ij_scale1,db_ij_scale2; - double b_ij; + t_3dvec force,temp; + t_3dvec *dist_ij; + double b,db,tmp; double f_c,df_c,f_a,df_a; - double chi,n,n_betan; + double chi,ni,betaini,nj,betajnj; double zeta; params=moldyn->pot2b_params; @@ -1264,44 +1268,53 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { if(!(exchange->run2bp_post)) return 0; - db_ij=&(exchange->db_ij); f_c=exchange->f_c; df_c=exchange->df_c; f_a=exchange->f_a; df_a=exchange->df_a; - n_betan=exchange->n_betan; - n=*(exchange->n); + betaini=exchange->betaini; + betajnj=exchange->betajnj; + ni=*(exchange->n_i); + nj=*(exchange->n_j); chi=exchange->chi; dist_ij=&(exchange->dist_ij); - zeta=exchange->zeta; - - db_ij_scale2=pow(zeta,n-1.0); -printf("DEBUG: %.15f %.15f\n",zeta,db_ij_scale2); - db_ij_scale1=db_ij_scale2*zeta; - db_ij_scale2*=n_betan; - db_ij_scale1=pow((1.0+n_betan*db_ij_scale1),-1.0/(2*n)-1); - b_ij=chi*db_ij_scale1*(1.0+n_betan*db_ij_scale1); - db_ij_scale1*=(-1.0*chi/(2*n)); - - /* db_ij part */ - v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2)); - v3_scale(db_ij,db_ij,f_a); - - /* df_a part */ - v3_scale(&temp,dist_ij,b_ij*df_a); - - /* db_ij + df_a part */ - v3_add(&force,&temp,db_ij); + + /* Vij and dVij */ + zeta=exchange->zeta_ij; + tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */ + b=(1+zeta*tmp); /* 1 + beta^n * zeta^n */ + db=chi*pow(b,-1.0/(2*ni)-1); /* chi * (...)^(-1/2n - 1) */ + b=db*b; /* b_ij */ + db*=-0.5*tmp; /* db_ij */ + v3_scale(&force,&(exchange->dzeta_ij),f_a*db); + v3_scale(&temp,dist_ij,df_a*b); + v3_add(&force,&force,&temp); v3_scale(&force,&force,f_c); - - /* df_c part */ - v3_scale(&temp,dist_ij,f_a*b_ij*df_c); + v3_scale(&temp,dist_ij,df_c*b*f_a); + v3_add(&force,&force,&temp); /* add energy of 3bp sum */ - moldyn->energy+=(0.5*f_c*b_ij*f_a); + moldyn->energy+=(0.5*f_c*b*f_a); + + /* add force (sub, as F = - dVij) */ + v3_sub(&(ai->f),&(ai->f),&force); + + /* dVji */ + zeta=exchange->zeta_ji; + tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */ + b=(1+zeta*tmp); /* 1 + beta^n * zeta^n */ + db=chi*pow(b,-1.0/(2*nj)-1); /* chi * (...)^(-1/2n - 1) */ + b=db*b; /* b_ij */ + db*=-0.5*tmp; /* db_ij */ + v3_scale(&force,&(exchange->dzeta_ji),f_a*db); + v3_scale(&temp,dist_ij,df_a*b); + v3_add(&force,&force,&temp); + v3_scale(&force,&force,f_c); + v3_scale(&temp,dist_ij,df_c*b*f_a); + v3_add(&force,&force,&temp); - /* add force of 3bp calculation (all three parts) */ - v3_add(&(ai->f),&temp,&force); + /* add force (sub, as F = - dVji) */ + v3_sub(&(ai->f),&(ai->f),&force); return 0; } @@ -1314,17 +1327,20 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_tersoff_exchange *exchange; t_3dvec dist_ij,dist_ik,dist_jk; t_3dvec temp1,temp2; + t_3dvec *dzeta; double R,S,s_r; + double B,mu; double d_ij,d_ik,d_jk; - double rxxryy,dxxdyy; - double f_c,df_c,f_a,df_a; + double rr,dd; + double f_c,df_c; double f_c_ik,df_c_ik,arg; + double f_c_jk; double n,c,d,h; double c2,d2,c2d2; double cos_theta,d_costheta1,d_costheta2; double h_cos,d2_h_cos2; - double frac; - double g; + double frac,g,zeta,chi; + double tmp; int num; params=moldyn->pot3b_params; @@ -1400,12 +1416,12 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { 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; + rr=v3_scalar_product(&dist_ij,&dist_ik); + dd=d_ij*d_ik; + cos_theta=rr/dd; /* d_costheta */ - tmp=1.0/dijdik; + tmp=1.0/dd; d_costheta1=cos_theta/(d_ij*d_ij)-tmp; d_costheta2=cos_theta/(d_ik*d_ik)-tmp; @@ -1424,6 +1440,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ /* f_c_ik & df_c_ik + {d,}zeta contribution */ + dzeta=&(exchange->dzeta_ij); if(d_ik f_c_ik=1.0; @@ -1433,7 +1450,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { exchange->zeta_ij+=g; /* dzeta_ij */ - v3_add(dzeta_ij,dzeta_ij,&temp1); + v3_add(dzeta,dzeta,&temp1); } else { /* {d,}f_c_ik */ @@ -1448,7 +1465,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* dzeta_ij */ v3_scale(&temp1,&temp1,f_c_ik); v3_scale(&temp2,&dist_ik,g*df_c_ik); - v3_add(dzeta_ij,&temp2,&temp1); + v3_add(&temp1,&temp1,&temp2); + v3_add(dzeta,dzeta,&temp1); } } @@ -1464,10 +1482,16 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { if(num==ak->bnum) { R=params->R[num]; S=params->S[num]; + B=params->B[num]; + mu=params->mu[num]; + chi=1.0; } else { R=params->Rmixed; S=params->Smixed; + B=params->Bmixed; + mu=params->mu_m; + chi=params->chi; } /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */ @@ -1483,9 +1507,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { 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; + rr=v3_scalar_product(&dist_ij,&dist_jk); /* times -1 */ + dd=d_ij*d_jk; + cos_theta=rr/dd; /* d_costheta */ d_costheta1=1.0/(d_jk*d_ij); @@ -1505,32 +1529,17 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { 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 !!! */ - } - /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */ + dzeta=&(exchange->dzeta_ji); if(d_jkzeta_ji+=g; - /* dzeta_ji */ - v3_add(dzeta_ji,dzeta_ji,&temp1); + v3_add(dzeta,dzeta,&temp1); } else { /* f_c_jk */ @@ -1543,8 +1552,25 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* dzeta_ij */ v3_scale(&temp1,&temp1,f_c_jk); - v3_add(dzeta_ji,dzeta_ji,&temp1); + v3_add(dzeta,dzeta,&temp1); + } + + /* dV_jk stuff | add force contribution on atom i immediately */ + if(exchange->d_ij_between_rs) { + zeta=f_c*g; + v3_scale(&temp1,&temp1,f_c); + v3_scale(&temp2,&dist_ij,df_c); + v3_add(&temp1,&temp1,&temp2); + } + else { + zeta=g; + // dzeta_jk is simply dg, which is temp1 } + /* betajnj * zeta_jk ^ nj-1 */ + tmp=exchange->betajnj*pow(zeta,(n-1.0)); + tmp=-chi/2.0*pow(1+tmp*zeta,-1.0/(2.0*n)-1)*tmp; + v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk); + v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */ } return 0; diff --git a/sic.c b/sic.c index 8073559..d141782 100644 --- a/sic.c +++ b/sic.c @@ -96,7 +96,7 @@ int main(int argc,char **argv) { /* set (initial) dimensions of simulation volume */ printf("[sic] setting dimensions\n"); - set_dim(&md,5*LC_SI,5*LC_SI,5*LC_SI,TRUE); + set_dim(&md,2*LC_SI,2*LC_SI,2*LC_SI,TRUE); /* set periodic boundary conditions in all directions */ printf("[sic] setting periodic boundary conditions\n"); @@ -107,7 +107,7 @@ int main(int argc,char **argv) { create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, //ATOM_ATTR_2BP|ATOM_ATTR_HB, - 0,5,5,5); + 0,2,2,2); /* setting a nearest neighbour distance for the moldyn checks */ set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */ -- 2.20.1