Water surface tension

Hi everyone,

I am trying to compute the surface tension of TIP4P water. I am not sure I am doing this right. I first created a 6 6 3 water box and performed EM, NVT equilibration and NPT equilibration. Then, I added vacuum in the z-axis using the following command:

gmx editconf -f *.gro -o with_vacuum.gro -box 6 6 18 -center 0 0 3

Then I performed an EM.

At the end, I have a box of 6 6 18. I then ran a 30 ns NPT simulation to compute surface tension. My .mdp file is:

integrator = md
nsteps = 30000000 ; 20 ns
dt = 0.001

cutoff-scheme = Verlet
nstxout = 1000
nstvout = 1000
nstenergy = 1000
nstlog = 1000
nstxout-compressed = 1000

rlist = 1.2
rcoulomb = 1.2
coulombtype = PME

; van der Waals
vdwtype = cutoff
vdw-modifier = potential-switch
rvdw-switch = 1.0
rvdw = 1.2

; Temperature coupling
tcoupl = Nose-Hoover
tc-grps = System
tau_t = 1.0
ref_t = 298

; Pressure coupling for surface tension
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
ref-p = 1.0 1.0
compressibility = 4.5e-5 4.5e-5

pbc = xyz
constraints = h-bonds

I got a value of 88 mN/m (with gmx energy), but, as I have 2 interfaces, I think I have to divide it by 2, which makes it incorrect. Also, interfaces remain throughout the simulation, but the z-axis size grows a lot. I’m thinking maybe I should do:

ref-p = 1.0 0.0
compressibility = 4.5e-5 0.0

instead of what I did, but I wanted to know your opinion on my simulation workflow (EM > NVT > NPT > add vacuum > EM > NPT production run); I don’t know if it’s best to create the vacuum before or after the first EM. Everytime I tried to create vacuum before the first EM, during the NPT the water molecules would just “invade” the vacuum, and I don’t know how to prevent this from happening. Thoughts?

Thanks!

You should not use pressure coupling when calculating the surface tension.

Thank you, hess. So the best method would be to build the box with interfaces, perform EM, then a NVT run and compute the surface tension of this run? I’m afraid if I do an EM with a box without interfaces > NVT > NPT > add vacuum > EM > NVT then the system is not equilibrated enough. I tried this and got a high potential energy (10^17).

A water box equilibrates very quickly. I usually take an equilibrated 3D-periodic water box, if necessary stack a few copies, add vacuum, run. No extra EM needed. I ignore the first 1 ns for equilibration.

Ok, so you add vacuum before equilibration and don’t run NPT? What if the system is more complex than water and takes more time to equilibrate / a NPT run to equilibrate?

I don’t understand. You can not run NPT when you have vacuum in the system, as the equilibrium state would be that the vaccuum disappears. You do want to allow some time for the interface to equilibrate, but I simply ignore some of the initial part of the production simulation.

I mean, if you add the interface before equilibrating the system, then you’re just supposed to equilibrate it with NVT? Or is it best to equilibrate a box without interfaces, add the interfaces and then run a NVT?

Equilibrate a bulk liquid (no interfaces) under NPT. Then extend the z dimension of the box such that there is vaccum on top and bottom and then run an NVT.
Arun