supply c2, d2, c2/d2

[physik/posic.git] / potentials / albe_fast.c

diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c
index 3b5f681..6331064 100644 (file)
--- 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) {