From 099c05fce45cacbbd7bbf6a7c5f1067e150d30ac Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 31 Jan 2008 23:56:26 +0100
Subject: [PATCH] added static lists (but: glibc double free prob (for both
 lists!))

---
 Makefile |   3 +-
 moldyn.c | 325 ++++++++++++++++++++++++++++++++++++++++++++++++++-----
 moldyn.h |  10 ++
 3 files changed, 311 insertions(+), 27 deletions(-)

diff --git a/Makefile b/Makefile
index 73645df..de46da4 100644
--- 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
--- 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
--- 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 *);
-- 
2.20.1