From 177cf8b5cb5a3c59e2330327b628937540f123ac Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 5 May 2006 14:13:58 +0000 Subject: [PATCH] linked cell list updates (not well working by now!) --- moldyn.c | 61 ++++++++++++++++++---- moldyn.h | 11 +++- posic.c | 150 +++++++++++++++++++++++++++++++------------------------ 3 files changed, 146 insertions(+), 76 deletions(-) diff --git a/moldyn.c b/moldyn.c index 4b58089..6e76bc0 100644 --- a/moldyn.c +++ b/moldyn.c @@ -35,6 +35,7 @@ int moldyn_usage(char **argv) { printf("--- physics options ---\n"); printf("-T [K] (%f)\n",MOLDYN_TEMP); printf("-t [s] (%.15f)\n",MOLDYN_TAU); + printf("-C [m] (%.15f)\n",MOLDYN_CUTOFF); printf("-R (%d)\n",MOLDYN_RUNS); printf(" -- integration algo --\n"); printf(" -I (%d)\n",MOLDYN_INTEGRATE_DEFAULT); @@ -98,6 +99,9 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) { case 't': moldyn->tau=atof(argv[++i]); break; + case 'C': + moldyn->cutoff=atof(argv[++i]); + break; case 'R': moldyn->time_steps=atoi(argv[++i]); break; @@ -153,12 +157,12 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) { return 0; } -int moldyn_log_init(t_moldyn *moldyn,void *v) { +int moldyn_log_init(t_moldyn *moldyn) { moldyn->lvstat=0; t_visual *vis; - vis=v; + vis=&(moldyn->vis); if(moldyn->ewrite) { moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC); @@ -204,7 +208,7 @@ int moldyn_log_init(t_moldyn *moldyn,void *v) { return 0; } -int moldyn_shutdown(t_moldyn *moldyn) { +int moldyn_log_shutdown(t_moldyn *moldyn) { if(moldyn->efd) close(moldyn->efd); if(moldyn->mfd) close(moldyn->efd); @@ -214,6 +218,33 @@ int moldyn_shutdown(t_moldyn *moldyn) { return 0; } +int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { + + int ret; + + ret=moldyn_parse_argv(moldyn,argc,argv); + if(ret<0) return ret; + + ret=moldyn_log_init(moldyn); + if(ret<0) return ret; + + rand_init(&(moldyn->random),NULL,1); + moldyn->random.status|=RAND_STAT_VERBOSE; + + moldyn->status=0; + + return 0; +} + +int moldyn_shutdown(t_moldyn *moldyn) { + + moldyn_log_shutdown(moldyn); + rand_close(&(moldyn->random)); + free(moldyn->atom); + + return 0; +} + int create_lattice(unsigned char type,int element,double mass,double lc, int a,int b,int c,t_atom **atom) { @@ -270,7 +301,7 @@ int destroy_lattice(t_atom *atom) { return 0; } -int thermal_init(t_moldyn *moldyn,t_random *random) { +int thermal_init(t_moldyn *moldyn) { /* * - gaussian distribution of velocities @@ -282,8 +313,10 @@ int thermal_init(t_moldyn *moldyn,t_random *random) { double v,sigma; t_3dvec p_total,delta; t_atom *atom; + t_random *random; atom=moldyn->atom; + random=&(moldyn->random); /* gaussian distribution of velocities */ v3_zero(&p_total); @@ -402,9 +435,13 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) { int link_cell_init(t_moldyn *moldyn) { t_linkcell *lc; + int i; lc=&(moldyn->lc); + /* list log fd */ + lc->listfd=open("/dev/null",O_WRONLY); + /* partitioning the md cell */ lc->nx=moldyn->dim.x/moldyn->cutoff; lc->x=moldyn->dim.x/lc->nx; @@ -413,7 +450,14 @@ int link_cell_init(t_moldyn *moldyn) { 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)); + lc->cells=lc->nx*lc->ny*lc->nz; + lc->subcell=malloc(lc->cells*sizeof(t_list)); + + printf("initializing linked cells (%d)\n",lc->cells); + + for(i=0;icells;i++) + //list_init(&(lc->subcell[i]),1); + list_init(&(lc->subcell[i]),lc->listfd); link_cell_update(moldyn); @@ -461,7 +505,7 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { lc=&(moldyn->lc); nx=lc->nx; ny=lc->ny; - nx=lc->nz; + nz=lc->nz; count1=1; count2=27; a=nx*ny; @@ -512,6 +556,8 @@ int link_cell_shutdown(t_moldyn *moldyn) { for(i=0;inx*lc->ny*lc->nz;i++) list_shutdown(&(moldyn->lc.subcell[i])); + if(lc->listfd) close(lc->listfd); + return 0; } @@ -548,9 +594,6 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; - /* create the neighbour list */ - link_cell_update(moldyn); - /* calculate initial forces */ moldyn->potential_force_function(moldyn); diff --git a/moldyn.h b/moldyn.h index e9f1495..0022208 100644 --- a/moldyn.h +++ b/moldyn.h @@ -23,7 +23,9 @@ typedef struct s_atom { } t_atom; typedef struct s_linkcell { + int listfd; int nx,ny,nz; + int cells; double x,y,z; t_list *subcell; } t_linkcell; @@ -53,6 +55,7 @@ typedef struct s_moldyn { /* energy */ double energy; /* logging & visualization */ + t_visual vis; unsigned char lvstat; unsigned int ewrite; char efb[64]; @@ -71,6 +74,8 @@ typedef struct s_moldyn { void *visual; /* moldyn general status */ unsigned char status; + /* random */ + t_random random; } t_moldyn; typedef struct s_ho_params { @@ -102,6 +107,7 @@ typedef struct s_lj_params { #define MOLDYN_TEMP 273.0 #define MOLDYN_TAU 1.0e-15 +#define MOLDYN_CUTOFF 10.0e-9 #define MOLDYN_RUNS 1000000 #define MOLDYN_INTEGRATE_VERLET 0x00 @@ -132,13 +138,14 @@ typedef struct s_lj_params { int moldyn_usage(char **argv); int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv); -int moldyn_log_init(t_moldyn *moldyn,void *v); +int moldyn_log_init(t_moldyn *moldyn); +int moldyn_init(t_moldyn *moldyn,int argc,char **argv); int moldyn_shutdown(t_moldyn *moldyn); int create_lattice(unsigned char type,int element,double mass,double lc, int a,int b,int c,t_atom **atom); int destroy_lattice(t_atom *atom); -int thermal_init(t_moldyn *moldyn,t_random *random); +int thermal_init(t_moldyn *moldyn); int scale_velocity(t_moldyn *moldyn); double get_e_kin(t_atom *atom,int count); double get_e_pot(t_moldyn *moldyn); diff --git a/posic.c b/posic.c index 236e208..2b18c79 100644 --- a/posic.c +++ b/posic.c @@ -17,9 +17,6 @@ int main(int argc,char **argv) { t_moldyn md; - t_atom *si; - t_visual vis; - t_random random; int a,b,c; double e; @@ -29,90 +26,115 @@ int main(int argc,char **argv) { t_lj_params lj; t_ho_params ho; - /* parse arguments */ - a=moldyn_parse_argv(&md,argc,argv); - if(a<0) return -1; - - /* init */ - moldyn_log_init(&md,&vis); - rand_init(&random,NULL,1); - random.status|=RAND_STAT_VERBOSE; - - /* testing random numbers */ - //for(a=0;a<1000000;a++) - // printf("%f %f\n",rand_get_gauss(&random), - // rand_get_gauss(&random)); - - a=LEN_X; - b=LEN_Y; - c=LEN_Z; + /* + * moldyn init + * + * - parsing argv + * - log init + * - random init + * + */ + moldyn_init(&md,argc,argv); - /* set for 'bounding atoms' */ - vis.dim.x=a*LC_SI; - vis.dim.y=b*LC_SI; - vis.dim.z=c*LC_SI; + /* + * the following overrides possibly set interaction methods by argv !! + */ - /* init lattice - printf("placing silicon atoms ... "); - md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si); - printf("(%d) ok!\n",md.count); - testing purpose */ - md.count=2; - si=malloc(2*sizeof(t_atom)); - si[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0; - si[0].r.y=0; - si[0].r.z=0; - si[0].element=SI; - si[0].mass=M_SI; - si[1].r.x=-si[0].r.x; - si[1].r.y=0; - si[1].r.z=0; - si[1].element=SI; - si[1].mass=M_SI; - /* */ - - /* moldyn init (now si is a valid address) */ - md.atom=si; + /* params */ + lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI; + help=lj.sigma6*lj.sigma6; + lj.sigma6*=help; + lj.sigma12=lj.sigma6*lj.sigma6; + lj.epsilon4=4.0*LJ_EPSILON_SI; + ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI; + ho.spring_constant=1.0; + /* assignement */ md.potential_force_function=lennard_jones; //md.potential_force_function=harmonic_oscillator; - md.cutoff=R_CUTOFF*LC_SI; md.pot_params=&lj; //md.pot_params=&ho; - md.status=0; - md.visual=&vis; - /* dimensions of the simulation cell */ + /* cutoff radius */ + md.cutoff=R_CUTOFF*LC_SI; + + /* + * testing random numbers + */ + +#ifdef DEBUG_RANDOM_NUMBER + for(a=0;a<1000000;a++) + printf("%f %f\n",rand_get_gauss(&(md.random)), + rand_get_gauss(&(md.random))); + return 0; +#endif + + /* + * geometry & particles + */ + + /* simulation cell volume in lattice constants */ + a=LEN_X; + b=LEN_Y; + c=LEN_Z; md.dim.x=a*LC_SI; md.dim.y=b*LC_SI; md.dim.z=c*LC_SI; + /* (un)set to (not) get visualized 'bounding atoms' */ + md.vis.dim.x=a*LC_SI; + md.vis.dim.y=b*LC_SI; + md.vis.dim.z=c*LC_SI; + + /* + * particles + */ + + /* lattice init */ + +#ifndef SIMPLE_TESTING + md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&(md.atom)); + printf("created silicon lattice (#atoms = %d)\n",md.count); +#else + md.count=2; + moldyn->atom=malloc(2*sizeof(t_atom)); + moldyn->atom[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0; + moldyn->atom[0].r.y=0; + moldyn->atom[0].r.z=0; + moldyn->atom[0].element=SI; + moldyn->atom[0].mass=M_SI; + moldyn->atom[1].r.x=-si[0].r.x; + moldyn->atom[1].r.y=0; + moldyn->atom[1].r.z=0; + moldyn->atom[1].element=SI; + moldyn->atom[1].mass=M_SI; +#endif + + /* initial thermal fluctuations of particles */ + +#ifndef SIMPLE_TESTING printf("setting thermal fluctuations (T=%f K)\n",md.t); - // thermal_init(&md,&random); + thermal_init(&md); +#else for(a=0;a