X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=4b580890f3eaab417494fa497130f53570130ce4;hb=e6f456c0fa807b86e1b25996e70efcdcfe390ea5;hp=17720675c4f5c0c76dea880ac40eb5787f772eb0;hpb=877139ec35cf357c96fdb314f690fd6075554c47;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 1772067..4b58089 100644 --- a/moldyn.c +++ b/moldyn.c @@ -36,6 +36,17 @@ int moldyn_usage(char **argv) { printf("-T [K] (%f)\n",MOLDYN_TEMP); printf("-t [s] (%.15f)\n",MOLDYN_TAU); printf("-R (%d)\n",MOLDYN_RUNS); + printf(" -- integration algo --\n"); + printf(" -I (%d)\n",MOLDYN_INTEGRATE_DEFAULT); + printf(" 0: velocity verlet\n"); + printf(" -- potential --\n"); + printf(" -P \n"); + printf(" 0: harmonic oscillator\n"); + printf(" param1: spring constant\n"); + printf(" param2: equilibrium distance\n"); + printf(" 1: lennard jones\n"); + printf(" param1: epsilon\n"); + printf(" param2: sigma\n"); printf("\n"); return 0; @@ -44,6 +55,9 @@ 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; + double s,e; memset(moldyn,0,sizeof(t_moldyn)); @@ -51,6 +65,8 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) { moldyn->t=MOLDYN_TEMP; moldyn->tau=MOLDYN_TAU; moldyn->time_steps=MOLDYN_RUNS; + moldyn->integrate=velocity_verlet; + moldyn->potential_force_function=lennard_jones; /* parse argv */ for(i=1;itime_steps=atoi(argv[++i]); break; + case 'I': + /* integration algorithm */ + switch(atoi(argv[++i])) { + case MOLDYN_INTEGRATE_VERLET: + moldyn->integrate=velocity_verlet; + break; + default: + printf("unknown integration algo %s\n",argv[i]); + moldyn_usage(argv); + return -1; + } + + case 'P': + /* potential + params */ + switch(atoi(argv[++i])) { + case MOLDYN_POTENTIAL_HO: + hop.spring_constant=atof(argv[++i]); + hop.equilibrium_distance=atof(argv[++i]); + moldyn->pot_params=malloc(sizeof(t_ho_params)); + memcpy(moldyn->pot_params,&hop,sizeof(t_ho_params)); + moldyn->potential_force_function=harmonic_oscillator; + break; + case MOLDYN_POTENTIAL_LJ: + e=atof(argv[++i]); + s=atof(argv[++i]); + ljp.epsilon4=4*e; + ljp.sigma6=s*s*s*s*s*s; + ljp.sigma12=ljp.sigma6*ljp.sigma6; + moldyn->pot_params=malloc(sizeof(t_lj_params)); + memcpy(moldyn->pot_params,&ljp,sizeof(t_lj_params)); + moldyn->potential_force_function=lennard_jones; + break; + default: + printf("unknown potential %s\n",argv[i]); + moldyn_usage(argv); + return -1; + } + default: printf("unknown option %s\n",argv[i]); moldyn_usage(argv); @@ -216,7 +270,7 @@ int destroy_lattice(t_atom *atom) { return 0; } -int thermal_init(t_moldyn *moldyn,t_random *random,int count) { +int thermal_init(t_moldyn *moldyn,t_random *random) { /* * - gaussian distribution of velocities @@ -233,7 +287,7 @@ int thermal_init(t_moldyn *moldyn,t_random *random,int count) { /* gaussian distribution of velocities */ v3_zero(&p_total); - for(i=0;icount;i++) { sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass); /* x direction */ v=sigma*rand_get_gauss(random); @@ -250,19 +304,19 @@ int thermal_init(t_moldyn *moldyn,t_random *random,int count) { } /* zero total momentum */ - v3_scale(&p_total,&p_total,1.0/count); - for(i=0;icount); + for(i=0;icount;i++) { v3_scale(&delta,&p_total,1.0/atom[i].mass); v3_sub(&(atom[i].v),&(atom[i].v),&delta); } /* velocity scaling */ - scale_velocity(moldyn,count); + scale_velocity(moldyn); return 0; } -int scale_velocity(t_moldyn *moldyn,int count) { +int scale_velocity(t_moldyn *moldyn) { int i; double e,c; @@ -274,10 +328,10 @@ int scale_velocity(t_moldyn *moldyn,int count) { * - velocity scaling (E = 3/2 N k T), E: kinetic energy */ e=0.0; - for(i=0;icount;i++) e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); - c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*moldyn->t)); - for(i=0;icount*K_BOLTZMANN*moldyn->t)); + for(i=0;icount;i++) v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c)); return 0; @@ -299,7 +353,7 @@ double get_e_kin(t_atom *atom,int count) { double get_e_pot(t_moldyn *moldyn) { - return(moldyn->potential(moldyn)); + return moldyn->energy; } double get_total_energy(t_moldyn *moldyn) { @@ -343,31 +397,7 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) { * numerical tricks */ -/* verlet list */ - -int verlet_list_init(t_moldyn *moldyn) { - - int i,fd; - - fd=open("/dev/null",O_WRONLY); - - for(i=0;icount;i++) - list_init(&(moldyn->atom[i].verlet),fd); - - moldyn->r_verlet=1.1*moldyn->cutoff; - /* +moldyn->tau*\ - sqrt(3.0*K_BOLTZMANN*moldyn->t/moldyn->atom[0].mass); */ - - printf("debug: r verlet = %.15f\n",moldyn->r_verlet); - printf(" r cutoff = %.15f\n",moldyn->cutoff); - printf(" dim = %.15f\n",moldyn->dim.x); - - /* make sure to update the list in the beginning */ - moldyn->dr_max1=moldyn->r_verlet; - moldyn->dr_max2=moldyn->r_verlet; - - return 0; -} +/* linked list / cell method */ int link_cell_init(t_moldyn *moldyn) { @@ -390,64 +420,27 @@ int link_cell_init(t_moldyn *moldyn) { return 0; } -int verlet_list_update(t_moldyn *moldyn) { - - int i,j; - t_3dvec d; - t_atom *atom; - - atom=moldyn->atom; - - puts("debug: list update start"); - - for(i=0;icount;i++) { - list_destroy(&(atom[i].verlet)); - for(j=0;jcount;j++) { - if(i!=j) { - v3_sub(&d,&(atom[i].r),&(atom[j].r)); - v3_per_bound(&d,&(moldyn->dim)); - if(v3_norm(&d)<=moldyn->r_verlet) - list_add_immediate_ptr(&(atom[i].verlet),&(atom[j])); - } - } - } - - moldyn->dr_max1=0.0; - moldyn->dr_max2=0.0; - - puts("debug: list update end"); - - return 0; -} - int link_cell_update(t_moldyn *moldyn) { int count,i,j,k; int nx,ny,nz; t_atom *atom; + t_linkcell *lc; atom=moldyn->atom; - nx=moldyn->lc.nx; ny=moldyn->lc.ny; nz=moldyn->lc.nz; + lc=&(moldyn->lc); + + nx=lc->nx; + ny=lc->ny; + nz=lc->nz; for(i=0;ilc.subcell[i])); for(count=0;countcount;count++) { - for(i=0;i=i*moldyn->dim.x) && \ - (atom[count].r.x<(i+1)*moldyn->dim.x)) - break; - } - for(j=0;j=j*moldyn->dim.y) && \ - (atom[count].r.y<(j+1)*moldyn->dim.y)) - break; - } - for(k=0;k=k*moldyn->dim.z) && \ - (atom[count].r.z<(k+1)*moldyn->dim.z)) - break; - } + i=atom[count].r.x/lc->x; + j=atom[count].r.y/lc->y; + k=atom[count].r.z/lc->z; list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]), &(atom[count])); } @@ -455,14 +448,58 @@ int link_cell_update(t_moldyn *moldyn) { return 0; } -int verlet_list_shutdown(t_moldyn *moldyn) { +int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { - int i; + t_linkcell *lc; + int a; + int count1,count2; + int ci,cj,ck; + int nx,ny,nz; + int x,y,z; + unsigned char bx,by,bz; - for(i=0;icount;i++) - list_shutdown(&(moldyn->atom[i].verlet)); + lc=&(moldyn->lc); + nx=lc->nx; + ny=lc->ny; + nx=lc->nz; + count1=1; + count2=27; + a=nx*ny; + + cell[0]=lc->subcell[i+j*nx+k*a]; + for(ci=-1;ci<=1;ci++) { + bx=0; + x=i+ci; + if((x<0)||(x>=nx)) { + x=(x+nx)%nx; + bx=1; + } + for(cj=-1;cj<=1;cj++) { + by=0; + y=j+cj; + if((y<0)||(y>=ny)) { + y=(y+ny)%ny; + by=1; + } + for(ck=-1;ck<=1;ck++) { + bz=0; + z=k+ck; + if((z<0)||(z>=nz)) { + z=(z+nz)%nz; + bz=1; + } + if(!(x|y|z)) continue; + if(bx|by|bz) { + cell[--count2]=lc->subcell[x+y*nx+z*a]; + } + else { + cell[count1++]=lc->subcell[x+y*nx+z*a]; + } + } + } + } - return 0; + return count2; } int link_cell_shutdown(t_moldyn *moldyn) { @@ -491,7 +528,6 @@ int moldyn_integrate(t_moldyn *moldyn) { int i; unsigned int e,m,s,d,v; t_3dvec p; - double rlc; int fd; char fb[128]; @@ -503,20 +539,20 @@ int moldyn_integrate(t_moldyn *moldyn) { d=moldyn->dwrite; v=moldyn->vwrite; - /* verlet list */ - rlc=moldyn->r_verlet-moldyn->cutoff; - if(!(moldyn->lvstat&MOLDYN_LVSTAT_INITIALIZED)) { printf("[moldyn] warning, lv system not initialized\n"); return -1; } + /* sqaure of some variables */ + moldyn->tau_square=moldyn->tau*moldyn->tau; + moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; + /* create the neighbour list */ - //verlet_list_update(moldyn); link_cell_update(moldyn); /* calculate initial forces */ - moldyn->force(moldyn); + moldyn->potential_force_function(moldyn); for(i=0;itime_steps;i++) { /* show runs */ @@ -524,15 +560,11 @@ int moldyn_integrate(t_moldyn *moldyn) { /* neighbour list update */ link_cell_update(moldyn); - //if(moldyn->dr_max1+moldyn->dr_max2>rlc) { - // printf("\n"); - // verlet_list_update(moldyn); - //} /* integration step */ moldyn->integrate(moldyn); - /* check for log & visualiziation */ + /* check for log & visualization */ if(e) { if(!(i%e)) dprintf(moldyn->efd, @@ -581,41 +613,30 @@ int moldyn_integrate(t_moldyn *moldyn) { int velocity_verlet(t_moldyn *moldyn) { int i,count; - double tau,tau_square,dr; + double tau,tau_square; t_3dvec delta; t_atom *atom; atom=moldyn->atom; count=moldyn->count; tau=moldyn->tau; - - tau_square=tau*tau; + tau_square=moldyn->tau_square; for(i=0;idim)); - /* set maximum dr (possible list update) [obsolete] */ - //dr=v3_norm(&(atom[i].dr)); - //if(dr>moldyn->dr_max1) { - // moldyn->dr_max2=moldyn->dr_max1; - // moldyn->dr_max1=dr; - //} - //else if(dr>moldyn->dr_max2) moldyn->dr_max2=dr; - /* velocities */ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass); v3_add(&(atom[i].v),&(atom[i].v),&delta); } /* forces depending on chosen potential */ - moldyn->force(moldyn); + moldyn->potential_force_function(moldyn); for(i=0;ipot_params; atom=moldyn->atom; + lc=&(moldyn->lc); sc=params->spring_constant; equi_dist=params->equilibrium_distance; count=moldyn->count; u=0.0; for(i=0;ix; + nj=atom[i].r.y/lc->y; + nk=atom[i].r.z/lc->z; + c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour); + + /* processing cell of atom i */ + this=&(neighbour[0]); + list_reset(this); /* list has 1 element at minimum */ + do { + btom=this->current->data; + if(btom==&(atom[i])) + continue; + v3_sub(&distance,&(atom[i].r),&(btom->r)); d=v3_norm(&distance); u+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); - } - } - - return u; -} - -int force_harmonic_oscillator(t_moldyn *moldyn) { - - t_ho_params *params; - int i,j,count; - t_atom *atom; - t_3dvec distance; - t_3dvec force; - double d; - double sc,equi_dist; + v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d))); + v3_add(&(atom[i].f),&(atom[i].f),&force); + } while(list_next(this)!=L_NO_NEXT_ELEMENT); + + /* neighbours not doing boundary condition overflow */ + for(j=1;jstart!=NULL) { + + do { + btom=this->current->data; + v3_sub(&distance,&(atom[i].r),&(btom->r)); + d=v3_norm(&distance); + if(d<=moldyn->cutoff) { + u+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); + v3_scale(&force,&distance, + -sc*(1.0-(equi_dist/d))); + v3_add(&(atom[i].f),&(atom[i].f), + &force); + } + } while(list_next(this)!=L_NO_NEXT_ELEMENT); - atom=moldyn->atom; - count=moldyn->count; - params=moldyn->pot_params; - sc=params->spring_constant; - equi_dist=params->equilibrium_distance; + } + } - for(i=0;istart!=NULL) { + + do { + btom=this->current->data; + v3_sub(&distance,&(atom[i].r),&(btom->r)); + v3_per_bound(&distance,&(moldyn->dim)); + d=v3_norm(&distance); + if(d<=moldyn->cutoff) { + u+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); + v3_scale(&force,&distance, + -sc*(1.0-(equi_dist/d))); + v3_add(&(atom[i].f),&(atom[i].f), + &force); + } + } while(list_next(this)!=L_NO_NEXT_ELEMENT); - for(i=0;idim)); - d=v3_norm(&distance); - if(d<=moldyn->cutoff) { - v3_scale(&force,&distance, - -sc*(1.0-(equi_dist/d))); - v3_add(&(atom[i].f),&(atom[i].f),&force); - v3_sub(&(atom[j].f),&(atom[j].f),&force); } } } + moldyn->energy=u; + return 0; } - /* lennard jones potential & force for one sort of atoms */ -double potential_lennard_jones(t_moldyn *moldyn) { +int lennard_jones(t_moldyn *moldyn) { t_lj_params *params; t_atom *atom,*btom; t_linkcell *lc; - int i; + t_list *this,neighbour[27]; + int i,j,c; int count; - t_3dvec distance; - double d,help; - double u; + t_3dvec force,distance; + double d,h1,h2,u; double eps,sig6,sig12; - int i,j,k; int ni,nj,nk; params=moldyn->pot_params; @@ -725,80 +774,101 @@ double potential_lennard_jones(t_moldyn *moldyn) { u=0.0; for(i=0;ix; nj=atom[i].r.y/lc->y; nk=atom[i].r.z/lc->z; - - while(1) { - btom=atom[i].verlet.current->data; + c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour); + + /* processing cell of atom i */ + this=&(neighbour[0]); + list_reset(this); /* list has 1 element at minimum */ + do { + btom=this->current->data; + if(btom==&(atom[i])) + continue; v3_sub(&distance,&(atom[i].r),&(btom->r)); - v3_per_bound(&distance,&(moldyn->dim)); - d=1.0/v3_absolute_square(&distance); /* 1/r^2 */ - help=d*d; /* 1/r^4 */ - help*=d; /* 1/r^6 */ - d=help*help; /* 1/r^12 */ - u+=eps*(sig12*d-sig6*help); - if(list_next(&(atom[i].verlet))==L_NO_NEXT_ELEMENT) - break; + d=1.0/v3_absolute_square(&distance); /* 1/r^2 */ + h1=d*d; /* 1/r^4 */ + h2*=d; /* 1/r^6 */ + h1=h2*h2; /* 1/r^12 */ + u+=eps*(sig12*h1-sig6*h2); + h2*=d; /* 1/r^8 */ + h1*=d; /* 1/r^14 */ + h2*=6*sig6; + h1*=12*sig12; + d=-h1+h2; + d*=eps; + v3_scale(&force,&distance,d); + v3_add(&(atom[i].f),&(atom[i].f),&force); + } while(list_next(this)!=L_NO_NEXT_ELEMENT); + + /* neighbours not doing boundary condition overflow */ + for(j=1;jstart!=NULL) { + + do { + btom=this->current->data; + v3_sub(&distance,&(atom[i].r),&(btom->r)); + d=v3_absolute_square(&distance); /* r^2 */ + if(d<=moldyn->cutoff_square) { + d=1.0/d; /* 1/r^2 */ + h1=d*d; /* 1/r^4 */ + h2*=d; /* 1/r^6 */ + h1=h2*h2; /* 1/r^12 */ + u+=eps*(sig12*h1-sig6*h2); + h2*=d; /* 1/r^8 */ + h1*=d; /* 1/r^14 */ + h2*=6*sig6; + h1*=12*sig12; + d=-h1+h2; + d*=eps; + v3_scale(&force,&distance,d); + v3_add(&(atom[i].f),&(atom[i].f), + &force); + } + } while(list_next(this)!=L_NO_NEXT_ELEMENT); + + } } - } - - return u; -} - -int force_lennard_jones(t_moldyn *moldyn) { - - t_lj_params *params; - int i,count; - t_atom *atom,*btom; - t_3dvec distance; - t_3dvec force; - double d,h1,h2; - double eps,sig6,sig12; - atom=moldyn->atom; - count=moldyn->count; - params=moldyn->pot_params; - eps=params->epsilon4; - sig6=6*params->sigma6; - sig12=12*params->sigma12; - - for(i=0;istart!=NULL) { + + do { + btom=this->current->data; + v3_sub(&distance,&(atom[i].r),&(btom->r)); + v3_per_bound(&distance,&(moldyn->dim)); + d=v3_absolute_square(&distance); /* r^2 */ + if(d<=moldyn->cutoff_square) { + d=1.0/d; /* 1/r^2 */ + h1=d*d; /* 1/r^4 */ + h2*=d; /* 1/r^6 */ + h1=h2*h2; /* 1/r^12 */ + u+=eps*(sig12*h1-sig6*h2); + h2*=d; /* 1/r^8 */ + h1*=d; /* 1/r^14 */ + h2*=6*sig6; + h1*=12*sig12; + d=-h1+h2; + d*=eps; + v3_scale(&force,&distance,d); + v3_add(&(atom[i].f),&(atom[i].f), + &force); + } + } while(list_next(this)!=L_NO_NEXT_ELEMENT); - for(i=0;idata; - v3_sub(&distance,&(atom[i].r),&(btom->r)); - v3_per_bound(&distance,&(moldyn->dim)); - d=v3_absolute_square(&distance); - if(d<=moldyn->cutoff_square) { - h1=1.0/d; /* 1/r^2 */ - d=h1*h1; /* 1/r^4 */ - h2=d*d; /* 1/r^8 */ - h1*=d; /* 1/r^6 */ - h1*=h2; /* 1/r^14 */ - h1*=sig12; - h2*=sig6; - /* actually there would be a '-', * - * but f=-d/dr potential */ - d=h1+h2; - d*=eps; - v3_scale(&force,&distance,d); - v3_add(&(atom[i].f),&(atom[i].f),&force); - //v3_sub(&(atom[j].f),&(atom[j].f),&force); } - if(list_next(&(atom[i].verlet))==L_NO_NEXT_ELEMENT) - break; } } + moldyn->energy=u; + return 0; }