X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=433be6824c9663f9bfa745ca0b82d08b566abd2b;hb=refs%2Fheads%2Forigin;hp=cd2a80309c6cf9de56e4b23132debbe9e46b02cd;hpb=0fbd532d8ddce9848df592ed586ffba439369284;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index cd2a803..433be68 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -218,6 +218,14 @@ int set_potential_params(t_moldyn *moldyn,void *params) {
 	return 0;
 }
 
+int set_avg_skip(t_moldyn *moldyn,int skip) {
+
+	printf("[moldyn] skip %d steps before starting average calc\n",skip);
+	moldyn->avg_skip=skip;
+
+	return 0;
+}
+
 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 
 	strncpy(moldyn->vlsdir,dir,127);
@@ -728,8 +736,6 @@ 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);
-	moldyn->t_sum+=moldyn->t;
-	moldyn->mean_t=moldyn->t_sum/moldyn->total_steps;
 
 	return moldyn->t;
 }
@@ -804,64 +810,88 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) {
 	return p;
 }
 
-double pressure_calc(t_moldyn *moldyn) {
+double virial_sum(t_moldyn *moldyn) {
 
 	int i;
 	double v;
 	t_virial *virial;
 
+	/* virial (sum over atom virials) */
+	v=0.0;
+	for(i=0;i<moldyn->count;i++) {
+		virial=&(moldyn->atom[i].virial);
+		v+=(virial->xx+virial->yy+virial->zz);
+	}
+	moldyn->virial=v;
+
+	/* global virial (absolute coordinates) */
+	virial=&(moldyn->gvir);
+	moldyn->gv=virial->xx+virial->yy+virial->zz;
+
+	return moldyn->virial;
+}
+
+double pressure_calc(t_moldyn *moldyn) {
+
 	/*
 	 * PV = NkT + <W>
-	 * W = 1/3 sum_i f_i r_i
+	 * with W = 1/3 sum_i f_i r_i (- skipped!)
 	 * virial = sum_i f_i r_i
 	 * 
 	 * => P = (2 Ekin + virial) / (3V)
 	 */
 
-	v=0.0;
-	for(i=0;i<moldyn->count;i++) {
-		virial=&(moldyn->atom[i].virial);
-		v+=(virial->xx+virial->yy+virial->zz);
-	}
+	/* assume up to date virial & up to date kinetic energy */
 
-	/* virial sum and mean virial */
-	moldyn->virial_sum+=v;
-	moldyn->mean_v=moldyn->virial_sum/moldyn->total_steps;
-
-	/* assume up to date kinetic energy */
-	moldyn->p=2.0*moldyn->ekin+moldyn->mean_v;
+	/* pressure (atom virials) */
+	moldyn->p=2.0*moldyn->ekin+moldyn->virial;
 	moldyn->p/=(3.0*moldyn->volume);
-	moldyn->p_sum+=moldyn->p;
-	moldyn->mean_p=moldyn->p_sum/moldyn->total_steps;
 
-	/* pressure from 'absolute coordinates' virial */
-	virial=&(moldyn->virial);
-	v=virial->xx+virial->yy+virial->zz;
-	moldyn->gp=2.0*moldyn->ekin+v;
+	/* pressure (absolute coordinates) */
+	moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
 	moldyn->gp/=(3.0*moldyn->volume);
-	moldyn->gp_sum+=moldyn->gp;
-	moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps;
 
 	return moldyn->p;
 }
 
-int energy_fluctuation_calc(t_moldyn *moldyn) {
+int average_and_fluctuation_calc(t_moldyn *moldyn) {
+
+	if(moldyn->total_steps<moldyn->avg_skip)
+		return 0;
 
-	/* assume up to date energies */
+	int denom=moldyn->total_steps+1-moldyn->avg_skip;
+
+	/* assume up to date energies, temperature, pressure etc */
 
 	/* kinetic energy */
 	moldyn->k_sum+=moldyn->ekin;
 	moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
-	moldyn->k_mean=moldyn->k_sum/moldyn->total_steps;
-	moldyn->k2_mean=moldyn->k2_sum/moldyn->total_steps;
-	moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean);
+	moldyn->k_avg=moldyn->k_sum/denom;
+	moldyn->k2_avg=moldyn->k2_sum/denom;
+	moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg);
 
 	/* potential energy */
 	moldyn->v_sum+=moldyn->energy;
 	moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
-	moldyn->v_mean=moldyn->v_sum/moldyn->total_steps;
-	moldyn->v2_mean=moldyn->v2_sum/moldyn->total_steps;
-	moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean);
+	moldyn->v_avg=moldyn->v_sum/denom;
+	moldyn->v2_avg=moldyn->v2_sum/denom;
+	moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg);
+
+	/* temperature */
+	moldyn->t_sum+=moldyn->t;
+	moldyn->t_avg=moldyn->t_sum/denom;
+
+	/* 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;
+
+	/* 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;
 
 	return 0;
 }
@@ -870,10 +900,14 @@ int get_heat_capacity(t_moldyn *moldyn) {
 
 	double temp2,ighc;
 
+	/* averages needed for heat capacity calc */
+	if(moldyn->total_steps<moldyn->avg_skip)
+		return 0;
+
 	/* (temperature average)^2 */
-	temp2=moldyn->mean_t*moldyn->mean_t;
+	temp2=moldyn->t_avg*moldyn->t_avg;
 	printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
-	       moldyn->mean_t);
+	       moldyn->t_avg);
 
 	/* ideal gas contribution */
 	ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
@@ -881,15 +915,16 @@ int get_heat_capacity(t_moldyn *moldyn) {
 	       ighc/moldyn->mass*KILOGRAM/JOULE);
 
 	/* specific heat for nvt ensemble */
-	moldyn->c_v_nvt=moldyn->dv2_mean/(K_BOLTZMANN*temp2)+ighc;
+	moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc;
 	moldyn->c_v_nvt/=moldyn->mass;
 
 	/* specific heat for nve ensemble */
-	moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_mean/(ighc*K_BOLTZMANN*temp2)));
+	moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2)));
 	moldyn->c_v_nve/=moldyn->mass;
 
 	printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
 	printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
+printf("  --> <dV2> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->count*K_B2*moldyn->t_avg*moldyn->t_avg*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM)));
 
 	return 0;
 }
@@ -1356,7 +1391,8 @@ return 0;
 	printf("[moldyn] integration start, go get a coffee ...\n");
 
 	/* executing the schedule */
-	for(sched->count=0;sched->count<sched->total_sched;sched->count++) {
+	sched->count=0;
+	while(sched->count<sched->total_sched) {
 
 		/* setting amount of runs and finite time step size */
 		moldyn->tau=sched->tau[sched->count];
@@ -1373,10 +1409,9 @@ return 0;
 		/* calculate kinetic energy, temperature and pressure */
 		e_kin_calc(moldyn);
 		temperature_calc(moldyn);
+		virial_sum(moldyn);
 		pressure_calc(moldyn);
-		energy_fluctuation_calc(moldyn);
-		//tp=thermodynamic_pressure_calc(moldyn);
-//printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR);
+		average_and_fluctuation_calc(moldyn);
 
 		/* p/t scaling */
 		if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
@@ -1406,15 +1441,15 @@ return 0;
 			if(!(i%p)) {
 				dprintf(moldyn->pfd,
 				        "%f %f %f %f %f\n",moldyn->time,
-				         moldyn->p/BAR,moldyn->mean_p/BAR,
-				         moldyn->gp/BAR,moldyn->mean_gp/BAR);
+				         moldyn->p/BAR,moldyn->p_avg/BAR,
+				         moldyn->gp/BAR,moldyn->gp_avg/BAR);
 			}
 		}
 		if(t) {
 			if(!(i%t)) {
 				dprintf(moldyn->tfd,
 				        "%f %f %f\n",
-				        moldyn->time,moldyn->t,moldyn->mean_t);
+				        moldyn->time,moldyn->t,moldyn->t_avg);
 			}
 		}
 		if(s) {
@@ -1440,15 +1475,12 @@ return 0;
 
 		/* display progress */
 		if(!(i%10)) {
-			printf("\rsched: %d, steps: %d, T: %f, P: %f %f V: %f",
-			       sched->count,i,
-			       moldyn->mean_t,
-			       moldyn->mean_p/BAR,
-			       moldyn->mean_gp/BAR,
-			       moldyn->volume);
-			fflush(stdout);
-printf("\n");
-get_heat_capacity(moldyn);
+	printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f",
+	       sched->count,i,
+	       moldyn->t,moldyn->t_avg,
+	       moldyn->p_avg/BAR,moldyn->p/BAR,
+	       moldyn->volume);
+	fflush(stdout);
 		}
 
 		/* increase absolute time */
@@ -1458,12 +1490,15 @@ get_heat_capacity(moldyn);
 	}
 
 		/* check for hooks */
-		if(sched->count+1<sched->total_sched)
-			if(sched->hook)
-				sched->hook(moldyn,sched->hook_params);
+		if(sched->hook) {
+			printf("\n ## schedule hook %d/%d start ##\n",
+			       sched->count+1,sched->total_sched-1);
+			sched->hook(moldyn,sched->hook_params);
+			printf(" ## schedule hook end ##\n");
+		}
 
-		/* get a new info line */
-		printf("\n");
+		/* increase the schedule counter */
+		sched->count+=1;
 
 	}
 
@@ -1542,7 +1577,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 	moldyn->energy=0.0;
 
 	/* reset global virial */
-	memset(&(moldyn->virial),0,sizeof(t_virial));
+	memset(&(moldyn->gvir),0,sizeof(t_virial));
 
 	/* reset force, site energy and virial of every atom */
 	for(i=0;i<count;i++) {
@@ -1768,12 +1803,12 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 	/* calculate global virial */
 	for(i=0;i<count;i++) {
-		moldyn->virial.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
-		moldyn->virial.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
-		moldyn->virial.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
-		moldyn->virial.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
-		moldyn->virial.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
-		moldyn->virial.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
+		moldyn->gvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
+		moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
+		moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
+		moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
+		moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
+		moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
 	}
 
 	return 0;