}
/* 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)
/* 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");
/* 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");
/* 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");
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");