# With free-energy=yes, angle energies are off **URL:** \ **Category:** User discussions\ **Tags:** mdrun\ **Created:** [December 8, 2021, 10:44pm UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272 "2021-12-08T22:44:40Z")\ **Posts on this page:** 8\ **Page:** 1
**Author:** ![michael.brunsteiner](https://avatars.discourse-cdn.com/v4/letter/m/7cd45c/32.png) [@michael.brunsteiner](https://gromacs.bioexcel.eu/u/michael.brunsteiner)\ **Post date:** [December 8, 2021, 10:44pm UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272/1 "2021-12-08T22:44:41Z")
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](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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [December 9, 2021, 2:32pm UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272/2 "2021-12-09T14:32:24Z")
> [@michael.brunsteiner](#): > > 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 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. ---
**Author:** ![michael.brunsteiner](https://avatars.discourse-cdn.com/v4/letter/m/7cd45c/32.png) [@michael.brunsteiner](https://gromacs.bioexcel.eu/u/michael.brunsteiner)\ **Post date:** [December 9, 2021, 9:55pm UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272/3 "2021-12-09T21:55:12Z")
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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/c/c80d20282d4bd30f7fe4b296b57b01f0053c5ff2.png) 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](https://doi.org/10.1021/jp981628n) and [https://doi.org/10.1021/jp981629f](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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [December 11, 2021, 12:48pm UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272/4 "2021-12-11T12:48:44Z")
> [@michael.brunsteiner](#): > > 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. 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. > [@michael.brunsteiner](#): > > 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. 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. ---
**Author:** ![michael.brunsteiner](https://avatars.discourse-cdn.com/v4/letter/m/7cd45c/32.png) [@michael.brunsteiner](https://gromacs.bioexcel.eu/u/michael.brunsteiner)\ **Post date:** [December 13, 2021, 4:52pm UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272/5 "2021-12-13T16:52:22Z")
> 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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [December 13, 2021, 6:46pm UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272/6 "2021-12-13T18:46:04Z")
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. > [@michael.brunsteiner](#): > > 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) 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). > [@michael.brunsteiner](#): > > 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 … 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. ---
**Author:** ![michael.brunsteiner](https://avatars.discourse-cdn.com/v4/letter/m/7cd45c/32.png) [@michael.brunsteiner](https://gromacs.bioexcel.eu/u/michael.brunsteiner)\ **Post date:** [December 23, 2021, 9:31am UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272/7 "2021-12-23T09:31:23Z")
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 ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [December 30, 2021, 1:54am UTC](https://gromacs.bioexcel.eu/t/with-free-energy-yes-angle-energies-are-off/3272/8 "2021-12-30T01:54:22Z")
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. > [@michael.brunsteiner](#): > > do you perhaps have any reference there? … perhaps this would help to > clarify this misunderstanding. 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.