From: hackbard <hackbard>
Date: Mon, 18 Dec 2006 09:30:22 +0000 (+0000)
Subject: foo
X-Git-Url: https://www.hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=commitdiff_plain;h=2d562a1946322ae5b70bdce180321f32adbeab62

foo
---

diff --git a/Makefile b/Makefile
index 42ff2bd..39457f7 100644
--- a/Makefile
+++ b/Makefile
@@ -2,8 +2,9 @@ CC=gcc-3.4
 #CC=gcc
 CFLAGS=-Wall
 CFLAGS+=-O3
-CFLAGS+=-g
 CFLAGS+=-ffloat-store
+CFLAGS+=-g
+CFLAGS+=-DDEBUG
 LDFLAGS=-lm
 
 OBJS=visual/visual.o random/random.o
diff --git a/moldyn.c b/moldyn.c
index db91859..faa6b7c 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -541,13 +541,13 @@ int scale_volume(t_moldyn *moldyn) {
 	}
 
 	/* just a guess so far ... */
-	v=sqrt(virial.xx*virial.xx+virial.yy*virial.yy+virial.zz+virial.zz);
+	v=virial.xx+virial.yy+virial.zz;
 
 printf("%f\n",v);
 	/* get pressure from virial */
-	moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*v;
+	moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+ONE_THIRD*v;
 	moldyn->p/=moldyn->volume;
-printf("%f\n",moldyn->p/(ATM));
+printf("%f | %f\n",moldyn->p/(ATM),moldyn->p_ref/ATM);
 
 	/* scale factor */
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN)
@@ -1427,6 +1427,15 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	/* add forces of 2bp (ij, ji) contribution
 	 * dVij = dVji and we sum up both: no 1/2) */
 	v3_add(&(ai->f),&(ai->f),&force);
+
+	/* virial */
+	ai->virial.xx-=force.x*dist_ij.x;
+	ai->virial.yy-=force.y*dist_ij.y;
+	ai->virial.zz-=force.z*dist_ij.z;
+	ai->virial.xy-=force.x*dist_ij.y;
+	ai->virial.xz-=force.x*dist_ij.z;
+	ai->virial.yz-=force.y*dist_ij.z;
+
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
 	printf("dVij, dVji (2bp) contrib:\n");
@@ -1526,6 +1535,15 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 	/* add force */
 	v3_add(&(ai->f),&(ai->f),&force);
+
+	/* virial */
+	ai->virial.xx-=force.x*dist_ij->x;
+	ai->virial.yy-=force.y*dist_ij->y;
+	ai->virial.zz-=force.z*dist_ij->z;
+	ai->virial.xy-=force.x*dist_ij->y;
+	ai->virial.xz-=force.x*dist_ij->z;
+	ai->virial.yz-=force.y*dist_ij->z;
+
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
 	printf("dVij (3bp) contrib:\n");
@@ -1562,6 +1580,15 @@ if(ai==&(moldyn->atom[0])) {
 
 	/* add force */
 	v3_add(&(ai->f),&(ai->f),&force);
+
+	/* virial - plus sign, as dist_ij = - dist_ji - (really??) */
+	ai->virial.xx+=force.x*dist_ij->x;
+	ai->virial.yy+=force.y*dist_ij->y;
+	ai->virial.zz+=force.z*dist_ij->z;
+	ai->virial.xy+=force.x*dist_ij->y;
+	ai->virial.xz+=force.x*dist_ij->z;
+	ai->virial.yz+=force.y*dist_ij->z;
+
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
 	printf("dVji (3bp) contrib:\n");
@@ -1830,6 +1857,15 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 		v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
 		v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
 		                                  /* scaled with 0.5 ^ */
+
+		/* virial */
+		ai->virial.xx-=temp2.x*dist_jk.x;
+		ai->virial.yy-=temp2.y*dist_jk.y;
+		ai->virial.zz-=temp2.z*dist_jk.z;
+		ai->virial.xy-=temp2.x*dist_jk.y;
+		ai->virial.xz-=temp2.x*dist_jk.z;
+		ai->virial.yz-=temp2.y*dist_jk.z;
+
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
 	printf("dVjk (3bp) contrib:\n");
diff --git a/sic.c b/sic.c
index c4355b3..aabb79c 100644
--- a/sic.c
+++ b/sic.c
@@ -112,26 +112,26 @@ int main(int argc,char **argv) {
 
 	/* create the lattice / place atoms */
 	printf("[sic] creating atoms\n");
-	//create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
-	//               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
-	//               0,5,5,5);
-	//moldyn_bc_check(&md);
+	create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+	               0,5,5,5);
+	moldyn_bc_check(&md);
 
 	/* testing configuration */
-	r.x=2.8/2;	v.x=0;
-	r.y=0;		v.y=0;
-	r.z=0;		v.z=0;
-	add_atom(&md,SI,M_SI,0,
-	           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,//|ATOM_ATTR_HB,
+	//r.x=2.8/2;	v.x=0;
+	//r.y=0;		v.y=0;
+	//r.z=0;		v.z=0;
+	//add_atom(&md,SI,M_SI,0,
+	//           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 	//           ATOM_ATTR_2BP,
-	           &r,&v);
-	r.x=-2.8/2;	v.x=0;
-	r.y=0;		v.y=0;
-	r.z=0;		v.z=0;
-	add_atom(&md,SI,M_SI,0,
-	           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,//|ATOM_ATTR_HB,
+	//           &r,&v);
+	//r.x=-2.8/2;	v.x=0;
+	//r.y=0;		v.y=0;
+	//r.z=0;		v.z=0;
+	//add_atom(&md,SI,M_SI,0,
+	//           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 	//           ATOM_ATTR_2BP,
-	           &r,&v);
+	//           &r,&v);
 
 	/* setting a nearest neighbour distance for the moldyn checks */
 	set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
@@ -153,7 +153,7 @@ int main(int argc,char **argv) {
 	printf("[sic] set p/t scaling\n");
 	//set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
 	//                 T_SCALE_BERENDSEN,100.0);
-	//set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
+	set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
 	
 	/* initial thermal fluctuations of particles (in equilibrium) */
 	printf("[sic] thermal init\n");
@@ -161,13 +161,13 @@ int main(int argc,char **argv) {
 
 	/* create the simulation schedule */
 	printf("[sic] adding schedule\n");
-	moldyn_add_schedule(&md,10000,.1);
+	moldyn_add_schedule(&md,100,1.0);
 
 	/* activate logging */
 	printf("[sic] activate logging\n");
 	moldyn_set_log_dir(&md,argv[1]);
 	moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
-	moldyn_set_log(&md,VISUAL_STEP,20);
+	moldyn_set_log(&md,VISUAL_STEP,1);
 
 	/*
 	 * let's do the actual md algorithm now