X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=3424d4946b18e4f74185554bb71cb9c7e58893c3;hb=2a902b8d09fd57bc20681448774e2f93894f10aa;hp=e02557c38f78d57e3336fff421697de55c571de1;hpb=dc70c570abec4596355df26ff19756658e33e762;p=physik%2Fposic.git 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; +} +