X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=94beb044eacc2211dc290e7e12729de4a7b18ec8;hb=2d95c5c11703731d2049b8a185f7cff5a6a5b329;hp=6d5f9e4ce39d5cf8124c615e99a330297803a71e;hpb=099c05fce45cacbbd7bbf6a7c5f1067e150d30ac;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 6d5f9e4..94beb04 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -595,7 +595,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 	int count;
 	
 	atom=moldyn->atom;
-	count=(moldyn->count)++;
+	count=(moldyn->count)++;	// asshole style!
 
 	ptr=realloc(atom,(count+1)*sizeof(t_atom));
 	if(!ptr) {
@@ -605,6 +605,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 	moldyn->atom=ptr;
 
 	atom=moldyn->atom;
+
+	/* initialize new atom */
+	memset(&(atom[count]),0,sizeof(t_atom));
 	atom[count].r=*r;
 	atom[count].v=*v;
 	atom[count].element=element;
@@ -954,12 +957,43 @@ double pressure_calc(t_moldyn *moldyn) {
 	return moldyn->p;
 }
 
+int average_reset(t_moldyn *moldyn) {
+
+	printf("[moldyn] average reset\n");
+
+	/* update skip value */
+	moldyn->avg_skip=moldyn->total_steps;
+
+	/* kinetic energy */
+	moldyn->k_sum=0.0;
+	moldyn->k2_sum=0.0;
+	
+	/* potential energy */
+	moldyn->v_sum=0.0;
+	moldyn->v2_sum=0.0;
+
+	/* temperature */
+	moldyn->t_sum=0.0;
+
+	/* virial */
+	moldyn->virial_sum=0.0;
+	moldyn->gv_sum=0.0;
+
+	/* pressure */
+	moldyn->p_sum=0.0;
+	moldyn->gp_sum=0.0;
+
+	return 0;
+}
+
 int average_and_fluctuation_calc(t_moldyn *moldyn) {
 
+	int denom;
+
 	if(moldyn->total_steps<moldyn->avg_skip)
 		return 0;
 
-	int denom=moldyn->total_steps+1-moldyn->avg_skip;
+	denom=moldyn->total_steps+1-moldyn->avg_skip;
 
 	/* assume up to date energies, temperature, pressure etc */
 
@@ -1031,7 +1065,8 @@ printf("  --> <dV2> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->coun
 
 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
-	t_3dvec dim,*tp;
+	t_3dvec dim;
+	//t_3dvec *tp;
 	double u_up,u_down,dv;
 	double scale,p;
 	t_atom *store;
@@ -1244,13 +1279,10 @@ 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);
 
@@ -1261,21 +1293,36 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 	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;
+
+#ifdef STATIC_LISTS
 	lc->subcell=malloc(lc->cells*sizeof(int*));
+#else
+	lc->subcell=malloc(lc->cells*sizeof(t_list));
+#endif
+
+	if(lc->subcell==NULL) {
+		perror("[moldyn] cell init (malloc)");
+		return -1;
+	}
 
 	if(lc->cells<27)
 		printf("[moldyn] FATAL: less then 27 subcells!\n");
 
 	if(vol) {
+#ifdef STATIC_LISTS
 		printf("[moldyn] initializing 'static' linked cells (%d)\n",
 		       lc->cells);
+#else
+		printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
+		       lc->cells);
+#endif
 		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);
 	}
 
+#ifdef STATIC_LISTS
 	/* list init */
 	for(i=0;i<lc->cells;i++) {
 		lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int));
@@ -1283,8 +1330,16 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 			perror("[moldyn] list init (malloc)");
 			return -1;
 		}
-//if(i==0) printf(" --- add one here! %d %p %p ----\n",i,lc->subcell,lc->subcell[0]);
+		/*
+		if(i==0)
+			printf(" ---> %d malloc %p (%p)\n",
+			       i,lc->subcell[0],lc->subcell);
+		*/
 	}
+#else
+	for(i=0;i<lc->cells;i++)
+		list_init_f(&(lc->subcell[i]));
+#endif
 
 	/* update the list */
 	link_cell_update(moldyn);
@@ -1298,7 +1353,9 @@ int link_cell_update(t_moldyn *moldyn) {
 	int nx,ny;
 	t_atom *atom;
 	t_linkcell *lc;
+#ifdef STATIC_LISTS
 	int p;
+#endif
 
 	atom=moldyn->atom;
 	lc=&(moldyn->lc);
@@ -1307,13 +1364,18 @@ int link_cell_update(t_moldyn *moldyn) {
 	ny=lc->ny;
 
 	for(i=0;i<lc->cells;i++)
+#ifdef STATIC_LISTS
 		memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+#else
+		list_destroy_f(&(lc->subcell[i]));
+#endif
 
 	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);
-
+	
+#ifdef STATIC_LISTS
 		p=0;
 		while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
 			p++;
@@ -1324,12 +1386,27 @@ int link_cell_update(t_moldyn *moldyn) {
 		}
 
 		lc->subcell[i+j*nx+k*nx*ny][p]=count;
+#else
+		list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
+		                     &(atom[count]));
+		/*
+		if(j==0&&k==0)
+			printf(" ---> %d %d malloc %p (%p)\n",
+			       i,count,lc->subcell[i].current,lc->subcell);
+		*/
+#endif
 	}
 
 	return 0;
 }
 
-int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int **cell) {
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
+#ifdef STATIC_LISTS
+                              int **cell
+#else
+                              t_list *cell
+#endif
+			     ) {
 
 	t_linkcell *lc;
 	int a;
@@ -1347,6 +1424,10 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int **cell) {
 	count2=27;
 	a=nx*ny;
 
+	if(i>=nx||j>=ny||k>=nz)
+		printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n",
+		       i,nx,j,ny,k,nz);
+
 	cell[0]=lc->subcell[i+j*nx+k*a];
 	for(ci=-1;ci<=1;ci++) {
 		bx=0;
@@ -1389,170 +1470,23 @@ 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);
+	for(i=0;i<lc->cells;i++) {
+#ifdef STATIC_LISTS
 		free(lc->subcell[i]);
-}
-
-	free(lc->subcell);
-
-	return 0;
-}
-
 #else
-
-int link_cell_init(t_moldyn *moldyn,u8 vol) {
-
-	t_linkcell *lc;
-	int i;
-
-	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(t_list));
-
-	if(lc->cells<27)
-		printf("[moldyn] FATAL: less then 27 subcells!\n");
-
-	if(vol) {
-		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);
-	}
-
-	for(i=0;i<lc->cells;i++)
-		list_init_f(&(lc->subcell[i]));
-
- 	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;
-	double x,y,z;
-
-	atom=moldyn->atom;
-	lc=&(moldyn->lc);
-
-	nx=lc->nx;
-	ny=lc->ny;
-
-	x=moldyn->dim.x/2;
-	y=moldyn->dim.y/2;
-	z=moldyn->dim.z/2;
-
-	for(i=0;i<lc->cells;i++)
+		//printf(" ---> %d free %p\n",i,lc->subcell[i].start);
 		list_destroy_f(&(lc->subcell[i]));
-	
-	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);
-		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;
-}
-
-int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *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];
-				}
-			}
-		}
+#endif
 	}
 
-	lc->dnlc=count1;
-
-	return count1;
-}
-
-int link_cell_shutdown(t_moldyn *moldyn) {
-
-	int i;
-	t_linkcell *lc;
-
-	lc=&(moldyn->lc);
-
-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;
@@ -1644,14 +1578,14 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
 	/* some stupid checks before we actually start calculating bullshit */
 	if(moldyn->cutoff>0.5*moldyn->dim.x)
-		printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
+		printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n");
 	if(moldyn->cutoff>0.5*moldyn->dim.y)
-		printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
+		printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
 	if(moldyn->cutoff>0.5*moldyn->dim.z)
-		printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
+		printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
 	ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
 	if(ds>0.05*moldyn->nnd)
-		printf("[moldyn] warning: forces too high / tau too small!\n");
+		printf("[moldyn] WARNING: forces too high / tau too small!\n");
 
 	/* zero absolute time */
 	moldyn->time=0.0;
@@ -1754,13 +1688,14 @@ int moldyn_integrate(t_moldyn *moldyn) {
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
-	printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
-	       sched->count,i,
-	       moldyn->t,moldyn->t_avg,
-	       moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
-	       moldyn->volume,
-	       (int)(t2.tv_sec-t1.tv_sec));
-	fflush(stdout);
+printf("\rsched:%d, steps:%d/%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
+       sched->count,i,moldyn->total_steps,
+       moldyn->t,moldyn->t_avg,
+       moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
+       moldyn->volume,
+       (int)(t2.tv_sec-t1.tv_sec));
+
+			fflush(stdout);
 
 			/* copy over time */
 			t1=t2;
@@ -1774,8 +1709,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
 		/* check for hooks */
 		if(sched->hook) {
-			printf("\n ## schedule hook %d/%d start ##\n",
-			       sched->count+1,sched->total_sched-1);
+			printf("\n ## schedule hook %d start ##\n",
+			       sched->count);
 			sched->hook(moldyn,sched->hook_params);
 			printf(" ## schedule hook end ##\n");
 		}
@@ -1816,6 +1751,9 @@ int velocity_verlet(t_moldyn *moldyn) {
 		v3_add(&(atom[i].v),&(atom[i].v),&delta);
 	}
 
+	/* criticial check */
+	moldyn_bc_check(moldyn);
+
 	/* neighbour list update */
 	link_cell_update(moldyn);
 
@@ -2156,14 +2094,21 @@ int potential_force_calc(t_moldyn *moldyn) {
 	}
 #endif
 
-	/* calculate global virial */
+	/* some postprocessing */
 	for(i=0;i<count;i++) {
-		moldyn->gvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
-		moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
-		moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
-		moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
-		moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
-		moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
+		/* calculate global virial */
+		moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
+		moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
+		moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
+		moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
+		moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
+		moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
+
+		/* check forces regarding the given timestep */
+		if(v3_norm(&(itom[i].f))>\
+		   0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
+			printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
+			       i);
 	}
 
 	return 0;
@@ -2286,10 +2231,14 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 	}
 
 	size=sizeof(t_moldyn);
-	cnt=read(fd,moldyn,size);
-	if(cnt!=size) {
-		perror("[moldyn] load save file read (moldyn)");
-		return cnt;
+
+	while(size) {
+		cnt=read(fd,moldyn,size);
+		if(cnt<0) {
+			perror("[moldyn] load save file read (moldyn)");
+			return cnt;
+		}
+		size-=cnt;
 	}
 
 	size=moldyn->count*sizeof(t_atom);
@@ -2300,13 +2249,16 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 		return -1;
 	}
 
-	cnt=read(fd,moldyn->atom,size);
-	if(cnt!=size) {
-		perror("[moldyn] load save file read (atoms)");
-		return cnt;
+	while(size) {
+		cnt=read(fd,moldyn->atom,size);
+		if(cnt<0) {
+			perror("[moldyn] load save file read (atoms)");
+			return cnt;
+		}
+		size-=cnt;
 	}
 
-	// hooks
+	// hooks etc ...
 
 	return 0;
 }
@@ -2362,15 +2314,19 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 	t_list *this;
 	unsigned char bc;
 	t_3dvec dist;
-	double d,norm;
+	double d;
+	//double norm;
 	int o,s;
 	unsigned char ibrand;
 
 	lc=&(moldyn->lc);
 
-	slots=(int)(moldyn->cutoff/dr);
+	slots=moldyn->cutoff/dr;
 	o=2*slots;
 
+	if(slots*dr<=moldyn->cutoff)
+		printf("[moldyn] WARNING: pcc (low #slots)\n");
+
 	printf("[moldyn] pair correlation calc info:\n");
 	printf("  time: %f\n",moldyn->time);
 	printf("  count: %d\n",moldyn->count);
@@ -2427,12 +2383,9 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 
 				jtom=this->current->data;
 #endif
-
-				if(jtom==&(itom[i]))
-					continue;
-
-				/* only count pairs once */
-				if(itom[i].tag>jtom->tag)
+				/* only count pairs once,
+				 * skip same atoms */
+				if(itom[i].tag>=jtom->tag)
 					continue;
 
 				/*
@@ -2444,14 +2397,23 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 				if(bc) check_per_bound(moldyn,&dist);
 				d=v3_absolute_square(&dist);
 
-				/* ignore if greater cutoff */
-				if(d>moldyn->cutoff_square)
+				/* ignore if greater or equal cutoff */
+				if(d>=moldyn->cutoff_square)
 					continue;
 
 				/* fill the slots */
 				d=sqrt(d);
 				s=(int)(d/dr);
 
+				/* should never happen but it does 8) -
+				 * related to -ffloat-store problem! */
+				if(s>=slots) {
+					printf("[moldyn] WARNING: pcc (%d/%d)",
+					       s,slots);
+					printf("\n");
+					s=slots-1;
+				}
+
 				if(ibrand!=jtom->brand) {
 					/* mixed */
 					stat[s]+=1;
@@ -2464,7 +2426,6 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 					/* type b - type b bonds */
 						stat[s+o]+=1;
 				}
-
 #ifdef STATIC_LISTS
 			}
 #else