bugfix dg of dV_ji,jk

diff --git a/moldyn.c b/moldyn.c
index faa6b7c..c21f59f 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -1129,6 +1129,9 @@ int potential_force_calc(t_moldyn *moldyn) {
        }
 #ifdef DEBUG
 printf("\n\n");
+#endif
+#ifdef VDEBUG
+printf("\n\n");
 #endif
 
        return 0;
@@ -1376,6 +1379,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* save for use in 3bp */
        exchange->d_ij=d_ij;
+       exchange->d_ij2=d_ij2;
        exchange->dist_ij=dist_ij;
 
        /* more constants */
@@ -1443,6 +1447,14 @@ if(ai==&(moldyn->atom[0])) {
        printf("%f | %f\n",force.y,ai->f.y);
        printf("%f | %f\n",force.z,ai->f.z);
 }
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+       printf("dVij, dVji (2bp) contrib:\n");
+       printf("%f | %f\n",force.x*dist_ij.x,ai->virial.xx);
+       printf("%f | %f\n",force.y*dist_ij.y,ai->virial.yy);
+       printf("%f | %f\n",force.z*dist_ij.z,ai->virial.zz);
+}
 #endif
 
        /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
@@ -1551,6 +1563,14 @@ if(ai==&(moldyn->atom[0])) {
        printf("%f | %f\n",force.y,ai->f.y);
        printf("%f | %f\n",force.z,ai->f.z);
 }
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+       printf("dVij (3bp) contrib:\n");
+       printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx);
+       printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy);
+       printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz);
+}
 #endif
 
        /* add energy of 3bp sum */
@@ -1596,6 +1616,14 @@ if(ai==&(moldyn->atom[0])) {
        printf("%f | %f\n",force.y,ai->f.y);
        printf("%f | %f\n",force.z,ai->f.z);
 }
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+       printf("dVji (3bp) contrib:\n");
+       printf("%f | %f\n",force.x*dist_ij->x,ai->virial.xx);
+       printf("%f | %f\n",force.y*dist_ij->y,ai->virial.yy);
+       printf("%f | %f\n",force.z*dist_ij->z,ai->virial.zz);
+}
 #endif
 
        return 0;
@@ -1610,9 +1638,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        t_3dvec dist_ij,dist_ik,dist_jk;
        t_3dvec temp1,temp2;
        t_3dvec *dzeta;
-       double R,S,s_r;
+       double R,S,S2,s_r;
        double B,mu;
-       double d_ij,d_ik,d_jk;
+       double d_ij,d_ik,d_jk,d_ij2,d_ik2,d_jk2;
        double rr,dd;
        double f_c,df_c;
        double f_c_ik,df_c_ik,arg;
@@ -1658,6 +1686,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        /* dist_ij, d_ij - this is < S_ij ! */
        dist_ij=exchange->dist_ij;
        d_ij=exchange->d_ij;
+       d_ij2=exchange->d_ij2;
 
        /* f_c_ij, df_c_ij (same for ji) */
        f_c=exchange->f_c;
@@ -1672,21 +1701,26 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        /* dist_ik, d_ik */
        v3_sub(&dist_ik,&(ak->r),&(ai->r));
        if(bc) check_per_bound(moldyn,&dist_ik);
-       d_ik=v3_norm(&dist_ik);
+       d_ik2=v3_absolute_square(&dist_ik);
 
        /* ik constants */
        brand=ai->brand;
        if(brand==ak->brand) {
                R=params->R[brand];
                S=params->S[brand];
+               S2=params->S2[brand];
        }
        else {
                R=params->Rmixed;
                S=params->Smixed;
+               S2=params->S2mixed;
        }
 
        /* zeta_ij/dzeta_ij contribution only for d_ik < S */
-       if(d_ik<S) {
+       if(d_ik2<S2) {
+
+               /* now we need d_ik */
+               d_ik=sqrt(d_ik2);
 
                /* get constants_i from exchange data */
                n=*(exchange->n_i);
@@ -1704,8 +1738,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
                /* d_costheta */
                tmp=1.0/dd;
-               d_costheta1=cos_theta/(d_ij*d_ij)-tmp;
-               d_costheta2=cos_theta/(d_ik*d_ik)-tmp;
+               d_costheta1=cos_theta/d_ij2-tmp;
+               d_costheta2=cos_theta/d_ik2-tmp;
 
                /* some usefull values */
                h_cos=(h-cos_theta);
@@ -1757,13 +1791,14 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        /* dist_jk, d_jk */
        v3_sub(&dist_jk,&(ak->r),&(aj->r));
        if(bc) check_per_bound(moldyn,&dist_jk);
-       d_jk=v3_norm(&dist_jk);
+       d_jk2=v3_absolute_square(&dist_jk);
 
        /* jk constants */
        brand=aj->brand;
        if(brand==ak->brand) {
                R=params->R[brand];
                S=params->S[brand];
+               S2=params->S2[brand];
                B=params->B[brand];
                mu=params->mu[brand];
                chi=1.0;
@@ -1771,13 +1806,17 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        else {
                R=params->Rmixed;
                S=params->Smixed;
+               S2=params->S2mixed;
                B=params->Bmixed;
                mu=params->mu_m;
                chi=params->chi;
        }
 
        /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */
-       if(d_jk<S) {
+       if(d_jk2<S2) {
+
+               /* now we need d_ik */
+               d_jk=sqrt(d_jk2);
 
                /* constants_j from exchange data */
                n=*(exchange->n_j);
@@ -1795,7 +1834,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
                /* d_costheta */
                d_costheta1=1.0/dd;
-               d_costheta2=cos_theta/(d_ij*d_ij);
+               d_costheta2=cos_theta/d_ij2;
 
                /* some usefull values */
                h_cos=(h-cos_theta);
@@ -1805,10 +1844,11 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                /* g(cos_theta) */
                g=1.0+c2d2-frac;
 
-               /* d_costheta_ij and dg(cos_theta) - needed in any case! */
+               /* d_costheta_jik and dg(cos_theta) - needed in any case! */
                v3_scale(&temp1,&dist_jk,d_costheta1);
                v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */
-               v3_add(&temp1,&temp1,&temp2);
+               //v3_add(&temp1,&temp1,&temp2);
+               v3_sub(&temp1,&temp1,&temp2); /* there is a minus! */
                v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
 
                /* store dg in temp2 and use it for dVjk later */
@@ -1873,6 +1913,14 @@ if(ai==&(moldyn->atom[0])) {
        printf("%f | %f\n",temp2.y,ai->f.y);
        printf("%f | %f\n",temp2.z,ai->f.z);
 }
+#endif
+#ifdef VDEBUG
+if(ai==&(moldyn->atom[0])) {
+       printf("dVjk (3bp) contrib:\n");
+       printf("%f | %f\n",temp2.x*dist_jk.x,ai->virial.xx);
+       printf("%f | %f\n",temp2.y*dist_jk.y,ai->virial.yy);
+       printf("%f | %f\n",temp2.z*dist_jk.z,ai->virial.zz);
+}
 #endif
 
        }

diff --git a/moldyn.h b/moldyn.h
index dd94b0e..23b5100 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -195,6 +195,7 @@ typedef struct s_tersoff_exchange {
        double f_a,df_a;
 
        t_3dvec dist_ij;
+       double d_ij2;
        double d_ij;
 
        double chi;