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

Thank you for your suggestion.

In the previous results, the temperatures extracted from the velocities of the atom groups needed to be corrected based on the no. of constraints imposed.

The command used for extracting the temperature values from the trajectory velocities: gmx_mpi traj -f nvt.trr -s nvt.tpr -ot temp_SOL1.xvg -w -b 0 -e 5000

Correction needed: T_real = T_calculated * [3natoms / (3natoms - nconstraints)] (Reference: How can I obtain the temperature of a specific group?)

Here

natoms = no. of atoms in the water molecule = 4
nconstraints = no. of constraints = 6 (3 for SETTLES in the topology file of water molecule + 3 for the constraining all the bonds by LINCS in MDP file)
T_real = T_calculated x 2

The corrected temperature profile is shown in the image below which matches with the set temperatures.

For the second simulation, as the temperature of both droplets are maintained at 300K. there is no evaporation and no fluctuations. However, the droplet splits into 2 over a longer time.

Kindly suggest as to what might be the cause of the droplet splitting.

Thanks and regards,
Tonmoy Sharma

Are you applying any type of COM motion removal? As this is essentially a system in the gas phase, I would expect comm-mode = angular to be appropriate here. It seems to be that SOL1 is spinning around SOL2 and I wonder if it is somehow acquiring some kind of artificial force that shears it apart.

Dear sir,

Thank you so much for your suggestion.

Previously I was using the following:

comm-mode = Linear
nstcomm = 1
comm-grps = System
This led to the movement of the droplets around the box.

Now, I have used the following:

comm-mode = angular
nstcomm = 1
comm-grps = SOL1 SOL2

Now the droplets are almost pinned to their initial positions and the droplet evaporation is also minimal. However, I faced too many LINCS warnings after ~ 4.2 ns which terminated the simulation. Is this due to the MDP options set in the new simulations?

Kindly also suggest if such MDP options (comm-mode) will hamper the free energy calculations later using WHAM analysis.

Thanks and regards,
Tonmoy Sharma

I don’t understand how you can use LINCS for constraining bonds in a water molecule that already uses SETTLE for constraints.

Dear sir,

I will be using other molecules in further simulations wherein the SETTLES is not inherently defined for the molecules. So I will have to use LINCS in such simulations, wherein the water molecules will also be treated the similar way. Also the SETTLES defines the constant distance between H atoms (dhh), which will not be constrained by the LINCS. Hence. I am using LINCS along with the SETTLES to constraint all bonds and distances.
Kindly suggest whether this method is appropriate or there are other simple alternate methods.

With regards
Tonmoy Sharma

So in the simulation where you got LINCS warnings you not only have water but also other types of molecules?

Dear sir,

This simulation only contains only water molecules. This simulation is the base for our further simulations.
I have figured out that the LINCS warning is given for an isolated water molecule moving in the box. I am not understanding why. Kindly refer the image below:

Thanks and regards
Tonmoy Sharma

Images are not useful.

I don’t understand how you can have LINCS on a water molecule when you are using SETTLE. Either you are using SETTLE or you are using LINCS, you can’t use both at the same time.

Dear sir,

Thank you for suggestion.

With regards
Tonmoy Sharma

But do you have only water, and both SETTLE and LINCS? Or do you have more than only water?

Dear sir,

In this simulation, I have only two types of water molecules (SOL1 and SOL2). Only difference between these are the resiude names (for simulation convenience), and one of these is defined using SETTLES and the other is defined using [constraints] same as the doh and dhh of SETTLES in the topology file. We needed to do this as SETTLES cannot be used for two moleculetypes (SOL1 & SOL2). The water model used is OPC water model. These are again further constrained using LINCS (in MDP file) during the mdrun.

I have simulated previously using the same conditions and I was able to validate the equilibrium state with experimental results such as surface tension and density. Even now, the LINCS warning is given when the angle is rotated more than 30 degree as specidified in the MDP file. Untill then there are no errors.

Sir, I need to constrain the water molecule O-H and H-H distances as defined for OPC water model, but how shall I specify these values without using SETTLES or [constraints] in topology file. Is there any other way?

Thanks and regards
Tonmoy Sharma

Sir, later I will be using organic molecules with these water molecules. But we havent done that yet. These organic molecules will require LINCS, but then I think we will not be able to exclude the water molecules from LINCS constraints in the same simulation. Hence, for now with only water molecules (SOL1 and SOL2), I am using both LINCS and constraints in the topology together.
I understand that using LINCS and constraints in the topology is creating redundancies, but I am a little confused as to how we can avoid this.

If the distances and masses are the same, you should be able to use SETTLE for both.