added pair correlation calc code

diff --git a/Makefile b/Makefile
index 0e7b122..73645df 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -17,7 +17,7 @@ DEPS = moldyn.o random/random.o list/list.o
 DEPS += potentials/lennard_jones.o potentials/harmonic_oscillator.o
 DEPS += potentials/tersoff.o potentials/albe.o
 
-all: posic sic fluctuation_calc postproc
+all: posic sic fluctuation_calc postproc pair_correlation_calc
 
 posic: $(DEPS)
 
@@ -25,6 +25,8 @@ sic: $(DEPS)
 
 postproc: $(DEPS)
 
+pair_correlation_calc: $(DEPS)
+
 .PHONY:clean
 clean:
        rm -vf sic postproc fluctuation_calc *.o */*.o

diff --git a/moldyn.c b/moldyn.c
index ad4f900..2e0e468 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -2057,6 +2057,13 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
        }
 
        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)");
@@ -2103,49 +2110,130 @@ int pair_correlation_init(t_moldyn *moldyn,double dr) {
        return 0;
 }
 
-int calculate_pair_correlation(t_moldyn *moldyn,double dr) {
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 
        int slots;
-       int *stat;
-       int i;
+       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;
 
-       stat=(int *)malloc(3*slots*sizeof(int));
-       if(stat==NULL) {
-               perror("[moldyn] pair correlation malloc");
-               return -1;
+       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);
-       
-       for(i=0;i<lc->cells;i++) {
-               /* check for atoms */
-               itom=lc->subcell[i].current->data;
-               if(itom==NULL)
-                       continue;
 
-               /* pick first atom and do neighbour indexing */
+       itom=moldyn->atom;
+       
+       for(i=0;i<moldyn->count;i++) {
+               /* neighbour indexing */
                link_cell_neighbour_index(moldyn,
-                                         (itom->r.x+moldyn->dim.x/2)/lc->x,
-                                         (itom->r.y+moldyn->dim.y/2)/lc->x,
-                                         (itom->r.z+moldyn->dim.z/2)/lc->x,
+                                         (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);
 
-               /* calculation of g(r) */
-               do {
-                       itom=lc->subcell[i].current->data;
-// GO ON HERE ...
-               } while(list_next_f(this)!=L_NO_NEXT_ELEMENT););
+               /* 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;

diff --git a/moldyn.h b/moldyn.h
index 2d4ca48..a7a81e9 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -475,8 +475,13 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a);
 
 int moldyn_bc_check(t_moldyn *moldyn);
 
+int moldyn_read_save_file(t_moldyn *moldyn,char *file);
+int moldyn_load(t_moldyn *moldyn);
 int get_line(int fd,char *line,int max);
 
+int pair_correlation_init(t_moldyn *moldyn,double dr);
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr);
+
 int visual_init(t_moldyn *moldyn,char *filebase);
 int visual_atoms(t_moldyn *moldyn);
 

diff --git a/pair_corr_calc_script b/pair_corr_calc_script
new file mode 100755 (executable)
index 0000000..27dc0ab
--- /dev/null
+++ b/pair_corr_calc_script
@@ -0,0 +1,29 @@
+#!/bin/sh
+
+#
+# calculate pair correlation functions for a whole simulation output
+# frank.zirkelbach@physik.uni-augsburg.de
+#
+
+do_it() {
+       echo "processing $1 ..."
+       ./pair_correlation_calc $1 $2
+       trgab=`echo $1 | sed 's%s-%pair_corr_ab-%' | sed 's%.save%%'`
+       trgaa=`echo $1 | sed 's%s-%pair_corr_aa-%' | sed 's%.save%%'`
+       trgbb=`echo $1 | sed 's%s-%pair_corr_bb-%' | sed 's%.save%%'`
+       mv pair_corr_func_ab.txt $trgab
+       mv pair_corr_func_aa.txt $trgaa
+       mv pair_corr_func_bb.txt $trgbb
+       echo "done"
+}
+
+if [ -d $1 ]; then
+       for file in $1/*.save; do
+               do_it $file $2
+       done
+fi
+
+if [ -f $1 ]; then
+       do_it $1 $2
+fi
+

diff --git a/pair_correlation_calc.c b/pair_correlation_calc.c
new file mode 100644 (file)
index 0000000..2f93ab7
--- /dev/null
+++ b/pair_correlation_calc.c
@@ -0,0 +1,83 @@
+/*
+ * calcultae pair correlation function
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#include "moldyn.h"
+
+int usage(char *prog) {
+
+       printf("\nusage:\n");
+       printf("  %s <save file> <dr>\n\n",prog);
+
+       return -1;
+}
+
+int main(int argc,char **argv) {
+
+       t_moldyn moldyn;
+       int ret;
+       double *stat;
+       int slots;
+       int i;
+       double dr;
+       int fd;
+
+       if(argc!=3) {
+               usage(argv[0]);
+               return -1;
+       }
+
+       ret=moldyn_read_save_file(&moldyn,argv[1]);
+       if(ret) {
+               printf("[pair corr calc] exit!\n");
+               return ret;
+       }
+
+       dr=atof(argv[2]);
+       slots=(int)(moldyn.cutoff/dr);
+
+       stat=(double *)malloc(3*slots*sizeof(double));
+       if(stat==NULL) {
+               perror("[pair corr calc] alloc mem");
+               return -1;
+       }
+
+       calculate_pair_correlation(&moldyn,dr,stat);
+
+       fd=open("pair_corr_func_ab.txt",
+               O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+       dprintf(fd,"# type a - type b bonds\n");
+       for(i=0;i<slots;i++)
+               dprintf(fd,"%f %f\n",i*dr,stat[i]);
+       close(fd);
+               
+       fd=open("pair_corr_func_aa.txt",
+               O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+       dprintf(fd,"# type a - type a bonds\n");
+       for(i=0;i<slots;i++)
+               dprintf(fd,"%f %f\n",i*dr,stat[slots+i]);
+       close(fd);
+               
+       fd=open("pair_corr_func_bb.txt",
+               O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+       dprintf(fd,"# type a - type b bonds\n");
+       for(i=0;i<slots;i++)
+               dprintf(fd,"%f %f\n",i*dr,stat[2*slots+i]);
+       close(fd);
+               
+       free(stat);
+
+       return 0;
+}