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!