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

[Joerg Sauter]

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)

Best,
Joerg


-- 
Joerg Sauter
Max Planck Institute of Colloids and Interfaces
Department of Theory & Bio-Systems
Am Mühlenberg 1
14476 Potsdam-Golm Germany

Fon: +49-(0)331-567-9622