From: hackbard <hackbard>
Date: Sun, 17 Dec 2006 01:29:39 +0000 (+0000)
Subject: testing foo
X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=150c5f804a5088a01ef2bb1e934b5e65d10ac23e;p=physik%2Fposic.git

testing foo
---

diff --git a/moldyn.c b/moldyn.c
index 5aa6d87..db91859 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -519,7 +519,8 @@ int scale_volume(t_moldyn *moldyn) {
 
 	t_atom *atom;
 	t_3dvec *dim,*vdim;
-	double virial,scale;
+	double scale,v;
+	t_virial virial;
 	t_linkcell *lc;
 	int i;
 
@@ -528,12 +529,23 @@ int scale_volume(t_moldyn *moldyn) {
 	vdim=&(moldyn->vis.dim);
 	lc=&(moldyn->lc);
 
-	for(i=0;i<moldyn->count;i++)
-		virial+=v3_norm(&(atom[i].virial));
+	memset(&virial,0,sizeof(t_virial));
+
+	for(i=0;i<moldyn->count;i++) {
+		virial.xx+=atom[i].virial.xx;
+		virial.yy+=atom[i].virial.yy;
+		virial.zz+=atom[i].virial.zz;
+		virial.xy+=atom[i].virial.xy;
+		virial.xz+=atom[i].virial.xz;
+		virial.yz+=atom[i].virial.yz;
+	}
+
+	/* just a guess so far ... */
+	v=sqrt(virial.xx*virial.xx+virial.yy*virial.yy+virial.zz+virial.zz);
 
-printf("%f\n",virial);
+printf("%f\n",v);
 	/* get pressure from virial */
-	moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial;
+	moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*v;
 	moldyn->p/=moldyn->volume;
 printf("%f\n",moldyn->p/(ATM));
 
@@ -984,6 +996,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 	int i,j,k,count;
 	t_atom *itom,*jtom,*ktom;
+	t_virial *virial;
 	t_linkcell *lc;
 	t_list neighbour_i[27];
 	t_list neighbour_i2[27];
@@ -1005,7 +1018,13 @@ int potential_force_calc(t_moldyn *moldyn) {
 		v3_zero(&(itom[i].f));
 
 		/* reset viral of atom i */
-		v3_zero(&(itom[i].virial));
+		virial=&(itom[i].virial);
+		virial->xx=0.0;
+		virial->yy=0.0;
+		virial->zz=0.0;
+		virial->xy=0.0;
+		virial->xz=0.0;
+		virial->yz=0.0;
 
 		/* reset site energy */
 		itom[i].e=0.0;
@@ -1224,7 +1243,10 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
 
 	printf("[moldyn] tersoff parameter completion\n");
+	p->S2[0]=p->S[0]*p->S[0];
+	p->S2[1]=p->S[1]*p->S[1];
 	p->Smixed=sqrt(p->S[0]*p->S[1]);
+	p->S2mixed=p->Smixed*p->Smixed;
 	p->Rmixed=sqrt(p->R[0]*p->R[1]);
 	p->Amixed=sqrt(p->A[0]*p->A[1]);
 	p->Bmixed=sqrt(p->B[0]*p->B[1]);
@@ -1285,8 +1307,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	t_tersoff_mult_params *params;
 	t_tersoff_exchange *exchange;
 	t_3dvec dist_ij,force;
-	double d_ij;
-	double A,B,R,S,lambda,mu;
+	double d_ij,d_ij2;
+	double A,B,R,S,S2,lambda,mu;
 	double f_r,df_r;
 	double f_c,df_c;
 	int brand;
@@ -1317,18 +1339,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	 *
 	 */
 
-	/* dist_ij, d_ij */
-	v3_sub(&dist_ij,&(aj->r),&(ai->r));
-	if(bc) check_per_bound(moldyn,&dist_ij);
-	d_ij=v3_norm(&dist_ij);
-
-	/* save for use in 3bp */
-	exchange->d_ij=d_ij;
-	exchange->dist_ij=dist_ij;
-
 	/* constants */
 	if(brand==ai->brand) {
 		S=params->S[brand];
+		S2=params->S2[brand];
 		R=params->R[brand];
 		A=params->A[brand];
 		B=params->B[brand];
@@ -1338,6 +1352,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	}
 	else {
 		S=params->Smixed;
+		S2=params->S2mixed;
 		R=params->Rmixed;
 		A=params->Amixed;
 		B=params->Bmixed;
@@ -1346,10 +1361,23 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		params->exchange.chi=params->chi;
 	}
 
-	/* if d_ij > S => no force & potential energy contribution */
-	if(d_ij>S)
+	/* dist_ij, d_ij */
+	v3_sub(&dist_ij,&(aj->r),&(ai->r));
+	if(bc) check_per_bound(moldyn,&dist_ij);
+	d_ij2=v3_absolute_square(&dist_ij);
+
+	/* if d_ij2 > S2 => no force & potential energy contribution */
+	if(d_ij2>S2)
 		return 0;
 
+	/* now we will need the distance */
+	//d_ij=v3_norm(&dist_ij);
+	d_ij=sqrt(d_ij2);
+
+	/* save for use in 3bp */
+	exchange->d_ij=d_ij;
+	exchange->dist_ij=dist_ij;
+
 	/* more constants */
 	exchange->beta_j=&(params->beta[brand]);
 	exchange->n_j=&(params->n[brand]);
@@ -1389,7 +1417,6 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		s_r=S-R;
 		arg=M_PI*(d_ij-R)/s_r;
 		f_c=0.5+0.5*cos(arg);
-		//df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));	/* MARK! */
 		df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 		/* two body contribution (ij, ji) */
 		v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
@@ -1685,7 +1712,6 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 			s_r=S-R;
 			arg=M_PI*(d_ik-R)/s_r;
 			f_c_ik=0.5+0.5*cos(arg);
-			//df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* MARK */
 			df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 
 			/* zeta_ij */
diff --git a/moldyn.h b/moldyn.h
index 38e4999..dd94b0e 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -22,12 +22,22 @@
 /* general */
 typedef unsigned char u8;
 
+/* virial */
+typedef struct s_virial {
+	double xx;	/*			| xx	 xy	xz |	*/
+	double yy;	/*	V	=	| yx	 yy	yz |	*/
+	double zz;	/*			| zx	 zy	zz |	*/
+	double xy;	/*						*/
+	double xz;	/*	with:	xy=yx, xz=zx, yz=zy		*/
+	double yz;	/*						*/
+} t_virial;
+
 /* the atom of the md simulation */
 typedef struct s_atom {
 	t_3dvec r;		/* position */
 	t_3dvec v;		/* velocity */
 	t_3dvec f;		/* force */
-	t_3dvec virial;		/* virial (v_xx, v_yy, v_zz) */
+	t_virial virial;	/* virial */
 	double e;		/* site energy */
 	int element;		/* number of element in pse */
 	double mass;		/* atom mass */
@@ -222,8 +232,10 @@ typedef struct s_tersoff_exchange {
 /* tersoff multi (2!) potential parameters */
 typedef struct s_tersoff_mult_params {
 	double S[2];		/* tersoff cutoff radii */
+	double S2[2];		/* tersoff cutoff radii squared */
 	double R[2];		/* tersoff cutoff radii */
 	double Smixed;		/* mixed S radius */
+	double S2mixed;		/* mixed S radius squared */
 	double Rmixed;		/* mixed R radius */
 	double A[2];		/* factor of tersoff attractive part */
 	double B[2];		/* factor of tersoff repulsive part */
@@ -310,6 +322,7 @@ typedef struct s_tersoff_mult_params {
 #define SI			0x0e
 #define LC_SI			(0.543105e-9*METER)		/* A */
 #define M_SI			28.08553			/* amu */
+
 #define LJ_SIGMA_SI		((0.25*sqrt(3.0)*LC_SI)/1.122462)	/* A */
 #define LJ_EPSILON_SI		(2.1678*EV)				/* NA */
 
diff --git a/sic.c b/sic.c
index 35931e0..c4355b3 100644
--- a/sic.c
+++ b/sic.c
@@ -112,24 +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.95/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.95/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.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);
 
 	/* setting a nearest neighbour distance for the moldyn checks */
 	set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
@@ -151,25 +153,21 @@ 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");
-	thermal_init(&md,TRUE);
+	//thermal_init(&md,TRUE);
 
 	/* create the simulation schedule */
 	printf("[sic] adding schedule\n");
-	moldyn_add_schedule(&md,20000,.1);
-	moldyn_add_schedule(&md,10000,.2);
-	moldyn_add_schedule(&md,6667,.3);
-	moldyn_add_schedule(&md,5000,.4);
-	moldyn_add_schedule(&md,4001,.5);
+	moldyn_add_schedule(&md,10000,.1);
 
 	/* activate logging */
 	printf("[sic] activate logging\n");
 	moldyn_set_log_dir(&md,argv[1]);
-	moldyn_set_log(&md,LOG_TOTAL_ENERGY,100);
-	moldyn_set_log(&md,VISUAL_STEP,100);
+	moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
+	moldyn_set_log(&md,VISUAL_STEP,20);
 
 	/*
 	 * let's do the actual md algorithm now