X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=716a087fa90f1b89b073411af22e2d55995aba8a;hb=15dd7ee084a5db24a3e1a15428cd47e8d61fda5b;hp=af0132dff157564f46f593c3280700b045a1810c;hpb=a9fbc66448c52bc4138176739b33d17ba86b7eae;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index af0132d..716a087 100644 --- a/moldyn.c +++ b/moldyn.c @@ -139,6 +139,14 @@ int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) { return 0; } +int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) { + + moldyn->func2b=func; + moldyn->pot2b_params=params; + + return 0; +} + int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) { moldyn->func3b=func; @@ -940,9 +948,9 @@ int potential_force_calc(t_moldyn *moldyn) { if(ktom==&(itom[i])) continue; -printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x); +//printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x); moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk); -printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x); +//printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x); } while(list_next(that)!=\ L_NO_NEXT_ELEMENT); @@ -950,6 +958,13 @@ printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[ } } while(list_next(this)!=L_NO_NEXT_ELEMENT); + + /* 2bp post function */ + if(moldyn->func2b_post) + moldyn->func2b_post(moldyn, + &(itom[i]), + jtom,bc_ij); + } } } @@ -1072,6 +1087,21 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) { p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]); p->mu_m=0.5*(p->mu[0]+p->mu[1]); + printf("[moldyn] tersoff mult parameter info:\n"); + printf(" S (m) | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed); + printf(" R (m) | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed); + printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV); + printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV); + printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1], + p->lambda_m); + printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m); + printf(" beta | %.10f | %.10f\n",p->beta[0],p->beta[1]); + printf(" n | %f | %f\n",p->n[0],p->n[1]); + printf(" c | %f | %f\n",p->c[0],p->c[1]); + printf(" d | %f | %f\n",p->d[0],p->d[1]); + printf(" h | %f | %f\n",p->h[0],p->h[1]); + printf(" chi | %f \n",p->chi); + return 0; } @@ -1098,6 +1128,7 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { 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; @@ -1139,18 +1170,17 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { d_ij=v3_norm(&dist_ij); /* save for use in 3bp */ - exchange->dist_ij=dist_ij; /* <- needed ? */ exchange->d_ij=d_ij; + exchange->dist_ij=dist_ij; + exchange->d_ij2=d_ij*d_ij; /* constants */ if(num==aj->bnum) { S=params->S[num]; R=params->R[num]; A=params->A[num]; + B=params->B[num]; lambda=params->lambda[num]; - /* more constants depending of atoms i and j, needed in 3bp */ - params->exchange.B=&(params->B[num]); - params->exchange.mu=&(params->mu[num]); mu=params->mu[num]; params->exchange.chi=1.0; } @@ -1158,10 +1188,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { S=params->Smixed; R=params->Rmixed; A=params->Amixed; + B=params->Bmixed; lambda=params->lambda_m; - /* more constants depending of atoms i and j, needed in 3bp */ - params->exchange.B=&(params->Bmixed); - params->exchange.mu=&(params->mu_m); mu=params->mu_m; params->exchange.chi=params->chi; } @@ -1203,6 +1231,64 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* enable the run of 3bp function */ exchange->run3bp=1; + /* reset 3bp sums */ + exchange->sum1_3bp=0.0; + exchange->sum2_3bp=0.0; + v3_zero(&(exchange->db_ij)); + + return 0; +} + +/* tersoff 2 body post part */ + +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 */ + + 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; + double f_c,df_c,f_a,df_a; + double chi,betan; + double help; + double n; + + params=moldyn->pot2b_params; + exchange=&(params->exchange); + + 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; + betan=exchange->betan; + n=*(exchange->n); + dist_ij=&(exchange->dist_ij); + + db_ij_scale1=(1+betan*exchange->sum1_3bp); + db_ij_scale2=(exchange->n_betan*exchange->sum2_3bp); + help=pow(db_ij_scale1,-1.0/(2*n)-1); + b_ij=chi*db_ij_scale1*help; + db_ij_scale1=-chi/(2*n)*help; + + v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2)); + v3_scale(db_ij,db_ij,f_a); + + v3_scale(&temp,dist_ij,b_ij*df_a); + + v3_add(&force,&temp,db_ij); + v3_scale(&force,&force,f_c); + + v3_scale(&temp,dist_ij,f_a*b_ij*df_c); + + /* add energy of 3bp sum */ + moldyn->energy+=(0.5*f_c*b_ij*f_a); + /* add force of 3bp calculation */ + v3_add(&(ai->f),&temp,&force); + return 0; } @@ -1216,18 +1302,16 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_3dvec temp,force; double R,S,s_r; double d_ij,d_ij2,d_ik,d_jk; - double f_c,df_c,b_ij,f_a,df_a; + double f_c,df_c,f_a,df_a; double f_c_ik,df_c_ik,arg; - double scale; - double chi; - double n,c,d,h,beta,betan; + double n,c,d,h; double c2,d2,c2d2; double numer,denom; double theta,cos_theta,sin_theta; double d_theta,d_theta1,d_theta2; - double h_cos,h_cos2,d2_h_cos2; - double frac1,bracket1,bracket2,bracket2_n_1,bracket2_n; - double bracket3,bracket3_pow_1,bracket3_pow; + double h_cos,d2_h_cos2; + double frac,bracket,bracket_n_1,bracket_n; + double g; int num; params=moldyn->pot3b_params; @@ -1245,9 +1329,13 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { d_ij=exchange->d_ij; d_ij2=exchange->d_ij2; + dist_ij=exchange->dist_ij; f_a=params->exchange.f_a; df_a=params->exchange.df_a; + + f_c=exchange->f_c; + df_c=exchange->df_c; /* d_ij is <= S, as we didn't return so far! */ @@ -1261,7 +1349,6 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { * */ - v3_sub(&dist_ik,&(ai->r),&(ak->r)); if(bc) check_per_bound(moldyn,&dist_ik); d_ik=v3_norm(&dist_ik); @@ -1277,11 +1364,11 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { } /* calc of f_c_ik */ - if(d_ik>S) - return 0; - - if(d_ikS) { + f_c_ik=0.0; + df_c_ik=0.0; + } + else if(d_ikbeta); - betan=exchange->betan; + /* get exchange data */ n=*(exchange->n); c=*(exchange->c); d=*(exchange->d); @@ -1309,6 +1395,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { numer=d_ij2+d_ik*d_ik-d_jk*d_jk; denom=2*d_ij*d_ik; cos_theta=numer/denom; + /* prefere law of cosines, dot product -> nan (often) */ + //cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); sin_theta=sqrt(1.0-(cos_theta*cos_theta)); theta=acos(cos_theta); d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom); @@ -1318,50 +1406,43 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { d_theta2*=d_theta; h_cos=(h-cos_theta); - h_cos2=h_cos*h_cos; - d2_h_cos2=d2-h_cos2; - - /* some usefull expressions */ - frac1=c2/(d2-h_cos2); - bracket1=1+c2d2-frac1; - bracket2=f_c_ik*bracket1; - bracket2_n_1=pow(bracket2,n-1.0); - bracket2_n=bracket2_n_1*bracket2; - bracket3=1+betan*bracket2_n; - bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0); - bracket3_pow=bracket3_pow_1*bracket3; - - /* now go on with calc of b_ij and derivation of b_ij */ - b_ij=chi*bracket3_pow; + d2_h_cos2=d2+(h_cos*h_cos); + + frac=c2/(d2_h_cos2); + g=1.0+c2d2-frac; + + if(f_c_ik==0.0) { + bracket=0.0; + bracket_n_1=0.0; + bracket_n=0.0; + printf("Foo -> 0: "); + } + else { + bracket=f_c_ik*g; + bracket_n_1=pow(bracket,n-1.0); + bracket_n=bracket_n_1*bracket; + printf("Foo -> 1: "); + } +//printf("%.15f %.15f %.15f\n",bracket_n_1,bracket_n,bracket); + + /* calc of db_ij and the 2 sums */ + exchange->sum1_3bp+=bracket_n; + exchange->sum2_3bp+=bracket_n_1; /* derivation of theta */ v3_scale(&force,&dist_ij,d_theta1); v3_scale(&temp,&dist_ik,d_theta2); v3_add(&force,&force,&temp); - /* part 1 of derivation of b_ij */ - v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac1); + /* part 1 of db_ij */ + v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac/d2_h_cos2); - /* part 2 of derivation of b_ij */ - v3_scale(&temp,&dist_ik,df_c_ik*bracket1); + /* part 2 of db_ij */ + v3_scale(&temp,&dist_ik,df_c_ik*g); - /* sum up and scale ... */ + /* sum up and add to db_ij */ v3_add(&temp,&temp,&force); - scale=bracket2_n_1*n*betan*(1+betan*bracket3_pow_1)*chi*(1.0/(2.0*n)); - v3_scale(&temp,&temp,scale); - - /* now construct an energy and a force out of that */ - v3_scale(&temp,&temp,f_a); - v3_scale(&force,&dist_ij,df_a*b_ij); - v3_add(&temp,&temp,&force); - v3_scale(&temp,&temp,f_c); - v3_scale(&force,&dist_ij,df_c*b_ij*f_a); - v3_add(&force,&force,&temp); - - /* add forces */ - v3_add(&(ai->f),&(ai->f),&force); - /* energy is 0.5 f_r f_c, but we will sum it up twice ... */ - moldyn->energy+=(0.25*f_a*b_ij*f_c); + v3_add(&(exchange->db_ij),&(exchange->db_ij),&temp); return 0; }