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

Tue Aug 6 01:46:13 CEST 2013

Hi, all-

That particular redmine was in the case of couple-intramol NOT being
specified.  So it may be different.

Can you upload samples files for this test case to redmine.gromacs.org
with specific instructions?  Please upload both the 'behaving
correctly' and 'not behaving correctly' versions.

Note that in may cases, you can get the correct behavior by changing
the topologies themselves, though this is more complicated.

On Mon, Aug 5, 2013 at 10:25 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> On 8/5/13 9:52 AM, Joerg Sauter wrote:
>>
>> On 08/05/2013 02:15 PM, Justin Lemkul wrote:
>>>
>>>
>>>
>>> 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
>>>
>> Visualization does not show any abnormalities beyond odd confomations in
>> the
>> free energy coupling case.
>>
>> Regarding the energies: I am not sure whether this is normal behaviour,
>> but, in
>> the free energy coupling case I do not see energy values for LJ 14 or
>> Coulomb 14
>> interactions. I thought that by selecting couple-intramol=0 all
>> intramolecular
>> interactions would be replaced by pairtype (14) interactions.
>>
>> In addition, for Coulomb (SR) I get in the free energy case
>> LJ (SR)                    -22.8084      0.076    5.90879 -0.379238
>> (kJ/mol)
>> in comparisson to
>> LJ (SR)                     16.8295       0.81    13.2303 4.51371
>> (kJ/mol)
>> in the regular vacuum simulation. So this looks like the 14 interactions
>> get
>> shifted to the LJ/Coulomb terms in the coupling simulation?! (However, the
>> numbers do not add up)
>>
>
> This is smelling buggy to me, and there have been issues with 1-4
> interactions in the free energy code before
> (http://redmine.gromacs.org/issues/1225), so there could be a lingering
> problem.  Hopefully Michael will chime in.  It might be worthwhile to open
> up a bug report on redmine.gromacs.org with all of the information you have
> provided (input files, output data, and the very nice description of the
> issue you have posted).
>
>
> -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
>
> ==================================================
> --
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the www
> interface or send it to gmx-users-request at gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists