X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=d74b391a72362cd33d5ed6e6e698de9e2bb359de;hb=0656efc2936da55ecf0b818fe33bb2acebc689c6;hp=1b52a5628b3e3fe4df32ec4a76483795a1491c15;hpb=d445878769b7fcc15959d2d10f2eeced803c5487;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 1b52a56..d74b391 100644 --- a/moldyn.c +++ b/moldyn.c @@ -58,6 +58,7 @@ 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)); @@ -81,10 +82,6 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) { moldyn->mwrite=atoi(argv[++i]); strncpy(moldyn->mfb,argv[++i],64); break; - case 'D': - moldyn->dwrite=atoi(argv[++i]); - strncpy(moldyn->dfb,argv[++i],64); - break; case 'S': moldyn->swrite=atoi(argv[++i]); strncpy(moldyn->sfb,argv[++i],64); @@ -187,16 +184,6 @@ int moldyn_log_init(t_moldyn *moldyn) { if(moldyn->swrite) moldyn->lvstat|=MOLDYN_LVSTAT_SAVE; - if(moldyn->dwrite) { - moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC); - if(moldyn->dfd<0) { - perror("[moldyn] dfd open"); - return moldyn->dfd; - } - write(moldyn->dfd,moldyn,sizeof(t_moldyn)); - moldyn->lvstat|=MOLDYN_LVSTAT_DUMP; - } - if((moldyn->vwrite)&&(vis)) { moldyn->visual=vis; visual_init(vis,moldyn->vfb); @@ -481,7 +468,7 @@ int link_cell_update(t_moldyn *moldyn) { for(i=0;icells;i++) list_destroy(&(moldyn->lc.subcell[i])); - for(count=0;countcount;count++) { + for(count=0;countcount;count++) { i=(atom[count].r.x+(moldyn->dim.x/2))/lc->x; j=(atom[count].r.y+(moldyn->dim.y/2))/lc->y; k=(atom[count].r.z+(moldyn->dim.z/2))/lc->z; @@ -544,6 +531,9 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { } } + lc->dnlc=count2; + lc->countn=27; + return count2; } @@ -632,14 +622,9 @@ int moldyn_integrate(t_moldyn *moldyn) { write(fd,moldyn->atom, moldyn->count*sizeof(t_atom)); } + close(fd); } } - if(d) { - if(!(i%d)) - write(moldyn->dfd,moldyn->atom, - moldyn->count*sizeof(t_atom)); - - } if(v) { if(!(i%v)) { visual_atoms(moldyn->visual,i*moldyn->tau, @@ -687,7 +672,8 @@ printf("done\n"); /* forces depending on chosen potential */ printf("calc potential/force ...\n"); - moldyn->potential_force_function(moldyn); + potential_force_calc(moldyn); + //moldyn->potential_force_function(moldyn); printf("done\n"); for(i=0;icount; + atom=moldyn->atom; + lc=&(moldyn->lc); + + /* reset energy */ + moldyn->energy=0.0; + + for(i=0;ipf_func1b(moldyn,&(atom[i])); + + /* 2 body pair potential/force */ + if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) { + + link_cell_neighbour_index(moldyn, + (atom[i].r.x+moldyn->dim.x/2)/lc->x, + (atom[i].r.y+moldyn->dim.y/2)/lc->y, + (atom[i].r.z+moldyn->dim.z/2)/lc->z, + neighbour); + + countn=lc->countn; + dnlc=lc->dnlc; + + for(j=0;jstart==NULL) + continue; + + bc=(jcurrent->data; + + if(btom==&(atom[i])) + continue; + + if((btom->attr&ATOM_ATTR_2BP)& + (atom[i].attr&ATOM_ATTR_2BP)) + moldyn->pf_func2b(moldyn, + &(atom[i]), + btom, + bc); + + /* 3 body potential/force */ + + if(!(atom[i].attr&ATOM_ATTR_3BP)|| + !(btom->attr&ATOM_ATTR_3BP)) + continue; + + link_cell_neighbour_index(moldyn, + (btom->r.x+moldyn->dim.x/2)/lc->x, + (btom->r.y+moldyn->dim.y/2)/lc->y, + (btom->r.z+moldyn->dim.z/2)/lc->z, + neighbourk); + + for(k=0;kcountn;k++) { + + thisk=&(neighbourk[k]); + list_reset(thisk); + + if(thisk->start==NULL) + continue; + + bck=(kdnlc)?0:1; + + do { + + ktom=thisk->current->data; + + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; + + if(ktom==btom) + continue; + + if(ktom==&(atom[i])) + continue; + + moldyn->pf_func3b(moldyn,&(atom[i]),btom,ktom,bck); + + } while(list_next(thisk)!=\ + L_NO_NEXT_ELEMENT); + + } while(list_next(this)!=L_NO_NEXT_ELEMENT); + } + } + } + + return 0; +} + +/* + * example potentials + */ + /* harmonic oscillator potential and force */ -int harmonic_oscillator(t_moldyn *moldyn) { +int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,unsigned char bc)) { t_ho_params *params; t_atom *atom,*btom; @@ -995,20 +1097,29 @@ int tersoff(t_moldyn *moldyn) { /* 2 body stuff */ + /* we need: f_c, df_c, f_r, df_r */ + v3_sub(&dist_ij,btom,&(atom[i])); d_ij=v3_norm(&dist_ij); if(d_ij<=S) { - S=; - R=; + /* determine the tersoff parameters */ + if(atom[i].element!=btom->element) { + S=sqrt(TERSOFF_S[e1]*TERSOFF_S[e2]); + R=R_m; A=; lambda=; B=; mu=; + chi=; + beta=; + betaN=; if(d_ij<=R) { - f_c=1.0; - df_c=0.0; + df_r=-lambda*A*exp(-lambda*d_ij)/d_ij; + v3_scale(&force,&dist_ij,df_r); + v3_add(&(atom[i].f),&(atom[i].f), + &force); } else { s_r=S-R; @@ -1016,7 +1127,11 @@ int tersoff(t_moldyn *moldyn) { f_c=0.5+0.5*cos(arg1); df_c=-0.5*sin(arg1)*(PI/(s_r*d_ij)); f_r=A*exp(-lambda*d_ij); - f_a=-B*exp(-mu*d_ij); + df_r=-lambda*f_r/d_ij; + scale=df_c*f_r+df_r*f_c; + v3_scale(&force,&dist_ij,scale); + v3_add(&(atom[i].f),&(atom[i].f), + &force); } } else @@ -1032,6 +1147,10 @@ int tersoff(t_moldyn *moldyn) { ck=link_cell_neighbour_index(moldyn,ki,kj,kk, neighbourk); + /* go for zeta - 3 body stuff! */ + zeta=0.0; + d_ij2=d_ij*d_ij; + /* cell of btom */ thisk=&(neighbourk[0]); list_reset(thisk); @@ -1044,12 +1163,47 @@ int tersoff(t_moldyn *moldyn) { /* 3 body stuff (1) */ - theta_ijk=; - sin_theta=; - cos_theta=; - hi_cos=; - hi_cos_square=; + v3_sub(&dist_ik,ktom,&(atom[i])); + d_ik=v3_norm(&dist_ik); + if(d_ik<=Sik) { + + Rik=; + Sik=; + Aik=; + lambda_ik=; + Bik=; + mu_ik=; + omega_ik=; + c_i=; + d_i=; + h_i=; + + if(d_ik<=Rik) { + f_cik=1.0; + df_cik=0.0; + } + else { + sik_rik=Sik-Rik; + arg1ik=PI*(d_ik-Rik)/sik_rik; + f_cik=0.5+0.5*cos(arg1ik); + df_cik=-0.5*sin(arg1ik)* \ + (PI/(sik_rik*d_ik)); + f_rik=Aik*exp(-lambda_ik*d_ik); + f_aik=-Bik*exp(-mu_ik*d_ik); + } + + v3_sub(&distance_jk,ktom,btom); + cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/\ + (2*d_ij*d_ik); + sin_theta=sqrt(1.0/\ + (cos_theta*cos_theta)); + theta=arccos(cos_theta); + + + } + else + continue; /* end 3 body stuff (1) */