X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=f09c93a27e10fb00dc7a4dc6087bb1cb3f921948;hb=0acbade9c8d1ed354c333955dc2c5575540a9236;hp=2631d6a675549e183a41fa91ca863dee215bed0f;hpb=746278fb5c566b621b70c78b38d8b6c0b4e2677f;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 2631d6a..f09c93a 100644 --- a/moldyn.c +++ b/moldyn.c @@ -141,7 +141,7 @@ int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) { int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) { - moldyn->func2b=func; + moldyn->func2b_post=func; moldyn->pot2b_params=params; return 0; @@ -576,9 +576,8 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { } lc->dnlc=count1; - lc->countn=27; - return count2; + return count1; } int link_cell_shutdown(t_moldyn *moldyn) { @@ -664,10 +663,11 @@ int moldyn_integrate(t_moldyn *moldyn) { /* sqaure of some variables */ moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; + /* calculate initial forces */ potential_force_calc(moldyn); - /* do some checks before we actually start calculating bullshit */ + /* some stupid checks before we actually start calculating bullshit */ if(moldyn->cutoff>0.5*moldyn->dim.x) printf("[moldyn] warning: cutoff > 0.5 x dim.x\n"); if(moldyn->cutoff>0.5*moldyn->dim.y) @@ -680,6 +680,11 @@ int moldyn_integrate(t_moldyn *moldyn) { /* zero absolute time */ moldyn->time=0.0; + + /* debugging, ignre */ + moldyn->debug=0; + + /* executing the schedule */ for(sched=0;schedschedule.content_count;sched++) { /* setting amount of runs and finite time step size */ @@ -692,6 +697,7 @@ int moldyn_integrate(t_moldyn *moldyn) { for(i=0;itime_steps;i++) { /* integration step */ +printf("MOVE\n"); moldyn->integrate(moldyn); /* p/t scaling */ @@ -736,7 +742,8 @@ int moldyn_integrate(t_moldyn *moldyn) { if(!(i%v)) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); - printf("\rsched: %d, steps: %d",sched,i); + printf("\rsched: %d, steps: %d, theta: %d", + sched,i,moldyn->debug); fflush(stdout); } } @@ -747,6 +754,9 @@ int moldyn_integrate(t_moldyn *moldyn) { if(schedule->hook) schedule->hook(moldyn,schedule->hook_params); + /* get a new info line */ + printf("\n"); + } return 0; @@ -784,7 +794,6 @@ int velocity_verlet(t_moldyn *moldyn) { /* forces depending on chosen potential */ potential_force_calc(moldyn); - //moldyn->potential_force_function(moldyn); for(i=0;icount; itom=moldyn->atom; @@ -823,6 +832,7 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset energy */ moldyn->energy=0.0; + /* get energy and force of every atom */ for(i=0;idim.z/2)/lc->z, neighbour_i); - countn=lc->countn; dnlc=lc->dnlc; - for(j=0;jattr&ATOM_ATTR_3BP)) continue; - /* neighbourhood of atom j is not needed! */ - - // link_cell_neighbour_index(moldyn, - // (jtom->r.x+moldyn->dim.x/2)/lc->x, - // (jtom->r.y+moldyn->dim.y/2)/lc->y, - // (jtom->r.z+moldyn->dim.z/2)/lc->z, - // neighbour_j); - -// /* neighbours of j */ -// for(k=0;kcountn;k++) { -// -// that=&(neighbour_j[k]); -// list_reset(that); -// -// if(that->start==NULL) -// continue; -// -// bc_ijk=(kdnlc)?0:1; -// -// do { -// -// ktom=that->current->data; -// -// if(!(ktom->attr&ATOM_ATTR_3BP)) -// continue; -// -// if(ktom==jtom) -// continue; -// -// if(ktom==&(itom[i])) -// continue; -// -// moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk); -// -/* } while(list_next(that)!=\ */ -// L_NO_NEXT_ELEMENT); -// -// } - /* copy the neighbour lists */ memcpy(neighbour_i2,neighbour_i, 27*sizeof(t_list)); /* get neighbours of i */ - for(k=0;kstart==NULL) continue; - bc_ijk=(kfunc3b(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); + moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ik|bc_ij); } while(list_next(that)!=\ L_NO_NEXT_ELEMENT); @@ -952,13 +920,15 @@ int potential_force_calc(t_moldyn *moldyn) { } while(list_next(this)!=L_NO_NEXT_ELEMENT); /* 2bp post function */ - if(moldyn->func2b_post) + if(moldyn->func2b_post) { moldyn->func2b_post(moldyn, &(itom[i]), jtom,bc_ij); + } } } +printf("debug atom %d: %.15f %.15f %.15f\n",i,itom[i].r.x,itom[i].v.x,itom[i].f.x); } return 0; @@ -1165,7 +1135,6 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* save for use in 3bp */ exchange->d_ij=d_ij; exchange->dist_ij=dist_ij; - exchange->d_ij2=d_ij*d_ij; /* constants */ if(num==aj->bnum) { @@ -1186,7 +1155,6 @@ 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) return 0; @@ -1226,8 +1194,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->run2bp_post=1; /* reset 3bp sums */ - exchange->sum1_3bp=0.0; - exchange->sum2_3bp=0.0; + exchange->zeta=0.0; v3_zero(&(exchange->db_ij)); return 0; @@ -1246,14 +1213,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { 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; + double chi,n,n_betan; + double zeta; params=moldyn->pot2b_params; exchange=&(params->exchange); - /* we do not run if f_c_ij was dtected to be 0! */ + /* we do not run if f_c_ij was detected to be 0! */ if(!(exchange->run2bp_post)) return 0; @@ -1262,16 +1228,19 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { df_c=exchange->df_c; f_a=exchange->f_a; df_a=exchange->df_a; - betan=exchange->betan; + n_betan=exchange->n_betan; n=*(exchange->n); chi=exchange->chi; dist_ij=&(exchange->dist_ij); + zeta=exchange->zeta; - 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; + 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)); @@ -1302,19 +1271,18 @@ 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,dist_jk; + t_3dvec dist_ij,dist_ik; t_3dvec temp,force; double R,S,s_r; - double d_ij,d_ij2,d_ik,d_jk; + double d_ij,d_ik; + double rijrik,dijdik; double f_c,df_c,f_a,df_a; double f_c_ik,df_c_ik,arg; 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 cos_theta,d_costheta1,d_costheta2; double h_cos,d2_h_cos2; - double frac,bracket,bracket_n_1,bracket_n; + double frac; double g; int num; @@ -1332,7 +1300,6 @@ 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; @@ -1348,8 +1315,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { * * - for b_ij: chi, beta, f_c_ik, w(=1), c, d, h, n, cos_theta * - * - for db_ij: d_theta, sin_theta, cos_theta, f_c_ik, df_c_ik, - * w_ik, + * - for db_ij: d_costheta, cos_theta, f_c_ik, df_c_ik, w_ik * */ @@ -1357,7 +1323,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { if(bc) check_per_bound(moldyn,&dist_ik); d_ik=v3_norm(&dist_ik); - /* constants for f_c_ik calc */ + /* constants */ if(num==ak->bnum) { R=params->R[num]; S=params->S[num]; @@ -1367,25 +1333,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { S=params->Smixed; } - /* calc of f_c_ik */ - if(d_ik>S) { - f_c_ik=0.0; - df_c_ik=0.0; - } - else if(d_ikr),&(ak->r)); - if(bc) check_per_bound(moldyn,&dist_jk); - d_jk=v3_norm(&dist_jk); + /* there is no contribution if f_c_ik = 0 */ + if(d_ik>S) + return 0; /* get exchange data */ n=*(exchange->n); @@ -1396,18 +1346,27 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { d2=exchange->d2; c2d2=exchange->c2d2; - 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); - d_theta1=2*denom-numer*2*d_ik/d_ij; - d_theta2=2*denom-numer*2*d_ij/d_ik; - d_theta1*=d_theta; - d_theta2*=d_theta; + /* 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; + } + else { + cos_theta=1.0; + } + } + + d_costheta1=dijdik-rijrik*d_ik/d_ij; + d_costheta2=dijdik-rijrik*d_ij/d_ik; h_cos=(h-cos_theta); d2_h_cos2=d2+(h_cos*h_cos); @@ -1415,38 +1374,42 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { 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); + /* 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 */ - /* calc of db_ij and the 2 sums */ - exchange->sum1_3bp+=bracket_n; - exchange->sum2_3bp+=bracket_n_1; + if(d_ik only d_costheta contrib to db_ij */ + // => do nothing ... - /* derivation of theta */ - v3_scale(&force,&dist_ij,d_theta1); - v3_scale(&temp,&dist_ik,d_theta2); - v3_add(&force,&force,&temp); + /* 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)); - /* part 1 of db_ij */ - v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac/d2_h_cos2); + /* scale d_costheta contrib with f_c_ik */ + v3_scale(&force,&force,f_c_ik); - /* part 2 of db_ij */ - v3_scale(&temp,&dist_ik,df_c_ik*g); + /* df_c_ik contrib to db_ij */ + v3_scale(&temp,&dist_ik,df_c_ik*g); - /* sum up and add to db_ij */ - v3_add(&temp,&temp,&force); - v3_add(&(exchange->db_ij),&(exchange->db_ij),&temp); + /* sum up both parts */ + v3_add(&force,&force,&temp); + + /* zeta */ + exchange->zeta+=f_c_ik*g; + } +printf("%.30f\n",exchange->zeta); + + /* add to db_ij */ + v3_add(&(exchange->db_ij),&(exchange->db_ij),&force); return 0; }