With free-energy=yes, angle energies are off

User discussions
1

GROMACS version:2021.4
GROMACS modification: No

Hi,
This is a little bit involved, please bear with me …
I was trying to emulate with gmx a kind of dual topology free-energy calculation,
using a simple ethane->methanol in vacuum test-case.
I first calculated single point energies for a particular conformer of each ethane, and methanol.
then did the same calculation for a hybrid structure with dual topology (CH3(-CH3)(-OH))
identical coordinates, free-energy = yes, delta-lambda = 0, and init-lambda-state
corresponding to either lamda=0 (ethane) or lambda=1 (methanol). So I have four cases:

  1. a normal simulation (single point energy) for ethane
  2. a normal simulation (single point energy) for methanol
  3. a simulation of the hybrid, free energy on, and lambda=0 (corresponding to ethane)
  4. a simulation of the hybrid, free energy on, and lambda=1 (corresponding to methanol)

Again, these are all only single point energies, and based on identical coordinates.
As expected the non-bonded energy terms are identical in cases 1 vs 3, as
well as in 2 vs 4. Also as expected, the bonded energy terms of case 3 and 4
are identical, and equal to the sum of the corresponding terms in case 1 and 2,
However the last statement is true only for the bond and the dihedral energies,
but not for the angle terms where the numbers clearly differ!
This came quite unexpected, I double and triple checked my input files, but found
nothing that would explain this …
I put the input files into an archive - which i cannot attach/upload here, therefore
I put the archive here: http://www.brunsteiner.net/angles.tgz from where it can be
downloaded - to reproduce the issue type:

case1

gmx grompp -f md.mdp -c e.pdb -p e.top -o mde.tpr -maxwarn 1
gmx mdrun -v -deffnm mde -c mde-out.pdb
#case 2
gmx grompp -f md.mdp -c m.pdb -p m.top -o mdm.tpr -maxwarn 1
gmx mdrun -v -deffnm mdm -c mdm-out.pdb
#case 3
gmx grompp -f md0.mdp -c ems.pdb -p fe.top -o fe0s.tpr -maxwarn 1
gmx mdrun -v -deffnm fe0s -c fe0s-out.pdb
#case 4
gmx grompp -f md1.mdp -c ems.pdb -p fe.top -o fe1s.tpr -maxwarn 1
gmx mdrun -v -deffnm fe1s -c fe1s-out.pdb

I’d be grateful for any suggestions here!
cheers,
Michael

2

This is not the correct way to perform a single-point energy evaluation. You should be using gmx mdrun -rerun to compute the energies:

gmx mdrun -v -deffnm fe1s -rerun ems.pdb

I noticed you are invoking -maxwarn 1, which is almost never appropriate. grompp is warning you that it is copying A and B state parameters, when this is almost certainly something that you don’t want. You’ve got bonds, angles, etc. that change between the A and B state, but you’re letting them share parameters. For instance, in the case of (what should be) ethane in the free energy topology, you have dummy atoms replacing O and H, but that have normal bonded parameters, so you get bond, angle, and dihedral contributions from them as if ethane has a hydroxyl group (in terms of the bonds, since they are dummy atoms this doesn’t show up in the LJ or electrostatics terms). When I do a single-point evaluation with mdrun -rerun, I see discrepancies in every internal term (bonds, angles, and dihedrals). This indicates your topologies are incorrect for doing the transformation.

3

Dear Justin,

Thanks for your reply!
You wrote:

This is not the correct way to perform a single-point energy evaluation. You should be
using gmx mdrun -rerun to compute the energies

There might be cases for which your point is valid, however, in this case it does not
matter. To be sure I tried your suggestion, and I do get exactly the same energies, if I
do it my way (as above) or in your way with a rerun.

I see discrepancies in every internal term (bonds, angles, and dihedrals).

If you really get different numbers there please tell me what you get, and which
gmx version you used, as this would be rather strange. i attached a screenshot
from a spreadsheet with my numbers.

e
e822×419 56.6 KB

you also say:

“I noticed you are invoking -maxwarn 1 , which is almost never
appropriate. grompp is warning you that it is copying A and B state
parameters, when this is almost certainly something that you dont want.”

I just tell the program that this particular parameter does not change/is the same in
both end-states/topologies (at lambda=0 and at lambda=1), which is exactly what I
want in those cases.

As for the rest of what you wrote, this might apply if I wanted to do a single
topology approach (which is what you used in the excellent tutorials you provided
for gromacs). What I try to do here is a dual topology free energy calculation.
Both approaches come with their respective pros and cons, and I have good reasons
for wanting dual topologies in this case. If you are not sure about what the difference
is (and if you have some time to spare;) you might want to take a look at:
https://doi.org/10.1021/jp981628n and https://doi.org/10.1021/jp981629f

Anyway, the energy differences I observed for the angle terms are in fact fine
this was a (somewhat stupid) misconception on my part - sorry for wasting your time!

cheers,
Michael

4

I get the same numbers you do - there is a difference between mde.log and fe0s.log when presumably there should not be, correct? I get the same values you are showing. If I understand your setup correctly, cases 1 and 3 should be equivalent for λ = 0, no? I have an idea of why this is but I’d like to confirm this is what should be the case.

But this confuses me. As I look at fe.top, I see you are explicitly copying A → B state parameters, e.g. for [bonds] that do not change (terminal methyl group) but then you do not have explicit B-state parameters for those interactions that are changing (hence grompp copies them and warns you that you have different end states that rely on copied bonded parameters, which is usually inappropriate). For instance, in ethane (A-state, λ = 0), there should not be an O-H bond. You have these assigned as dummy atoms with zero charge (therefore not contributing to LJ or Coulombic energies) but you do have bonded parameters listed in the A-state between these dummies, so naturally the energies between the standard MD topology and the free energy topology disagree. You’re assigning different states, carbons with 5 bonds, etc.

5

cases 1 and 3 should be equivalent for λ = 0, no?

correct

for [bonds] that do not change (terminal methyl group) but then you do not have
explicit B-state parameters for those interactions that are changing (hence grompp
copies them and warns you that you have different end states that rely on copied
bonded parameters, which is usually inappropriate)

all the parameters that change do have different numbers in A and B, those that
do not change have only paramters in A, or in both (the same) - the latter is
an artifact due to the algorithm i used for making the merged topology,
but it should not matter as having only a parameter in A, or the same number
in both A and B has the same effect.

For instance, in ethane (A-state, λ = 0), there should not be an O-H bond.

well, this is the point about dual topology approach, ALL atoms are there in BOTH
end-states - but in each of the two end-states different atoms are made dummies.
the advantage of this is, at least as far as I see it, that you have no problems when
converting constraints in harmonic bonds, or vice versa, also making or breaking
rings is straight forward. With the single topology approach you took in your tutorial
both these cases would be problematic, to say the least - and both occur quite
frequently when you do typical alchemical transformations of organic molecules.
For more details see the publications above. I guess I should stop here, this is
becoming somewhat off-topic …
best,
michael

6

For cases 1 and 3, here is what I get with GROMACS 2021.3:

Case 1:

   Energies (kJ/mol)
           Bond          Angle    Proper Dih.          LJ-14     Coulomb-14
    9.10096e-04    4.90315e-01    1.08333e-05    1.55416e-01    2.00272e+00
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential
    0.00000e+00   -5.57615e-02    5.42466e-02    2.64785e+00

Case 3:

           Bond          Angle    Proper Dih.          LJ-14     Coulomb-14
    2.14947e-03    5.46054e-01    1.00363e-04    1.55416e-01    2.00272e+00
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    dVremain/dl
    0.00000e+00   -5.57610e-02    5.42466e-02    2.70492e+00    1.03866e+01

As is clear from these results, the Bond, Angle, and Proper Dih. energies are all larger when using the free energy topology, implying there are additional interactions being considered that shouldn’t. The fact that the LJ and Coulombic terms are all identical means your topology setup is a frankenmolecule of bonded ethane+some methanol and entirely nonbonded interactions of ethane.

But this is not what you are doing; you are in fact doing the opposite. Consider your [bonds] directive:

[ bonds ]
       1       2   1   1.0970e-01   3.1455e+05   1.0970e-01   3.1455e+05 ; both 1-2 1-2
       1       3   1   1.0970e-01   3.1455e+05   1.0970e-01   3.1455e+05 ; both 1-3 1-3
       1       4   1   1.0970e-01   3.1455e+05   1.0970e-01   3.1455e+05 ; both 1-4 1-4

Here, the terminal methyl group is the only thing that is not changing, so there is no need to explicitly define the A- and B-states. Everything else below does change when transforming as a function of λ but you’re only ever assigning an A-state to the bonded parameters. So these get copied, and you have all the ethane bonds and all the methanol bonds present in the A-state (along with angles and dihedrals, but I won’t list those here).

       1       5   1   1.5380e-01   1.9456e+05 ; only 0 1-5
       5       6   1   1.0970e-01   3.1455e+05 ; only 0 5-6
       5       7   1   1.0970e-01   3.1455e+05 ; only 0 5-7
       5       8   1   1.0970e-01   3.1455e+05 ; only 0 5-8
       1       9   1   1.4230e-01   2.4552e+05 ; only 1 1-5
       9      10   1   9.7300e-02   4.7154e+05 ; only 1 5-6

Everything above needs to have altered A- and B-state parameters explicitly defined. That’s what grompp checks for - it tells you when an atom type has changed but the associated bonded parameters have not, because this is unexpected (may be correct in some cases but not here).

I am not saying one should use my tutorial approach for your system; indeed that would be entirely inappropriate for the type of calculation here (relative free energy difference) vs. what is in the tutorial (absolute free energy difference). I’m quite well versed in doing dual-topology calculations in GROMACS, hence why I am offering some insight here without ever once referring to the tutorial. The key aspect of the papers you are mentioning is that the bonded connectivity must be defined for both states simultaneously, but not that the parameters should be the same. Indeed, they should be different as a function of the two end states.

7

Dear Justin,

Thanks your comments - I believe that there is a misunderstanding here.
you wrote:

As is clear from these results, the Bond, Angle, and Proper Dih. energies are all
larger when using the free energy topology, implying there are additional interactions
being considered that shouldn’t. The fact that the LJ and Coulombic terms are all
identical means your topology setup is a frankenmolecule of bonded ethane+some
methanol and entirely nonbonded interactions of ethane.

there are indeed additional interactions, and they SHOULD be considered, otherwise
the molecules would fall apart/disintegrate. In a dual topology approach you need
to keep all bonds, and the corresponding interactions, the crucial point is:
if you do a double annihilation, i.e., calculate a DeltaDeltaG (relative free energy
of solvation, or relative binding free energy) as opposed to a single free energy
(one leg in a thermodynamic cycle) then these additional contributions cancel out,
which is why a double topology approach works.
If you changed the bonded interactions, as you suggested, then in one end-state
some of these interactions would be gone altogether … consider the the O-H bond,
in my hybrid topology the force constant and the equilibrium bond length are identical
throughout - if I wanted to change these parameters, as you suggest, which values
would you use at lambda=0 (ethane)? … If I follow your logic, then at lambda=0
the force constant of this interaction would be ZERO … that is the bond would be gone
and I’d end up with two molecules instead of one, and I cannot see how you’d
account for this in terms of the calculated free energy…

you also wrote:

I’m quite well versed in doing dual-topology calculations in GROMACS

do you perhaps have any reference there? … perhaps this would help to
clarify this misunderstanding.

cheers
Michael

8

Perhaps my focus on bonds got things a bit derailed, because you can certainly define arbitrary parameters for the structure of the dummies that then become the B-state, but you still have a large number of interactions that are actually perturbed that you’re not transforming. Consider your pair interactions, which include H-H in the B-state. You’re not defining anything there, so the A-state (dummy, zero LJ) gets copied to the B-state (which are H, but then have zero LJ). These need to be accounted for because they are actual transformations.

I worked with the CHARMM-GUI team to get things working in their new(ish) free energy setup. I gave them examples of working topologies that matched published values for simple transformations like these. I haven’t checked the status of whether or not they’ve got the GROMACS input generation publicly available, but that’s where I’ve done a bunch of this work.