From e6f456c0fa807b86e1b25996e70efcdcfe390ea5 Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 28 Apr 2006 15:02:29 +0000 Subject: [PATCH] implemented link cell method, segfaulting by now! --- moldyn.c | 506 +++++++++++++++++++++++++++++++------------------------ moldyn.h | 34 ++-- posic.c | 41 ++--- 3 files changed, 319 insertions(+), 262 deletions(-) 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; } diff --git a/moldyn.h b/moldyn.h index 53a22fd..e9f1495 100644 --- a/moldyn.h +++ b/moldyn.h @@ -18,10 +18,8 @@ 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 verlet; /* verlet neighbour list */ } t_atom; typedef struct s_linkcell { @@ -38,15 +36,11 @@ typedef struct s_moldyn { t_atom *atom; t_3dvec dim; /* potential, force & parameters */ - double (*potential)(struct s_moldyn *moldyn); - int (*force)(struct s_moldyn *moldyn); + int (*potential_force_function)(struct s_moldyn *moldyn); void *pot_params; - /* cut off radius, verlet list & co */ + /* cut off radius */ double cutoff; double cutoff_square; - double r_verlet; - double dr_max1; - double dr_max2; /* linked list / cell method */ t_linkcell lc; /* temperature */ @@ -55,6 +49,9 @@ typedef struct s_moldyn { int (*integrate)(struct s_moldyn *moldyn); int time_steps; double tau; + double tau_square; + /* energy */ + double energy; /* logging & visualization */ unsigned char lvstat; unsigned int ewrite; @@ -107,6 +104,13 @@ typedef struct s_lj_params { #define MOLDYN_TAU 1.0e-15 #define MOLDYN_RUNS 1000000 +#define MOLDYN_INTEGRATE_VERLET 0x00 +#define MOLDYN_INTEGRATE_DEFAULT MOLDYN_INTEGRATE_VERLET + +#define MOLDYN_POTENTIAL_HO 0x00 +#define MOLDYN_POTENTIAL_LJ 0x01 +#define MOLDYN_POTENTIAL_DEFAULT MOLDYN_POTENTIAL_LJ + /* phsical values */ #define K_BOLTZMANN 1.3807e-27 /* Nm/K */ @@ -134,8 +138,8 @@ 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 count); -int scale_velocity(t_moldyn *moldyn,int count); +int thermal_init(t_moldyn *moldyn,t_random *random); +int scale_velocity(t_moldyn *moldyn); double get_e_kin(t_atom *atom,int count); double get_e_pot(t_moldyn *moldyn); double get_total_energy(t_moldyn *moldyn); @@ -143,19 +147,15 @@ 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_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell); int link_cell_shutdown(t_moldyn *moldyn); int moldyn_integrate(t_moldyn *moldyn); int velocity_verlet(t_moldyn *moldyn); -double potential_harmonic_oscillator(t_moldyn *moldyn); -int force_harmonic_oscillator(t_moldyn *moldyn); -double potential_lennard_jones(t_moldyn *moldyn); -int force_lennard_jones(t_moldyn *moldyn); +int harmonic_oscillator(t_moldyn *moldyn); +int lennard_jones(t_moldyn *moldyn); #endif diff --git a/posic.c b/posic.c index 7cde6e2..236e208 100644 --- a/posic.c +++ b/posic.c @@ -25,7 +25,6 @@ int main(int argc,char **argv) { double e; double help; t_3dvec p; - int count; t_lj_params lj; t_ho_params ho; @@ -53,12 +52,12 @@ int main(int argc,char **argv) { vis.dim.y=b*LC_SI; vis.dim.z=c*LC_SI; - /* init lattice */ + /* init lattice printf("placing silicon atoms ... "); - count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si); - printf("(%d) ok!\n",count); - /* testing purpose - count=2; + 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; @@ -70,22 +69,15 @@ int main(int argc,char **argv) { si[1].r.z=0; si[1].element=SI; si[1].mass=M_SI; - */ + /* */ /* moldyn init (now si is a valid address) */ - md.count=count; md.atom=si; - md.potential=potential_lennard_jones; - md.force=force_lennard_jones; - //md.potential=potential_harmonic_oscillator; - //md.force=force_harmonic_oscillator; + md.potential_force_function=lennard_jones; + //md.potential_force_function=harmonic_oscillator; 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; - //md.time_steps=RUNS; - //md.tau=TAU; md.status=0; md.visual=&vis; /* dimensions of the simulation cell */ @@ -93,23 +85,19 @@ int main(int argc,char **argv) { 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); - //for(a=0;a