Project

General

Profile

Bug #3074

box volume is affected by decoupled molecule when using couple-moltype mdp option

Added by Vytautas Gapsys about 1 month ago. Updated about 2 hours ago.

Status:
Rejected
Priority:
Normal
Assignee:
Category:
mdrun
Target version:
-
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

Alchemically decoupled molecule (via .mdp option couple-moltype) affects the box volume, i.e. at constant pressure the volume of a water box with a decoupled ligand is significantly different from the volume of a pure water box (without a decoupled ligand) with the same number of water molecules.

This effect is not present when decoupling is performed via .itp by explicitly defining the B-state. In turn, the solvation free energy is significantly different for both decoupling methods (.mdp and .itp based).

Attached is the figure illustrating differences in volumes and free energies. Also, attached are the input files for decoupling vdw interactions of anthracene in a water box of 913 molecules: .mdp based decoupling and .itp based decoupling.

volume_vdw_dG.png (280 KB) volume_vdw_dG.png volume and free energy differences for different decoupling schemes Vytautas Gapsys, 09/04/2019 07:20 PM
itp_decoupling.tar.gz (44.7 KB) itp_decoupling.tar.gz input files for .itp based decoupling Vytautas Gapsys, 09/04/2019 07:21 PM
mdp_decoupling.tar.gz (44.6 KB) mdp_decoupling.tar.gz input files for .mdp based decoupling Vytautas Gapsys, 09/04/2019 07:21 PM
itp_decoupling_vacuum.tar.gz (6.87 KB) itp_decoupling_vacuum.tar.gz Vytautas Gapsys, 10/14/2019 06:20 PM
mdp_decoupling_after15ns_vacuum.png (178 KB) mdp_decoupling_after15ns_vacuum.png Vytautas Gapsys, 10/14/2019 09:08 PM
volume_per_water.png (296 KB) volume_per_water.png Vytautas Gapsys, 10/14/2019 09:15 PM
volume_per_water_v2.png (333 KB) volume_per_water_v2.png Vytautas Gapsys, 10/14/2019 09:52 PM
average_energies.png (357 KB) average_energies.png Vytautas Gapsys, 10/15/2019 06:58 PM
dg.png (244 KB) dg.png Vytautas Gapsys, 10/17/2019 09:27 AM
dg_updated_lincs.png (306 KB) dg_updated_lincs.png Vytautas Gapsys, 10/18/2019 09:47 AM
volume_per_water_lincs.png (385 KB) volume_per_water_lincs.png Vytautas Gapsys, 10/18/2019 10:32 AM

History

#1 Updated by Berk Hess 4 days ago

  • Status changed from New to Accepted
  • Assignee set to Berk Hess

I think the decoupled molecule still contributes to the virial and the kinetic energy which goes into the pressure calculation.
There is a technical as well as a user interaction challenge to removing these contributions. If you use couple_moltype is might be clear when you want what to re removed. But the user can also manually modify the A/B states in the topology.

#2 Updated by Berk Hess 4 days ago

  • Status changed from Accepted to Feedback wanted

There is only partial information here. It looks like the difference in this leg is that the itp decoupling decoupled the 1-4 interactions whereas the mdp decoupling does not.
How did you compute the vacuum decoupling parts? Are the 1-4 coupled the same way as in the solvent runs?

#3 Updated by Vytautas Gapsys 4 days ago

For the itp_decoupling, vacuum-leg of the cycle uses the same topology as water-leg (attached).

#4 Updated by Berk Hess 4 days ago

Yes. And for the mdp decoupling in water?

#5 Updated by Vytautas Gapsys 4 days ago

mdp_decoupling uses only one leg of the cycle. Since couple-intramol flag is set to "no", the decoupled state corresponds to the proper vacuum state.

#6 Updated by Berk Hess 4 days ago

Isn't the free-energy issue simply because the molecule folds onto itself and doesn't sample when you do not decouple the intramolecular interactions?

#7 Updated by Vytautas Gapsys 4 days ago

True, this might happen depending on a molecule. However, in the current situation this is not the case: attached is a snapshot of the simulated molecule after 15 ns in the decoupled state using mdp_decoupling, couple-intramol=no. The molecule is rigid, thus the aforementioned problem does not occur.

I think that the first idea about the decoupled molecule still contributing to the virial is likely to be correct. This is further supported by another test where I simulated anthracene in water boxes of different size (attachment volume_per_water.png). In the figure "Decoupled anthracene" denotes mdp_decoupling with anthracene set to be fully decoupled from the system. "Anthracene without nonbondeds" corresponds to the itp_decoupling, where the molecule is also fully decoupled. There is a box size dependence for the mdp_decoupling. Similarly, free energy becomes box size dependent as well, when using mdp_decoupling.

#8 Updated by Berk Hess 4 days ago

If the issue is really caused by this pressure issue, it should go away when you turn on couple-intramol, or?

#9 Updated by Vytautas Gapsys 4 days ago

I considered this option as well, however, with the couple-intramol=yes the effect remains: attached is again volume dependence on the box size with an additional couple-intramol=yes.

From these observations it appears that the issue is triggered by the mdp_decoupling only. Although effectively both decoupling variants should achieve the same decoupled state (when couple-intramol=yes), could it be that the contribution of the mdp_decoupled molecule to the virial has lambda dependence missing?

#10 Updated by Berk Hess 3 days ago

Vytautas Gapsys wrote:

I considered this option as well, however, with the couple-intramol=yes the effect remains: attached is again volume dependence on the box size with an additional couple-intramol=yes.

From these observations it appears that the issue is triggered by the mdp_decoupling only. Although effectively both decoupling variants should achieve the same decoupled state (when couple-intramol=yes), could it be that the contribution of the mdp_decoupled molecule to the virial has lambda dependence missing?

No, the virial is computed from the forces.

Do you see any differences when you compare the average energies from the itp decoupling vs mdp decoupling with couple_intermol=yes? Could you attach log files for these two cases with lambda=1?

#11 Updated by Vytautas Gapsys 3 days ago

To my understanding, the energy analysis shows the expected results (attached average_energies.png).

The overall difference between the itp_decoupling and mdp_decoupling is explained by the intramolecular (MOL-MOL) interactions.

The difference between the couple-intramol=no and couple-intramol=yes is explained by the SR MOL-MOL interactions, whereas 1-4 interactions are retained for both of these schemes.

#12 Updated by Berk Hess 3 days ago

Indeed, the 1-4 is never decoupled with the couple mdp option. But I still don't see how this would lead to a different pressure and/or an incorrect solvation free-energy.

What free-energy difference does the mdp coupling scheme with couple-intramol=yes give?

#13 Updated by Vytautas Gapsys 1 day ago

The free energy summary is now attached (dg.png).
mdp decoupling with couple-intramol=yes has the same behavior as the couple-intramol=no (a techical note: for the couple-intramol=yes, I needed to use the second leg of the cycle in vacuum to make the dG values directly comparable to those from couple-intramol=no).

Considering the results from free energy and potential energy calculations, it seems that the box size dependence could be due to pressure contribution. I would like to revisit the virial question. The intra-molecular forces in "itp decoupling" scheme are not present, thus they don't contribute to the virial. For the "mdp decoupling" schemes there are intra-molecular forces that do contribute to the virial. However, if a molecule is decoupled from the rest of the system, these intra-molecular forces should not have any contribution to pressure. Therefore, the force contribution to the virial ought to carry lambda dependence. Of course, forces themselves do depend on lambda and in the "itp decoupling" scheme this dependence directly translates to the virial dependence on lambda. For the "mdp decoupling", this is not the case for the intra-molecular forces.

All in all, I think that the intra-molecular force contribution for the "mdp decoupling" scheme needs to be explicitly scaled with lambda when calculating the virial.

#14 Updated by Berk Hess 1 day ago

Thanks, now it's clear to me what the results are.

I understand one part of the issue, but not another.

Having a decoupled molecule contribute to the virial is wrong, but the effect should be small. All internal kinetic energy and virial contributions should cancel and only the Ekin of the 3 degrees of freedom of the COM should contribute. This seems consistent with the volume change which I computed from your results as 0.036 nm^3, which at 1 bar gives a free-energy difference of 0.02 kJ/mol. This is a spurious contribution, but is negligible for all practical purposes.

But the difference you get in DG is much larger. I don't understand where this comes from. The 1-4 interactions should have no net contribution to the virial. They should be canceled out by other bonded interactions and constraints. There could be other sources of error, but the box dependence seems to indicate there is a pressure contribution.
Is the volume difference between itp and mdp runs independent of the number of water?

#15 Updated by Berk Hess 1 day ago

I found the volume plot and indeed the volume change seems to be independent of the number of water molecules.

But if the compressibility is 4.6e-5 bar, I get a volume change of 0.006 nm^3 for the Ekin of the COM, which is 6 times smaller than the observed difference. I don't know where this difference comes from.
But as I said, even the 0.036 nm^2 would only give 0.02 kJ/mol free-energy difference.

#16 Updated by Berk Hess 1 day ago

There is an issue with your constraint accuracy.
I ran your molecule in vacuum and get a 7 times larger pressure than what it should, which explains the factor 6 I mentioned before. Lincs with order 4 and 2 iterations is not sufficiently accurate when constraining all bonds. You need at least order=6 and iter=2 or change to h-bond only constraints.
The question is now how this constraint error affects the free-energy difference. The factor 6 is not enough to explain this, but the error might affect dG/dl in other ways as well.

#17 Updated by Vytautas Gapsys about 2 hours ago

Indeed, higher order LINCS setting removed the box size dependence for free energies (dg_updated_lincs.png) and volume (volume_per_water_lincs.png).

Would it make sense to add a suggestion to the manual about the higher lincs order for free energy calculations? Currently, the manual suggests to increase lincs-order only for larger time-steps with virtual sites, BD and accurate energy minimization.

Also, the issue of decoupled molecule contributing to the virial (mdp_decoupling scheme) is still there, isn't it? Could it manifest again even with the higher order lincs, if I were to decouple a much larger molecule, i.e. smaller water/molecule atom ratio?

#18 Updated by Berk Hess about 2 hours ago

  • Status changed from Feedback wanted to Rejected

Vytautas Gapsys wrote:

Indeed, higher order LINCS setting removed the box size dependence for free energies (dg_updated_lincs.png) and volume (volume_per_water_lincs.png).

Great!

Would it make sense to add a suggestion to the manual about the higher lincs order for free energy calculations? Currently, the manual suggests to increase lincs-order only for larger time-steps with virtual sites, BD and accurate energy minimization.

I don't know. I would say that in most cases the issue is using constraints=all-bonds. Most force fields are not parametrized with this, so you should use constraints=h-bonds. Also I think this is only an issue in very rigid system such as you have here. Ideally LINCS would automatically choose the number of iterations based on the eigenvalues and some tolerance, but that's not easy to accomplish.

Also, the issue of decoupled molecule contributing to the virial (mdp_decoupling scheme) is still there, isn't it? Could it manifest again even with the higher order lincs, if I were to decouple a much larger molecule, i.e. smaller water/molecule atom ratio?

The virial contribution is zero on average. The kinetic energy contribution to the pressure is 3kT/V, which is practically negligible.

Also available in: Atom PDF