increased relaxation steps + introduced average reset + added critical

check in verlet routine + modified sic delta t action + cleanups

diff --git a/config.h b/config.h
index 8783dc6..f0919b7 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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 (file)
--- 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);