From 881d130efd1104bd886507e3d90f262f7a82cf2d Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 25 Apr 2006 15:15:34 +0000 Subject: [PATCH] basic linked list / cell method support implemented --- Makefile | 4 +- moldyn.c | 197 +++++++++++++++++++++++++++++++++++++++++++++++++++---- moldyn.h | 23 ++++++- posic.c | 31 ++++----- run | 1 + 5 files changed, 228 insertions(+), 28 deletions(-) diff --git a/Makefile b/Makefile index 5f25159..b696aae 100644 --- a/Makefile +++ b/Makefile @@ -1,10 +1,12 @@ CC=gcc CFLAGS=-Wall -OBJS=init/init.o visual/visual.o math/math.o random/random.o moldyn.o +OBJS=init/init.o visual/visual.o math/math.o random/random.o list/list.o +OBJS+=moldyn.o all: posic +.PHONY:posic posic: $(OBJS) moldyn.o $(CC) $(CFLAGS) -lm -o $@ $(OBJS) $(LIBS) posic.c diff --git a/moldyn.c b/moldyn.c index 7f9745a..dbb3e6c 100644 --- a/moldyn.c +++ b/moldyn.c @@ -21,6 +21,7 @@ #include "init/init.h" #include "random/random.h" #include "visual/visual.h" +#include "list/list.h" int moldyn_usage(char **argv) { @@ -338,6 +339,144 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) { return tau; } +/* + * 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; +} + +int link_cell_init(t_moldyn *moldyn) { + + t_linkcell *lc; + + lc=&(moldyn->lc); + + /* partitioning the md cell */ + lc->nx=moldyn->dim.x/moldyn->cutoff; + lc->x=moldyn->dim.x/lc->nx; + lc->ny=moldyn->dim.y/moldyn->cutoff; + lc->y=moldyn->dim.y/lc->ny; + lc->nz=moldyn->dim.z/moldyn->cutoff; + lc->z=moldyn->dim.z/lc->nz; + + lc->subcell=malloc(lc->nx*lc->ny*lc->nz*sizeof(t_list)); + + link_cell_update(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; + + atom=moldyn->atom; + nx=moldyn->lc.nx; ny=moldyn->lc.ny; nz=moldyn->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; + } + list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]), + &(atom[count])); + } + + return 0; +} + +int verlet_list_shutdown(t_moldyn *moldyn) { + + int i; + + for(i=0;icount;i++) + list_shutdown(&(moldyn->atom[i].verlet)); + + return 0; +} + +int link_cell_shutdown(t_moldyn *moldyn) { + + int i; + t_linkcell *lc; + + lc=&(moldyn->lc); + + for(i=0;inx*lc->ny*lc->nz;i++) + list_shutdown(&(moldyn->lc.subcell[i])); + + return 0; +} /* * @@ -352,25 +491,40 @@ 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]; + /* logging & visualization */ e=moldyn->ewrite; m=moldyn->mwrite; s=moldyn->swrite; 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; } + /* create the verlet list */ + verlet_list_update(moldyn); + /* calculate initial forces */ moldyn->force(moldyn); for(i=0;itime_steps;i++) { + /* check for velet list update */ + printf("."); + if(moldyn->dr_max1+moldyn->dr_max2>rlc) { + printf("\n"); + verlet_list_update(moldyn); + } + /* integration step */ moldyn->integrate(moldyn); @@ -423,7 +577,7 @@ int moldyn_integrate(t_moldyn *moldyn) { int velocity_verlet(t_moldyn *moldyn) { int i,count; - double tau,tau_square; + double tau,tau_square,dr; t_3dvec delta; t_atom *atom; @@ -437,10 +591,20 @@ int velocity_verlet(t_moldyn *moldyn) { /* new positions */ v3_scale(&delta,&(atom[i].v),tau); v3_add(&(atom[i].r),&(atom[i].r),&delta); + v3_add(&(atom[i].dr),&(atom[i].dr),&delta); v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass); v3_add(&(atom[i].r),&(atom[i].r),&delta); + v3_add(&(atom[i].dr),&(atom[i].dr),&delta); v3_per_bound(&(atom[i].r),&(moldyn->dim)); + /* set maximum dr (possible list update) */ + 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); @@ -536,8 +700,8 @@ int force_harmonic_oscillator(t_moldyn *moldyn) { double potential_lennard_jones(t_moldyn *moldyn) { t_lj_params *params; - t_atom *atom; - int i,j; + t_atom *atom,*btom; + int i; int count; t_3dvec distance; double d,help; @@ -553,13 +717,19 @@ double potential_lennard_jones(t_moldyn *moldyn) { u=0.0; for(i=0;idata; + 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*(sig6*help-sig12*d); + u+=eps*(sig12*d-sig6*help); + if(list_next(&(atom[i].verlet))==L_NO_NEXT_ELEMENT) + break; } } @@ -569,8 +739,8 @@ double potential_lennard_jones(t_moldyn *moldyn) { int force_lennard_jones(t_moldyn *moldyn) { t_lj_params *params; - int i,j,count; - t_atom *atom; + int i,count; + t_atom *atom,*btom; t_3dvec distance; t_3dvec force; double d,h1,h2; @@ -586,8 +756,11 @@ int force_lennard_jones(t_moldyn *moldyn) { 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) { @@ -604,8 +777,10 @@ int force_lennard_jones(t_moldyn *moldyn) { 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); + //v3_sub(&(atom[j].f),&(atom[j].f),&force); } + if(list_next(&(atom[i].verlet))==L_NO_NEXT_ELEMENT) + break; } } diff --git a/moldyn.h b/moldyn.h index 266436b..53a22fd 100644 --- a/moldyn.h +++ b/moldyn.h @@ -10,6 +10,7 @@ #include "math/math.h" #include "random/random.h" +#include "list/list.h" /* datatypes */ @@ -17,11 +18,18 @@ typedef struct s_atom { t_3dvec r; /* positions */ t_3dvec v; /* velocities */ t_3dvec f; /* forces */ + t_3dvec dr; /* delta r for verlet list updates */ int element; /* number of element in pse */ double mass; /* atom mass */ - //t_list vicinity /* verlet neighbour list */ + t_list verlet; /* verlet neighbour list */ } t_atom; +typedef struct s_linkcell { + int nx,ny,nz; + double x,y,z; + t_list *subcell; +} t_linkcell; + #include "visual/visual.h" typedef struct s_moldyn { @@ -33,8 +41,14 @@ typedef struct s_moldyn { double (*potential)(struct s_moldyn *moldyn); int (*force)(struct s_moldyn *moldyn); void *pot_params; + /* cut off radius, verlet list & co */ double cutoff; double cutoff_square; + double r_verlet; + double dr_max1; + double dr_max2; + /* linked list / cell method */ + t_linkcell lc; /* temperature */ double t; /* integration of newtons equations */ @@ -129,6 +143,13 @@ t_3dvec get_total_p(t_atom *atom,int count); double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t); +int verlet_list_init(t_moldyn *moldyn); +int link_cell_init(t_moldyn *moldyn); +int verlet_list_update(t_moldyn *moldyn); +int link_cell_update(t_moldyn *moldyn); +int verlet_list_shutdown(t_moldyn *moldyn); +int link_cell_shutdown(t_moldyn *moldyn); + int moldyn_integrate(t_moldyn *moldyn); int velocity_verlet(t_moldyn *moldyn); diff --git a/posic.c b/posic.c index 8bbfe0e..7cde6e2 100644 --- a/posic.c +++ b/posic.c @@ -22,7 +22,7 @@ int main(int argc,char **argv) { t_random random; int a,b,c; - double e,u; + double e; double help; t_3dvec p; int count; @@ -79,8 +79,8 @@ int main(int argc,char **argv) { md.force=force_lennard_jones; //md.potential=potential_harmonic_oscillator; //md.force=force_harmonic_oscillator; - md.cutoff=R_CUTOFF; - md.cutoff_square=(R_CUTOFF*R_CUTOFF); + md.cutoff=R_CUTOFF*LC_SI; + md.cutoff_square=md.cutoff*md.cutoff; md.pot_params=&lj; //md.pot_params=&ho; md.integrate=velocity_verlet; @@ -88,6 +88,14 @@ int main(int argc,char **argv) { //md.tau=TAU; md.status=0; md.visual=&vis; + /* dimensions of the simulation cell */ + md.dim.x=a*LC_SI; + md.dim.y=b*LC_SI; + md.dim.z=c*LC_SI; + + /* verlet list init */ + // later integrated in moldyn_init function! + verlet_list_init(&md); printf("setting thermal fluctuations (T=%f K)\n",md.t); thermal_init(&md,&random,count); @@ -97,13 +105,14 @@ int main(int argc,char **argv) { e=get_e_kin(si,count); printf("kinetic energy: %.40f [J]\n",e); - printf("3/2 N k T = %.40f [J]\n",1.5*count*K_BOLTZMANN*md.t); + printf("3/2 N k T = %.40f [J] (T=%f [K])\n", + 1.5*count*K_BOLTZMANN*md.t,md.t); /* check total momentum */ p=get_total_p(si,count); printf("total momentum: %.30f [Ns]\n",v3_norm(&p)); - /* check potential energy */ + /* potential paramters */ lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI; help=lj.sigma6*lj.sigma6; lj.sigma6*=help; @@ -113,16 +122,6 @@ int main(int argc,char **argv) { ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; ho.spring_constant=1.0; - u=get_e_pot(&md); - - printf("potential energy: %.40f [J]\n",u); - printf("total energy (1): %.40f [J]\n",e+u); - printf("total energy (2): %.40f [J]\n",get_total_energy(&md)); - - md.dim.x=a*LC_SI; - md.dim.y=b*LC_SI; - md.dim.z=c*LC_SI; - printf("estimated accurate time step: %.30f [s]\n", estimate_time_step(&md,3.0,md.t)); @@ -140,6 +139,8 @@ int main(int argc,char **argv) { /* close */ + verlet_list_shutdown(&md); + rand_close(&random); moldyn_shutdown(&md); diff --git a/run b/run index 4c753fe..584a446 100755 --- a/run +++ b/run @@ -1,3 +1,4 @@ +mkdir -p saves video ./clean ./posic $@ if [ "$?" == "0" ]; then -- 2.39.2