From 6a467846e1bd9212bcfd4c5ac29f8954305aaff5 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Wed, 6 Feb 2008 18:45:36 +0100
Subject: [PATCH] increased relaxation steps + introduced average reset + added
 critical check in verlet routine + modified sic delta t action + cleanups

---
 config.h |  4 ++--
 moldyn.c | 36 +++++++++++++++++++++++++++++++++++-
 moldyn.h |  3 ++-
 sic.c    | 36 ++++++++++++++++++++++++++----------
 4 files changed, 65 insertions(+), 14 deletions(-)

diff --git a/config.h b/config.h
index 8783dc6..f0919b7 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	10
+#define INS_RELAX	50
 #define INS_TAU		1.0
 
 // postrun
@@ -49,7 +49,7 @@
 #define POST_RUNS	430
 #define POST_DELTA_TC	1.0
 #define POST_DT		1.0
-#define POST_RELAX	10
+#define POST_RELAX	50
 #define POST_TAU	1.0
 
 // log
diff --git a/moldyn.c b/moldyn.c
index edda007..ec96712 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -954,6 +954,33 @@ double pressure_calc(t_moldyn *moldyn) {
 	return moldyn->p;
 }
 
+int average_reset(t_moldyn *moldyn) {
+
+	/* 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;
@@ -1715,6 +1742,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);
 
@@ -2358,7 +2388,11 @@ 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) s=slots-1;
+				if(s>=slots) {
+					printf("[moldyn] WARNING pcc (%d/%d)\n",
+					       s,slots);
+					s=slots-1;
+				}
 
 				if(ibrand!=jtom->brand) {
 					/* mixed */
diff --git a/moldyn.h b/moldyn.h
index 2290743..5a10197 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -444,7 +444,8 @@ double get_temperature(t_moldyn *moldyn);
 int scale_velocity(t_moldyn *moldyn,u8 equi_init);
 double virial_sum(t_moldyn *moldyn);
 double pressure_calc(t_moldyn *moldyn);
-int energy_fluctuation_calc(t_moldyn *moldyn);
+int average_reset(t_moldyn *moldyn);
+int average_and_fluctuation_calc(t_moldyn *moldyn);
 int get_heat_capacity(t_moldyn *moldyn);
 double thermodynamic_pressure_calc(t_moldyn *moldyn);
 double get_pressure(t_moldyn *moldyn);
diff --git a/sic.c b/sic.c
index e1103bd..a7613e2 100644
--- a/sic.c
+++ b/sic.c
@@ -40,7 +40,7 @@ int insert_atoms(t_moldyn *moldyn) {
 	int i,j;
 	u8 run;
 	t_3dvec r,v,dist;
-	double d;
+	double d,dmin;
 
 	t_atom *atom;
 
@@ -83,6 +83,7 @@ int insert_atoms(t_moldyn *moldyn) {
 			r.z+=INS_OFFSET;
 			/* assume valid coordinates */
 			run=0;
+			dmin=10000000000.0;		// for sure too high!
 			for(i=0;i<moldyn->count;i++) {
 				atom=&(moldyn->atom[i]);
 				v3_sub(&dist,&(atom->r),&r);
@@ -94,6 +95,8 @@ int insert_atoms(t_moldyn *moldyn) {
 					run=1;
 					break;
 				}
+				if(d<dmin)
+					dmin=d;
 			}
 		}
 		add_atom(moldyn,INS_TYPE,INS_MASS,INS_BRAND,
@@ -101,6 +104,8 @@ int insert_atoms(t_moldyn *moldyn) {
 			 //ATOM_ATTR_HB|ATOM_ATTR_VB,
 			 ATOM_ATTR_HB,
 		         &r,&v);
+		printf(" %02d: atom %d | %f %f %f | %f\n",
+		       j,moldyn->count-1,r.x,r.y,r.z,dmin);
 	}
 
 	return 0;
@@ -112,6 +117,7 @@ int sic_hook(void *moldyn,void *hook_params) {
 	t_moldyn *md;
 	int steps;
 	double tau;
+	double dt;
 
 	hp=hook_params;
 	md=moldyn;
@@ -129,12 +135,15 @@ int sic_hook(void *moldyn,void *hook_params) {
 	/* act according to state */
 	switch(hp->state) {
 		case STATE_INSERT:
+			/* assigne values */
+			steps=INS_RELAX;
+			tau=INS_TAU;
 			/* check temperature */
-			if(md->t_avg-md->t_ref>INS_DELTA_TC) {
-				steps=INS_RELAX;
-				tau=INS_TAU;
+			dt=md->t_avg-md->t_ref;
+			if(dt<0)
+				dt=-dt;
+			if(dt>INS_DELTA_TC)
 				break;
-			}
 			/* insert atoms */
 			hp->insert_count+=1;
 			printf("   ### insert atoms (%d/%d) ###\n",
@@ -145,11 +154,15 @@ int sic_hook(void *moldyn,void *hook_params) {
 				hp->state=STATE_POSTRUN;
 			break;
 		case STATE_POSTRUN:
-			/* settings */
-			if(md->t-md->t_ref>POST_DELTA_TC) {
-				steps=POST_RELAX;
-				tau=POST_TAU;
-			}
+			/* assigne values */
+			steps=POST_RELAX;
+			tau=POST_TAU;
+			/* check temperature */
+			dt=md->t_avg-md->t_ref;
+			if(dt<0)
+				dt=-dt;
+			if(dt>INS_DELTA_TC)
+				break;
 			/* decrease temperature */
 			hp->postrun_count+=1;
 			printf(" ### postrun (%d/%d) ###\n",
@@ -163,6 +176,9 @@ int sic_hook(void *moldyn,void *hook_params) {
 			break;
 	}
 
+	/* reset the average counters */
+	average_reset(md);
+
 	/* add schedule */
 	moldyn_add_schedule(md,steps,tau);
 
-- 
2.20.1