X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=16c811ed5e0001d905fb7760239fa0b71a055305;hb=25c22fe95e80365056c6a7fadc548119360ca8ce;hp=cd2a80309c6cf9de56e4b23132debbe9e46b02cd;hpb=0fbd532d8ddce9848df592ed586ffba439369284;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index cd2a803..16c811e 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -218,6 +218,14 @@ int set_potential_params(t_moldyn *moldyn,void *params) {
 	return 0;
 }
 
+int set_mean_skip(t_moldyn *moldyn,int skip) {
+
+	printf("[moldyn] skip %d steps before starting average calc\n",skip);
+	moldyn->mean_skip=skip;
+
+	return 0;
+}
+
 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 
 	strncpy(moldyn->vlsdir,dir,127);
@@ -728,8 +736,12 @@ 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->total_steps<moldyn->mean_skip)
+		return 0;
+
 	moldyn->t_sum+=moldyn->t;
-	moldyn->mean_t=moldyn->t_sum/moldyn->total_steps;
+	moldyn->mean_t=moldyn->t_sum/(moldyn->total_steps+1-moldyn->mean_skip);
 
 	return moldyn->t;
 }
@@ -826,41 +838,54 @@ double pressure_calc(t_moldyn *moldyn) {
 
 	/* virial sum and mean virial */
 	moldyn->virial_sum+=v;
-	moldyn->mean_v=moldyn->virial_sum/moldyn->total_steps;
+	if(moldyn->total_steps>=moldyn->mean_skip)
+		moldyn->mean_v=moldyn->virial_sum/
+		               (moldyn->total_steps+1-moldyn->mean_skip);
 
 	/* assume up to date kinetic energy */
 	moldyn->p=2.0*moldyn->ekin+moldyn->mean_v;
 	moldyn->p/=(3.0*moldyn->volume);
-	moldyn->p_sum+=moldyn->p;
-	moldyn->mean_p=moldyn->p_sum/moldyn->total_steps;
+	if(moldyn->total_steps>=moldyn->mean_skip) {
+		moldyn->p_sum+=moldyn->p;
+		moldyn->mean_p=moldyn->p_sum/
+		               (moldyn->total_steps+1-moldyn->mean_skip);
+	}
 
 	/* pressure from 'absolute coordinates' virial */
 	virial=&(moldyn->virial);
 	v=virial->xx+virial->yy+virial->zz;
 	moldyn->gp=2.0*moldyn->ekin+v;
 	moldyn->gp/=(3.0*moldyn->volume);
-	moldyn->gp_sum+=moldyn->gp;
-	moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps;
+	if(moldyn->total_steps>=moldyn->mean_skip) {
+		moldyn->gp_sum+=moldyn->gp;
+		moldyn->mean_gp=moldyn->gp_sum/
+		               (moldyn->total_steps+1-moldyn->mean_skip);
+	}
 
 	return moldyn->p;
 }
 
 int energy_fluctuation_calc(t_moldyn *moldyn) {
 
+	if(moldyn->total_steps<moldyn->mean_skip)
+		return 0;
+
 	/* assume up to date energies */
 
 	/* 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->k_mean=moldyn->k_sum/(moldyn->total_steps+1-moldyn->mean_skip);
+	moldyn->k2_mean=moldyn->k2_sum/
+	                (moldyn->total_steps+1-moldyn->mean_skip);
 	moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean);
 
 	/* 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->v_mean=moldyn->v_sum/(moldyn->total_steps+1-moldyn->mean_skip);
+	moldyn->v2_mean=moldyn->v2_sum/
+	                (moldyn->total_steps+1-moldyn->mean_skip);
 	moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean);
 
 	return 0;
@@ -870,6 +895,10 @@ int get_heat_capacity(t_moldyn *moldyn) {
 
 	double temp2,ighc;
 
+	/* averages needed for heat capacity calc */
+	if(moldyn->total_steps<moldyn->mean_skip)
+		return 0;
+
 	/* (temperature average)^2 */
 	temp2=moldyn->mean_t*moldyn->mean_t;
 	printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
@@ -890,6 +919,7 @@ int get_heat_capacity(t_moldyn *moldyn) {
 
 	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_mean,1.5*moldyn->count*K_B2*moldyn->mean_t*moldyn->mean_t*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM)));
 
 	return 0;
 }