Rigid TIP4P energy and force computations

User discussions
,
1

GROMACS version: 2024.3 (double precision)
GROMACS modification: No

Hello,

My goal is to compute the Hessian for the TIP4P2005 water model, that has LJ interactions (between oxygens), and coulomb interactions (between H and dummy particles). Another feature that’s important here is that the model has rigid constraints.

My current understanding is that GROMACS cannot minimize with rigid constraints (right?), so my strategy has been to export all the GROMACS simulation output files into MDAnalysis, and compute in python all the related quantities and perform my minimization there (the questions below do not have anything to do with MDAnalysis).

My procedure is as follows: I equilibrate the water molecules at high temperature, and then cool it down. After cooling it down I follow with a super damped Langevin dynamics (that allows me to maintain the pressure at its wanted level, while taking away the kinetic energy). I do this with the following command (mdp and top files attached):

gmx_mpi_d grompp -f "damped_quench_p=1000_T=80_s=1.mdp" -c "Tmp_p=1000_T=80_s=1.gro" -p tip4p2005-1000_tmp.top -o "DQuenched_p=1000_T=80_s=1.tpr"
gmx_mpi_d mdrun -s "DQuenched_p=1000_T=80_s=1.tpr" -mp tip4p2005-1000.top -c "DQuenched_p=1000_T=80_s=1.gro" -e "DQuenched_p=1000_T=80_s=1.edr" -o "DQuenched_p=1000_T=80_s=1.trr" -pin on -ntomp 1 -v -deffnm em

This gives me .edr, .gro and .trr files which I then load in python.

I now have two questions. The first part has to do with energy computation (mostly the Coulomb (SR) part), and the second question has to do with the resulting forces and how to extract/interpret those.

Part 1:
My first goal is to reproduce the energy values stored after the dynamics in the .edr file. As I remove most of the kinetic energy, I aim first at recovering the potential energy. The potential energies GROMACS gives are divided into three components - LJ (SR), Coulomb (SR) and Coul. recip.. I always compare the last values stored (e.g., I use edrFile.data_dict[‘Coulomb (SR)’][-1]). I would like to be able to recover these three values from the .gro and .trr files (together with my model parameters of course). To be clear I compare the computed values to those obtained by GROMACS using this: \Delta = \frac{E_\text{edr} - E_\text{computed}}{E_\text{edr}}.

I start with the LJ component as it only operates between the oxygens. I am able to compute it in python with \Delta_{LJ} = 0.000010 - so to a fairly good extent. I would like to achieve a similar value with the Coulomb interactions.

My issue starts with the short-range coulomb interactions, which I obtain to \Delta_{C(SR)} \simeq 0.24 - so something is off. I want to understand what I am missing in my Coulomb energy calculations (we will keep the Ewald/PME contributions aside for now).

Here is what I compute - first, I define the charges for the H’s (Q = 0.5564), and Dummy particles (-2Q), and set the constants \alpha = 2.6035 and F = 138.935458. Then, I run over all pairs that do not belong to the same molecule, and compute the coulomb energy as

E_{C(SR)} = F q_1 q_2 \frac{\text{erfc}(\alpha r)}{r}.

I also tried computing with the cutoff as E_{C(SR)} = F q_1 q_2 \left[\frac{\text{erfc}(\alpha r)}{r} - \frac{\text{erfc}(\alpha r_{cut})}{r_{cut}}\right] but this does not improve the situation. I also add the “self” energy E_{C(SR)} = -F q_1 q_2 \frac{\text{erfc}(\alpha r_{cut})}{r_{cut}} for pairs within the same molecule, but this contribution is neglegible with respect to the first one.

This causes me to believe that I am missing something in the first Coulomb computation, but I am not sure what is missing (is this something to do with the PBCs, an aspect I am missing from the potential calculations, or other) - any help would be much appreciated.

Part 2:
After getting the energy the next stage would be to compute gradients and perform a rigid-body minimization in python. To do this I want to compare the forces in my GROMACS to the ones I compute from taking the energy derivatives.

The .trr file above gives me the forces on all the particles. I observe the forces at the end of the simulation by using MDAnalysis - if u1 is the universe, then I set u1.trajectory[-1] to load the last frame, and get the forces by u1.atoms.forces.
Now I have another run using this last frame as an initial state. Say I try and minimize (whatever GROMACS does in this case) using (quench1.mdp file attached):

gmx_mpi_d grompp -f "quench1_p=1000_T=80_s=1.mdp" -c "DQuenched_p=1000_T=80_s=1.gro" -p tip4p2005-1000_tmp.top -o "Min_p=1000_T=80_s=1.tpr"
gmx_mpi_d mdrun -s "Min_p=1000_T=80_s=1.tpr" -mp tip4p2005-1000.top -c "Min_p=1000_T=80_s=1.gro" -e "Min_p=1000_T=80_s=1.edr" -o "Min_p=1000_T=80_s=1.trr" -pin on -ntomp 1 -v -deffnm em

This also produces a .trr file. I now want to compare the forces from the last frame in the first trr file - denoted here by u1 - , and the first frame in the second .trr file, denoted here as u2.
When I load these to MDAnalysis, the forces from u1 are several orders of magnitude larger than in u2 (both different from the values I compute explicitly with taking an explicit gradient, even when I consider only the oxygens, which “feel” only the LJ part of the forces).
The forces are different in both magnitude and in direction (so it seems like there is no simple ``conversion factor’’ that can explain this discrepancy).
Again - I seem to be overlooking something. Would love to have some input here.

I have attached mdp files and the top file for anyone interested in running these (couldn’t upload the .gro files or the resuling .edr and .trr files).

Many thanks in advance
damped_quench_p=1000_T=80_s=1.mdp|attachment (1.7 KB)
quench1_p=1000_T=80_s=1.mdp|attachment (773 Bytes)
tip4p2005-1000.top (1.3 KB)

2

GROMACS can minimize with rigid constraints.

For the normal mode analysis you would need to do some algebra to project out the constraints dimensions.

For plain Coulomb, GROMACS adds a correct to excluded pairs. I have recently added a note on this to the mdp option in the user guide.

3

Hi @hess - many thanks for the prompt answer.

So do I understand correctly that GROMACS minimizes properly the system using the minimization code I provided?
When I run this minimization procedure I get that “Steepest Descents converged to machine precision in 258 steps, but did not reach the requested Fmax < 0.0001.”

For my purposes I want to get the forces as close to zero as possible, but the resulting forces after GROMACS minimization are still fairly large (max force of 1.2+02). This probably introduces some issues in the normal modes. Is there anything I can do to improve this?

For the Coulomb - are you referring to the fact that the excluded pairs contribute by the value of the potential at the cutoff?
If that’s what you had in mind, at least what I see is that this contribution is negligible, and does not fix the issue above (e.g., I see that the pairs contribute ~ -64K but the excluded pairs contribute 1). Could you be more specific about the corrections you refer to?

Many thanks again

4

What do you mean with code? All GROMACS minimizers support constraints.

I don’t know from what configuration you started, but minimizing liquid water to that small a force is practically impossible. Also note that the global energy minimum is ice.

5

Well,

I meant that I provided an mdp file in the original message, that has minimization using the steepest descent algorithm. I understand that the implementation of this does support constraints (but for example conjugate-gradient does not support constraints - if I understood this).

Regarding minimization - a (mathematically) “real” minimum would have zero net forces on each of the particles. I am interested in local energy minima as obtained after cooling sufficiently fast (so amorphous ice). Extracting the Hessian when there are non-vanishing (or non-negligible) forces is not informative for me - I really need to get those forces to be vanishingly-small.

As my steepest-descent minimization starts from a state that has been obtained by super-damped Langevin dynamics, I suspect that the initial state is not to blame here. I just wonder if there is a way for GROMACS to minimize so that the resulting forces are really negligible (as should occur in an energetic minimum).

Why are you saying its “practically impossible”?
Is it a precision issue, or a GROMACS implementation issue, or a water model issue?

6

All minimizers work with constraints. That comment in the manual is (very) out of date. Steepest-descents in not sufficient for minimization for normal modes.

The problem with minimizing liquids, and water in particular, is that a liquid is, by definition, liquid. This is the effect of thermal motion. Minimizing means going to 0 K, in which case all liquids become a solid. Iterative energy minimization is a very inefficient procedure to transform a liquid into a solid. A question is if you want the NM of ice or of a supercooled state of water. In the first case you should start from a crystal. In the second case minimization might take too long.

7

“All minimizers work with constraints. That comment in the manual is (very) out of date. Steepest-descents in not sufficient for minimization for normal modes” - OK, so I will now use CG, thanks.

As mentioned above, my goal is to look at amorphous states of ice. Practically, I take water (at room T values, say), and cool it rapidly (I am doing this by recurring equilibrations at decreasing temperature values). Once I reach a sufficiently low temperature, I use overdamped Langevin dynamics to damp out most/all the kinetic energy (as implemented in damped_quench…mdp code attached in the original message). After running this I reduce the kinetic energy substantially (from ~2K to 10). Just for comparison, the potential energy here is ~-70K, so the kinetic part is completely negligible.

From this state - which should already have very small net forces on the particles - I want to find the closest minimum, that should result in ~zero net forces. So, indeed I want to examine an (amourphous) solid here, and the minimization procedure should converge rather rapidly as I have already taken away most of the kinetic energy. I am just not sure why its so difficult to converge better to an energetic minimum given that the system is already pretty close to it (this is the reason I had to export the GROMACS data to python in the first place).

Is there a way to improve the minimization procedure in GROMACS, so that the resulting forces will be truly negligible?

Additionally, I still do not understand what was wrong in my energy computation of the Coulomb SR interactions (could it be the pbcs? any additional term that I am missing from GROMACS computation, etc), and I also wonder about the forces - how come when I finish the first simulation (damped_quench) I get forces that are orders of magnitude larger (and in different directions) from the first frame in the following simulation (quench1) - am I doing something wrong in loading the data back to GROMACS?

Thanks again for all the information and help

8

I have never tried myself, but I still think it will be difficult to minimize amorphous ice. The hydrogen bond network is very connected, which means that slightly changing the orientation of one water molecule will affect many others. Maybe L-BFGS works better.

My guess is that the energy mismatch is due to the subtraction of the energy between excluded pairs. Not shifting the potential gets rid of this effect. Buy for you case a plain Coulomb cut-off is completely unsuitable. For minimization and NM you need continuous forces.

My guess of you last issue is that you might have used the gro file for the next grompp step, which has only 3 digits of precision. Use a trr file instead.

9

I also tried with L-BFGS, and seems to be working as well, however still all the algorithms converge to the same order of magnitude of forces, which still implies that the forces are not vanishingly small (and by definition this is not a minimum).

I also tried loading the .trr file instead of the .gro but this didn’t improve the minimization process.

I wonder, is it possible to compile GROMACS to quad precision? Maybe this would help?

I am still not sure I completely understand your comment about the potential shift in the excluded pairs - I computed this contribution (i.e., I added the value of the potential at its cutoff distance to all the excluded pairs - those that are within the same molecules), and it was completely negligible with respect to the total energy. Is this what you refer to?