From 2a902b8d09fd57bc20681448774e2f93894f10aa Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 3 Aug 2006 19:58:16 +0000 Subject: [PATCH] lj fixes + first lines of tersoff potential --- moldyn.c | 285 +++++++++++++++++++++++++++++++++++++++++++++++++++++-- moldyn.h | 11 ++- posic.c | 15 +-- 3 files changed, 296 insertions(+), 15 deletions(-) diff --git a/moldyn.c b/moldyn.c index e02557c..3424d49 100644 --- a/moldyn.c +++ b/moldyn.c @@ -457,7 +457,7 @@ int link_cell_init(t_moldyn *moldyn) { for(i=0;icells;i++) //list_init(&(lc->subcell[i]),1); - list_init(&(lc->subcell[i]),lc->listfd); + list_init(&(lc->subcell[i])); link_cell_update(moldyn); @@ -681,10 +681,14 @@ int velocity_verlet(t_moldyn *moldyn) { } /* neighbour list update */ +printf("list update ...\n"); link_cell_update(moldyn); +printf("done\n"); /* forces depending on chosen potential */ +printf("calc potential/force ...\n"); moldyn->potential_force_function(moldyn); +printf("done\n"); for(i=0;icutoff_square) { d=1.0/d; /* 1/r^2 */ - h1=d*d; /* 1/r^4 */ + h2=d*d; /* 1/r^4 */ h2*=d; /* 1/r^6 */ h1=h2*h2; /* 1/r^12 */ u+=eps*(sig12*h1-sig6*h2); @@ -869,7 +873,7 @@ int lennard_jones(t_moldyn *moldyn) { h1*=d; /* 1/r^14 */ h2*=6*sig6; h1*=12*sig12; - d=-h1+h2; + d=+h1-h2; d*=eps; v3_scale(&force,&distance,d); v3_add(&(atom[i].f),&(atom[i].f),&force); @@ -888,7 +892,7 @@ int lennard_jones(t_moldyn *moldyn) { 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*d; /* 1/r^4 */ h2*=d; /* 1/r^6 */ h1=h2*h2; /* 1/r^12 */ u+=eps*(sig12*h1-sig6*h2); @@ -896,7 +900,7 @@ int lennard_jones(t_moldyn *moldyn) { h1*=d; /* 1/r^14 */ h2*=6*sig6; h1*=12*sig12; - d=-h1+h2; + d=+h1-h2; d*=eps; v3_scale(&force,&distance,d); v3_add(&(atom[i].f),&(atom[i].f), @@ -920,7 +924,7 @@ int lennard_jones(t_moldyn *moldyn) { 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*d; /* 1/r^4 */ h2*=d; /* 1/r^6 */ h1=h2*h2; /* 1/r^12 */ u+=eps*(sig12*h1-sig6*h2); @@ -928,7 +932,7 @@ int lennard_jones(t_moldyn *moldyn) { h1*=d; /* 1/r^14 */ h2*=6*sig6; h1*=12*sig12; - d=-h1+h2; + d=+h1-h2; d*=eps; v3_scale(&force,&distance,d); v3_add(&(atom[i].f),&(atom[i].f), @@ -945,3 +949,270 @@ int lennard_jones(t_moldyn *moldyn) { return 0; } +/* tersoff potential & force for 2 sorts of atoms */ + +int tersoff(t_moldyn *moldyn) { + + t_tersoff_params *params; + t_atom *atom,*btom,*ktom; + t_linkcell *lc; + t_list *this,*thisk,neighbour[27],neighbourk[27]; + int i,j,k,c,ck; + int count; + double u; + int ni,nj,nk; + int ki,kj,kk; + + + params=moldyn->pot_params; + atom=moldyn->atom; + lc=&(moldyn->lc); + count=moldyn->count; + + /* reset enrgy 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; + c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour); + + /* + * processing cell of atom i + * => no need to check for empty list (1 element at minimum) + */ + this=&(neighbour[0]); + list_reset(this); + do { + btom=this->current->data; + if(btom==&(atom[i])) + continue; + + /* 2 body stuff */ + + /* determine cell neighbours of btom */ + ki=(btom->r.x+(moldyn->dim.x/2))/lc->x; + kj=(btom->r.y+(moldyn->dim.y/2))/lc->y; + kk=(btom->r.z+(moldyn->dim.z/2))/lc->z; + ck=link_cell_neighbour_index(moldyn,ki,kj,kk, + neighbourk); + + /* cell of btom */ + thisk=&(neighbourk[0]); + list_reset(thisk); + do { + ktom=thisk->current->data; + if(ktom==btom) + continue; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff (1) */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + /* direct neighbours of btom cell */ + for(k=1;kstart!=NULL) { + + do { + ktom=thisk->current->data; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff (2) */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + } + } + + /* indirect neighbours of btom cell */ + for(k=ck;k<27;k++) { + thisk=&(neighbourk[k]); + list_reset(thisk); + if(thisk->start!=NULL) { + + do { + ktom=thisk->current->data; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + } + } + + + } while(list_next(this)!=L_NO_NEXT_ELEMENT); + + /* + * direct neighbour cells of atom i + */ + for(j=1;jstart!=NULL) { + + do { + btom=this->current->data; + + /* 2 body stuff */ + + + /* determine cell neighbours of btom */ + ki=(btom->r.x+(moldyn->dim.x/2))/lc->x; + kj=(btom->r.y+(moldyn->dim.y/2))/lc->y; + kk=(btom->r.z+(moldyn->dim.z/2))/lc->z; + ck=link_cell_neighbour_index(moldyn,ki,kj,kk, + neighbourk); + + /* cell of btom */ + thisk=&(neighbourk[0]); + list_reset(thisk); + do { + ktom=thisk->current->data; + if(ktom==btom) + continue; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff (1) */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + /* direct neighbours of btom cell */ + for(k=1;kstart!=NULL) { + + do { + ktom=thisk->current->data; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff (2) */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + } + } + + /* indirect neighbours of btom cell */ + for(k=ck;k<27;k++) { + thisk=&(neighbourk[k]); + list_reset(thisk); + if(thisk->start!=NULL) { + + do { + ktom=thisk->current->data; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff (3) */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + } + } + + + } while(list_next(this)!=L_NO_NEXT_ELEMENT); + + } + } + + /* + * indirect neighbour cells of atom i + */ + for(j=c;j<27;j++) { + this=&(neighbour[j]); + list_reset(this); + if(this->start!=NULL) { + + do { + btom=this->current->data; + + /* 2 body stuff */ + + + /* determine cell neighbours of btom */ + ki=(btom->r.x+(moldyn->dim.x/2))/lc->x; + kj=(btom->r.y+(moldyn->dim.y/2))/lc->y; + kk=(btom->r.z+(moldyn->dim.z/2))/lc->z; + ck=link_cell_neighbour_index(moldyn,ki,kj,kk, + neighbourk); + + /* cell of btom */ + thisk=&(neighbourk[0]); + list_reset(thisk); + do { + ktom=thisk->current->data; + if(ktom==btom) + continue; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff (1) */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + /* direct neighbours of btom cell */ + for(k=1;kstart!=NULL) { + + do { + ktom=thisk->current->data; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff (2) */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + } + } + + /* indirect neighbours of btom cell */ + for(k=ck;k<27;k++) { + thisk=&(neighbourk[k]); + list_reset(thisk); + if(thisk->start!=NULL) { + + do { + ktom=thisk->current->data; + if(ktom==&(atom[i])) + continue; + + /* 3 body stuff (3) */ + + } while(list_next(thisk)!=L_NO_NEXT_ELEMENT); + + } + } + + + } while(list_next(this)!=L_NO_NEXT_ELEMENT); + + } + } + + } + + moldyn->energy=0.5*u; + + return 0; +} + diff --git a/moldyn.h b/moldyn.h index 3016bc4..49ca0b2 100644 --- a/moldyn.h +++ b/moldyn.h @@ -89,6 +89,15 @@ typedef struct s_lj_params { double epsilon4; } t_lj_params; +typedef struct s_tersoff_params { + double l_1,l_2; + double m_1,m_2; + double a_1,a_2; + double b_1,b_2; + double r_1,r_2; + double s_1,s_2; +} t_tersoff_params; + /* * defines */ @@ -131,7 +140,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.20*sqrt(3.0)*LC_SI)/1.122462) /* m */ +#define LJ_SIGMA_SI ((0.25*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 fba11d2..d53af6f 100644 --- a/posic.c +++ b/posic.c @@ -25,6 +25,7 @@ int main(int argc,char **argv) { t_lj_params lj; t_ho_params ho; + t_tersoff_params tp; /* * moldyn init @@ -95,9 +96,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=3; + md.count=2; 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.x=0.23*sqrt(3.0)*LC_SI/2.0; md.atom[0].r.y=0; md.atom[0].r.z=0; md.atom[0].element=SI; @@ -108,11 +109,11 @@ int main(int argc,char **argv) { 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[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; -- 2.20.1