Constraining two water droplets at constant distance between their COMs using pull code

GROMACS version: 2021.3
GROMACS modification: No

Hello all,

I need to constrain two water droplets, each having 5000 molecules with different residue names (SOL1 and SOL2) and atom names. I have recently tried pull code with the following MDP options:

pull = yes
pull_ngroups = 2
pull_ncoords = 1
pull_group1_name = SOL1
pull_group2_name = SOL2
pull_coord1_type = umbrella ; harmonic biasing force
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_groups = 1 2
pull_coord1_dim = Y Y Y
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_coord1_start = yes ; define initial COM distance > 0

The initial configuration is shown below at the left side:

The water droplet on the right side is evaporating fast and the distance between the COMs are not being constrained as the initial configuration. The final configuration is shown at right side after 1380 ps.

Kindly help me constrain the COMs of both the droplets at their initial configurations. I plan to find the free energy at such constratined distances between their COMs.

Thanks and regards,
Tonmoy Sharma

What is the deviation of the distance? I can’t see that from your picture.
What distance does the pull coord report in pullx.xvg?

It could be that your force constant is too low.

Hello hess,

Thank you for your reply.

Please see the image below which is after 6.5 ns. Here, the smaller droplet is droplet which was towards the right side. It has evaporated and all the water molecules have been scattered inside the box. I have checked the pullx.xvg file, the distances are almost constant at the start but increases towards the end. The main problem is the evaporation. Kindly suggest regarding this.

I have also checked for pull force values, i.e., 10, 1000, 10000 kJ mol^-1 nm^-2. But the behaviour remains almost the same for all.

Thanks and regards
Tonmoy Sharma

If a droplet evaporates, its COM is no longer well defined as the distance of evaporated molecules jumps up and down with periodic boundary conditions.

So, how can I prevent this undesired evaporation? Only the right side droplet is evaporating with huge fluctuations in positions of water molecules, while the left side droplet is moving about the right droplet with comparatively much less fluctuations and no evaporation. I am using OPC water model for the simulations. Kindly suggest any changes to be made in the MDP file.

I don’t know why your droplet is evaporating. Did you equilibrate the droplets properly?

I remember that we also experienced problems with evaporating droplets when system sizes got very large. If a box had a long extension in, say, the x-direction, then water droplets near the origin would behave as expected, while droplets with high x positions would heat up, start to rotate and evaporate. Also, if you take an exact same droplet that runs fine in a small box, and you put it into a far corner of a much (factor 100) larger box, it would start to evaporate. This all happens in single precision, but not in double precision. I think single precision might not be enough to calculate differences of positions (as e.g. needed to fulfil the constraints within each water molecule) accurately enough for molecules with large x-, y-, or z coordinates when using floating point. So you could try to check whether running in double precision solves the issue. Are you running with or without thermostat and barostat?

Hi, if you are running single precision and setting tc-grps=system, you could also look at the temperature of each group (SOL1 and SOL2) separately (e.g. by using gmx traj). Another ‘signature’ of the single-precision issue that Carsten highlighted is the unphisical transfer of thermal energy from one group to the other, i.e. you should see a droplet heating up and the other cooling down.

Yes, the droplets were created individually by NVT for 50ns and later merged together into the same GRO file. The radial density of the droplets was around 1 g/cc.

Dear Carsten,

Thank you for your suggestion.

Yes, I am using Nose-Hoover thermostat but no barostat. Would you kindly guide me how to execute mdrun in double precision.

Thanks an regards
Tonmoy Sharma

Hi Tonmoy,

you have to use the GROMACS version that is compiled in double precision. If it is provided on the cluster you use, it is probably named gmx_d instead of just gmx. Or did you compile it yourself? Then you will have to make an extra double precision build.

Best,
Carsten

You could also try to thermostat each group separately. Maybe that’s enough to kill unphysical evaporation.

I think a bit of evaporation is physical at 300K :)

Exactly. This might be an unavoidable issue.

But check the actual temperature of the droplets.

Sure, but I had expereinced quite extreme heating/cooling of water droplets when running single-precision in large domain. So I wouldn’t be surpised if the temperature of one of the two droplets is much greated than 300K (and consequently the other may be even frozen).

That happens at distance larger than 100 nm. I guess the distances are smaller here, but I can’t really say from the images.

Dear Carsten

We have installed the GROMACS version 2021.3 in single precision. Would you kindly share how to make an extra double precision build or share any available source in internet. It will be really helpful for me.

With regards,
Tonmoy

Dear Michele,

I have simulated two systems based on your suggestions:

  1. With the mdp option “tc-grps=System” and “ref_t = 300”
  2. With the mdp option “tc-grps=SOL1 SOL2” & “ref_t = 300 300”, i.e., separate for both droplets.

Kindly observe the results of the temperature of both droplets with time. For case 1, as the temperature of one droplet is rising, the temperature of the other droping. Whereas, the temperature of both remains almost same (about 150K) thoguhout the simulation for case 2. However, the sum of the temperatures of both the droplets is about 300K.

Kindly shed some light as to why both the droplets are not maintained at 300K even though I have set the reference temperature as 300K for both.

Thanks and regards,
Tonmoy Sharma

Dear Hess,

I have encountered recently that as I use separate tc-grps for SOL1 and SOL2, the droplet with SOL1 residues breaks into 2 droplets after 4.3 ns of the simulation. The splitting takes place within 1.2 ns. I have used the following mdp options:

;TEMPERATURE coupling
Tcoupl = nose-hoover
tc-grps = SOL1 SOL2 ; Groups to couple separately
nsttcouple = -1
tau_t = 0.5 0.5 ; Time constant (ps) and
ref_t = 300.0 300.0 ; reference temperature (K)
; PULL CODE
pull = yes
pull_ngroups = 2
pull_ncoords = 1
pull_group1_name = SOL1
pull_group2_name = SOL2
pull_coord1_type = umbrella ; harmonic biasing force
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_groups = 1 2
pull_coord1_dim = Y Y Y
;pull_coord1_rate = 0 ; nm/ps (-ve means in the -ve Z direction)
pull_coord1_k = 10 ; kJ mol^-1 nm^-2
pull_coord1_start = yes ; define initial COM distance > 0

Kindly suggest what might be the problem here.

Thanks and regards,
Tonmoy Sharma

You add the -DGMX_DOUBLE=on cmake flag to build GROMACS in double precision, see the installation instructions at https://manual.gromacs.org/current/install-guide/index.html. There will however be no GPU support in double precision.