From 000e702c884e337f27ae66b6506ff9f9b613f8a8 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Sun, 3 Dec 2006 01:54:29 +0000
Subject: [PATCH] tests ...

---
 math/math.c |  5 +++++
 math/math.h |  1 +
 moldyn.c    | 49 +++++++++++++++++++++++++++++--------------------
 moldyn.h    |  3 ---
 4 files changed, 35 insertions(+), 23 deletions(-)

diff --git a/math/math.c b/math/math.c
index 6f4345e..54d79fd 100644
--- a/math/math.c
+++ b/math/math.c
@@ -66,6 +66,11 @@ int v3_cmp(t_3dvec *a,t_3dvec *b) {
 	return(memcmp(a,b,sizeof(t_3dvec)));
 }
 
+double v3_scalar_product(t_3dvec *a,t_3dvec *b) {
+
+	return(a->x*b->x+a->y*b->y+a->z*b->z);
+}
+
 double v3_absolute_square(t_3dvec *a) {
 
 	return(a->x*a->x+a->y*a->y+a->z*a->z);
diff --git a/math/math.h b/math/math.h
index 77100f3..8bbcc6c 100644
--- a/math/math.h
+++ b/math/math.h
@@ -25,6 +25,7 @@ int v3_zero(t_3dvec *vec);
 int v3_set(t_3dvec *vec,double *ptr);
 int v3_copy(t_3dvec *trg,t_3dvec *src);
 int v3_cmp(t_3dvec *a,t_3dvec *b);
+double v3_scalar_product(t_3dvec *a,t_3dvec *b);
 double v3_absolute_square(t_3dvec *a);
 double v3_norm(t_3dvec *a);
 int v3_per_bound(t_3dvec *a,t_3dvec *dim);
diff --git a/moldyn.c b/moldyn.c
index af0132d..f504e68 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -1139,18 +1139,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 +1157,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;
 	}
@@ -1261,7 +1258,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 +1273,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,13 +1298,15 @@ 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;
 
 	numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
 	denom=2*d_ij*d_ik;
-	cos_theta=numer/denom;
+	//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);
@@ -1324,11 +1322,22 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	/* some usefull expressions */
 	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);
+printf("THETA: %.15f %.15f\n",cos_theta,theta*180/(2*M_PI));
 	bracket3_pow=bracket3_pow_1*bracket3;
 
 	/* now go on with calc of b_ij and derivation of b_ij */
@@ -1360,8 +1369,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;
 }
diff --git a/moldyn.h b/moldyn.h
index 8b40d4c..813af68 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -182,9 +182,6 @@ typedef struct s_tersoff_exchange {
 
 	double chi;
 
-	double *B;
-	double *mu;
-
 	double *beta;
 	double *n;
 	double *c;
-- 
2.39.2