bugfix dg of dV_ji,jk

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index db91859..c21f59f 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -541,13 +541,13 @@ int scale_volume(t_moldyn *moldyn) {
        }
 
        /* just a guess so far ... */
-       v=sqrt(virial.xx*virial.xx+virial.yy*virial.yy+virial.zz+virial.zz);
+       v=virial.xx+virial.yy+virial.zz;
 
 printf("%f\n",v);
        /* get pressure from virial */
-       moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*v;
+       moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+ONE_THIRD*v;
        moldyn->p/=moldyn->volume;
-printf("%f\n",moldyn->p/(ATM));
+printf("%f | %f\n",moldyn->p/(ATM),moldyn->p_ref/ATM);
 
        /* scale factor */
        if(moldyn->pt_scale&P_SCALE_BERENDSEN)
@@ -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 */
@@ -1427,6 +1431,15 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        /* add forces of 2bp (ij, ji) contribution
         * dVij = dVji and we sum up both: no 1/2) */
        v3_add(&(ai->f),&(ai->f),&force);
+
+       /* virial */
+       ai->virial.xx-=force.x*dist_ij.x;
+       ai->virial.yy-=force.y*dist_ij.y;
+       ai->virial.zz-=force.z*dist_ij.z;
+       ai->virial.xy-=force.x*dist_ij.y;
+       ai->virial.xz-=force.x*dist_ij.z;
+       ai->virial.yz-=force.y*dist_ij.z;
+
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
        printf("dVij, dVji (2bp) contrib:\n");
@@ -1434,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 ... */
@@ -1526,6 +1547,15 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* add force */
        v3_add(&(ai->f),&(ai->f),&force);
+
+       /* virial */
+       ai->virial.xx-=force.x*dist_ij->x;
+       ai->virial.yy-=force.y*dist_ij->y;
+       ai->virial.zz-=force.z*dist_ij->z;
+       ai->virial.xy-=force.x*dist_ij->y;
+       ai->virial.xz-=force.x*dist_ij->z;
+       ai->virial.yz-=force.y*dist_ij->z;
+
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
        printf("dVij (3bp) contrib:\n");
@@ -1533,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 */
@@ -1562,6 +1600,15 @@ if(ai==&(moldyn->atom[0])) {
 
        /* add force */
        v3_add(&(ai->f),&(ai->f),&force);
+
+       /* virial - plus sign, as dist_ij = - dist_ji - (really??) */
+       ai->virial.xx+=force.x*dist_ij->x;
+       ai->virial.yy+=force.y*dist_ij->y;
+       ai->virial.zz+=force.z*dist_ij->z;
+       ai->virial.xy+=force.x*dist_ij->y;
+       ai->virial.xz+=force.x*dist_ij->z;
+       ai->virial.yz+=force.y*dist_ij->z;
+
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
        printf("dVji (3bp) contrib:\n");
@@ -1569,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;
@@ -1583,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;
@@ -1631,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;
@@ -1645,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);
@@ -1677,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);
@@ -1730,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;
@@ -1744,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);
@@ -1768,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);
@@ -1778,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 */
@@ -1830,6 +1897,15 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
                v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
                                                  /* scaled with 0.5 ^ */
+
+               /* virial */
+               ai->virial.xx-=temp2.x*dist_jk.x;
+               ai->virial.yy-=temp2.y*dist_jk.y;
+               ai->virial.zz-=temp2.z*dist_jk.z;
+               ai->virial.xy-=temp2.x*dist_jk.y;
+               ai->virial.xz-=temp2.x*dist_jk.z;
+               ai->virial.yz-=temp2.y*dist_jk.z;
+
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
        printf("dVjk (3bp) contrib:\n");
@@ -1837,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
 
        }