From dc70c570abec4596355df26ff19756658e33e762 Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 13 Jun 2006 15:21:37 +0000 Subject: [PATCH] fixed several stupid mistakes --- moldyn.c | 99 +++++++++++++++++++++++++++---------------------- moldyn.h | 2 +- posic.c | 25 +++++++++---- posic.h | 3 +- visual/visual.c | 1 + 5 files changed, 77 insertions(+), 53 deletions(-) diff --git a/moldyn.c b/moldyn.c index 707fe03..e02557c 100644 --- a/moldyn.c +++ b/moldyn.c @@ -456,8 +456,8 @@ int link_cell_init(t_moldyn *moldyn) { 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); + //list_init(&(lc->subcell[i]),1); + list_init(&(lc->subcell[i]),lc->listfd); link_cell_update(moldyn); @@ -512,7 +512,6 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { cell[0]=lc->subcell[i+j*nx+k*a]; - printf("%d\n",i+j*nx+k*a); for(ci=-1;ci<=1;ci++) { bx=0; x=i+ci; @@ -535,16 +534,11 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { bz=1; } if(!(ci|cj|ck)) continue; - printf(" %d %d %d \n",x,y,z); if(bx|by|bz) { cell[--count2]=lc->subcell[x+y*nx+z*a]; - printf("%d\n",x+y*nx+z*a); - printf("--- %d\n",count2); } else { cell[count1++]=lc->subcell[x+y*nx+z*a]; - printf("%d\n",x+y*nx+z*a); - printf("--- %d\n",count1); } } } @@ -585,6 +579,9 @@ int moldyn_integrate(t_moldyn *moldyn) { int fd; char fb[128]; + /* initialize linked cell method */ + link_cell_init(moldyn); + /* logging & visualization */ e=moldyn->ewrite; m=moldyn->mwrite; @@ -606,9 +603,6 @@ int moldyn_integrate(t_moldyn *moldyn) { for(i=0;itime_steps;i++) { - /* neighbour list update */ - link_cell_update(moldyn); - /* integration step */ moldyn->integrate(moldyn); @@ -686,6 +680,9 @@ int velocity_verlet(t_moldyn *moldyn) { v3_add(&(atom[i].v),&(atom[i].v),&delta); } + /* neighbour list update */ + link_cell_update(moldyn); + /* forces depending on chosen potential */ moldyn->potential_force_function(moldyn); @@ -727,30 +724,43 @@ int harmonic_oscillator(t_moldyn *moldyn) { equi_dist=params->equilibrium_distance; count=moldyn->count; + /* reset energy counter */ u=0.0; + for(i=0;idim.x/2))/lc->x; nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y; nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z; - printf("%d %d %d\n",ni,nj,nk); c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour); - /* processing cell of atom i */ + /* + * processing cell of atom i + * => no need to check for empty list (1 element at minimum) + */ this=&(neighbour[0]); - list_reset(this); /* list has 1 element at minimum */ + list_reset(this); 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)); - v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d))); - v3_add(&(atom[i].f),&(atom[i].f),&force); + 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); - /* neighbours not doing boundary condition overflow */ + /* + * direct neighbour cells + * => no boundary condition check necessary + */ for(j=1;j check boundary conditions + */ for(j=c;j<27;j++) { this=&(neighbour[j]); list_reset(this); /* check boundary conditions */ @@ -796,7 +809,7 @@ int harmonic_oscillator(t_moldyn *moldyn) { } } - moldyn->energy=u; + moldyn->energy=0.5*u; return 0; } @@ -824,18 +837,18 @@ int lennard_jones(t_moldyn *moldyn) { sig6=params->sigma6; sig12=params->sigma12; + /* reset energy counter */ u=0.0; + for(i=0;idim.x/2))/lc->x; nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y; nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z; - printf("hier atom = %08x\n",&(atom[i])); c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour); - printf("da atom = %08x\n",&(atom[i])); - printf("da atom = %08x\n",&(moldyn->atom[i])); - - printf("c = %d (%d %d %d)\n",c,ni,nj,nk); /* processing cell of atom i */ this=&(neighbour[0]); @@ -844,23 +857,23 @@ int lennard_jones(t_moldyn *moldyn) { btom=this->current->data; if(btom==&(atom[i])) continue; - puts("foo"); v3_sub(&distance,&(atom[i].r),&(btom->r)); - puts("foo"); - 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); - printf("test!!\n"); + d=v3_absolute_square(&distance); /* 1/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); /* neighbours not doing boundary condition overflow */ @@ -870,7 +883,6 @@ int lennard_jones(t_moldyn *moldyn) { if(this->start!=NULL) { do { - printf("in bound: %d\n",j); btom=this->current->data; v3_sub(&distance,&(atom[i].r),&(btom->r)); d=v3_absolute_square(&distance); /* r^2 */ @@ -902,7 +914,6 @@ int lennard_jones(t_moldyn *moldyn) { if(this->start!=NULL) { do { - printf("out bound: %d\n",j); btom=this->current->data; v3_sub(&distance,&(atom[i].r),&(btom->r)); v3_per_bound(&distance,&(moldyn->dim)); @@ -929,7 +940,7 @@ int lennard_jones(t_moldyn *moldyn) { } } - moldyn->energy=u; + moldyn->energy=0.5*u; return 0; } diff --git a/moldyn.h b/moldyn.h index 0022208..3016bc4 100644 --- a/moldyn.h +++ b/moldyn.h @@ -131,7 +131,7 @@ typedef struct s_lj_params { #define SI 0x0e #define LC_SI 0.543105e-9 /* m */ #define M_SI (28.085*AMU) /* kg */ -#define LJ_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) /* m */ +#define LJ_SIGMA_SI ((0.20*sqrt(3.0)*LC_SI)/1.122462) /* m */ #define LJ_EPSILON_SI (2.1678*1.60e-19) /* Nm */ /* function prototypes */ diff --git a/posic.c b/posic.c index c185d16..fba11d2 100644 --- a/posic.c +++ b/posic.c @@ -48,7 +48,7 @@ int main(int argc,char **argv) { 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; + ho.spring_constant=1; /* assignement */ md.potential_force_function=lennard_jones; //md.potential_force_function=harmonic_oscillator; @@ -95,9 +95,9 @@ int main(int argc,char **argv) { 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; - md.atom=malloc(2*sizeof(t_atom)); - md.atom[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0; + md.count=3; + md.atom=malloc(md.count*sizeof(t_atom)); + md.atom[0].r.x=0.21*sqrt(3.0)*LC_SI/2.0; md.atom[0].r.y=0; md.atom[0].r.z=0; md.atom[0].element=SI; @@ -107,6 +107,18 @@ int main(int argc,char **argv) { md.atom[1].r.z=0; md.atom[1].element=SI; md.atom[1].mass=M_SI; + + md.atom[2].r.x=0.5*(a-1)*LC_SI; + md.atom[2].r.y=0.5*(b-1)*LC_SI; + md.atom[2].r.z=0; + md.atom[2].element=C; + md.atom[2].mass=M_C; + + //md.atom[3].r.x=0.5*(a-1)*LC_SI; + //md.atom[3].r.y=0; + //md.atom[3].r.z=0; + //md.atom[3].element=SI; + //md.atom[3].mass=M_SI; #endif /* initial thermal fluctuations of particles */ @@ -116,6 +128,8 @@ int main(int argc,char **argv) { thermal_init(&md); #else for(a=0;afd,"rotate y 10\n"); dprintf(v->fd,"set ambient 20\n"); dprintf(v->fd,"set specular on\n"); + dprintf(v->fd,"label\n"); sprintf(file,"%s-%.15f.ppm",v->fb,time); dprintf(v->fd,"write ppm %s\n",file); dprintf(v->fd,"zap\n"); -- 2.20.1