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

diff --git a/moldyn.c b/moldyn.c
index aebd604..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);
@@ -404,11 +412,14 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
 	if(moldyn->rfd) {
 		dprintf(moldyn->rfd,report_end);
 		close(moldyn->rfd);
-		snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
+		snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
+		         moldyn->vlsdir);
 		system(sc);
-		snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
+		snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
+		         moldyn->vlsdir);
 		system(sc);
-		snprintf(sc,255,"cd %s && dvipdf report",moldyn->vlsdir);
+		snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1",
+		         moldyn->vlsdir);
 		system(sc);
 	}
 	if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
@@ -501,6 +512,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 		check_per_bound(moldyn,&(atom[ret].r));
 	}
 
+	/* update total system mass */
+	total_mass_calc(moldyn);
+
 	return ret;
 }
 
@@ -642,6 +656,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 	atom[count].tag=count;
 	atom[count].attr=attr;
 
+	/* update total system mass */
+	total_mass_calc(moldyn);
+
 	return 0;
 }
 
@@ -702,13 +719,29 @@ int thermal_init(t_moldyn *moldyn,u8 equi_init) {
 	return 0;
 }
 
+double total_mass_calc(t_moldyn *moldyn) {
+
+	int i;
+
+	moldyn->mass=0.0;
+
+	for(i=0;i<moldyn->count;i++)
+		moldyn->mass+=moldyn->atom[i].mass;
+
+	return moldyn->mass;
+}
+
 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;
 }
@@ -805,24 +838,91 @@ 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+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+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;
+}
+
+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",
+	       moldyn->mean_t);
+
+	/* ideal gas contribution */
+	ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
+	printf("  ideal gas contribution: %f\n",
+	       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->mass;
+
+	/* specific heat for nve ensemble */
+	moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_mean/(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_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;
+}
 
 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
@@ -1304,6 +1404,7 @@ return 0;
 		e_kin_calc(moldyn);
 		temperature_calc(moldyn);
 		pressure_calc(moldyn);
+		energy_fluctuation_calc(moldyn);
 		//tp=thermodynamic_pressure_calc(moldyn);
 //printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR);
 
@@ -1376,6 +1477,8 @@ return 0;
 			       moldyn->mean_gp/BAR,
 			       moldyn->volume);
 			fflush(stdout);
+printf("\n");
+get_heat_capacity(moldyn);
 		}
 
 		/* increase absolute time */
@@ -1385,8 +1488,9 @@ return 0;
 	}
 
 		/* check for hooks */
-		if(sched->hook)
-			sched->hook(moldyn,sched->hook_params);
+		if(sched->count+1<sched->total_sched)
+			if(sched->hook)
+				sched->hook(moldyn,sched->hook_params);
 
 		/* get a new info line */
 		printf("\n");
@@ -1807,19 +1911,19 @@ int moldyn_bc_check(t_moldyn *moldyn) {
 }
 
 /*
- * postprocessing functions
+ * post processing functions
  */
-#define LINE_MAX 128
-int read_line(int fd,char *line) {
+
+int get_line(int fd,char *line,int max) {
 
 	int count,ret;
 
 	count=0;
 
 	while(1) {
-		if(count==LINE_MAX) return count;
+		if(count==max) return count;
 		ret=read(fd,line+count,1);
-		if(ret<0) return ret;
+		if(ret<=0) return ret;
 		if(line[count]=='\n') {
 			line[count]='\0';
 			return count+1;
@@ -1828,79 +1932,3 @@ int read_line(int fd,char *line) {
 	}
 }
 
-int calc_fluctuations(double start,double end,char *file) {
-
-	int fd;
-	int count,ret;
-	double time,pot,kin,tot;
-	double p_m,k_m,t_m;
-	double p2_m,k2_m,t2_m;
-	double p_sum,k_sum,t_sum;
-	char buf[LINE_MAX];
-
-	fd=open(file,O_RDONLY);
-	if(fd<0) {
-		perror("[moldyn] post proc open");
-		return fd;
-	}
-
-	/* first calc the averages */
-	p_sum=0.0;
-	k_sum=0.0;
-	t_sum=0.0;
-	count=0;
-	while(1) {
-		ret=read_line(fd,buf);
-		if(ret<0) break;
-printf("%d\n",ret);
-		if(buf[0]=='#') continue;
-		sscanf(buf,"%lf %lf %lf %lf",&time,&kin,&pot,&tot);
-printf("%f %f %f %f\n",time,pot,kin,tot);
-		//if(time>end) break;
-		if((time>=start)&(time<=end)) {
-			p_sum+=pot;
-			k_sum+=kin;
-			t_sum+=tot;
-			count+=1;
-		}
-	}
-
-	p_m=p_sum/count;
-	k_m=k_sum/count;
-	t_m=t_sum/count;
-
-	/* mean square fluctuations */
-	if(lseek(fd,SEEK_SET,0)<0) {
-		perror("[moldyn] lseek");
-		return -1;
-	}
-	count=0;
-	p_sum=0.0;
-	k_sum=0.0;
-	t_sum=0.0;
-	while(1) {
-		ret=read_line(fd,buf);
-		if(ret<0) break;
-		if(buf[0]=='#') continue;
-		sscanf(buf,"%lf %lf %lf %lf",&time,&kin,&pot,&tot);
-		if(time>end) break;
-		if((time>=start)&(time<=end)) {
-			p_sum+=((pot-p_m)*(pot-p_m));
-			k_sum+=((kin-k_m)*(kin-k_m));
-			t_sum+=((tot-t_m)*(tot-t_m));
-			count+=1;
-		}
-	}
-
-	p2_m=p_sum/count;
-	k2_m=k_sum/count;
-	t2_m=t_sum/count;
-
-	printf("[moldyn] fluctuations (%f - %f):\n",start,end);
-	printf("  - averages   : %f %f %f\n",k_m,p_m,t_m);
-	printf("  - mean square: %f %f %f\n",k2_m,p2_m,t2_m);
-
-	close(fd);
-
-	return 0;
-}