From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Wed, 9 Apr 2008 12:47:25 +0000 (+0200)
Subject: introduce a 2body bond function with callback, modified pair corr calc
X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2346d076f3d0955fb3a7bfb47f2b17b9eae6dd5b;p=physik%2Fposic.git

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
--- 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
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
--- 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
--- 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
--- 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
--- 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;
 }