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

diff --git a/moldyn.c b/moldyn.c
index 7496223..91fead7 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -191,15 +191,12 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 		moldyn->vis.dim.z=z;
 	}
 
-	moldyn->dv=0.000001*moldyn->volume;
-
 	printf("[moldyn] dimensions in A and A^3 respectively:\n");
 	printf("  x: %f\n",moldyn->dim.x);
 	printf("  y: %f\n",moldyn->dim.y);
 	printf("  z: %f\n",moldyn->dim.z);
 	printf("  volume: %f\n",moldyn->volume);
 	printf("  visualize simulation box: %s\n",visualize?"yes":"no");
-	printf("  delta volume (pressure calc): %f\n",moldyn->dv);
 
 	return 0;
 }
@@ -370,12 +367,25 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 			dprintf(moldyn->tfd,"# temperature log file\n");
 			printf("temperature (%d)\n",timer);
 			break;
+		case LOG_VOLUME:
+			moldyn->vwrite=timer;
+			snprintf(filename,127,"%s/volume",moldyn->vlsdir);
+			moldyn->vfd=open(filename,
+			                 O_WRONLY|O_CREAT|O_EXCL,
+			                 S_IRUSR|S_IWUSR);
+			if(moldyn->vfd<0) {
+				perror("[moldyn] volume log file\n");
+				return moldyn->vfd;
+			}
+			dprintf(moldyn->vfd,"# volume log file\n");
+			printf("volume (%d)\n",timer);
+			break;
 		case SAVE_STEP:
 			moldyn->swrite=timer;
 			printf("save file (%d)\n",timer);
 			break;
 		case VISUAL_STEP:
-			moldyn->vwrite=timer;
+			moldyn->awrite=timer;
 			ret=visual_init(moldyn,moldyn->vlsdir);
 			if(ret<0) {
 				printf("[moldyn] visual init failure\n");
@@ -916,16 +926,26 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) {
 double virial_sum(t_moldyn *moldyn) {
 
 	int i;
-	double v;
 	t_virial *virial;
 
 	/* virial (sum over atom virials) */
-	v=0.0;
+	moldyn->virial=0.0;
+	moldyn->vir.xx=0.0;
+	moldyn->vir.yy=0.0;
+	moldyn->vir.zz=0.0;
+	moldyn->vir.xy=0.0;
+	moldyn->vir.xz=0.0;
+	moldyn->vir.yz=0.0;
 	for(i=0;i<moldyn->count;i++) {
 		virial=&(moldyn->atom[i].virial);
-		v+=(virial->xx+virial->yy+virial->zz);
+		moldyn->virial+=(virial->xx+virial->yy+virial->zz);
+		moldyn->vir.xx+=virial->xx;
+		moldyn->vir.yy+=virial->yy;
+		moldyn->vir.zz+=virial->zz;
+		moldyn->vir.xy+=virial->xy;
+		moldyn->vir.xz+=virial->xz;
+		moldyn->vir.yz+=virial->yz;
 	}
-	moldyn->virial=v;
 
 	/* global virial (absolute coordinates) */
 	virial=&(moldyn->gvir);
@@ -982,6 +1002,7 @@ int average_reset(t_moldyn *moldyn) {
 	/* pressure */
 	moldyn->p_sum=0.0;
 	moldyn->gp_sum=0.0;
+	moldyn->tp_sum=0.0;
 
 	return 0;
 }
@@ -1026,6 +1047,8 @@ int average_and_fluctuation_calc(t_moldyn *moldyn) {
 	moldyn->p_avg=moldyn->p_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;
 
 	return 0;
 }
@@ -1067,8 +1090,9 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
 	t_3dvec dim;
 	//t_3dvec *tp;
-	double u_up,u_down,dv;
-	double scale,p;
+	double h,dv;
+	double y0,y1;
+	double su,sd;
 	t_atom *store;
 
 	/*
@@ -1078,54 +1102,56 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	 *
 	 */
 
-	scale=0.00001;
-	dv=8*scale*scale*scale*moldyn->volume;
-
+	/* store atomic configuration + dimension */
 	store=malloc(moldyn->count*sizeof(t_atom));
 	if(store==NULL) {
 		printf("[moldyn] allocating store mem failed\n");
 		return -1;
 	}
-
-	/* save unscaled potential energy + atom/dim configuration */
 	memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
 	dim=moldyn->dim;
 
+	/* x1, y1 */
+	sd=0.00001;
+	h=(1.0-sd)*(1.0-sd)*(1.0-sd);
+	su=pow(2.0-h,ONE_THIRD)-1.0;
+	dv=(1.0-h)*moldyn->volume;
+
 	/* scale up dimension and atom positions */
-	scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
-	scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
+	scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
+	scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	u_up=moldyn->energy;
+	y1=moldyn->energy;
 
 	/* restore atomic configuration + dim */
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 	moldyn->dim=dim;
 
 	/* scale down dimension and atom positions */
-	scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
-	scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
+	scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
+	scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	u_down=moldyn->energy;
+	y0=moldyn->energy;
 	
 	/* calculate pressure */
-	p=-(u_up-u_down)/dv;
-printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR);
+	moldyn->tp=-(y1-y0)/(2.0*dv);
 
-	/* restore atomic configuration + dim */
+	/* restore atomic configuration */
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 	moldyn->dim=dim;
-
-	/* restore energy */
-	potential_force_calc(moldyn);
-
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
+	//potential_force_calc(moldyn);
 
-	return p;
+	/* free store buffer */
+	if(store)
+		free(store);
+
+	return moldyn->tp;
 }
 
 double get_pressure(t_moldyn *moldyn) {
@@ -1186,13 +1212,12 @@ int scale_volume(t_moldyn *moldyn) {
 
 	/* scaling factor */
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
-		scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc;
+		scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc;
 		scale=pow(scale,ONE_THIRD);
 	}
 	else {
 		scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
 	}
-moldyn->debug=scale;
 
 	/* scale the atoms and dimensions */
 	scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
@@ -1538,7 +1563,7 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
 int moldyn_integrate(t_moldyn *moldyn) {
 
 	int i;
-	unsigned int e,m,s,v,p,t;
+	unsigned int e,m,s,v,p,t,a;
 	t_3dvec momentum;
 	t_moldyn_schedule *sched;
 	t_atom *atom;
@@ -1560,6 +1585,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	m=moldyn->mwrite;
 	s=moldyn->swrite;
 	v=moldyn->vwrite;
+	a=moldyn->awrite;
 	p=moldyn->pwrite;
 	t=moldyn->twrite;
 
@@ -1621,6 +1647,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		temperature_calc(moldyn);
 		virial_sum(moldyn);
 		pressure_calc(moldyn);
+		//thermodynamic_pressure_calc(moldyn);
+
+		/* calculate fluctuations + averages */
 		average_and_fluctuation_calc(moldyn);
 
 		/* p/t scaling */
@@ -1650,9 +1679,10 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		if(p) {
 			if(!(moldyn->total_steps%p)) {
 				dprintf(moldyn->pfd,
-				        "%f %f %f %f %f\n",moldyn->time,
+				        "%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->gp/BAR,moldyn->gp_avg/BAR,
+					 moldyn->tp/BAR,moldyn->tp_avg/BAR);
 			}
 		}
 		if(t) {
@@ -1662,6 +1692,12 @@ int moldyn_integrate(t_moldyn *moldyn) {
 				        moldyn->time,moldyn->t,moldyn->t_avg);
 			}
 		}
+		if(v) {
+			if(!(moldyn->total_steps%v)) {
+				dprintf(moldyn->vfd,
+				        "%f %f\n",moldyn->time,moldyn->volume);
+			}
+		}
 		if(s) {
 			if(!(moldyn->total_steps%s)) {
 				snprintf(dir,128,"%s/s-%07.f.save",
@@ -1677,8 +1713,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
 				close(fd);
 			}	
 		}
-		if(v) {
-			if(!(moldyn->total_steps%v)) {
+		if(a) {
+			if(!(moldyn->total_steps%a)) {
 				visual_atoms(moldyn);
 			}
 		}
@@ -1688,10 +1724,10 @@ int moldyn_integrate(t_moldyn *moldyn) {
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
-printf("\rsched:%d, steps:%d/%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
+printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.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->p/BAR,moldyn->p_avg/BAR,
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));