[gmx-developers] pbc and virial problem

[Alexander Herz]

Hi,

I've added a force similar to the walls. Now I'm trying to calculate
surface tensions using this new correction via the virial.
We used to do these calculations with pbc=xyz with an older gromacs
version so I'm trying to do the same in order to be able compare the
results before going to pbc=xy.

The problem is that I get non sense for the virial xx and yy components.
To figure out why, I compared the averaged virials for the same sim with
and without the new correction.

The virial contribution from the non-shifted particles is pretty much
identical but the virial from the shifted
particles:
  (calc_virial in sim_util.c):
  clear_mat(vir_part);
 
calc_vir(fplog,SHIFTS,fr->shift_vec,fr->fshift,vir_part,ePBC==epbcSCREW,box);
  inc_nrnb(nrnb,eNR_VIRIAL,SHIFTS);

differs by almost exactly the amount which I'd expect.

The correction I apply only adds forces in z direction and the
simulation is setup such that no pbc in z is applied because of the cut
off (there is a water slab of ca. 6 nm thickness in the middle of the
box (z direction) surrounded by ca. 3 nm of nothing in both z
directions). Debugging through the code I can see that no images for
z!=0 are created.
The correction force is quite small compared to the forces applied by
gromacs non_bonded
(correction: avrg. force=0.03, max. force=8.0, gromacs: avrg.:456.0 ,
max.: 4831.0).

The force is calculated just after the walls right before do_nonbonded.

Finally, the question:
Any ideas how the small additional force in z direction can have such a
strong effect on the virial xx and yy components
(apparently fshifted[XX] and fshifted[YY] are influenced significantly
whereas the virial ZZ value remains virtually identical
no matter if the correction is applied or not, as fshifted[ZZ]==0 in
both cases).

Thx,
Alex