added static lists (but: glibc double free prob (for both lists!))

diff --git a/Makefile b/Makefile
index 73645df..de46da4 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -6,6 +6,7 @@ CFLAGS += -g
 #CFLAGS += -ffloat-store
 
 CFLAGS += -DALBE
+CFLAGS += -DSTATIC_LISTS
 #CFLAGS += -DDEBUG
 CFLAGS += -DDSTART=400 -DDEND=600
 #CFLAGS += -DVDEBUG
@@ -29,4 +30,4 @@ pair_correlation_calc: $(DEPS)
 
 .PHONY:clean
 clean:
-       rm -vf sic postproc fluctuation_calc *.o */*.o
+       rm -vf sic postproc fluctuation_calc pair_correlation_calc *.o */*.o

diff --git a/moldyn.c b/moldyn.c
index 2e0e468..6d5f9e4 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -1244,6 +1244,168 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 
 /* linked list / cell method */
 
+#ifdef STATIC_LISTS
+
+int link_cell_init(t_moldyn *moldyn,u8 vol) {
+
+       t_linkcell *lc;
+       int i;
+       int *foo;
+
+       lc=&(moldyn->lc);
+
+       /* partitioning the md cell */
+       lc->nx=moldyn->dim.x/moldyn->cutoff;
+       lc->x=moldyn->dim.x/lc->nx;
+       lc->ny=moldyn->dim.y/moldyn->cutoff;
+       lc->y=moldyn->dim.y/lc->ny;
+       lc->nz=moldyn->dim.z/moldyn->cutoff;
+       lc->z=moldyn->dim.z/lc->nz;
+
+       lc->cells=lc->nx*lc->ny*lc->nz;
+       lc->subcell=malloc(lc->cells*sizeof(int*));
+
+       if(lc->cells<27)
+               printf("[moldyn] FATAL: less then 27 subcells!\n");
+
+       if(vol) {
+               printf("[moldyn] initializing 'static' linked cells (%d)\n",
+                      lc->cells);
+               printf("  x: %d x %f A\n",lc->nx,lc->x);
+               printf("  y: %d x %f A\n",lc->ny,lc->y);
+               printf("  z: %d x %f A\n",lc->nz,lc->z);
+       }
+
+       /* list init */
+       for(i=0;i<lc->cells;i++) {
+               lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int));
+               if(lc->subcell[i]==NULL) {
+                       perror("[moldyn] list init (malloc)");
+                       return -1;
+               }
+//if(i==0) printf(" --- add one here! %d %p %p ----\n",i,lc->subcell,lc->subcell[0]);
+       }
+
+       /* update the list */
+       link_cell_update(moldyn);
+
+       return 0;
+}
+
+int link_cell_update(t_moldyn *moldyn) {
+
+       int count,i,j,k;
+       int nx,ny;
+       t_atom *atom;
+       t_linkcell *lc;
+       int p;
+
+       atom=moldyn->atom;
+       lc=&(moldyn->lc);
+
+       nx=lc->nx;
+       ny=lc->ny;
+
+       for(i=0;i<lc->cells;i++)
+               memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+
+       for(count=0;count<moldyn->count;count++) {
+               i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
+               j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
+               k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
+
+               p=0;
+               while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
+                       p++;
+
+               if(p>=MAX_ATOMS_PER_LIST) {
+                       printf("[moldyn] FATAL: amount of atoms too high!\n");
+                       return -1;
+               }
+
+               lc->subcell[i+j*nx+k*nx*ny][p]=count;
+       }
+
+       return 0;
+}
+
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int **cell) {
+
+       t_linkcell *lc;
+       int a;
+       int count1,count2;
+       int ci,cj,ck;
+       int nx,ny,nz;
+       int x,y,z;
+       u8 bx,by,bz;
+
+       lc=&(moldyn->lc);
+       nx=lc->nx;
+       ny=lc->ny;
+       nz=lc->nz;
+       count1=1;
+       count2=27;
+       a=nx*ny;
+
+       cell[0]=lc->subcell[i+j*nx+k*a];
+       for(ci=-1;ci<=1;ci++) {
+               bx=0;
+               x=i+ci;
+               if((x<0)||(x>=nx)) {
+                       x=(x+nx)%nx;
+                       bx=1;
+               }
+               for(cj=-1;cj<=1;cj++) {
+                       by=0;
+                       y=j+cj;
+                       if((y<0)||(y>=ny)) {
+                               y=(y+ny)%ny;
+                               by=1;
+                       }
+                       for(ck=-1;ck<=1;ck++) {
+                               bz=0;
+                               z=k+ck;
+                               if((z<0)||(z>=nz)) {
+                                       z=(z+nz)%nz;
+                                       bz=1;
+                               }
+                               if(!(ci|cj|ck)) continue;
+                               if(bx|by|bz) {
+                                       cell[--count2]=lc->subcell[x+y*nx+z*a];
+                               }
+                               else {
+                                       cell[count1++]=lc->subcell[x+y*nx+z*a];
+                               }
+                       }
+               }
+       }
+
+       lc->dnlc=count1;
+
+       return count1;
+}
+
+int link_cell_shutdown(t_moldyn *moldyn) {
+
+       int i;
+       t_linkcell *lc;
+       int *foo;
+
+       lc=&(moldyn->lc);
+
+       for(i=0;i<lc->cells;i++)
+{
+//printf(" --- free %p , %d\n",lc->subcell[i],i);
+               free(lc->subcell[i]);
+}
+
+       free(lc->subcell);
+
+       return 0;
+}
+
+#else
+
 int link_cell_init(t_moldyn *moldyn,u8 vol) {
 
        t_linkcell *lc;
@@ -1266,7 +1428,8 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
                printf("[moldyn] FATAL: less then 27 subcells!\n");
 
        if(vol) {
-               printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
+               printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
+                      lc->cells);
                printf("  x: %d x %f A\n",lc->nx,lc->x);
                printf("  y: %d x %f A\n",lc->ny,lc->y);
                printf("  z: %d x %f A\n",lc->nz,lc->z);
@@ -1307,6 +1470,7 @@ int link_cell_update(t_moldyn *moldyn) {
                k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
                list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
                                     &(atom[count]));
+//if(i==0&&j==0&&k==0) printf(" --- add one here! %d %p ----\n",count,lc->subcell[0].current);
        }
 
        return 0;
@@ -1375,14 +1539,20 @@ int link_cell_shutdown(t_moldyn *moldyn) {
 
        lc=&(moldyn->lc);
 
-       for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
+printf("FOO:\n");
+       for(i=0;i<lc->nx*lc->ny*lc->nz;i++) {
+printf(" %d\n",i);
                list_destroy_f(&(moldyn->lc.subcell[i]));
+printf(" %d!\n",i);
+}
 
        free(lc->subcell);
 
        return 0;
 }
 
+#endif
+
 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
 
        int count;
@@ -1676,15 +1846,24 @@ int potential_force_calc(t_moldyn *moldyn) {
        t_atom *itom,*jtom,*ktom;
        t_virial *virial;
        t_linkcell *lc;
+#ifdef STATIC_LISTS
+       int *neighbour_i[27];
+       int p,q;
+       t_atom *atom;
+#else
        t_list neighbour_i[27];
        t_list neighbour_i2[27];
        t_list *this,*that;
+#endif
        u8 bc_ij,bc_ik;
        int dnlc;
 
        count=moldyn->count;
        itom=moldyn->atom;
        lc=&(moldyn->lc);
+#ifdef STATIC_LISTS
+       atom=moldyn->atom;
+#endif
 
        /* reset energy */
        moldyn->energy=0.0;
@@ -1739,14 +1918,33 @@ int potential_force_calc(t_moldyn *moldyn) {
                if(moldyn->func2b) {
                        for(j=0;j<27;j++) {
 
+                               bc_ij=(j<dnlc)?0:1;
+#ifdef STATIC_LISTS
+                               p=0;
+
+                               while(neighbour_i[j][p]!=0) {
+
+                                       jtom=&(atom[neighbour_i[j][p]]);
+                                       p++;
+
+                                       if(jtom==&(itom[i]))
+                                               continue;
+
+                                       if((jtom->attr&ATOM_ATTR_2BP)&
+                                          (itom[i].attr&ATOM_ATTR_2BP)) {
+                                               moldyn->func2b(moldyn,
+                                                              &(itom[i]),
+                                                              jtom,
+                                                              bc_ij);
+                                       }
+                               }
+#else
                                this=&(neighbour_i[j]);
                                list_reset_f(this);
 
                                if(this->start==NULL)
                                        continue;
 
-                               bc_ij=(j<dnlc)?0:1;
-
                                do {
                                        jtom=this->current->data;
 
@@ -1761,6 +1959,7 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                               bc_ij);
                                        }
                                } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
 
                        }
                }
@@ -1771,21 +1970,34 @@ int potential_force_calc(t_moldyn *moldyn) {
                        continue;
 
                /* copy the neighbour lists */
+#ifdef STATIC_LISTS
+               /* no copy needed for static lists */
+#else
                memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
+#endif
 
                /* second loop over atoms j */
                for(j=0;j<27;j++) {
 
+                       bc_ij=(j<dnlc)?0:1;
+#ifdef STATIC_LISTS
+                       p=0;
+
+                       while(neighbour_i[j][p]!=0) {
+
+                               jtom=&(atom[neighbour_i[j][p]]);
+                               p++;
+#else
                        this=&(neighbour_i[j]);
                        list_reset_f(this);
 
                        if(this->start==NULL)
                                continue;
 
-                       bc_ij=(j<dnlc)?0:1;
-
                        do {
+
                                jtom=this->current->data;
+#endif
 
                                if(jtom==&(itom[i]))
                                        continue;
@@ -1811,17 +2023,24 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                for(k=0;k<27;k++) {
 
+                                       bc_ik=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+                                       q=0;
+
+                                       while(neighbour_i[j][q]!=0) {
+
+                                               ktom=&(atom[neighbour_i[k][q]]);
+                                               q++;
+#else
                                        that=&(neighbour_i2[k]);
                                        list_reset_f(that);
                                        
                                        if(that->start==NULL)
                                                continue;
 
-                                       bc_ik=(k<dnlc)?0:1;
-
                                        do {
-
                                                ktom=that->current->data;
+#endif
 
                                                if(!(ktom->attr&ATOM_ATTR_3BP))
                                                        continue;
@@ -1837,9 +2056,12 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                                  jtom,
                                                                  ktom,
                                                                  bc_ik|bc_ij);
-
+#ifdef STATIC_LISTS
+                                       }
+#else
                                        } while(list_next_f(that)!=\
                                                L_NO_NEXT_ELEMENT);
+#endif
 
                                }
 
@@ -1856,17 +2078,24 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                for(k=0;k<27;k++) {
 
+                                       bc_ik=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+                                       q=0;
+
+                                       while(neighbour_i[j][q]!=0) {
+
+                                               ktom=&(atom[neighbour_i[k][q]]);
+                                               q++;
+#else
                                        that=&(neighbour_i2[k]);
                                        list_reset_f(that);
                                        
                                        if(that->start==NULL)
                                                continue;
 
-                                       bc_ik=(k<dnlc)?0:1;
-
                                        do {
-
                                                ktom=that->current->data;
+#endif
 
                                                if(!(ktom->attr&ATOM_ATTR_3BP))
                                                        continue;
@@ -1883,8 +2112,12 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                                  ktom,
                                                                  bc_ik|bc_ij);
 
+#ifdef STATIC_LISTS
+                                       }
+#else
                                        } while(list_next_f(that)!=\
                                                L_NO_NEXT_ELEMENT);
+#endif
 
                                }
                                
@@ -1896,8 +2129,11 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                          &(itom[i]),
                                                          jtom,bc_ij);
                                }
-                                       
+#ifdef STATIC_LISTS
+                       }
+#else
                        } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
                
                }
                
@@ -2116,7 +2352,12 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
        double *stat;
        int i,j;
        t_linkcell *lc;
+#ifdef STATIC_LISTS
+       int *neighbour[27];
+       int p;
+#else
        t_list neighbour[27];
+#endif
        t_atom *itom,*jtom;
        t_list *this;
        unsigned char bc;
@@ -2130,6 +2371,12 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
        slots=(int)(moldyn->cutoff/dr);
        o=2*slots;
 
+       printf("[moldyn] pair correlation calc info:\n");
+       printf("  time: %f\n",moldyn->time);
+       printf("  count: %d\n",moldyn->count);
+       printf("  cutoff: %f\n",moldyn->cutoff);
+       printf("  temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
+
        if(ptr!=NULL) {
                stat=(double *)ptr;
        }
@@ -2159,20 +2406,27 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
                ibrand=itom[i].brand;
        
                for(j=0;j<27;j++) {
-                       /* prepare the neighbour cell list */
+
+                       bc=(j<lc->dnlc)?0:1;
+
+#ifdef STATIC_LISTS
+                       p=0;
+
+                       while(neighbour[j][p]!=0) {
+
+                               jtom=&(moldyn->atom[neighbour[j][p]]);
+                               p++;
+#else
                        this=&(neighbour[j]);
                        list_reset_f(this);
 
-                       /* check for atoms */
                        if(this->start==NULL)
                                continue;
 
-                       /* boundary check */    
-                       bc=(j<lc->dnlc)?0:1;
-
                        do {
-                               jtom=this->current->data;
 
+                               jtom=this->current->data;
+#endif
 
                                if(jtom==&(itom[i]))
                                        continue;
@@ -2211,21 +2465,24 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
                                                stat[s+o]+=1;
                                }
 
+#ifdef STATIC_LISTS
+                       }
+#else
                        } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
                }
        }
 
-       /* normalization */
+       /* normalization 
        for(i=1;i<slots;i++) {
-               /*
-                * normalization: 4 pi r r dr
-                * here: not double counting pairs -> 2 pi r r dr
-                */
+                // normalization: 4 pi r r dr
+                // here: not double counting pairs -> 2 pi r r dr
                norm=2*M_PI*moldyn->count*(i*dr*i*dr)*dr;
                stat[i]/=norm;
                stat[slots+i]/=norm;
                stat[o+i]/=norm;
        }
+       */
 
        if(ptr==NULL) {
                /* todo: store/print pair correlation function */
@@ -2268,7 +2525,12 @@ int visual_atoms(t_moldyn *moldyn) {
        t_atom *atom;
        t_atom *btom;
        t_linkcell *lc;
+#ifdef STATIC_LISTS
+       int *neighbour[27];
+       int p;
+#else
        t_list neighbour[27];
+#endif
        u8 bc;
        t_3dvec dist;
        double d2;
@@ -2322,12 +2584,19 @@ int visual_atoms(t_moldyn *moldyn) {
                                          (atom[i].r.z+moldyn->dim.z/2)/lc->z,
                                          neighbour);
                for(j=0;j<27;j++) {
+                       bc=j<lc->dnlc?0:1;
+#ifdef STATIC_LISTS
+                       p=0;
+                       while(neighbour[j][p]!=0) {
+                               btom=&(atom[neighbour[j][p]]);
+                               p++;
+#else
                        list_reset_f(&neighbour[j]);
                        if(neighbour[j].start==NULL)
                                continue;
-                       bc=j<lc->dnlc?0:1;
                        do {
                                btom=neighbour[j].current->data;
+#endif
                                if(btom==&atom[i])      // skip identical atoms
                                        continue;
                                //if(btom<&atom[i])     // skip half of them
@@ -2347,7 +2616,11 @@ int visual_atoms(t_moldyn *moldyn) {
                                dprintf(fd,"# [B] %f %f %f %f %f %f\n",
                                        atom[i].r.x,atom[i].r.y,atom[i].r.z,
                                        btom->r.x,btom->r.y,btom->r.z);
+#ifdef STATIC_LISTS
+                       }
+#else
                        } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT);
+#endif
                }
        }
 

diff --git a/moldyn.h b/moldyn.h
index a7a81e9..ee45228 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -62,10 +62,16 @@ typedef struct s_linkcell {
        int cells;              /* total amount of cells */
        double len;             /* prefered cell edge length */
        double x,y,z;           /* the actual cell lengthes */
+#ifdef STATIC_LISTS
+       int **subcell;          /* pointer to the cell lists */
+#else
        t_list *subcell;        /* pointer to the cell lists */
+#endif
        int dnlc;               /* direct neighbour lists counter */
 } t_linkcell;
 
+#define MAX_ATOMS_PER_LIST     10
+
 /* moldyn schedule structure */
 typedef struct s_moldyn_schedule {
        int count;
@@ -454,7 +460,11 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist);
 
 int link_cell_init(t_moldyn *moldyn,u8 vol);
 int link_cell_update(t_moldyn *moldyn);
+#ifdef STATIC_LISTS
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int **cell);
+#else
 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell);
+#endif
 int link_cell_shutdown(t_moldyn *moldyn);
 
 typedef int (*set_hook)(void *,void *);