X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=fe54330f002f414d345295ff73dd4d1e6b6b557d;hb=1097ea2efb575855edd332c28a9cd8f807cd0716;hp=830b8f044cbab93357598b57f59180b17f4c5d45;hpb=47bc224440051c430cefe4b82a7fb320ba76eba4;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 830b8f0..fe54330 100644 --- a/moldyn.c +++ b/moldyn.c @@ -56,10 +56,6 @@ int moldyn_usage(char **argv) { int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) { int i; - t_ho_params hop; - t_lj_params ljp; - t_tersoff_params tp; - double s,e; memset(moldyn,0,sizeof(t_moldyn)); @@ -68,7 +64,6 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) { moldyn->tau=MOLDYN_TAU; moldyn->time_steps=MOLDYN_RUNS; moldyn->integrate=velocity_verlet; - moldyn->potential_force_function=lennard_jones; /* parse argv */ for(i=1;iatom; if(type==FCC) count*=4; if(type==DIAMOND) count*=8; - *atom=malloc(count*sizeof(t_atom)); - if(*atom==NULL) { + atom=malloc(count*sizeof(t_atom)); + if(atom==NULL) { perror("malloc (atoms)"); return -1; } @@ -254,10 +251,10 @@ int create_lattice(u8 type,int element,double mass,double lc, switch(type) { case FCC: - ret=fcc_init(a,b,c,lc,*atom,&origin); + ret=fcc_init(a,b,c,lc,atom,&origin); break; case DIAMOND: - ret=diamond_init(a,b,c,lc,*atom,&origin); + ret=diamond_init(a,b,c,lc,atom,&origin); break; default: printf("unknown lattice type (%02x)\n",type); @@ -273,17 +270,45 @@ int create_lattice(u8 type,int element,double mass,double lc, } while(count) { - (*atom)[count-1].element=element; - (*atom)[count-1].mass=mass; + atom[count-1].element=element; + atom[count-1].mass=mass; + atom[count-1].attr=attr; + atom[count-1].bnum=bnum; count-=1; } return ret; } -int destroy_lattice(t_atom *atom) { +int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr, + t_3dvec r,t_3dvec v) { - if(atom) free(atom); + t_atom *atom; + void *ptr; + int count; + + atom=moldyn->atom; + count=++(moldyn->count); + + ptr=realloc(atom,count*sizeof(t_atom)); + if(!ptr) { + perror("[moldyn] realloc (add atom)"); + return -1; + } + + atom=ptr; + atom->r=r; + atom->v=v; + atom->element=element; + atom->bnum=bnum; + atom->attr=attr; + + return 0; +} + +int destroy_atoms(t_moldyn *moldyn) { + + if(moldyn->atom) free(moldyn->atom); return 0; } @@ -552,6 +577,40 @@ int link_cell_shutdown(t_moldyn *moldyn) { return 0; } +int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau ) { + + int count; + void *ptr; + t_moldyn_schedule *schedule; + + schedule=moldyn->schedule; + count=++(schedule->content_count); + + ptr=realloc(moldyn->schedule.runs,count*sizeof(int)); + if(!ptr) { + perror("[moldyn] realloc (runs)"); + return -1; + } + moldyn->schedule.runs[count-1]=runs; + + ptr=realloc(schedule->tau,count*sizeof(double)); + if(!ptr) { + perror("[moldyn] realloc (tau)"); + return -1; + } + moldyn->schedule.tau[count-1]=tau; + + return 0; +} + +int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) { + + moldyn->schedule.hook=hook; + moldyn->schedule.hook_params=hook_params; + + return 0; +} + /* * * 'integration of newtons equation' - algorithms @@ -562,7 +621,7 @@ int link_cell_shutdown(t_moldyn *moldyn) { int moldyn_integrate(t_moldyn *moldyn) { - int i; + int i,sched; unsigned int e,m,s,d,v; t_3dvec p; @@ -591,6 +650,15 @@ int moldyn_integrate(t_moldyn *moldyn) { /* calculate initial forces */ moldyn->potential_force_function(moldyn); + for(sched=0;schedschedule.content_count;sched++) { + + /* setting amont of runs and finite time step size */ + moldyn->tau=schedule->tau[sched]; + moldyn->tau_square=moldyn->tau*moldyn->tau; + moldyn->timesteps=schedule->runs[sched]; + + /* integration according to schedule */ + for(i=0;itime_steps;i++) { /* integration step */ @@ -635,6 +703,10 @@ int moldyn_integrate(t_moldyn *moldyn) { } } + /* check for hooks */ + if(schedule->hook) + schedule->hook(moldyn,schedule->hook_params); + return 0; } @@ -925,6 +997,7 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) { exchange->d=&(params->d[num]); exchange->h=&(params->h[num]); + exchange->betan=pow(*(exchange->beta),*(exchange->n)); exchange->c2=params->c[num]*params->c[num]; exchange->d2=params->d[num]*params->d[num]; exchange->c2d2=exchange->c2/exchange->d2; @@ -1035,11 +1108,12 @@ 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 db_ij,temp,force; + t_3dvec temp,force; double R,S,s_r; double d_ij,d_ik,d_jk; double f_c,df_c,b_ij,f_a,df_a; - double B,mu; + double n,c,d,h,neta,betan,betan_1; + double theta,cos_theta,sin_theta; int num; params=moldyn->pot_params; @@ -1058,9 +1132,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; - B=*(params->exchange.B); - mu=*(params->exchange.mu); - f_a=params->exchange.f_a; df_a=params->exchange.df_a; @@ -1111,17 +1182,68 @@ 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_jk); d_jk=v3_norm(&dist_jk); - - // GO ON HERE !!! - - cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/(2*d_ij*d_ik); + beta=*(exchange->beta); + betan=exchange->betan; + n=*(exchange->n); + c=*(exchange->c); + d=*(exchange->d); + h=*(exchange->h); + c2=exchange->c2; + 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; sin_theta=sqrt(1.0-(cos_theta*cos_theta)); theta=arccos(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; 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; + + /* derivation of theta */ + v3_scale(&force,&dist_ij,d1_theta); + v3_scale(&temp,&dist_ik,d_theta2); + v3_add(&force,&force,&temp); + + /* part 1 of derivation of b_ij */ + v3_scale(&force,sin_theta*2*h_cos*f_c_ik*frac1); + + /* part 2 of derivation of b_ij */ + v3_scale(&temp,&dist_ik,df_c_ik*bracket1); + + /* sum up and scale ... */ + 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 ... */