From: hackbard <hackbard>
Date: Mon, 11 Dec 2006 16:03:33 +0000 (+0000)
Subject: bc check probs
X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=bb3e6f4ab3dfea9083e0d5b7d403ee6197f4041d;p=physik%2Fposic.git

bc check probs
---

diff --git a/moldyn.c b/moldyn.c
index 1f6391f..4ddc42e 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -781,12 +781,15 @@ int velocity_verlet(t_moldyn *moldyn) {
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
 		v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
+//if(i==0) printf("und dann %.20f\n",atom[i].r.x*1e10);
 		check_per_bound(moldyn,&(atom[i].r));
+//if(i==0) printf("und danach %.20f\n",atom[i].r.x);
 
 		/* velocities */
 		v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 		v3_add(&(atom[i].v),&(atom[i].v),&delta);
 	}
+moldyn_bc_check(moldyn);
 
 	/* neighbour list update */
 	link_cell_update(moldyn);
@@ -1189,7 +1192,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 	/* f_r_ij = f_r_ji, df_r_ij = df_r_ji */
 	f_r=A*exp(-lambda*d_ij);
-	df_r=-lambda*f_r/d_ij;
+	df_r=lambda*f_r/d_ij;
 
 	/* f_a, df_a calc (again, same for ij and ji) | save for later use! */
 	exchange->f_a=-B*exp(-mu*d_ij);
@@ -1289,6 +1292,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		moldyn->debug++;		/* just for debugging ... */
 		db=0.0;
 		b=chi;
+		v3_scale(&force,dist_ij,df_a*b*f_c);
 	}
 	else {
 		tmp=betaini*pow(zeta,ni-1.0);		/* beta^n * zeta^n-1 */
@@ -1300,13 +1304,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		v3_scale(&temp,dist_ij,df_a*b);
 		v3_add(&force,&force,&temp);
 		v3_scale(&force,&force,f_c);
-		v3_scale(&temp,dist_ij,df_c*b*f_a);
-		v3_add(&force,&force,&temp);
-		v3_scale(&force,&force,-0.5);
-
-		/* add force */
-		v3_add(&(ai->f),&(ai->f),&force);
 	}
+	v3_scale(&temp,dist_ij,df_c*b*f_a);
+	v3_add(&force,&force,&temp);
+	v3_scale(&force,&force,-0.5);
+
+	/* add force */
+	v3_add(&(ai->f),&(ai->f),&force);
 
 	/* add energy of 3bp sum */
 	moldyn->energy+=(0.5*f_c*b*f_a);
@@ -1314,8 +1318,9 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	/* dVji */
 	zeta=exchange->zeta_ji;
 	if(zeta==0.0) {
-		// do nothing ... (db=0, b=chi)
 		moldyn->debug++;
+		b=chi;
+		v3_scale(&force,dist_ij,df_a*b*f_c);
 	}
 	else {
 		tmp=betajnj*pow(zeta,nj-1.0);		/* beta^n * zeta^n-1 */
@@ -1327,13 +1332,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		v3_scale(&temp,dist_ij,df_a*b);
 		v3_add(&force,&force,&temp);
 		v3_scale(&force,&force,f_c);
-		v3_scale(&temp,dist_ij,0.5*df_c*b*f_a);
-		v3_add(&force,&force,&temp);
-		v3_scale(&force,&force,-0.5);
-
-		/* add force */
-		v3_sub(&(ai->f),&(ai->f),&force);
 	}
+	v3_scale(&temp,dist_ij,df_c*b*f_a);
+	v3_add(&force,&force,&temp);
+	v3_scale(&force,&force,-0.5);
+
+	/* add force */
+	v3_sub(&(ai->f),&(ai->f),&force);
 
 	return 0;
 }
@@ -1596,3 +1601,32 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	return 0;
 }
 
+
+/*
+ * debugging / critical check functions
+ */
+
+int moldyn_bc_check(t_moldyn *moldyn) {
+
+	t_atom *atom;
+	t_3dvec *dim;
+	int i;
+
+	atom=moldyn->atom;
+	dim=&(moldyn->dim);
+
+	for(i=0;i<moldyn->count;i++) {
+		if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
+			printf("FATAL: atom %d: x: %.20f (%.20f)\n",
+			       i,atom[i].r.x*1e10,dim->x/2*1e10);
+		if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
+			printf("FATAL: atom %d: y: %.20f (%.20f)\n",
+			       i,atom[i].r.y*1e10,dim->y/2*1e10);
+		if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
+			printf("FATAL: atom %d: z: %.20f (%.20f)\n",
+			       i,atom[i].r.z*1e10,dim->z/2*1e10);
+	}
+
+	return 0;
+}
+
diff --git a/moldyn.h b/moldyn.h
index d62ec47..534287b 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -388,4 +388,6 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
 
+int moldyn_bc_check(t_moldyn *moldyn);
+
 #endif
diff --git a/sic.c b/sic.c
index cbb81a3..0d3b798 100644
--- a/sic.c
+++ b/sic.c
@@ -108,20 +108,19 @@ 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,
-	  //             //ATOM_ATTR_2BP|ATOM_ATTR_HB,
-	   //            0,5,5,5);
+	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);
 
 	/* testing configuration */
-	r.x=-0.55*0.25*sqrt(3.0)*LC_SI;	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=+0.55*0.25*sqrt(3.0)*LC_SI;	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=-0.55*0.25*sqrt(3.0)*LC_SI;	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=+0.55*0.25*sqrt(3.0)*LC_SI;	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);
 
 	/* setting a nearest neighbour distance for the moldyn checks */
 	set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
@@ -141,12 +140,12 @@ int main(int argc,char **argv) {
 
 	/* create the simulation schedule */
 	printf("[sic] adding schedule\n");
-	moldyn_add_schedule(&md,1000,1.0e-15);
+	moldyn_add_schedule(&md,100,1.0e-15);
 
 	/* activate logging */
 	printf("[sic] activate logging\n");
-	moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",10);
-	moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",10);
+	moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",1);
+	moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",1);
 
 	/*
 	 * let's do the actual md algorithm now