From b776d78d5fe14e550e2b4a51a2b837742ed6f850 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Fri, 29 Jun 2007 12:08:29 +0000
Subject: [PATCH] still playing around with heat capacity ...

---
 moldyn.c | 38 +++++++++++---------------------------
 sic.c    | 19 +++++++++++++------
 2 files changed, 24 insertions(+), 33 deletions(-)

diff --git a/moldyn.c b/moldyn.c
index c130ef4..4c43f2e 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -1838,6 +1838,7 @@ int calc_fluctuations(double start,double end,t_moldyn *moldyn) {
 	int count,ret;
 	double time,pot,kin,tot;
 	double p_sum,k_sum,t_sum;
+	double p2_sum,k2_sum,t2_sum;
 	char buf[64];
 	char file[128+7];
 
@@ -1850,7 +1851,7 @@ int calc_fluctuations(double start,double end,t_moldyn *moldyn) {
 		return fd;
 	}
 
-	/* first calc the averages */
+	/* calc the averages of A and A^2 */
 	p_sum=0.0;
 	k_sum=0.0;
 	t_sum=0.0;
@@ -1865,38 +1866,21 @@ int calc_fluctuations(double start,double end,t_moldyn *moldyn) {
 		p_sum+=pot;
 		k_sum+=kin;
 		t_sum+=tot;
+		p2_sum+=(pot*pot);
+		k2_sum+=(kin*kin);
+		t2_sum+=(tot*tot);
 		count+=1;
 	}
 
-	moldyn->p_m=p_sum/count;
+	/* averages */
 	moldyn->k_m=k_sum/count;
+	moldyn->p_m=p_sum/count;
 	moldyn->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=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;
-		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;
-	}
-
-	moldyn->dp2_m=p_sum/count;
-	moldyn->dk2_m=k_sum/count;
-	moldyn->dt2_m=t_sum/count;
+	/* rms */
+	moldyn->dk2_m=k2_sum/count-moldyn->k_m*moldyn->k_m;
+	moldyn->dp2_m=p2_sum/count-moldyn->p_m*moldyn->p_m;
+	moldyn->dt2_m=t2_sum/count-moldyn->t_m*moldyn->t_m;
 
 	printf("  averages   : %f %f %f\n",moldyn->k_m,
 	                                moldyn->p_m,
diff --git a/sic.c b/sic.c
index d50c35e..07f32bc 100644
--- a/sic.c
+++ b/sic.c
@@ -209,19 +209,25 @@ int main(int argc,char **argv) {
 	albe_mult_complete_params(&ap);
 
 	/* set (initial) dimensions of simulation volume */
+#ifdef ALBLE
 	set_dim(&md,6*LC_SI_ALBE,6*LC_SI_ALBE,6*LC_SI_ALBE,TRUE);
-	//set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
 	//set_dim(&md,6*LC_C_ALBE,6*LC_C_ALBE,6*LC_C_ALBE,TRUE);
-	//set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE);
 	//set_dim(&md,6*LC_SIC_ALBE,6*LC_SIC_ALBE,6*LC_SIC_ALBE,TRUE);
+#else
+	set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+	//set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE);
+#endif
 
 	/* set periodic boundary conditions in all directions */
 	set_pbc(&md,TRUE,TRUE,TRUE);
 
 	/* create the lattice / place atoms */
-	//create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
-	//create_lattice(&md,DIAMOND,LC_C_ALBE,C,M_C,
+#ifdef ALBLE
 	create_lattice(&md,DIAMOND,LC_SI_ALBE,SI,M_SI,
+#else
+	create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+	//create_lattice(&md,DIAMOND,LC_C_ALBE,C,M_C,
+#endif
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 	//               ATOM_ATTR_2BP|ATOM_ATTR_HB,
 	               0,6,6,6,NULL);
@@ -278,6 +284,7 @@ int main(int argc,char **argv) {
 	set_pressure(&md,BAR);
 
 	/* set p/t scaling */
+	//set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
 	//set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,
 	//                 T_SCALE_BERENDSEN,100.0);
 	//set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
@@ -288,7 +295,7 @@ int main(int argc,char **argv) {
 
 	/* create the simulation schedule */
 	/* initial configuration */
-	moldyn_add_schedule(&md,10000,1.0);
+	moldyn_add_schedule(&md,3000,1.0);
 	//moldyn_add_schedule(&md,1000,1.0);
 	//moldyn_add_schedule(&md,1000,1.0);
 	//moldyn_add_schedule(&md,1000,1.0);
@@ -331,7 +338,7 @@ return 0;
 	 */
 
 	/* response functions expressed by energy fluctuations */
-	calc_fluctuations(1000.0,9999.0,&md);
+	calc_fluctuations(1000.0,2999.0,&md);
 	get_heat_capacity(&md);
 
 	/* close */
-- 
2.20.1