float-store bugfix (unstatisfying!) + slightly new design of sic.c

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 5681779..edda007 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -956,10 +956,12 @@ double pressure_calc(t_moldyn *moldyn) {
 
 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 +1033,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 +1439,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 +1652,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 +1673,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");
                }
@@ -2182,10 +2185,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 +2203,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 +2268,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 +2334,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 +2348,18 @@ 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) s=slots-1;
+
                                if(ibrand!=jtom->brand) {
                                        /* mixed */
                                        stat[s]+=1;
@@ -2360,7 +2372,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