X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=52a716e2bdfb6552a1ee13bb22f0e5eb60e29380;hb=dfbe75140f5a113f898809c529a882034507f6d1;hp=af0132dff157564f46f593c3280700b045a1810c;hpb=a9fbc66448c52bc4138176739b33d17ba86b7eae;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index af0132d..52a716e 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -950,6 +950,13 @@ printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[
 					}
 					
 				} while(list_next(this)!=L_NO_NEXT_ELEMENT);
+		
+				/* 2bp post function */
+				if(moldyn->func2b_post)
+					mlodyn->func2b_post(moldyn,
+					                    &(itom[i]),
+					                    jtom,bc_ij);
+
 			}
 		}
 	}
@@ -1072,6 +1079,21 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
 	p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
 	p->mu_m=0.5*(p->mu[0]+p->mu[1]);
 
+	printf("[moldyn] tersoff mult parameter info:\n");
+	printf("  S (m)  | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed);
+	printf("  R (m)  | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed);
+	printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
+	printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
+	printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
+	                                  p->lambda_m);
+	printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
+	printf("  beta   | %.10f | %.10f\n",p->beta[0],p->beta[1]);
+	printf("  n      | %f | %f\n",p->n[0],p->n[1]);
+	printf("  c      | %f | %f\n",p->c[0],p->c[1]);
+	printf("  d      | %f | %f\n",p->d[0],p->d[1]);
+	printf("  h      | %f | %f\n",p->h[0],p->h[1]);
+	printf("  chi    | %f \n",p->chi);
+
 	return 0;
 }
 
@@ -1139,18 +1161,17 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	d_ij=v3_norm(&dist_ij);
 
 	/* save for use in 3bp */
-	exchange->dist_ij=dist_ij; /* <- needed ? */
 	exchange->d_ij=d_ij;
+	exchange->dist_ij=dist_ij;
+	exchange->d_ij2=d_ij*d_ij;
 
 	/* constants */
 	if(num==aj->bnum) {
 		S=params->S[num];
 		R=params->R[num];
 		A=params->A[num];
+		B=params->B[num];
 		lambda=params->lambda[num];
-		/* more constants depending of atoms i and j, needed in 3bp */
-		params->exchange.B=&(params->B[num]);
-		params->exchange.mu=&(params->mu[num]);
 		mu=params->mu[num];
 		params->exchange.chi=1.0;
 	}
@@ -1158,10 +1179,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		S=params->Smixed;
 		R=params->Rmixed;
 		A=params->Amixed;
+		B=params->Bmixed;
 		lambda=params->lambda_m;
-		/* more constants depending of atoms i and j, needed in 3bp */
-		params->exchange.B=&(params->Bmixed);
-		params->exchange.mu=&(params->mu_m);
 		mu=params->mu_m;
 		params->exchange.chi=params->chi;
 	}
@@ -1203,6 +1222,57 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	/* enable the run of 3bp function */
 	exchange->run3bp=1;
 
+	/* reset 3bp sums */
+	exchange->3bp_sum1=0.0;
+	exchange->3bp_sum2=0.0;
+
+	return 0;
+}
+
+/* tersoff 2 body post part */
+
+int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+
+	/* here we have to allow for the 3bp sums */
+
+	t_tersoff_mult_params *params;
+	t_tersoff_exchange *exchange;
+
+	t_3dvec force,temp,*db_ij;
+	double db_ij_scale1,db_ij_scale2;
+	double b_ij;
+	double f_c,df_c,f_a,df_a;
+
+	params=moldyn->pot2b_params;
+	exchange=&(moldyn->exchange);
+
+	db_ij=&(exchange->db_ij);
+	f_c=exchange->f_c;
+	df_c=exchange->df_c;
+	f_a=exchange->f_a;
+	df_a=exchange->df_a;
+
+	db_ij_scale1=(1+betan*3bp_sum1);
+	db_ij_scale2=(n*betan*3bp_sum2);
+	help=pow(db_ij_scale1,-1.0/(2*n)-1);
+	b_ij=chi*db_ij_scale1*help;
+	db_ij_scale1=-chi/(2*n)*help;
+
+	v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
+	v3_scale(db_ij,db_ij,f_a);
+
+	v3_scale(&temp,dist_ij,b_ij*df_a);
+
+	v3_add(&force,&temp,db_ij);
+	v3_scale(&force,&force,f_c);
+
+	v3_scale(&temp,&dist_ij,f_a*b_ij*df_c);
+
+	/* add energy of 3bp sum */
+	moldyn->energy+=(0.5*f_c*b_ij*f_a);
+	/* add force of 3bp calculation */	
+	v3_add(&(ai->f),&temp,&force);
+
 	return 0;
 }
 
@@ -1245,9 +1315,13 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
 	d_ij=exchange->d_ij;
 	d_ij2=exchange->d_ij2;
+	dist_ij=exchange->dist_ij;
 
 	f_a=params->exchange.f_a;
 	df_a=params->exchange.df_a;
+
+	f_c=exchange->f_c;
+	df_c=exchange->df_c;
 	
 	/* d_ij is <= S, as we didn't return so far! */
 
@@ -1261,7 +1335,6 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	 *
 	 */
 
-	
 	v3_sub(&dist_ik,&(ai->r),&(ak->r));
 	if(bc) check_per_bound(moldyn,&dist_ik);
 	d_ik=v3_norm(&dist_ik);
@@ -1277,11 +1350,11 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	}
 
 	/* calc of f_c_ik */
-	if(d_ik>S)
-		return 0;
-
-	if(d_ik<R) {
-		/* f_c_ik = 1, df_c_ik = 0 */
+	if(d_ik>S) {
+		f_c_ik=0.0;
+		df_c_ik=0.0;
+	}
+	else if(d_ik<R) {
 		f_c_ik=1.0;
 		df_c_ik=0.0;
 	}
@@ -1302,6 +1375,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	c=*(exchange->c);
 	d=*(exchange->d);
 	h=*(exchange->h);
+	chi=exchange->chi;
 	c2=exchange->c2;
 	d2=exchange->d2;
 	c2d2=exchange->c2d2;
@@ -1309,6 +1383,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
 	denom=2*d_ij*d_ik;
 	cos_theta=numer/denom;
+	//cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 	sin_theta=sqrt(1.0-(cos_theta*cos_theta));
 	theta=acos(cos_theta);
 	d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom);
@@ -1319,17 +1394,28 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
 	h_cos=(h-cos_theta);
 	h_cos2=h_cos*h_cos;
-	d2_h_cos2=d2-h_cos2;
+	d2_h_cos2=d2+h_cos2;
 
 	/* some usefull expressions */
-	frac1=c2/(d2-h_cos2);
+	frac1=c2/(d2_h_cos2);
 	bracket1=1+c2d2-frac1;
-	bracket2=f_c_ik*bracket1;
-	bracket2_n_1=pow(bracket2,n-1.0);
-	bracket2_n=bracket2_n_1*bracket2;
-	bracket3=1+betan*bracket2_n;
+	if(f_c_ik==0.0) {
+		bracket2=0.0;
+		bracket2_n_1=0.0;
+		bracket2_n=0.0;
+		bracket3=1.0;
+		printf("Foo -> 0: ");
+	}
+	else {
+		bracket2=f_c_ik*bracket1;
+		bracket2_n_1=pow(bracket2,n-1.0);
+		bracket2_n=bracket2_n_1*bracket2;
+		bracket3=1.0+betan*bracket2_n;
+		printf("Foo -> 1: ");
+	}
 	bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
 	bracket3_pow=bracket3_pow_1*bracket3;
+printf("%.15f %.15f %.15f\n",bracket2_n_1,bracket2_n);
 
 	/* now go on with calc of b_ij and derivation of b_ij */
 	b_ij=chi*bracket3_pow;
@@ -1360,8 +1446,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
 	/* add forces */
 	v3_add(&(ai->f),&(ai->f),&force);
-	/* energy is 0.5 f_r f_c, but we will sum it up twice ... */
-	moldyn->energy+=(0.25*f_a*b_ij*f_c);
+	/* energy is 0.5 f_r f_c */
+	moldyn->energy+=(0.5*f_a*b_ij*f_c);
 				
 	return 0;
 }