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

diff --git a/moldyn.c b/moldyn.c
index aebd604..c130ef4 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -404,11 +404,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));
@@ -1385,8 +1388,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");
@@ -1809,17 +1813,17 @@ int moldyn_bc_check(t_moldyn *moldyn) {
 /*
  * postprocessing 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,19 +1832,21 @@ int read_line(int fd,char *line) {
 	}
 }
 
-int calc_fluctuations(double start,double end,char *file) {
+int calc_fluctuations(double start,double end,t_moldyn *moldyn) {
 
 	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];
+	char buf[64];
+	char file[128+7];
+
+	printf("[moldyn] calculating energy fluctuations [eV]:\n");
 
+	snprintf(file,128+7,"%s/energy",moldyn->vlsdir);
 	fd=open(file,O_RDONLY);
 	if(fd<0) {
-		perror("[moldyn] post proc open");
+		perror("[moldyn] post proc energy open");
 		return fd;
 	}
 
@@ -1850,24 +1856,21 @@ int calc_fluctuations(double start,double end,char *file) {
 	t_sum=0.0;
 	count=0;
 	while(1) {
-		ret=read_line(fd,buf);
-		if(ret<0) break;
-printf("%d\n",ret);
+		ret=get_line(fd,buf,63);
+		if(ret<=0) break;
 		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;
-		}
+		if(time<start) continue;
+		if(time>end) break;
+		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;
+	moldyn->p_m=p_sum/count;
+	moldyn->k_m=k_sum/count;
+	moldyn->t_m=t_sum/count;
 
 	/* mean square fluctuations */
 	if(lseek(fd,SEEK_SET,0)<0) {
@@ -1879,28 +1882,60 @@ printf("%f %f %f %f\n",time,pot,kin,tot);
 	k_sum=0.0;
 	t_sum=0.0;
 	while(1) {
-		ret=read_line(fd,buf);
-		if(ret<0) break;
+		ret=get_line(fd,buf,63);
+		if(ret<=0) break;
 		if(buf[0]=='#') continue;
 		sscanf(buf,"%lf %lf %lf %lf",&time,&kin,&pot,&tot);
+		if(time<start) continue;
 		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;
-		}
+		k_sum+=((kin-moldyn->k_m)*(kin-moldyn->k_m));
+		p_sum+=((pot-moldyn->p_m)*(pot-moldyn->p_m));
+		t_sum+=((tot-moldyn->t_m)*(tot-moldyn->t_m));
+		count+=1;
 	}
 
-	p2_m=p_sum/count;
-	k2_m=k_sum/count;
-	t2_m=t_sum/count;
+	moldyn->dp2_m=p_sum/count;
+	moldyn->dk2_m=k_sum/count;
+	moldyn->dt2_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);
+	printf("  averages   : %f %f %f\n",moldyn->k_m,
+	                                moldyn->p_m,
+	                                moldyn->t_m);
+	printf("  mean square: %f %f %f\n",moldyn->dk2_m,
+	                                   moldyn->dp2_m,
+	                                   moldyn->dt2_m);
 
 	close(fd);
 
 	return 0;
 }
+
+int get_heat_capacity(t_moldyn *moldyn) {
+
+	double temp2,mass,ighc;
+	int i;
+
+	/* (temperature average)^2 */
+	temp2=2.0*moldyn->k_m*EV/(3.0*K_BOLTZMANN);
+	printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",temp2);
+	temp2*=temp2;
+
+	/* total mass */
+	mass=0.0;
+	for(i=0;i<moldyn->count;i++)
+		mass+=moldyn->atom[i].mass;
+
+	/* ideal gas contribution */
+	ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
+	printf("  ideal gas contribution: %f\n",ighc/mass*KILOGRAM/JOULE);
+
+	moldyn->c_v_nvt=moldyn->dp2_m*moldyn->count*moldyn->count*EV/(K_BOLTZMANN*temp2)+ighc;
+	moldyn->c_v_nvt/=mass;
+	moldyn->c_v_nve=ighc/(1.0-(moldyn->dp2_m*moldyn->count*moldyn->count*EV/(ighc*K_BOLTZMANN*temp2)));
+	moldyn->c_v_nve/=mass;
+
+	printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
+	printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
+
+	return 0;
+}