From a5a3025f6e2767714c6488a5bf1eec2ac9fac848 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Fri, 24 Nov 2006 12:11:43 +0000
Subject: [PATCH] finished basic tersoff stuff, now start fixing ... !

---
 moldyn.c | 68 ++++++++++++++++++++++++++++++++++++++++++++++++--------
 moldyn.h |  1 +
 2 files changed, 60 insertions(+), 9 deletions(-)

diff --git a/moldyn.c b/moldyn.c
index 830b8f0..14ae320 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -925,6 +925,7 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
 	exchange->d=&(params->d[num]);
 	exchange->h=&(params->h[num]);
 
+	exchange->betan=pow(*(exchange->beta),*(exchange->n));
 	exchange->c2=params->c[num]*params->c[num];
 	exchange->d2=params->d[num]*params->d[num];
 	exchange->c2d2=exchange->c2/exchange->d2;
@@ -1035,11 +1036,12 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	t_tersoff_mult_params *params;
 	t_tersoff_exchange *exchange;
 	t_3dvec dist_ij,dist_ik,dist_jk;
-	t_3dvec db_ij,temp,force;
+	t_3dvec temp,force;
 	double R,S,s_r;
 	double d_ij,d_ik,d_jk;
 	double f_c,df_c,b_ij,f_a,df_a;
-	double B,mu;
+	double n,c,d,h,neta,betan,betan_1;
+	double theta,cos_theta,sin_theta;
 	int num;
 
 	params=moldyn->pot_params;
@@ -1058,9 +1060,6 @@ 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;
 
-	B=*(params->exchange.B);
-	mu=*(params->exchange.mu);
-
 	f_a=params->exchange.f_a;
 	df_a=params->exchange.df_a;
 	
@@ -1111,17 +1110,68 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 	if(bc) check_per_bound(moldyn,&dist_jk);
 	d_jk=v3_norm(&dist_jk);
 
-
-	// GO ON HERE !!!
-
-	cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/(2*d_ij*d_ik);
+	beta=*(exchange->beta);
+	betan=exchange->betan;
+	n=*(exchange->n);
+	c=*(exchange->c);
+	d=*(exchange->d);
+	h=*(exchange->h);
+	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;
 	sin_theta=sqrt(1.0-(cos_theta*cos_theta));
 	theta=arccos(cos_theta);
+	d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom);
+	d_theta1=2*denom-numer*2*d_ik/d_ij;
+	d_theta2=2*denom-numer*2*d_ij/d_ik;
+	d_theta1*=d_theta;
+	d_theta2*=d_theta;
 
 	h_cos=(h-cos_theta);
 	h_cos2=h_cos*h_cos;
 	d2_h_cos2=d2-h_cos2;
 
+	/* 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;
+	bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
+	bracket3_pow=bracket3_pow_1*bracket3;
+
+	/* now go on with calc of b_ij and derivation of b_ij */
+	b_ij=chi*bracket3_pow;
+
+	/* derivation of theta */
+	v3_scale(&force,&dist_ij,d1_theta);
+	v3_scale(&temp,&dist_ik,d_theta2);
+	v3_add(&force,&force,&temp);
+
+	/* part 1 of derivation of b_ij */
+	v3_scale(&force,sin_theta*2*h_cos*f_c_ik*frac1);
+
+	/* part 2 of derivation of b_ij */
+	v3_scale(&temp,&dist_ik,df_c_ik*bracket1);
+
+	/* sum up and scale ... */
+	v3_add(&temp,&temp,&force);
+	scale=bracket2_n_1*n*betan*(1+betan*bracket3_pow_1)*chi*(1.0/(2.0*n));
+	v3_scale(&temp,&temp,scale);
+
+	/* now construct an energy and a force out of that */
+	v3_scale(&temp,&temp,f_a);
+	v3_scale(&force,&dist_ij,df_a*b_ij);
+	v3_add(&temp,&temp,&force);
+	v3_scale(&temp,&temp,f_c);
+	v3_scale(&force,&dist_ij,df_c*b_ij*f_a);
+	v3_add(&force,&force,&temp);
+
 	/* add forces */
 	v3_add(&(ai->f),&(ai->f),&force);
 	/* energy is 0.5 f_r f_c, but we will sum it up twice ... */
diff --git a/moldyn.h b/moldyn.h
index 9c63307..0a6b4f2 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -145,6 +145,7 @@ typedef struct s_tersoff_exchange {
 	double c2;
 	double d2;
 	double c2d2;
+	double betan;
 
 	u8 run3bp;
 } t_tersoff_exchange;
-- 
2.20.1