From a198ee2eb93c0f14bdaea14722f4d98e5890935f Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 1 Sep 2008 14:28:15 +0200
Subject: [PATCH] supply c2, d2, c2/d2

---
 potentials/albe.c      |  9 +++++++++
 potentials/albe.h      |  6 ++++++
 potentials/albe_fast.c | 42 ++++++++++++++++++++++++++++++------------
 3 files changed, 45 insertions(+), 12 deletions(-)

diff --git a/potentials/albe.c b/potentials/albe.c
index 02c254f..08ce4fa 100644
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -96,6 +96,15 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 	p->S2[0]=p->S[0]*p->S[0];
 	p->S2[1]=p->S[1]*p->S[1];
 	p->S2mixed=p->Smixed*p->Smixed;
+	p->c2[0]=p->c[0]*p->c[0];
+	p->c2[1]=p->c[1]*p->c[1];
+	p->c2_mixed=p->c_mixed*p->c_mixed;
+	p->d2[0]=p->d[0]*p->d[0];
+	p->d2[1]=p->d[1]*p->d[1];
+	p->d2_mixed=p->d_mixed*p->d_mixed;
+	p->c2d2[0]=p->c2[0]/p->d2[0];
+	p->c2d2[1]=p->c2[1]/p->d2[1];
+	p->c2d2_m=p->c2_mixed/p->d2_mixed;
 
 	printf("[albe] mult parameter info:\n");
 	printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
diff --git a/potentials/albe.h b/potentials/albe.h
index 15386e7..87f45a0 100644
--- a/potentials/albe.h
+++ b/potentials/albe.h
@@ -66,10 +66,16 @@ typedef struct s_albe_mult_params {
 	double gamma_m;
 	double c[2];
 	double c_mixed;
+	double c2[2];
+	double c2_mixed;
 	double d[2];
 	double d_mixed;
+	double d2[2];
+	double d2_mixed;
 	double h[2];
 	double h_mixed;
+	double c2d2[2];
+	double c2d2_m;
 
 	t_albe_exchange exchange;	/* exchange between 2bp and 3bp calc */
 } t_albe_mult_params;
diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c
index 3b5f681..6331064 100644
--- a/potentials/albe_fast.c
+++ b/potentials/albe_fast.c
@@ -90,6 +90,14 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
 	t_3dvec dcosdrj,dcosdrk;
 	t_3dvec tmp;
 
+	// even more ...
+	double gamma_i;
+	double c_i;
+	double d_i;
+	double h_i;
+	double ci2;
+	double di2;
+	double ci2di2;
 
 	count=moldyn->count;
 	itom=moldyn->atom;
@@ -267,24 +275,27 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
 		Sk=params->S[brand_i];
 		Sk2=params->S2[brand_i];
 		/* albe needs i,k depending c,d,h and gamma values */
-		exchange->gamma_i=&(params->gamma[brand_i]);
-		exchange->c_i=&(params->c[brand_i]);
-		exchange->d_i=&(params->d[brand_i]);
-		exchange->h_i=&(params->h[brand_i]);
+		gamma_i=params->gamma[brand_i];
+		c_i=params->c[brand_i];
+		d_i=params->d[brand_i];
+		h_i=params->h[brand_i];
+		ci2=params->c2[brand_i];
+		di2=params->d2[brand_i];
+		ci2di2=params->c2d2[brand_i];
 	}
 	else {
 		Rk=params->Rmixed;
 		Sk=params->Smixed;
 		Sk2=params->S2mixed;
 		/* albe needs i,k depending c,d,h and gamma values */
-		exchange->gamma_i=&(params->gamma_m);
-		exchange->c_i=&(params->c_mixed);
-		exchange->d_i=&(params->d_mixed);
-		exchange->h_i=&(params->h_mixed);
+		gamma_i=params->gamma_m;
+		c_i=params->c_mixed;
+		d_i=params->d_mixed;
+		h_i=params->h_mixed;
+		ci2=params->c2_mixed;
+		di2=params->d2_mixed;
+		ci2di2=params->c2d2_m;
 	}
-	exchange->ci2=*(exchange->c_i)**(exchange->c_i);
-	exchange->di2=*(exchange->d_i)**(exchange->d_i);
-	exchange->ci2di2=exchange->ci2/exchange->di2;
 
 	/* dist_ik, d_ik2 */
 	v3_sub(&dist_ik,&(ktom->r),&(ai->r));
@@ -307,12 +318,19 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
 	/* cos theta */
 	cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 
-	/* g_ijk */
+	/* g_ijk 
 	h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
 	d2_h_cos2=exchange->di2+(h_cos*h_cos);
 	frac=exchange->ci2/d2_h_cos2;
 	g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 	dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
+	*/
+
+	h_cos=h_i+cos_theta; // + in albe formalism
+	d2_h_cos2=di2+(h_cos*h_cos);
+	frac=ci2/d2_h_cos2;
+	g=gamma_i*(1.0+ci2di2-frac);
+	dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
 
 	/* zeta sum += f_c_ik * g_ijk */
 	if(d_ik<=Rk) {
-- 
2.20.1