From: hackbard <hackbard>
Date: Mon, 11 Dec 2006 10:19:02 +0000 (+0000)
Subject: security checkin
X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=1d32be0a3b47add316f22f506b8595d56ac09a34;p=physik%2Fposic.git

security checkin
---

diff --git a/moldyn.c b/moldyn.c
index 6736e45..1f6391f 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -681,7 +681,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	/* zero absolute time */
 	moldyn->time=0.0;
 
-	/* debugging, ignre */
+	/* debugging, ignore */
 	moldyn->debug=0;
 
 	/* executing the schedule */
@@ -741,7 +741,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 			if(!(i%v)) {
 				visual_atoms(&(moldyn->vis),moldyn->time,
 				             moldyn->atom,moldyn->count);
-				printf("\rsched: %d, steps: %d, theta: %d",
+				printf("\rsched: %d, steps: %d, debug: %d",
 				       sched,i,moldyn->debug);
 				fflush(stdout);
 			}
@@ -920,15 +920,14 @@ int potential_force_calc(t_moldyn *moldyn) {
 		
 				/* 2bp post function */
 				if(moldyn->func2b_post) {
-printf("DEBUG: vor 2bp post\n");
 					moldyn->func2b_post(moldyn,
 					                    &(itom[i]),
 					                    jtom,bc_ij);
-printf("DEBUG: nach 2bp post\n");
 				}
 
 			}
 		}
+//printf("DEBUG: %.15f \n",itom[i].f.x);
 	}
 
 	return 0;
@@ -1119,6 +1118,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	/* clear 3bp and 2bp post run */
 	exchange->run3bp=0;
 	exchange->run2bp_post=0;
+
+	/* reset S > r > R mark */
+	exchange->d_ij_between_rs=0;
 	
 	/*
 	 * calc of 2bp contribution of V_ij and dV_ij/ji
@@ -1208,6 +1210,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		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);
+		/* tell 3bp that S > r > R */
+		exchange->d_ij_between_rs=1;
 	}
 
 	/* add forces of 2bp (ij, ji) contribution
@@ -1281,40 +1285,55 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	
 	/* Vij and dVij */
 	zeta=exchange->zeta_ij;
-	tmp=betaini*pow(zeta,ni-1.0);		/* beta^n * zeta^n-1 */
-	b=(1+zeta*tmp);				/* 1 + beta^n * zeta^n */
-	db=chi*pow(b,-1.0/(2*ni)-1);		/* chi * (...)^(-1/2n - 1) */
-	b=db*b;					/* b_ij */
-	db*=-0.5*tmp;				/* db_ij */
-	v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
-	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);
+	if(zeta==0.0) {
+		moldyn->debug++;		/* just for debugging ... */
+		db=0.0;
+		b=chi;
+	}
+	else {
+		tmp=betaini*pow(zeta,ni-1.0);		/* beta^n * zeta^n-1 */
+		b=(1+zeta*tmp);				/* 1 + beta^n zeta^n */
+		db=chi*pow(b,-1.0/(2*ni)-1);		/* x(...)^(-1/2n - 1) */
+		b=db*b;					/* b_ij */
+		db*=-0.5*tmp;				/* db_ij */
+		v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
+		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);
+	}
 
 	/* add energy of 3bp sum */
 	moldyn->energy+=(0.5*f_c*b*f_a);
 
-	/* add force (sub, as F = - dVij) */
-	v3_sub(&(ai->f),&(ai->f),&force);
-
 	/* dVji */
 	zeta=exchange->zeta_ji;
-	tmp=betajnj*pow(zeta,nj-1.0);		/* beta^n * zeta^n-1 */
-	b=(1+zeta*tmp);				/* 1 + beta^n * zeta^n */
-	db=chi*pow(b,-1.0/(2*nj)-1);		/* chi * (...)^(-1/2n - 1) */
-	b=db*b;					/* b_ij */
-	db*=-0.5*tmp;				/* db_ij */
-	v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
-	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);
-
-	/* add force (sub, as F = - dVji) */
-	v3_sub(&(ai->f),&(ai->f),&force);
+	if(zeta==0.0) {
+		// do nothing ... (db=0, b=chi)
+		moldyn->debug++;
+	}
+	else {
+		tmp=betajnj*pow(zeta,nj-1.0);		/* beta^n * zeta^n-1 */
+		b=(1+zeta*tmp);				/* 1 + beta^n zeta^n */
+		db=chi*pow(b,-1.0/(2*nj)-1);		/* x(...)^(-1/2n - 1) */
+		b=db*b;					/* b_ij */
+		db*=-0.5*tmp;				/* db_ij */
+		v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
+		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);
+	}
 
 	return 0;
 }
@@ -1569,8 +1588,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 		/* betajnj * zeta_jk ^ nj-1 */
 		tmp=exchange->betajnj*pow(zeta,(n-1.0));
 		tmp=-chi/2.0*pow(1+tmp*zeta,-1.0/(2.0*n)-1)*tmp;
-		v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk);
+		v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
 		v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */
+		                                  /* scaled with 0.5 ^ */
 	}
 
 	return 0;
diff --git a/sic.c b/sic.c
index d141782..cbb81a3 100644
--- a/sic.c
+++ b/sic.c
@@ -26,6 +26,9 @@ int main(int argc,char **argv) {
 	/* misc parameters */
 	double tau;
 
+	/* testing location & velocity vector */
+	t_3dvec r,v;
+
 	/* values */
 	tau=1.0e-15;	/* delta t = 1 fs */
 
@@ -93,10 +96,11 @@ int main(int argc,char **argv) {
 	/* cutoff radius */
 	printf("[sic] setting cutoff radius\n");
 	set_cutoff(&md,TM_S_SI);
+	//set_cutoff(&md,2*LC_SI);
 
 	/* set (initial) dimensions of simulation volume */
 	printf("[sic] setting dimensions\n");
-	set_dim(&md,2*LC_SI,2*LC_SI,2*LC_SI,TRUE);
+	set_dim(&md,5*LC_SI,5*LC_SI,5*LC_SI,TRUE);
 
 	/* set periodic boundary conditions in all directions */
 	printf("[sic] setting periodic boundary conditions\n");
@@ -104,10 +108,20 @@ 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,2,2,2);
+	//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);
+
+	/* 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);
 
 	/* setting a nearest neighbour distance for the moldyn checks */
 	set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
@@ -127,12 +141,12 @@ int main(int argc,char **argv) {
 
 	/* create the simulation schedule */
 	printf("[sic] adding schedule\n");
-	moldyn_add_schedule(&md,100,1.0e-15);
+	moldyn_add_schedule(&md,1000,1.0e-15);
 
 	/* activate logging */
 	printf("[sic] activate logging\n");
-	moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",1);
-	moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",1);
+	moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",10);
+	moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",10);
 
 	/*
 	 * let's do the actual md algorithm now
diff --git a/visual/visual.c b/visual/visual.c
index 88bc6e6..687c561 100644
--- a/visual/visual.c
+++ b/visual/visual.c
@@ -93,7 +93,7 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
 	dprintf(v->fd,"set ambient 20\n");
 	dprintf(v->fd,"set specular on\n");
 	dprintf(v->fd,"set boundbox on\n");
-	//dprintf(v->fd,"label\n");
+	dprintf(v->fd,"label\n");
 	sprintf(file,"%s-%.15f.ppm",v->fb,time);
 	dprintf(v->fd,"write ppm %s\n",file);
 	dprintf(v->fd,"zap\n");