How to set up a water-vapor interface calculation?

User discussions
1

GROMACS version: 2020.3
GROMACS modification: No
Hello everyone. I am very new to GROMACS and molecular dynamics in general, but I want to simulate a small ion in the water-vapor (or vacuum) interface. I am curious about how exactly to set up this calculation.
What I did,

  1. Simulated the ion in bulk following the lysozyme tutorial
  2. I took the gro file generated in the previous step (after minimization) and I doubled the z coordinate of the box.
  3. I performed minimization calculation with with emtol 100.0, which was successful but on visualizing the em.gro file I discovered hydrogens on top of the box.

Considering PBC, that’s not abnormal but the following NVT calculation takes the em.gro as input, so I’m wondering if I have to wrap the system to fix the broken bonds?

Please let me know if my procedure makes sense and if not, how should I proceed?

Thank you.

2

Hi RonitS,

With PBC you cannot have only one interface. A liquid slab obtained by doubling the box size in one direction will have two interfaces, so it’s perfectly natural for some atoms to ‘stick-out’ on the other side.

If energy minimization has gone through well then probably there is no ‘broken’ bond (you can always check if the final energy is ‘negative enough’).

In case there is a broken bond, you can initialize the liquid box before doubling its volume by running gmx editconf with the option -nopbc, so that molecules are made whole.

Michele

3

Hello Michele,

Thank you for your answer. You’re absolutely correct in asserting that there will be two interfaces with PBC, so I will triple the size of my box. While looking through the manual for editconf and solvate I found that they can be used to set box sizes. I assume using those I’ll be able to generate the two interfaces.

As for the broken bonds, I used the editconf with -nopbc option and the bonds were indeed not broken but I was curious if proceeding with the original .gro file (with the seemingly ‘broken’ bonds) to the NVT equilibration would produce any artifacts instead of modifying the .gro file to make the molecules whole.

Regardless, I will look into how to triple the box while placing my solvated system at the center.

4

Hi again,

You are welcome! There is no need to triple the size of the simulation box. If you find more more convenient to place the liquid slab in the center you can just run gmx editconf with -c option. However as far as I understand this should be irrelevat when it comes to computing surface tension, you should get the same result regardless if the slab is centered or not.