[gmx-users] Unphysical conformations in decoupled free energy simulation

[Justin Lemkul]

On 8/5/13 8:12 AM, Joerg Sauter wrote:
> On 08/05/2013 01:54 PM, Justin Lemkul wrote:
>>
>>
>> On 8/5/13 7:44 AM, Joerg Sauter wrote:
>>> Dear all,
>>>
>>> I am trying to compute the free energy of hydration for cellobiose (a beta (1-4)
>>> glucose dimer) using BAR in Gromacs 4.6.1. However, I encounter a problem.
>>>
>>> I find that the vacuum conformations of the molecule in a regular vacuum
>>> simulation differ from the conformations in the decoupled simulation in the free
>>> energy case i.e., with the additional mdp entries:
>>>
>>> free-energy              = yes
>>> init-lambda                 = 0
>>> couple-lambda0           = none
>>> couple-moltype           = solute
>>> couple-intramol          = no
>>>
>>> Here is a histogram of the dihedral angles of the glycosidic linkage in
>>> the vacuum case
>>> https://dl.dropboxusercontent.com/u/70358077/reg.pdf (it stays in the global
>>> free energy minimum)
>>> and this is the decoupled free energy case (same starting conformation in the
>>> global free energy minimum)
>>> https://dl.dropboxusercontent.com/u/70358077/fe.pdf
>>>
>>> I understand that Gromacs replaces the intramolecular interactions with explicit
>>> pair
>>> interactions. Therefore, I had to increase table-extension but that did not
>>> change much. I was thinking that maybe this could be a problem specific to this
>>> topology (a GLYCAM06h conversion from Amber), however, I do not see how this can
>>> occur.
>>> I hope someone has an idea what is going wrong.
>>>
>>> mdp file:    https://dl.dropboxusercontent.com/u/70358077/md.mdp (same behaviour
>>> for longer runs)
>>> top file: https://dl.dropboxusercontent.com/u/70358077/cellob.top
>>> gro file: https://dl.dropboxusercontent.com/u/70358077/cellob.gro
>>>
>>> This is sort of a minimal example. The same problem occurs
>>> in a regular simulation when decoupling from water using multiple lambdas.
>>>
>>
>> Have you tried running with rlist=rcoulomb=rvdw = 0 rather than 99? If nothing
>> else, I would try again with version 4.6.3.  There have been tons of fixes to
>> the free energy code since 4.6.1, though I don't know offhand which, if any,
>> would be affecting your results. Regardless, the current version should always
>> be the "best" version ;)
>>
>> -Justin
>>
> Hi Justin,
>
> I just tried it with version 4.6.3. and rlist=rcoulomb=rvdw = 0. Unfortunately,
> I get the same results.
>

OK, good to know.  We can rule out 4.6.1-era bugs, then.  Does visualization or 
an analysis of energy terms stored in the .edr file provide any clues as to the 
source of the problem?

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441

==================================================