increased relaxation steps + introduced average reset + added critical

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 5681779..ec96712 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -954,12 +954,41 @@ 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;
+
        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 +1060,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;
@@ -1436,7 +1466,6 @@ int link_cell_shutdown(t_moldyn *moldyn) {
 
        for(i=0;i<lc->cells;i++) {
 #ifdef STATIC_LISTS
-               //printf(" ---> %d free %p\n",i,lc->subcell[i]);
                free(lc->subcell[i]);
 #else
                //printf(" ---> %d free %p\n",i,lc->subcell[i].start);
@@ -1650,13 +1679,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;
@@ -1670,8 +1700,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");
                }
@@ -1712,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);
 
@@ -2182,10 +2215,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);
@@ -2196,10 +2233,13 @@ 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 etc ...
@@ -2258,13 +2298,14 @@ 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;
 
        printf("[moldyn] pair correlation calc info:\n");
@@ -2323,12 +2364,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;
 
                                /*
@@ -2340,14 +2378,22 @@ 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)\n",
+                                              s,slots);
+                                       s=slots-1;
+                               }
+
                                if(ibrand!=jtom->brand) {
                                        /* mixed */
                                        stat[s]+=1;
@@ -2360,7 +2406,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