lj fixes + first lines of tersoff potential

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;
+}
+

diff --git a/moldyn.h b/moldyn.h
index 3016bc4..49ca0b2 100644 (file)
--- 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 (file)
--- 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;