X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=edda007d9290ff6cfe1001d4105c3f8a9ab12863;hb=e1080fc0dd66b0cf5b7715c5e99e7a34ac04a8cf;hp=997c32ccb24d44607d746d1e6aed87027e9de6cd;hpb=b5b47daaa3718c4dec2056fe5147668023575b8e;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 997c32c..edda007 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -956,10 +956,12 @@ double pressure_calc(t_moldyn *moldyn) {
 
 int average_and_fluctuation_calc(t_moldyn *moldyn) {
 
+	int denom;
+
 	if(moldyn->total_steps<moldyn->avg_skip)
 		return 0;
 
-	int denom=moldyn->total_steps+1-moldyn->avg_skip;
+	denom=moldyn->total_steps+1-moldyn->avg_skip;
 
 	/* assume up to date energies, temperature, pressure etc */
 
@@ -1031,7 +1033,8 @@ printf("  --> <dV2> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->coun
 
 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
-	t_3dvec dim,*tp;
+	t_3dvec dim;
+	//t_3dvec *tp;
 	double u_up,u_down,dv;
 	double scale,p;
 	t_atom *store;
@@ -1258,25 +1261,57 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 	lc->y=moldyn->dim.y/lc->ny;
 	lc->nz=moldyn->dim.z/moldyn->cutoff;
 	lc->z=moldyn->dim.z/lc->nz;
-
 	lc->cells=lc->nx*lc->ny*lc->nz;
+
+#ifdef STATIC_LISTS
+	lc->subcell=malloc(lc->cells*sizeof(int*));
+#else
 	lc->subcell=malloc(lc->cells*sizeof(t_list));
+#endif
+
+	if(lc->subcell==NULL) {
+		perror("[moldyn] cell init (malloc)");
+		return -1;
+	}
 
 	if(lc->cells<27)
 		printf("[moldyn] FATAL: less then 27 subcells!\n");
 
 	if(vol) {
-		printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
+#ifdef STATIC_LISTS
+		printf("[moldyn] initializing 'static' linked cells (%d)\n",
+		       lc->cells);
+#else
+		printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
+		       lc->cells);
+#endif
 		printf("  x: %d x %f A\n",lc->nx,lc->x);
 		printf("  y: %d x %f A\n",lc->ny,lc->y);
 		printf("  z: %d x %f A\n",lc->nz,lc->z);
 	}
 
+#ifdef STATIC_LISTS
+	/* list init */
+	for(i=0;i<lc->cells;i++) {
+		lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int));
+		if(lc->subcell[i]==NULL) {
+			perror("[moldyn] list init (malloc)");
+			return -1;
+		}
+		/*
+		if(i==0)
+			printf(" ---> %d malloc %p (%p)\n",
+			       i,lc->subcell[0],lc->subcell);
+		*/
+	}
+#else
 	for(i=0;i<lc->cells;i++)
 		list_init_f(&(lc->subcell[i]));
+#endif
+
+	/* update the list */
+	link_cell_update(moldyn);
 
- 	link_cell_update(moldyn);
-	
 	return 0;
 }
 
@@ -1286,7 +1321,9 @@ int link_cell_update(t_moldyn *moldyn) {
 	int nx,ny;
 	t_atom *atom;
 	t_linkcell *lc;
-	double x,y,z;
+#ifdef STATIC_LISTS
+	int p;
+#endif
 
 	atom=moldyn->atom;
 	lc=&(moldyn->lc);
@@ -1294,25 +1331,50 @@ int link_cell_update(t_moldyn *moldyn) {
 	nx=lc->nx;
 	ny=lc->ny;
 
-	x=moldyn->dim.x/2;
-	y=moldyn->dim.y/2;
-	z=moldyn->dim.z/2;
-
 	for(i=0;i<lc->cells;i++)
+#ifdef STATIC_LISTS
+		memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+#else
 		list_destroy_f(&(lc->subcell[i]));
-	
+#endif
+
 	for(count=0;count<moldyn->count;count++) {
 		i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
 		j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
 		k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
+
+#ifdef STATIC_LISTS
+		p=0;
+		while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
+			p++;
+
+		if(p>=MAX_ATOMS_PER_LIST) {
+			printf("[moldyn] FATAL: amount of atoms too high!\n");
+			return -1;
+		}
+
+		lc->subcell[i+j*nx+k*nx*ny][p]=count;
+#else
 		list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
 		                     &(atom[count]));
+		/*
+		if(j==0&&k==0)
+			printf(" ---> %d %d malloc %p (%p)\n",
+			       i,count,lc->subcell[i].current,lc->subcell);
+		*/
+#endif
 	}
 
 	return 0;
 }
 
-int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
+#ifdef STATIC_LISTS
+                              int **cell
+#else
+                              t_list *cell
+#endif
+			     ) {
 
 	t_linkcell *lc;
 	int a;
@@ -1375,8 +1437,14 @@ int link_cell_shutdown(t_moldyn *moldyn) {
 
 	lc=&(moldyn->lc);
 
-	for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
-		list_destroy_f(&(moldyn->lc.subcell[i]));
+	for(i=0;i<lc->cells;i++) {
+#ifdef STATIC_LISTS
+		free(lc->subcell[i]);
+#else
+		//printf(" ---> %d free %p\n",i,lc->subcell[i].start);
+		list_destroy_f(&(lc->subcell[i]));
+#endif
+	}
 
 	free(lc->subcell);
 
@@ -1580,21 +1648,22 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		}
 
 		/* display progress */
-		if(!(moldyn->total_steps%10)) {
+		//if(!(moldyn->total_steps%10)) {
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
-	printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
-	       sched->count,i,
-	       moldyn->t,moldyn->t_avg,
-	       moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
-	       moldyn->volume,
-	       (int)(t2.tv_sec-t1.tv_sec));
-	fflush(stdout);
+printf("\rsched:%d, steps:%d/%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
+       sched->count,i,moldyn->total_steps,
+       moldyn->t,moldyn->t_avg,
+       moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
+       moldyn->volume,
+       (int)(t2.tv_sec-t1.tv_sec));
+
+			fflush(stdout);
 
 			/* copy over time */
 			t1=t2;
-		}
+		//}
 
 		/* increase absolute time */
 		moldyn->time+=moldyn->tau;
@@ -1604,8 +1673,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
 		/* check for hooks */
 		if(sched->hook) {
-			printf("\n ## schedule hook %d/%d start ##\n",
-			       sched->count+1,sched->total_sched-1);
+			printf("\n ## schedule hook %d start ##\n",
+			       sched->count);
 			sched->hook(moldyn,sched->hook_params);
 			printf(" ## schedule hook end ##\n");
 		}
@@ -1676,15 +1745,24 @@ int potential_force_calc(t_moldyn *moldyn) {
 	t_atom *itom,*jtom,*ktom;
 	t_virial *virial;
 	t_linkcell *lc;
+#ifdef STATIC_LISTS
+	int *neighbour_i[27];
+	int p,q;
+	t_atom *atom;
+#else
 	t_list neighbour_i[27];
 	t_list neighbour_i2[27];
 	t_list *this,*that;
+#endif
 	u8 bc_ij,bc_ik;
 	int dnlc;
 
 	count=moldyn->count;
 	itom=moldyn->atom;
 	lc=&(moldyn->lc);
+#ifdef STATIC_LISTS
+	atom=moldyn->atom;
+#endif
 
 	/* reset energy */
 	moldyn->energy=0.0;
@@ -1739,14 +1817,33 @@ int potential_force_calc(t_moldyn *moldyn) {
 		if(moldyn->func2b) {
 			for(j=0;j<27;j++) {
 
+				bc_ij=(j<dnlc)?0:1;
+#ifdef STATIC_LISTS
+				p=0;
+
+				while(neighbour_i[j][p]!=0) {
+
+					jtom=&(atom[neighbour_i[j][p]]);
+					p++;
+
+					if(jtom==&(itom[i]))
+						continue;
+
+					if((jtom->attr&ATOM_ATTR_2BP)&
+					   (itom[i].attr&ATOM_ATTR_2BP)) {
+						moldyn->func2b(moldyn,
+						               &(itom[i]),
+						               jtom,
+					        	       bc_ij);
+					}
+				}
+#else
 				this=&(neighbour_i[j]);
 				list_reset_f(this);
 
 				if(this->start==NULL)
 					continue;
 
-				bc_ij=(j<dnlc)?0:1;
-
 				do {
 					jtom=this->current->data;
 
@@ -1761,6 +1858,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 					        	       bc_ij);
 					}
 				} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
 
 			}
 		}
@@ -1771,21 +1869,34 @@ int potential_force_calc(t_moldyn *moldyn) {
 			continue;
 
 		/* copy the neighbour lists */
+#ifdef STATIC_LISTS
+		/* no copy needed for static lists */
+#else
 		memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
+#endif
 
 		/* second loop over atoms j */
 		for(j=0;j<27;j++) {
 
+			bc_ij=(j<dnlc)?0:1;
+#ifdef STATIC_LISTS
+			p=0;
+
+			while(neighbour_i[j][p]!=0) {
+
+				jtom=&(atom[neighbour_i[j][p]]);
+				p++;
+#else
 			this=&(neighbour_i[j]);
 			list_reset_f(this);
 
 			if(this->start==NULL)
 				continue;
 
-			bc_ij=(j<dnlc)?0:1;
-
 			do {
+
 				jtom=this->current->data;
+#endif
 
 				if(jtom==&(itom[i]))
 					continue;
@@ -1811,17 +1922,24 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 				for(k=0;k<27;k++) {
 
+					bc_ik=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+					q=0;
+
+					while(neighbour_i[j][q]!=0) {
+
+						ktom=&(atom[neighbour_i[k][q]]);
+						q++;
+#else
 					that=&(neighbour_i2[k]);
 					list_reset_f(that);
 					
 					if(that->start==NULL)
 						continue;
 
-					bc_ik=(k<dnlc)?0:1;
-
 					do {
-
 						ktom=that->current->data;
+#endif
 
 						if(!(ktom->attr&ATOM_ATTR_3BP))
 							continue;
@@ -1837,9 +1955,12 @@ int potential_force_calc(t_moldyn *moldyn) {
 						                  jtom,
 						                  ktom,
 						                  bc_ik|bc_ij);
-
+#ifdef STATIC_LISTS
+					}
+#else
 					} while(list_next_f(that)!=\
 					        L_NO_NEXT_ELEMENT);
+#endif
 
 				}
 
@@ -1856,17 +1977,24 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 				for(k=0;k<27;k++) {
 
+					bc_ik=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+					q=0;
+
+					while(neighbour_i[j][q]!=0) {
+
+						ktom=&(atom[neighbour_i[k][q]]);
+						q++;
+#else
 					that=&(neighbour_i2[k]);
 					list_reset_f(that);
 					
 					if(that->start==NULL)
 						continue;
 
-					bc_ik=(k<dnlc)?0:1;
-
 					do {
-
 						ktom=that->current->data;
+#endif
 
 						if(!(ktom->attr&ATOM_ATTR_3BP))
 							continue;
@@ -1883,8 +2011,12 @@ int potential_force_calc(t_moldyn *moldyn) {
 						                  ktom,
 						                  bc_ik|bc_ij);
 
+#ifdef STATIC_LISTS
+					}
+#else
 					} while(list_next_f(that)!=\
 					        L_NO_NEXT_ELEMENT);
+#endif
 
 				}
 				
@@ -1896,8 +2028,11 @@ int potential_force_calc(t_moldyn *moldyn) {
 					                  &(itom[i]),
 					                  jtom,bc_ij);
 				}
-					
+#ifdef STATIC_LISTS
+			}
+#else
 			} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
 		
 		}
 		
@@ -2038,6 +2173,50 @@ int moldyn_bc_check(t_moldyn *moldyn) {
  * restore function
  */
 
+int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
+
+	int fd;
+	int cnt,size;
+
+	fd=open(file,O_RDONLY);
+	if(fd<0) {
+		perror("[moldyn] load save file open");
+		return fd;
+	}
+
+	size=sizeof(t_moldyn);
+
+	while(size) {
+		cnt=read(fd,moldyn,size);
+		if(cnt<0) {
+			perror("[moldyn] load save file read (moldyn)");
+			return cnt;
+		}
+		size-=cnt;
+	}
+
+	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;
+	}
+
+	while(size) {
+		cnt=read(fd,moldyn->atom,size);
+		if(cnt<0) {
+			perror("[moldyn] load save file read (atoms)");
+			return cnt;
+		}
+		size-=cnt;
+	}
+
+	// hooks etc ...
+
+	return 0;
+}
+
 int moldyn_load(t_moldyn *moldyn) {
 
 	// later ...
@@ -2067,6 +2246,163 @@ int get_line(int fd,char *line,int max) {
 	}
 }
 
+int pair_correlation_init(t_moldyn *moldyn,double dr) {
+
+	
+	return 0;
+}
+
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+
+	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;
+
+	lc=&(moldyn->lc);
+
+	slots=moldyn->cutoff/dr;
+	o=2*slots;
+
+	printf("[moldyn] pair correlation calc info:\n");
+	printf("  time: %f\n",moldyn->time);
+	printf("  count: %d\n",moldyn->count);
+	printf("  cutoff: %f\n",moldyn->cutoff);
+	printf("  temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
+
+	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);
+
+	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);
+
+		/* brand of atom i */
+		ibrand=itom[i].brand;
+	
+		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
+				/* only count pairs once,
+				 * skip same atoms */
+				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 or equal cutoff */
+				if(d>=moldyn->cutoff_square)
+					continue;
+
+				/* 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) s=slots-1;
+
+				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
+		}
+	}
+
+	/* 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;
+}
+
 int analyze_bonds(t_moldyn *moldyn) {
 
 	
@@ -2096,7 +2432,12 @@ int visual_atoms(t_moldyn *moldyn) {
 	t_atom *atom;
 	t_atom *btom;
 	t_linkcell *lc;
+#ifdef STATIC_LISTS
+	int *neighbour[27];
+	int p;
+#else
 	t_list neighbour[27];
+#endif
 	u8 bc;
 	t_3dvec dist;
 	double d2;
@@ -2150,12 +2491,19 @@ int visual_atoms(t_moldyn *moldyn) {
 		                          (atom[i].r.z+moldyn->dim.z/2)/lc->z,
 		                          neighbour);
 		for(j=0;j<27;j++) {
+			bc=j<lc->dnlc?0:1;
+#ifdef STATIC_LISTS
+			p=0;
+			while(neighbour[j][p]!=0) {
+				btom=&(atom[neighbour[j][p]]);
+				p++;
+#else
 			list_reset_f(&neighbour[j]);
 			if(neighbour[j].start==NULL)
 				continue;
-			bc=j<lc->dnlc?0:1;
 			do {
 				btom=neighbour[j].current->data;
+#endif
 				if(btom==&atom[i])	// skip identical atoms
 					continue;
 				//if(btom<&atom[i])	// skip half of them
@@ -2175,7 +2523,11 @@ int visual_atoms(t_moldyn *moldyn) {
 				dprintf(fd,"# [B] %f %f %f %f %f %f\n",
 				        atom[i].r.x,atom[i].r.y,atom[i].r.z,
 		       		        btom->r.x,btom->r.y,btom->r.z);
+#ifdef STATIC_LISTS
+			}
+#else
 			} while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT);
+#endif
 		}
 	}