X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=moldyn.c;h=36086581e666b082c6a3d554706609800beafbb2;hb=ade81aa2afb15f22e98ed9595ff303d4fedfe122;hp=716a087fa90f1b89b073411af22e2d55995aba8a;hpb=15dd7ee084a5db24a3e1a15428cd47e8d61fda5b;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 716a087..3608658 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 */ @@ -736,7 +741,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 +753,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 +793,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 +831,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; - /* - * according to mr. nordlund, we dont need to take the - * sum over all atoms now, as 'this is centered' around - * atom i ... - * i am not quite sure though! there is a not vanishing - * part even if f_c_ik is zero ... - * this analytical potentials suck! - * switching from mc to md to dft soon! - */ - - // 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); @@ -960,10 +919,13 @@ 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) { +printf("pre(%d): %.15f %.15f %.15f\n",i,itom[i].f.x,itom[i].r.x,itom[i].v.x); moldyn->func2b_post(moldyn, &(itom[i]), jtom,bc_ij); +printf("post(%d): %.15f %.15f %.15f\n",i,itom[i].f.x,itom[i].r.x,itom[i].v.x); + } } } @@ -1156,6 +1118,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange=&(params->exchange); exchange->run3bp=0; + exchange->run2bp_post=0; /* * we need: f_c, df_c, f_r, df_r @@ -1193,14 +1156,13 @@ 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; f_r=A*exp(-lambda*d_ij); df_r=-lambda*f_r/d_ij; - /* f_a, df_a calc + save for 3bp use */ + /* f_a, df_a calc + save for later use */ exchange->f_a=-B*exp(-mu*d_ij); exchange->df_a=-mu*exchange->f_a/d_ij; @@ -1228,8 +1190,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->f_c=f_c; exchange->df_c=df_c; - /* enable the run of 3bp function */ + /* enable the run of 3bp function and 2bp post processing */ exchange->run3bp=1; + exchange->run2bp_post=1; /* reset 3bp sums */ exchange->sum1_3bp=0.0; @@ -1259,6 +1222,10 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { params=moldyn->pot2b_params; exchange=&(params->exchange); + /* we do not run if f_c_ij was detected to be 0! */ + if(!(exchange->run2bp_post)) + return 0; + db_ij=&(exchange->db_ij); f_c=exchange->f_c; df_c=exchange->df_c; @@ -1266,6 +1233,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { df_a=exchange->df_a; betan=exchange->betan; n=*(exchange->n); + chi=exchange->chi; dist_ij=&(exchange->dist_ij); db_ij_scale1=(1+betan*exchange->sum1_3bp); @@ -1273,20 +1241,26 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { help=pow(db_ij_scale1,-1.0/(2*n)-1); b_ij=chi*db_ij_scale1*help; db_ij_scale1=-chi/(2*n)*help; +printf("debug: %.20f %.20f %.20f\n",db_ij->x,exchange->sum1_3bp,exchange->sum2_3bp); + /* 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); + /* df_c part */ 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 */ + + /* add force of 3bp calculation (all three parts) */ v3_add(&(ai->f),&temp,&force); return 0; @@ -1392,11 +1366,26 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { d2=exchange->d2; c2d2=exchange->c2d2; + /* cosine of theta by scalaproduct, * + * derivation of theta by law of cosines! */ 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); + cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); +printf("cos theta: %.25f\n",cos_theta); + + /* hack - cos theta machine accuracy problems! */ + if(cos_theta>1.0||cos_theta<-1.0) { + 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; + printf("THETA CORRECTION\n"); + } + 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); @@ -1404,6 +1393,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { d_theta2=2*denom-numer*2*d_ij/d_ik; d_theta1*=d_theta; d_theta2*=d_theta; +printf("FOO %.15f %.15f\n",sin_theta,cos_theta); h_cos=(h-cos_theta); d2_h_cos2=d2+(h_cos*h_cos); @@ -1415,15 +1405,12 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { 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; @@ -1434,6 +1421,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { v3_scale(&temp,&dist_ik,d_theta2); v3_add(&force,&force,&temp); +printf("DA:%.20f %.20f %.20f\n",d_theta1,force.x,temp.x); /* part 1 of db_ij */ v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac/d2_h_cos2);