From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Fri, 8 Feb 2008 00:30:25 +0000 (+0100)
Subject: fixed missing init for new atoms (especially force!) + increased insert
X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b3f72740276662c09ff43c294b04409a5bb2e013;p=physik%2Fposic.git

fixed missing init for new atoms (especially force!) + increased insert
relaxation time + added some more critical debugging checks
---

diff --git a/config.h b/config.h
index f0919b7..f43f006 100644
--- a/config.h
+++ b/config.h
@@ -41,7 +41,7 @@
 
 #define INS_R_C		1.5
 #define INS_DELTA_TC	1.0
-#define INS_RELAX	50
+#define INS_RELAX	100
 #define INS_TAU		1.0
 
 // postrun
diff --git a/moldyn.c b/moldyn.c
index ec96712..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;
@@ -956,6 +959,8 @@ double pressure_calc(t_moldyn *moldyn) {
 
 int average_reset(t_moldyn *moldyn) {
 
+	printf("[moldyn] average reset\n");
+
 	/* update skip value */
 	moldyn->avg_skip=moldyn->total_steps;
 
@@ -1369,7 +1374,7 @@ int link_cell_update(t_moldyn *moldyn) {
 		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)
@@ -1419,6 +1424,10 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
 	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;
@@ -1569,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;
@@ -2085,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;
@@ -2308,6 +2324,9 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 	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);
@@ -2389,8 +2408,9 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 				/* should never happen but it does 8) -
 				 * related to -ffloat-store problem! */
 				if(s>=slots) {
-					printf("[moldyn] WARNING pcc (%d/%d)\n",
+					printf("[moldyn] WARNING: pcc (%d/%d)",
 					       s,slots);
+					printf("\n");
 					s=slots-1;
 				}
 
diff --git a/sic.c b/sic.c
index a7613e2..fb20a56 100644
--- a/sic.c
+++ b/sic.c
@@ -122,6 +122,9 @@ int sic_hook(void *moldyn,void *hook_params) {
 	hp=hook_params;
 	md=moldyn;
 
+	tau=1.0;
+	steps=0;
+
 	/* switch on t scaling */
 	if(md->schedule.count==0)
 		set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0);
@@ -438,6 +441,9 @@ int main(int argc,char **argv) {
 	moldyn_set_log(&md,SAVE_STEP,LOG_S);
 	moldyn_set_log(&md,CREATE_REPORT,0);
 
+	/* next neighbour distance for critical checking */
+	set_nn_dist(&md,0.25*ALBE_LC_SI*sqrt(3.0));
+
 	/*
 	 * let's do the actual md algorithm now
 	 *