X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;fp=moldyn.c;h=a9341941909c452a613ff62877e6b016bbf10571;hb=d843bfcc38ac4ad1ceffbb82a1aba8096f791b41;hp=8451c4436f43db2213bec83677c359fb4062f133;hpb=787848dc9bf3daa2ccb76296de904688cc504a45;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 8451c44..a934194 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -596,29 +596,48 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 
 	/* debug */
 	if(ret!=new) {
-		printf("[moldyn] creating lattice failed\n");
+		printf("[moldyn] creating %s lattice (lc=%f) incomplete\n",
+		       name,lc);
+		printf("  (ignore in case of partial lattice creation)\n");
 		printf("  amount of atoms\n");
 		printf("  - expected: %d\n",new);
 		printf("  - created: %d\n",ret);
-		return -1;
 	}
 
-	moldyn->count+=new;
-	printf("[moldyn] created %s lattice with %d atoms\n",name,new);
+	moldyn->count+=ret;
+	if(ret==new)
+		printf("[moldyn] created %s lattice with %d atoms\n",name,ret);
 
-	for(ret=0;ret<new;ret++) {
-		atom[ret].element=element;
-		atom[ret].mass=mass;
-		atom[ret].attr=attr;
-		atom[ret].brand=brand;
-		atom[ret].tag=count+ret;
-		check_per_bound(moldyn,&(atom[ret].r));
-		atom[ret].r_0=atom[ret].r;
+	for(new=0;new<ret;new++) {
+		atom[new].element=element;
+		atom[new].mass=mass;
+		atom[new].attr=attr;
+		atom[new].brand=brand;
+		atom[new].tag=count+new;
+		check_per_bound(moldyn,&(atom[new].r));
+		atom[new].r_0=atom[new].r;
 #ifdef PTHREADS
-		pthread_mutex_init(&(mutex[ret]),NULL);
+		pthread_mutex_init(&(mutex[new]),NULL);
 #endif
 	}
 
+	/* fix allocation */
+	ptr=realloc(moldyn->atom,moldyn->count*sizeof(t_atom));
+	if(!ptr) {
+		perror("[moldyn] realloc (create lattice - alloc fix)");
+		return -1;
+	}
+	moldyn->atom=ptr;
+
+#ifdef LOWMEM_LISTS
+	ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int));
+	if(!ptr) {
+		perror("[moldyn] list realloc (create lattice)");
+		return -1;
+	}
+	moldyn->lc.subcell->list=ptr;
+#endif
+
 	/* update total system mass */
 	total_mass_calc(moldyn);
 
@@ -719,6 +738,9 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 	int i,j,k;
 	t_3dvec o;
 	t_3dvec dist;
+	t_3dvec p;
+
+	p.x=0; p.y=0; p.z=0;
 
 	count=0;
 	if(origin)
@@ -728,9 +750,9 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 
 	/* shift partition values */
 	if(p_type) {
-		p_vals->p.x+=(a*lc)/2.0;
-		p_vals->p.y+=(b*lc)/2.0;
-		p_vals->p.z+=(c*lc)/2.0;
+		p.x=p_vals->p.x+(a*lc)/2.0;
+		p.y=p_vals->p.y+(b*lc)/2.0;
+		p.z=p_vals->p.z+(c*lc)/2.0;
 	}
 
 	r.x=o.x;
@@ -741,17 +763,35 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 			for(k=0;k<c;k++) {
 				switch(p_type) {
 					case PART_INSIDE_R:
-						v3_sub(&dist,&r,&(p_vals->p));
+						v3_sub(&dist,&r,&p);
 			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
 				v3_copy(&(atom[count].r),&r);
 				count+=1;
 			}
 						break;
 					case PART_OUTSIDE_R:
-						v3_sub(&dist,&r,&(p_vals->p));
+						v3_sub(&dist,&r,&p);
 			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
 				v3_copy(&(atom[count].r),&r);
 				count+=1;
+			}
+						break;
+					case PART_INSIDE_D:
+						v3_sub(&dist,&r,&p);
+			if((fabs(dist.x)<p_vals->d.x)&&
+			   (fabs(dist.y)<p_vals->d.y)&&
+			   (fabs(dist.z)<p_vals->d.z)) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
+			}
+						break;
+					case PART_OUTSIDE_D:
+						v3_sub(&dist,&r,&p);
+			if((fabs(dist.x)>=p_vals->d.x)||
+			   (fabs(dist.y)>=p_vals->d.y)||
+			   (fabs(dist.z)>=p_vals->d.z)) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
 			}
 						break;
 					default:	
@@ -784,6 +824,9 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 	t_3dvec o,r,n;
 	t_3dvec basis[3];
 	t_3dvec dist;
+	t_3dvec p;
+
+	p.x=0; p.y=0; p.z=0;
 
 	count=0;
 	if(origin)
@@ -802,9 +845,9 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 
 	/* shift partition values */
 	if(p_type) {
-		p_vals->p.x+=(a*lc)/2.0;
-		p_vals->p.y+=(b*lc)/2.0;
-		p_vals->p.z+=(c*lc)/2.0;
+		p.x=p_vals->p.x+(a*lc)/2.0;
+		p.y=p_vals->p.y+(b*lc)/2.0;
+		p.z=p_vals->p.z+(c*lc)/2.0;
 	}
 
 	/* fill up the room */
@@ -817,17 +860,35 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 				/* first atom */
 				switch(p_type) {
 					case PART_INSIDE_R:
-						v3_sub(&dist,&r,&(p_vals->p));
+						v3_sub(&dist,&r,&p);
 			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
 				v3_copy(&(atom[count].r),&r);
 				count+=1;
 			}
 						break;
 					case PART_OUTSIDE_R:
-						v3_sub(&dist,&r,&(p_vals->p));
+						v3_sub(&dist,&r,&p);
 			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
 				v3_copy(&(atom[count].r),&r);
 				count+=1;
+			}
+						break;
+					case PART_INSIDE_D:
+						v3_sub(&dist,&r,&p);
+			if((fabs(dist.x)<p_vals->d.x)&&
+			   (fabs(dist.y)<p_vals->d.y)&&
+			   (fabs(dist.z)<p_vals->d.z)) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
+			}
+						break;
+					case PART_OUTSIDE_D:
+						v3_sub(&dist,&r,&p);
+			if((fabs(dist.x)>=p_vals->d.x)||
+			   (fabs(dist.y)>=p_vals->d.y)||
+			   (fabs(dist.z)>=p_vals->d.z)) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
 			}
 						break;
 					default:
@@ -840,19 +901,37 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 					v3_add(&n,&r,&basis[l]);
 					switch(p_type) {
 						case PART_INSIDE_R:
-			v3_sub(&dist,&n,&(p_vals->p));
+			v3_sub(&dist,&n,&p);
 			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
-				v3_copy(&(atom[count].r),&r);
+				v3_copy(&(atom[count].r),&n);
 				count+=1;
 			}
 							break;
 						case PART_OUTSIDE_R:
-			v3_sub(&dist,&n,&(p_vals->p));
+			v3_sub(&dist,&n,&p);
 			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
-				v3_copy(&(atom[count].r),&r);
+				v3_copy(&(atom[count].r),&n);
 				count+=1;
 			}
 							break;
+					case PART_INSIDE_D:
+						v3_sub(&dist,&n,&p);
+			if((fabs(dist.x)<p_vals->d.x)&&
+			   (fabs(dist.y)<p_vals->d.y)&&
+			   (fabs(dist.z)<p_vals->d.z)) {
+				v3_copy(&(atom[count].r),&n);
+				count+=1;
+			}
+						break;
+					case PART_OUTSIDE_D:
+						v3_sub(&dist,&n,&p);
+			if((fabs(dist.x)>=p_vals->d.x)||
+			   (fabs(dist.y)>=p_vals->d.y)||
+			   (fabs(dist.z)>=p_vals->d.z)) {
+				v3_copy(&(atom[count].r),&n);
+				count+=1;
+			}
+						break;
 						default:
 					v3_copy(&(atom[count].r),&n);
 					count+=1;
@@ -968,7 +1047,9 @@ double temperature_calc(t_moldyn *moldyn) {
 
 	/* assume up to date kinetic energy, which is 3/2 N k_B T */
 
-	moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
+	if(moldyn->count)
+		moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
+	else moldyn->t=0.0;
 
 	return moldyn->t;
 }
@@ -1068,8 +1149,8 @@ double virial_sum(t_moldyn *moldyn) {
 	}
 
 	/* global virial (absolute coordinates) */
-	virial=&(moldyn->gvir);
-	moldyn->gv=virial->xx+virial->yy+virial->zz;
+	//virial=&(moldyn->gvir);
+	//moldyn->gv=virial->xx+virial->yy+virial->zz;
 
 	return moldyn->virial;
 }
@@ -1090,9 +1171,16 @@ double pressure_calc(t_moldyn *moldyn) {
 	moldyn->p=2.0*moldyn->ekin+moldyn->virial;
 	moldyn->p/=(3.0*moldyn->volume);
 
+	moldyn->px=2.0*moldyn->ekinx+moldyn->vir.xx;
+	moldyn->px/=moldyn->volume;
+	moldyn->py=2.0*moldyn->ekiny+moldyn->vir.yy;
+	moldyn->py/=moldyn->volume;
+	moldyn->pz=2.0*moldyn->ekinz+moldyn->vir.zz;
+	moldyn->pz/=moldyn->volume;
+
 	/* pressure (absolute coordinates) */
-	moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
-	moldyn->gp/=(3.0*moldyn->volume);
+	//moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
+	//moldyn->gp/=(3.0*moldyn->volume);
 
 	return moldyn->p;
 }
@@ -1117,11 +1205,11 @@ int average_reset(t_moldyn *moldyn) {
 
 	/* virial */
 	moldyn->virial_sum=0.0;
-	moldyn->gv_sum=0.0;
+	//moldyn->gv_sum=0.0;
 
 	/* pressure */
 	moldyn->p_sum=0.0;
-	moldyn->gp_sum=0.0;
+	//moldyn->gp_sum=0.0;
 	moldyn->tp_sum=0.0;
 
 	return 0;
@@ -1159,14 +1247,14 @@ int average_and_fluctuation_calc(t_moldyn *moldyn) {
 	/* virial */
 	moldyn->virial_sum+=moldyn->virial;
 	moldyn->virial_avg=moldyn->virial_sum/denom;
-	moldyn->gv_sum+=moldyn->gv;
-	moldyn->gv_avg=moldyn->gv_sum/denom;
+	//moldyn->gv_sum+=moldyn->gv;
+	//moldyn->gv_avg=moldyn->gv_sum/denom;
 
 	/* pressure */
 	moldyn->p_sum+=moldyn->p;
 	moldyn->p_avg=moldyn->p_sum/denom;
-	moldyn->gp_sum+=moldyn->gp;
-	moldyn->gp_avg=moldyn->gp_sum/denom;
+	//moldyn->gp_sum+=moldyn->gp;
+	//moldyn->gp_avg=moldyn->gp_sum/denom;
 	moldyn->tp_sum+=moldyn->tp;
 	moldyn->tp_avg=moldyn->tp_sum/denom;
 
@@ -1320,17 +1408,47 @@ int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
 	return 0;
 }
 
+int scale_atoms_ind(t_moldyn *moldyn,double x,double y,double z) {
+
+	int i;
+	t_3dvec *r;
+
+	for(i=0;i<moldyn->count;i++) {
+		r=&(moldyn->atom[i].r);
+		r->x*=x;
+		r->y*=y;
+		r->z*=z;
+	}
+
+	return 0;
+}
+
+int scale_dim_ind(t_moldyn *moldyn,double x,double y,double z) {
+
+	t_3dvec *dim;
+
+	dim=&(moldyn->dim);
+
+	dim->x*=x;
+	dim->y*=y;
+	dim->z*=z;
+
+	return 0;
+}
+
 int scale_volume(t_moldyn *moldyn) {
 
 	t_3dvec *dim,*vdim;
-	double scale;
+	//double scale;
 	t_linkcell *lc;
+	double sx,sy,sz;
 
 	vdim=&(moldyn->vis.dim);
 	dim=&(moldyn->dim);
 	lc=&(moldyn->lc);
 
 	/* scaling factor */
+	/*
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
 		scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau;
 		scale=pow(scale,ONE_THIRD);
@@ -1338,10 +1456,20 @@ int scale_volume(t_moldyn *moldyn) {
 	else {
 		scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
 	}
+	*/
+
+	sx=1.0-(moldyn->p_ref-moldyn->px)*moldyn->p_tc*moldyn->tau;
+	sy=1.0-(moldyn->p_ref-moldyn->py)*moldyn->p_tc*moldyn->tau;
+	sz=1.0-(moldyn->p_ref-moldyn->pz)*moldyn->p_tc*moldyn->tau;
+	sx=pow(sx,ONE_THIRD);
+	sy=pow(sy,ONE_THIRD);
+	sz=pow(sz,ONE_THIRD);
 
 	/* scale the atoms and dimensions */
-	scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
-	scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+	//scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+	//scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+	scale_atoms_ind(moldyn,sx,sy,sz);
+	scale_dim_ind(moldyn,sx,sy,sz);
 
 	/* visualize dimensions */
 	if(vdim->x!=0) {
@@ -1360,9 +1488,12 @@ int scale_volume(t_moldyn *moldyn) {
 		link_cell_shutdown(moldyn);
 		link_cell_init(moldyn,QUIET);
 	} else {
-		lc->x*=scale;
-		lc->y*=scale;
-		lc->z*=scale;
+		//lc->x*=scale;
+		//lc->y*=scale;
+		//lc->z*=scale;
+		lc->x*=sx;
+		lc->y*=sx;
+		lc->z*=sy;
 	}
 
 	return 0;
@@ -1376,10 +1507,16 @@ double e_kin_calc(t_moldyn *moldyn) {
 
 	atom=moldyn->atom;
 	moldyn->ekin=0.0;
+	moldyn->ekinx=0.0;
+	moldyn->ekiny=0.0;
+	moldyn->ekinz=0.0;
 
 	for(i=0;i<moldyn->count;i++) {
 		atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 		moldyn->ekin+=atom[i].ekin;
+		moldyn->ekinx+=0.5*atom[i].mass*atom[i].v.x*atom[i].v.x;
+		moldyn->ekiny+=0.5*atom[i].mass*atom[i].v.y*atom[i].v.y;
+		moldyn->ekinz+=0.5*atom[i].mass*atom[i].v.z*atom[i].v.z;
 	}
 
 	return moldyn->ekin;
@@ -1884,7 +2021,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 				dprintf(moldyn->pfd,
 				        "%f %f %f %f %f %f %f\n",moldyn->time,
 				         moldyn->p/BAR,moldyn->p_avg/BAR,
-				         moldyn->gp/BAR,moldyn->gp_avg/BAR,
+				         moldyn->p/BAR,moldyn->p_avg/BAR,
 					 moldyn->tp/BAR,moldyn->tp_avg/BAR);
 			}
 		}
@@ -1898,7 +2035,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		if(v) {
 			if(!(moldyn->total_steps%v)) {
 				dprintf(moldyn->vfd,
-				        "%f %f\n",moldyn->time,moldyn->volume);
+				        "%f %f %f %f %f\n",moldyn->time,
+				                           moldyn->volume,
+				                           moldyn->dim.x,
+				                           moldyn->dim.y,
+				                           moldyn->dim.z);
 			}
 		}
 		if(s) {
@@ -2596,6 +2737,16 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 		size-=cnt;
 	}
 
+#ifdef PTHREADS
+	amutex=malloc(moldyn->count*sizeof(pthread_mutex_t));
+	if(amutex==NULL) {
+		perror("[moldyn] load save file (mutexes)");
+		return -1;
+	}
+	for(cnt=0;cnt<moldyn->count;cnt++)
+		pthread_mutex_init(&(amutex[cnt]),NULL);
+#endif
+
 	// hooks etc ...
 
 	return 0;