testing foo

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 5aa6d87..db91859 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -519,7 +519,8 @@ int scale_volume(t_moldyn *moldyn) {
 
        t_atom *atom;
        t_3dvec *dim,*vdim;
-       double virial,scale;
+       double scale,v;
+       t_virial virial;
        t_linkcell *lc;
        int i;
 
@@ -528,12 +529,23 @@ int scale_volume(t_moldyn *moldyn) {
        vdim=&(moldyn->vis.dim);
        lc=&(moldyn->lc);
 
-       for(i=0;i<moldyn->count;i++)
-               virial+=v3_norm(&(atom[i].virial));
+       memset(&virial,0,sizeof(t_virial));
+
+       for(i=0;i<moldyn->count;i++) {
+               virial.xx+=atom[i].virial.xx;
+               virial.yy+=atom[i].virial.yy;
+               virial.zz+=atom[i].virial.zz;
+               virial.xy+=atom[i].virial.xy;
+               virial.xz+=atom[i].virial.xz;
+               virial.yz+=atom[i].virial.yz;
+       }
+
+       /* just a guess so far ... */
+       v=sqrt(virial.xx*virial.xx+virial.yy*virial.yy+virial.zz+virial.zz);
 
-printf("%f\n",virial);
+printf("%f\n",v);
        /* get pressure from virial */
-       moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial;
+       moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*v;
        moldyn->p/=moldyn->volume;
 printf("%f\n",moldyn->p/(ATM));
 
@@ -984,6 +996,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
        int i,j,k,count;
        t_atom *itom,*jtom,*ktom;
+       t_virial *virial;
        t_linkcell *lc;
        t_list neighbour_i[27];
        t_list neighbour_i2[27];
@@ -1005,7 +1018,13 @@ int potential_force_calc(t_moldyn *moldyn) {
                v3_zero(&(itom[i].f));
 
                /* reset viral of atom i */
-               v3_zero(&(itom[i].virial));
+               virial=&(itom[i].virial);
+               virial->xx=0.0;
+               virial->yy=0.0;
+               virial->zz=0.0;
+               virial->xy=0.0;
+               virial->xz=0.0;
+               virial->yz=0.0;
 
                /* reset site energy */
                itom[i].e=0.0;
@@ -1224,7 +1243,10 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
 
        printf("[moldyn] tersoff parameter completion\n");
+       p->S2[0]=p->S[0]*p->S[0];
+       p->S2[1]=p->S[1]*p->S[1];
        p->Smixed=sqrt(p->S[0]*p->S[1]);
+       p->S2mixed=p->Smixed*p->Smixed;
        p->Rmixed=sqrt(p->R[0]*p->R[1]);
        p->Amixed=sqrt(p->A[0]*p->A[1]);
        p->Bmixed=sqrt(p->B[0]*p->B[1]);
@@ -1285,8 +1307,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        t_tersoff_mult_params *params;
        t_tersoff_exchange *exchange;
        t_3dvec dist_ij,force;
-       double d_ij;
-       double A,B,R,S,lambda,mu;
+       double d_ij,d_ij2;
+       double A,B,R,S,S2,lambda,mu;
        double f_r,df_r;
        double f_c,df_c;
        int brand;
@@ -1317,18 +1339,10 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
         *
         */
 
-       /* dist_ij, d_ij */
-       v3_sub(&dist_ij,&(aj->r),&(ai->r));
-       if(bc) check_per_bound(moldyn,&dist_ij);
-       d_ij=v3_norm(&dist_ij);
-
-       /* save for use in 3bp */
-       exchange->d_ij=d_ij;
-       exchange->dist_ij=dist_ij;
-
        /* constants */
        if(brand==ai->brand) {
                S=params->S[brand];
+               S2=params->S2[brand];
                R=params->R[brand];
                A=params->A[brand];
                B=params->B[brand];
@@ -1338,6 +1352,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        }
        else {
                S=params->Smixed;
+               S2=params->S2mixed;
                R=params->Rmixed;
                A=params->Amixed;
                B=params->Bmixed;
@@ -1346,10 +1361,23 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                params->exchange.chi=params->chi;
        }
 
-       /* if d_ij > S => no force & potential energy contribution */
-       if(d_ij>S)
+       /* dist_ij, d_ij */
+       v3_sub(&dist_ij,&(aj->r),&(ai->r));
+       if(bc) check_per_bound(moldyn,&dist_ij);
+       d_ij2=v3_absolute_square(&dist_ij);
+
+       /* if d_ij2 > S2 => no force & potential energy contribution */
+       if(d_ij2>S2)
                return 0;
 
+       /* now we will need the distance */
+       //d_ij=v3_norm(&dist_ij);
+       d_ij=sqrt(d_ij2);
+
+       /* save for use in 3bp */
+       exchange->d_ij=d_ij;
+       exchange->dist_ij=dist_ij;
+
        /* more constants */
        exchange->beta_j=&(params->beta[brand]);
        exchange->n_j=&(params->n[brand]);
@@ -1389,7 +1417,6 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                s_r=S-R;
                arg=M_PI*(d_ij-R)/s_r;
                f_c=0.5+0.5*cos(arg);
-               //df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij)); /* MARK! */
                df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
                /* two body contribution (ij, ji) */
                v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
@@ -1685,7 +1712,6 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                        s_r=S-R;
                        arg=M_PI*(d_ik-R)/s_r;
                        f_c_ik=0.5+0.5*cos(arg);
-                       //df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* MARK */
                        df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 
                        /* zeta_ij */