tests ...

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index af0132d..f504e68 100644 (file)
--- 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;
 }