lj fixes + first lines of tersoff potential

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index e02557c..3424d49 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -457,7 +457,7 @@ int link_cell_init(t_moldyn *moldyn) {
 
        for(i=0;i<lc->cells;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;i<count;i++) {
                /* again velocities */
@@ -861,7 +865,7 @@ int lennard_jones(t_moldyn *moldyn) {
                        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*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;i<count;i++) {
+               /* reset force */
+               v3_zero(&(atom[i].f));
+
+               /* determin cell neighbours */
+               ni=(atom[i].r.x+(moldyn->dim.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;k<ck;k++) {
+                               thisk=&(neighbourk[k]);
+                               list_reset(thisk);
+                               if(thisk->start!=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;j<c;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;k<ck;k++) {
+                               thisk=&(neighbourk[k]);
+                               list_reset(thisk);
+                               if(thisk->start!=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;k<ck;k++) {
+                               thisk=&(neighbourk[k]);
+                               list_reset(thisk);
+                               if(thisk->start!=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;
+}
+