X-Git-Url: https://www.hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.c;h=1f6391fdaa23a1a03c424f7b9f4da2c113b8620c;hp=6736e457b28b9b97fe1b5dd89c9093a15d7e1c30;hb=1d32be0a3b47add316f22f506b8595d56ac09a34;hpb=b30d3d558354102cf5d0c06392961032b21bf9a9

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;