introduce a 2body bond function with callback, modified pair corr calc

approp + introduced minimal bond analyze code

diff --git a/Makefile b/Makefile
index 1b875fe..c839759 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -22,7 +22,10 @@ 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 pair_correlation_calc diffusion_calc
+ALL = posic sic fluctuation_calc postproc pair_correlation_calc diffusion_calc
+ALL += bond_analyze
+
+all: $(ALL)
 
 posic: $(DEPS)
 
@@ -34,6 +37,8 @@ pair_correlation_calc: $(DEPS)
 
 diffusion_calc: $(DEPS)
 
+bond_analyze: $(DEPS)
+
 .PHONY:clean
 clean:
-       rm -vf sic postproc fluctuation_calc pair_correlation_calc *.o */*.o
+       rm -vf $(ALL) *.o */*.o

diff --git a/bond_analyze.c b/bond_analyze.c
new file mode 100644 (file)
index 0000000..0feffef
--- /dev/null
+++ b/bond_analyze.c
@@ -0,0 +1,54 @@
+/*
+ * bonding analyzation code
+ *
+ * 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>\n\n",prog);
+
+       return -1;
+}
+
+int main(int argc,char **argv) {
+
+       t_moldyn moldyn;
+       int ret;
+       double quality;
+
+       if(argc!=2) {
+               usage(argv[0]);
+               return -1;
+       }
+
+       memset(&moldyn,0,sizeof(t_moldyn));
+
+       printf("[bond analyze] reading save file ...\n");
+       ret=moldyn_read_save_file(&moldyn,argv[1]);
+       if(ret) {
+               printf("[bond analyze] exit!\n");
+               return ret;
+       }
+
+       bond_analyze(&moldyn,&quality);
+
+       printf("[bond analyze] quality = %f\n",quality);
+
+       moldyn_free_save_file(&moldyn);
+
+       return 0;
+}

diff --git a/moldyn.c b/moldyn.c
index 91fead7..84cf700 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -2323,6 +2323,13 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
        return 0;
 }
 
+int moldyn_free_save_file(t_moldyn *moldyn) {
+
+       free(moldyn->atom);
+
+       return 0;
+}
+
 int moldyn_load(t_moldyn *moldyn) {
 
        // later ...
@@ -2330,6 +2337,78 @@ int moldyn_load(t_moldyn *moldyn) {
        return 0;
 }
 
+/*
+ * function to find/callback all combinations of 2 body bonds
+ */
+
+int process_2b_bonds(t_moldyn *moldyn,void *data,
+                     int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+                                   void *data,u8 bc)) {
+
+       t_linkcell *lc;
+#ifdef STATIC_LISTS
+       int *neighbour[27];
+       int p;
+#else
+       t_list neighbour[27];
+#endif
+       u8 bc;
+       t_atom *itom,*jtom;
+       int i,j;
+       t_list *this;
+
+       lc=&(moldyn->lc);
+
+       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);
+
+               for(j=0;j<27;j++) {
+
+                       bc=(j<lc->dnlc)?0:1;
+
+#ifdef STATIC_LISTS
+                       p=0;
+
+                       while(neighbour[j][p]!=0) {
+
+                               jtom=&(moldyn->atom[neighbour[j][p]]);
+                               p++;
+#else
+                       this=&(neighbour[j]);
+                       list_reset_f(this);
+
+                       if(this->start==NULL)
+                               continue;
+
+                       do {
+
+                               jtom=this->current->data;
+#endif
+
+                               /* process bond */
+                               process(moldyn,&(itom[i]),jtom,data,bc);
+
+#ifdef STATIC_LISTS
+                       }
+#else
+                       } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
+               }
+       }
+
+       return 0;
+
+}
+
 /*
  * post processing functions
  */
@@ -2399,33 +2478,80 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
        return 0;
 }
 
-int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+int bonding_analyze(t_moldyn *moldyn,double *cnt) {
+
+       return 0;
+}
+
+int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
+                                       t_atom *jtom,void *data,u8 bc) {
 
-       int slots;
-       double *stat;
-       int i,j;
-       t_linkcell *lc;
-#ifdef STATIC_LISTS
-       int *neighbour[27];
-       int p;
-#else
-       t_list neighbour[27];
-#endif
-       t_atom *itom,*jtom;
-       t_list *this;
-       unsigned char bc;
        t_3dvec dist;
        double d;
-       double norm;
-       int o,s;
-       unsigned char ibrand;
+       int s;
+       t_pcc *pcc;
 
-       lc=&(moldyn->lc);
+       /* only count pairs once,
+        * skip same atoms */
+       if(itom->tag>=jtom->tag)
+               return 0;
+
+       /*
+        * pair correlation calc
+        */
+
+       /* get pcc data */
+       pcc=data;
+
+       /* distance */
+       v3_sub(&dist,&(jtom->r),&(itom->r));
+       if(bc) check_per_bound(moldyn,&dist);
+       d=v3_absolute_square(&dist);
+
+       /* ignore if greater or equal 2 times cutoff */
+       if(d>=4.0*moldyn->cutoff_square)
+               return 0;
+
+       /* fill the slots */
+       d=sqrt(d);
+       s=(int)(d/pcc->dr);
+
+       /* should never happen but it does 8) -
+        * related to -ffloat-store problem! */
+       if(s>=pcc->o1) {
+               printf("[moldyn] WARNING: pcc (%d/%d)",
+                      s,pcc->o1);
+               printf("\n");
+               s=pcc->o1-1;
+       }
+
+       if(itom->brand!=jtom->brand) {
+               /* mixed */
+               pcc->stat[s]+=1;
+       }
+       else {
+               /* type a - type a bonds */
+               if(itom->brand==0)
+                       pcc->stat[s+pcc->o1]+=1;
+               else
+               /* type b - type b bonds */
+                       pcc->stat[s+pcc->o2]+=1;
+       }
 
-       slots=2.0*moldyn->cutoff/dr;
-       o=2*slots;
+       return 0;
+}
 
-       if(slots*dr<=moldyn->cutoff)
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+
+       t_pcc pcc;
+       double norm;
+       int i;
+
+       pcc.dr=dr;
+       pcc.o1=2.0*moldyn->cutoff/dr;
+       pcc.o2=2*pcc.o1;
+
+       if(pcc.o1*dr<=moldyn->cutoff)
                printf("[moldyn] WARNING: pcc (low #slots)\n");
 
        printf("[moldyn] pair correlation calc info:\n");
@@ -2435,134 +2561,129 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
        printf("  temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
 
        if(ptr!=NULL) {
-               stat=(double *)ptr;
+               pcc.stat=(double *)ptr;
        }
        else {
-               stat=(double *)malloc(3*slots*sizeof(double));
-               if(stat==NULL) {
+               pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double));
+               if(pcc.stat==NULL) {
                        perror("[moldyn] pair correlation malloc");
                        return -1;
                }
        }
 
-       memset(stat,0,3*slots*sizeof(double));
+       memset(pcc.stat,0,3*pcc.o1*sizeof(double));
 
-       link_cell_init(moldyn,VERBOSE);
+       /* process */
+       process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process);
 
-       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);
+       /* normalization */
+       for(i=1;i<pcc.o1;i++) {
+                // normalization: 4 pi r^2 dr
+                // here: not double counting pairs -> 2 pi r r dr
+                // ... and actually it's a constant times r^2
+               norm=i*i*dr*dr;
+               pcc.stat[i]/=norm;
+               pcc.stat[pcc.o1+i]/=norm;
+               pcc.stat[pcc.o2+i]/=norm;
+       }
+       /* */
 
-               /* brand of atom i */
-               ibrand=itom[i].brand;
-       
-               for(j=0;j<27;j++) {
+       if(ptr==NULL) {
+               /* todo: store/print pair correlation function */
+               free(pcc.stat);
+       }
 
-                       bc=(j<lc->dnlc)?0:1;
+       return 0;
+}
 
-#ifdef STATIC_LISTS
-                       p=0;
+int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+                         void *data,u8 bc) {
 
-                       while(neighbour[j][p]!=0) {
+       t_ba *ba;
+       t_3dvec dist;
+       double d;
 
-                               jtom=&(moldyn->atom[neighbour[j][p]]);
-                               p++;
-#else
-                       this=&(neighbour[j]);
-                       list_reset_f(this);
+       if(itom->tag>=jtom->tag)
+               return 0;
 
-                       if(this->start==NULL)
-                               continue;
+       /* distance */
+       v3_sub(&dist,&(jtom->r),&(itom->r));
+       if(bc) check_per_bound(moldyn,&dist);
+       d=v3_absolute_square(&dist);
 
-                       do {
+       /* ignore if greater or equal cutoff */
+       if(d>=moldyn->cutoff_square)
+               return 0;
 
-                               jtom=this->current->data;
-#endif
-                               /* only count pairs once,
-                                * skip same atoms */
-                               if(itom[i].tag>=jtom->tag)
-                                       continue;
+       /* now count this bonding ... */
+       ba=data;
 
-                               /*
-                                * pair correlation calc
-                                */
+       /* increase total bond counter
+        * ... double counting!
+        */
+       ba->tcnt+=2;
 
-                               /* distance */
-                               v3_sub(&dist,&(jtom->r),&(itom[i].r));
-                               if(bc) check_per_bound(moldyn,&dist);
-                               d=v3_absolute_square(&dist);
+       if(itom->brand==0)
+               ba->acnt[jtom->tag]+=1;
+       else
+               ba->bcnt[jtom->tag]+=1;
+       
+       if(jtom->brand==0)
+               ba->acnt[itom->tag]+=1;
+       else
+               ba->bcnt[itom->tag]+=1;
 
-                               /* ignore if greater or equal cutoff */
-                               if(d>=4.0*moldyn->cutoff_square)
-                                       continue;
+       return 0;
+}
 
-                               /* 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)",
-                                              s,slots);
-                                       printf("\n");
-                                       s=slots-1;
-                               }
+int bond_analyze(t_moldyn *moldyn,double *quality) {
 
-                               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;
-                               }
-#ifdef STATIC_LISTS
-                       }
-#else
-                       } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
-#endif
-               }
-       }
+       // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
 
-       /* normalization */
-       for(i=1;i<slots;i++) {
-                // normalization: 4 pi r^2 dr
-                // here: not double counting pairs -> 2 pi r r dr
-                // ... and actually it's a constant times r^2
-               norm=i*i*dr*dr;
-               stat[i]/=norm;
-               stat[slots+i]/=norm;
-               stat[o+i]/=norm;
+       int qcnt;
+       t_ba ba;
+       int i;
+       t_atom *atom;
+
+       ba.acnt=malloc(moldyn->count*sizeof(int));
+       if(ba.acnt==NULL) {
+               perror("[moldyn] bond analyze malloc (a)");
+               return -1;
        }
-       /* */
+       memset(ba.acnt,0,moldyn->count*sizeof(int));
 
-       if(ptr==NULL) {
-               /* todo: store/print pair correlation function */
-               free(stat);
+       ba.bcnt=malloc(moldyn->count*sizeof(int));
+       if(ba.bcnt==NULL) {
+               perror("[moldyn] bond analyze malloc (b)");
+               return -1;
        }
+       memset(ba.bcnt,0,moldyn->count*sizeof(int));
 
-       free(moldyn->atom);
+       ba.tcnt=0;
 
-       link_cell_shutdown(moldyn);
+       qcnt=0;
 
-       return 0;
-}
+       atom=moldyn->atom;
 
-int analyze_bonds(t_moldyn *moldyn) {
+       process_2b_bonds(moldyn,&ba,bond_analyze_process);
 
-       
-       
+       for(i=0;i<moldyn->count;i++) {
+               if(atom[i].brand==0) {
+                       if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
+                               qcnt+=4;
+               }
+               else {
+                       if((ba.acnt[i]==4)&(ba.bcnt[i]==0))
+                               qcnt+=4;
+               }
+       }
+
+printf("%d %d\n",qcnt,ba.tcnt);
+       if(quality)
+               *quality=1.0*qcnt/ba.tcnt;
+       else
+               printf("[moldyn] bond analyze: quality = %f\n",
+                      1.0*qcnt/ba.tcnt);
 
        return 0;
 }

diff --git a/moldyn.h b/moldyn.h
index 83b3a86..1647c1c 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -219,6 +219,19 @@ typedef struct s_moldyn {
        double debug;           /* debugging stuff, ignore */
 } t_moldyn;
 
+typedef struct s_pcc {
+       int o1;
+       int o2;
+       double dr;
+       double *stat;
+} t_pcc;
+
+typedef struct s_ba {
+       int *acnt;
+       int *bcnt;
+       int tcnt;
+} t_ba;
+
 /*
  *
  *  defines
@@ -494,12 +507,21 @@ 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_free_save_file(t_moldyn *moldyn);
 int moldyn_load(t_moldyn *moldyn);
+int process_2b_bonds(t_moldyn *moldyn,void *data,
+                     int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+                                    void *data,u8 bc));
 int get_line(int fd,char *line,int max);
 
 int pair_correlation_init(t_moldyn *moldyn,double dr);
 int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc);
+int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
+                                       t_atom *jtom,void *data,u8 bc);
 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr);
+int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+                         void *data,u8 bc);
+int bond_analyze(t_moldyn *moldyn,double *quality);
 
 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
index 7f25a74..b005960 100755 (executable)
--- a/pair_corr_calc_script
+++ b/pair_corr_calc_script
@@ -27,7 +27,11 @@ fi
 
 if [ "$3" = "g" ]; then
 
-pdir=`dirname $1`
+if [ -d $1 ]; then
+       pdir=$1
+else
+       pdir=`dirname $1`
+fi
 pfile=$pdir/pair_corr.scr
 
 cat > $pfile <<-EOF

diff --git a/pair_correlation_calc.c b/pair_correlation_calc.c
index e5903aa..c7716a3 100644 (file)
--- a/pair_correlation_calc.c
+++ b/pair_correlation_calc.c
@@ -76,5 +76,7 @@ int main(int argc,char **argv) {
 
        free(stat);
 
+       moldyn_free_save_file(&moldyn);
+
        return 0;
 }