I am trying to calculate the surface tension between a ‘surrogate’ solid substrate I obtain by restraining atomic positions and either a fluid or a gas phase. The results I get are a bit strange, in the sense that I obtain a negative estimate for the surface tension at the solid/something interface, in both cases.
By definition, surface tension is computed as the difference between normal and parallel components (w.r.t. the interface) of the pressure tensor; my understanding is that the pressure tensor is in turn computed from the virial, which may be ill-defined when position restrains are present.
When doing NPT one should select a proper refcoord-scaling, but what about NVT calculation? Is the virial always well-defined? Or maybe never well defined, and so one should never do a NVT simulation to estimate surface tension when position restreaints are present? Is it possible to extract the pos. res. contribution to surface tension alone?
For the calculating the pressure, it is not relevant whether you are using an NVT or NPT ensemble. For the definition of the pressure, it is relevant how the system scales when the volume scales. So it matters how the position restraint reference coordinates would scale. When a whole, periodic substrate is restrained, the only option is to scale the coordinates along. This can lead to very large virial contributions. But they should cancel in most physically meaningful quantities, as one usually computes a difference for two systems.
I see. The problem is that the contribution does not cancel out when taking differences (or at least I cannot find a way to cancel it without it resulting in a under-detemined system), so I was thinking about whether it is possible to estimate the contribution due to the restraints alone.
Is there any reference on how the position restraints contribution specifically is accounted for in the virial calculation (I cannot find it on the GROMACS reference manual)? I guess the most general formula (e.g. Tuckermann) still holds, but maybe there is some ‘trick’ implemented in the code to make it easier.
That is what I thought at first, however it’s not as simple.
Imagine comparing a substrate-gas system to a substrate-fluid system when I compare the two results I get on one hand 2surf_tens_SV + restraints_contribution, while on the other hand I get 2surf_tens_SL + restraints_contribution. Even assuming that the contribution of position restraints is the same in both cases, the resulting system is under-determined.
What I did initally was to compare a solid-gas (‘vacuum’ really) configuration with a soild-fluid-vapour one and that also lead to an underdetermined system (even provided one has an estimate of the equilibrium contact angle between the three phases). Doing a pure solid-liquid NVT simulation is not so trivial on the other hand because one would need to make sure liquid density exacly matches the one of the film in the S-L-V configuration.
Anyway, I believe these are now only technical issue; the most fundamental problem was if the surface tension calculation itself made sense when using position restraints, which now I believe it does. Thanks.