testing

diff --git a/moldyn.c b/moldyn.c
index 815fc876717415cb17f1464ed6db5a04b50b0c7e..26f298fbb58f2488f048d04dd4fbaad550e9563e 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -635,7 +635,6 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
        /* zero absolute time */
        moldyn->time=0.0;
-
        for(sched=0;sched<moldyn->schedule.content_count;sched++) {
 
                /* setting amount of runs and finite time step size */
@@ -759,15 +758,15 @@ int velocity_verlet(t_moldyn *moldyn) {
 int potential_force_calc(t_moldyn *moldyn) {
 
        int i,j,k,count;
-       t_atom *atom,*btom,*ktom;
+       t_atom *itom,*jtom,*ktom;
        t_linkcell *lc;
-       t_list neighbour[27],neighbourk[27];
-       t_list *this,*thisk;
-       u8 bc,bck;
+       t_list neighbour_i[27],neighbour_j[27];
+       t_list *this,*that;
+       u8 bc_ij,bc_ijk;
        int countn,dnlc;
 
        count=moldyn->count;
-       atom=moldyn->atom;
+       itom=moldyn->atom;
        lc=&(moldyn->lc);
 
        /* reset energy */
@@ -776,85 +775,117 @@ int potential_force_calc(t_moldyn *moldyn) {
        for(i=0;i<count;i++) {
 
                /* reset force */
-               v3_zero(&(atom[i].f));
+               v3_zero(&(itom[i].f));
 
                /* single particle potential/force */
-               if(atom[i].attr&ATOM_ATTR_1BP)
-                       moldyn->func1b(moldyn,&(atom[i]));
+               if(itom[i].attr&ATOM_ATTR_1BP)
+                       moldyn->func1b(moldyn,&(itom[i]));
 
                /* 2 body pair potential/force */
-               if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
+               if(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
        
                        link_cell_neighbour_index(moldyn,
-                               (atom[i].r.x+moldyn->dim.x/2)/lc->x,
-                               (atom[i].r.y+moldyn->dim.y/2)/lc->y,
-                               (atom[i].r.z+moldyn->dim.z/2)/lc->z,
-                               neighbour);
+                               (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+                               (itom[i].r.y+moldyn->dim.y/2)/lc->y,
+                               (itom[i].r.z+moldyn->dim.z/2)/lc->z,
+                               neighbour_i);
 
                        countn=lc->countn;
                        dnlc=lc->dnlc;
 
                        for(j=0;j<countn;j++) {
 
-                               this=&(neighbour[j]);
+                               this=&(neighbour_i[j]);
                                list_reset(this);
 
                                if(this->start==NULL)
                                        continue;
 
-                               bc=(j<dnlc)?0:1;
+                               bc_ij=(j<dnlc)?0:1;
 
                                do {
-                                       btom=this->current->data;
+                                       jtom=this->current->data;
 
-                                       if(btom==&(atom[i]))
+                                       if(jtom==&(itom[i]))
                                                continue;
 
-                                       if((btom->attr&ATOM_ATTR_2BP)&
-                                          (atom[i].attr&ATOM_ATTR_2BP))
+                                       if((jtom->attr&ATOM_ATTR_2BP)&
+                                          (itom[i].attr&ATOM_ATTR_2BP))
                                                moldyn->func2b(moldyn,
-                                                              &(atom[i]),
-                                                              btom,
-                                                              bc);
+                                                              &(itom[i]),
+                                                              jtom,
+                                                              bc_ij);
 
                                        /* 3 body potential/force */
 
-                                       if(!(atom[i].attr&ATOM_ATTR_3BP)||
-                                          !(btom->attr&ATOM_ATTR_3BP))
+                                       if(!(itom[i].attr&ATOM_ATTR_3BP)||
+                                          !(jtom->attr&ATOM_ATTR_3BP))
                                                continue;
 
                                        link_cell_neighbour_index(moldyn,
-                                          (btom->r.x+moldyn->dim.x/2)/lc->x,
-                                          (btom->r.y+moldyn->dim.y/2)/lc->y,
-                                          (btom->r.z+moldyn->dim.z/2)/lc->z,
-                                          neighbourk);
+                                          (jtom->r.x+moldyn->dim.x/2)/lc->x,
+                                          (jtom->r.y+moldyn->dim.y/2)/lc->y,
+                                          (jtom->r.z+moldyn->dim.z/2)/lc->z,
+                                          neighbour_j);
 
+                                       /* neighbours of j */
                                        for(k=0;k<lc->countn;k++) {
 
-                                               thisk=&(neighbourk[k]);
-                                               list_reset(thisk);
+                                               that=&(neighbour_j[k]);
+                                               list_reset(that);
+                                       
+                                               if(that->start==NULL)
+                                                       continue;
+
+                                               bc_ijk=(k<lc->dnlc)?0:1;
+
+                                               do {
+
+                       ktom=that->current->data;
+
+                       if(!(ktom->attr&ATOM_ATTR_3BP))
+                               continue;
+
+                       if(ktom==jtom)
+                               continue;
+
+                       if(ktom==&(itom[i]))
+                               continue;
+
+                       moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+
+                                               } while(list_next(that)!=\
+                                                       L_NO_NEXT_ELEMENT);
+
+                                       }
+                                       
+                                       /* neighbours of i */
+                                       for(k=0;k<countn;k++) {
+
+                                               that=&(neighbour_i[k]);
+                                               list_reset(that);
                                        
-                                               if(thisk->start==NULL)
+                                               if(that->start==NULL)
                                                        continue;
 
-                                               bck=(k<lc->dnlc)?0:1;
+                                               bc_ijk=(k<dnlc)?0:1;
 
                                                do {
 
-                       ktom=thisk->current->data;
+                       ktom=that->current->data;
 
                        if(!(ktom->attr&ATOM_ATTR_3BP))
                                continue;
 
-                       if(ktom==btom)
+                       if(ktom==jtom)
                                continue;
 
-                       if(ktom==&(atom[i]))
+                       if(ktom==&(itom[i]))
                                continue;
 
-                       moldyn->func3b(moldyn,&(atom[i]),btom,ktom,bck);
+                       moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
 
-                                               } while(list_next(thisk)!=\
+                                               } while(list_next(that)!=\
                                                        L_NO_NEXT_ELEMENT);
 
                                        }

diff --git a/sic.c b/sic.c
index 2cd32eb04853a4c14518561999f9d6a429d9e5a1..82b371018cfc11fa99bef055030b541b079b03d2 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -86,6 +86,7 @@ printf("%d\n",sizeof(t_atom));
        /* cutoff radius */
        printf("[sic] setting cutoff radius\n");
        set_cutoff(&md,TM_S_SI);
+       //set_cutoff(&md,1.0*LC_SI);
 
        /* set (initial) dimensions of simulation volume */
        printf("[sic] setting dimensions\n");
@@ -99,6 +100,7 @@ printf("%d\n",sizeof(t_atom));
        printf("[sic] creating atoms\n");
        create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,
+                      //ATOM_ATTR_2BP,
                       0,4,4,4);
 
        /* setting a nearest neighbour distance for the moldyn checks */
@@ -106,7 +108,7 @@ printf("%d\n",sizeof(t_atom));
 
        /* set temperature */
        printf("[sic] setting temperature\n");
-       set_temperature(&md,273.0+450.0);
+       set_temperature(&md,10.0);
        
        /* initial thermal fluctuations of particles */
        printf("[sic] thermal init\n");
@@ -114,12 +116,12 @@ printf("%d\n",sizeof(t_atom));
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");
-       moldyn_add_schedule(&md,10000,1.0e-15);
+       moldyn_add_schedule(&md,1000,1.0e-15);
 
        /* activate logging */
        printf("[sic] activate logging\n");
-       moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",100);
-       moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",100);
+       moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",10);
+       moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",10);
 
        /*
         * let's do the actual md algorithm now