added pair correlation calc code

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 258985e..2e0e468 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -1469,7 +1469,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
        /* calculate initial forces */
        potential_force_calc(moldyn);
 #ifdef DEBUG
-return 0;
+//return 0;
 #endif
 
        /* some stupid checks before we actually start calculating bullshit */
@@ -1580,7 +1580,7 @@ return 0;
                }
 
                /* display progress */
-               if(!(moldyn->total_steps%10)) {
+               //if(!(moldyn->total_steps%10)) {
                        /* get current time */
                        gettimeofday(&t2,NULL);
 
@@ -1594,7 +1594,7 @@ return 0;
 
                        /* copy over time */
                        t1=t2;
-               }
+               //}
 
                /* increase absolute time */
                moldyn->time+=moldyn->tau;
@@ -1911,7 +1911,13 @@ int potential_force_calc(t_moldyn *moldyn) {
        }
 
 #ifdef DEBUG
-       printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
+       //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
+       if(moldyn->time>DSTART&&moldyn->time<DEND) {
+               printf("force:\n");
+               printf("  x: %0.40f\n",moldyn->atom[5832].f.x);
+               printf("  y: %0.40f\n",moldyn->atom[5832].f.y);
+               printf("  z: %0.40f\n",moldyn->atom[5832].f.z);
+       }
 #endif
 
        /* calculate global virial */
@@ -2032,6 +2038,43 @@ int moldyn_bc_check(t_moldyn *moldyn) {
  * restore function
  */
 
+int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
+
+       int fd;
+       int cnt,size;
+
+       fd=open(file,O_RDONLY);
+       if(fd<0) {
+               perror("[moldyn] load save file open");
+               return fd;
+       }
+
+       size=sizeof(t_moldyn);
+       cnt=read(fd,moldyn,size);
+       if(cnt!=size) {
+               perror("[moldyn] load save file read (moldyn)");
+               return cnt;
+       }
+
+       size=moldyn->count*sizeof(t_atom);
+
+       moldyn->atom=(t_atom *)malloc(size);
+       if(moldyn->atom==NULL) {
+               perror("[moldyn] load save file malloc (atoms)");
+               return -1;
+       }
+
+       cnt=read(fd,moldyn->atom,size);
+       if(cnt!=size) {
+               perror("[moldyn] load save file read (atoms)");
+               return cnt;
+       }
+
+       // hooks
+
+       return 0;
+}
+
 int moldyn_load(t_moldyn *moldyn) {
 
        // later ...
@@ -2061,6 +2104,141 @@ int get_line(int fd,char *line,int max) {
        }
 }
 
+int pair_correlation_init(t_moldyn *moldyn,double dr) {
+
+       
+       return 0;
+}
+
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+
+       int slots;
+       double *stat;
+       int i,j;
+       t_linkcell *lc;
+       t_list neighbour[27];
+       t_atom *itom,*jtom;
+       t_list *this;
+       unsigned char bc;
+       t_3dvec dist;
+       double d,norm;
+       int o,s;
+       unsigned char ibrand;
+
+       lc=&(moldyn->lc);
+
+       slots=(int)(moldyn->cutoff/dr);
+       o=2*slots;
+
+       if(ptr!=NULL) {
+               stat=(double *)ptr;
+       }
+       else {
+               stat=(double *)malloc(3*slots*sizeof(double));
+               if(stat==NULL) {
+                       perror("[moldyn] pair correlation malloc");
+                       return -1;
+               }
+       }
+
+       memset(stat,0,3*slots*sizeof(double));
+
+       link_cell_init(moldyn,VERBOSE);
+
+       itom=moldyn->atom;
+       
+       for(i=0;i<moldyn->count;i++) {
+               /* neighbour indexing */
+               link_cell_neighbour_index(moldyn,
+                                         (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+                                         (itom[i].r.y+moldyn->dim.y/2)/lc->x,
+                                         (itom[i].r.z+moldyn->dim.z/2)/lc->x,
+                                         neighbour);
+
+               /* brand of atom i */
+               ibrand=itom[i].brand;
+       
+               for(j=0;j<27;j++) {
+                       /* prepare the neighbour cell list */
+                       this=&(neighbour[j]);
+                       list_reset_f(this);
+
+                       /* check for atoms */
+                       if(this->start==NULL)
+                               continue;
+
+                       /* boundary check */    
+                       bc=(j<lc->dnlc)?0:1;
+
+                       do {
+                               jtom=this->current->data;
+
+
+                               if(jtom==&(itom[i]))
+                                       continue;
+
+                               /* only count pairs once */
+                               if(itom[i].tag>jtom->tag)
+                                       continue;
+
+                               /*
+                                * pair correlation calc
+                                */
+
+                               /* distance */
+                               v3_sub(&dist,&(jtom->r),&(itom[i].r));
+                               if(bc) check_per_bound(moldyn,&dist);
+                               d=v3_absolute_square(&dist);
+
+                               /* ignore if greater cutoff */
+                               if(d>moldyn->cutoff_square)
+                                       continue;
+
+                               /* fill the slots */
+                               d=sqrt(d);
+                               s=(int)(d/dr);
+
+                               if(ibrand!=jtom->brand) {
+                                       /* mixed */
+                                       stat[s]+=1;
+                               }
+                               else {
+                                       /* type a - type a bonds */
+                                       if(ibrand==0)
+                                               stat[s+slots]+=1;
+                                       else
+                                       /* type b - type b bonds */
+                                               stat[s+o]+=1;
+                               }
+
+                       } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+               }
+       }
+
+       /* normalization */
+       for(i=1;i<slots;i++) {
+               /*
+                * normalization: 4 pi r r dr
+                * here: not double counting pairs -> 2 pi r r dr
+                */
+               norm=2*M_PI*moldyn->count*(i*dr*i*dr)*dr;
+               stat[i]/=norm;
+               stat[slots+i]/=norm;
+               stat[o+i]/=norm;
+       }
+
+       if(ptr==NULL) {
+               /* todo: store/print pair correlation function */
+               free(stat);
+       }
+
+       free(moldyn->atom);
+
+       link_cell_shutdown(moldyn);
+
+       return 0;
+}
+
 int analyze_bonds(t_moldyn *moldyn) {