X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=2e0e46810f4ba065ee8d384d0fd556749d53db93;hb=dff2917d22ed07707d222bc10fab7370356699dc;hp=ad4f9003fcc9cdb6d6ccf44b59f58f18b418d346;hpb=799fb13ad0a987bacfd259bfc34ba3c4703dd402;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index ad4f900..2e0e468 100644
--- 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;