X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=84cf7003c13f70cd6e7180a7a23b53ddda6722da;hb=2346d076f3d0955fb3a7bfb47f2b17b9eae6dd5b;hp=91fead7739a0aeb4d45165ccfd48ea1427e56bc1;hpb=a32468230b319b32819f1b20fd28aa9659574d45;p=physik%2Fposic.git

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