X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=moldyn.c;h=4ddc42e98eb8bf1e6d2fa735ad4fd614de5f5f9b;hb=bb3e6f4ab3dfea9083e0d5b7d403ee6197f4041d;hp=8d1cb2334395528873284387debca19a6e210d2a;hpb=8bb8358a7dca3d77f5c504f9da57002d2cf30303;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 8d1cb23..4ddc42e 100644 --- a/moldyn.c +++ b/moldyn.c @@ -681,7 +681,7 @@ int moldyn_integrate(t_moldyn *moldyn) { /* zero absolute time */ moldyn->time=0.0; - /* debugging, ignre */ + /* debugging, ignore */ moldyn->debug=0; /* executing the schedule */ @@ -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 */ @@ -742,7 +741,7 @@ printf("MOVE\n"); if(!(i%v)) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); - printf("\rsched: %d, steps: %d, theta: %d", + printf("\rsched: %d, steps: %d, debug: %d", sched,i,moldyn->debug); fflush(stdout); } @@ -782,12 +781,15 @@ int velocity_verlet(t_moldyn *moldyn) { v3_add(&(atom[i].r),&(atom[i].r),&delta); v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass); v3_add(&(atom[i].r),&(atom[i].r),&delta); +//if(i==0) printf("und dann %.20f\n",atom[i].r.x*1e10); check_per_bound(moldyn,&(atom[i].r)); +//if(i==0) printf("und danach %.20f\n",atom[i].r.x); /* velocities */ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass); v3_add(&(atom[i].v),&(atom[i].v),&delta); } +moldyn_bc_check(moldyn); /* neighbour list update */ link_cell_update(moldyn); @@ -928,7 +930,7 @@ int potential_force_calc(t_moldyn *moldyn) { } } -printf("debug atom %d: %.15f %.15f %.15f\n",i,itom[i].r.x,itom[i].v.x,itom[i].f.x); +//printf("DEBUG: %.15f \n",itom[i].f.x); } return 0; @@ -1111,19 +1113,17 @@ 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; + + /* reset S > r > R mark */ + exchange->d_ij_between_rs=0; /* * calc of 2bp contribution of V_ij and dV_ij/ji @@ -1192,7 +1192,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* 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; + df_r=lambda*f_r/d_ij; /* f_a, df_a calc (again, same for ij and ji) | save for later use! */ exchange->f_a=-B*exp(-mu*d_ij); @@ -1213,6 +1213,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij)); /* two body contribution (ij, ji) */ v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c); + /* tell 3bp that S > r > R */ + exchange->d_ij_between_rs=1; } /* add forces of 2bp (ij, ji) contribution @@ -1233,10 +1235,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 +1245,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 +1275,70 @@ 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); - v3_scale(&force,&force,f_c); + + /* Vij and dVij */ + zeta=exchange->zeta_ij; + if(zeta==0.0) { + moldyn->debug++; /* just for debugging ... */ + db=0.0; + b=chi; + v3_scale(&force,dist_ij,df_a*b*f_c); + } + else { + 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); /* x(...)^(-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); + } + v3_scale(&temp,dist_ij,df_c*b*f_a); + v3_add(&force,&force,&temp); + v3_scale(&force,&force,-0.5); - /* df_c part */ - v3_scale(&temp,dist_ij,f_a*b_ij*df_c); + /* add force */ + v3_add(&(ai->f),&(ai->f),&force); /* add energy of 3bp sum */ - moldyn->energy+=(0.5*f_c*b_ij*f_a); + moldyn->energy+=(0.5*f_c*b*f_a); + + /* dVji */ + zeta=exchange->zeta_ji; + if(zeta==0.0) { + moldyn->debug++; + b=chi; + v3_scale(&force,dist_ij,df_a*b*f_c); + } + else { + 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); /* x(...)^(-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); + v3_scale(&force,&force,-0.5); - /* add force of 3bp calculation (all three parts) */ - v3_add(&(ai->f),&temp,&force); + /* add force */ + v3_sub(&(ai->f),&(ai->f),&force); return 0; } @@ -1314,17 +1351,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 +1440,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 +1464,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 +1474,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 +1489,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 +1506,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 +1531,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 +1553,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 +1576,55 @@ 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*0.5); + v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */ + /* scaled with 0.5 ^ */ + } + + return 0; +} + + +/* + * debugging / critical check functions + */ + +int moldyn_bc_check(t_moldyn *moldyn) { + + t_atom *atom; + t_3dvec *dim; + int i; + + atom=moldyn->atom; + dim=&(moldyn->dim); + + for(i=0;icount;i++) { + if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) + printf("FATAL: atom %d: x: %.20f (%.20f)\n", + i,atom[i].r.x*1e10,dim->x/2*1e10); + if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2) + printf("FATAL: atom %d: y: %.20f (%.20f)\n", + i,atom[i].r.y*1e10,dim->y/2*1e10); + if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2) + printf("FATAL: atom %d: z: %.20f (%.20f)\n", + i,atom[i].r.z*1e10,dim->z/2*1e10); } return 0;