improved log/report subsystem, playing around w/ pressure, sic hook

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

diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index 972075d..16f771c 100644 (file)
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -95,6 +95,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        double s_r;
        double arg;
 
+       /* use newtons third law */
+       //if(ai<aj) return 0;
+
        params=moldyn->pot_params;
        brand=aj->brand;
        exchange=&(params->exchange);
@@ -119,27 +122,11 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
         *
         */
 
-       /* constants */
-       if(brand==ai->brand) {
-               S=params->S[brand];
+       /* determine cutoff square */
+       if(brand==ai->brand)
                S2=params->S2[brand];
-               R=params->R[brand];
-               A=params->A[brand];
-               B=params->B[brand];
-               lambda=params->lambda[brand];
-               mu=params->mu[brand];
-               exchange->chi=1.0;
-       }
-       else {
-               S=params->Smixed;
+       else
                S2=params->S2mixed;
-               R=params->Rmixed;
-               A=params->Amixed;
-               B=params->Bmixed;
-               lambda=params->lambda_m;
-               mu=params->mu_m;
-               params->exchange.chi=params->chi;
-       }
 
        /* dist_ij, d_ij */
        v3_sub(&dist_ij,&(aj->r),&(ai->r));
@@ -166,12 +153,26 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        exchange->d_j=&(params->d[brand]);
        exchange->h_j=&(params->h[brand]);
        if(brand==ai->brand) {
+               S=params->S[brand];
+               R=params->R[brand];
+               A=params->A[brand];
+               B=params->B[brand];
+               lambda=params->lambda[brand];
+               mu=params->mu[brand];
+               exchange->chi=1.0;
                exchange->betajnj=exchange->betaini;
                exchange->cj2=exchange->ci2;
                exchange->dj2=exchange->di2;
                exchange->cj2dj2=exchange->ci2di2;
        }
        else {
+               S=params->Smixed;
+               R=params->Rmixed;
+               A=params->Amixed;
+               B=params->Bmixed;
+               lambda=params->lambda_m;
+               mu=params->mu_m;
+               params->exchange.chi=params->chi;
                exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
                exchange->cj2=params->c[brand]*params->c[brand];
                exchange->dj2=params->d[brand]*params->d[brand];
@@ -210,12 +211,13 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        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;
+       virial_calc(ai,&force,&dist_ij);
+       //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])) {
@@ -325,12 +327,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        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;
+       virial_calc(ai,&force,dist_ij);
+       //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])) {
@@ -379,12 +382,13 @@ if(ai==&(moldyn->atom[0])) {
 
        /* virial - plus sign, as dist_ij = - dist_ji - (really??) */
 // TEST ... with a minus instead
-       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;
+       virial_calc(ai,&force,dist_ij);
+       //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])) {