[gmx-users] FEP

[Michael Shirts]

> So BAR is only
> referring to the mathematical code used to calculate the overall free
> energy for the FEP, correct?

Yes.  The information required is the same, with the exception that
exponential averaging requires energy differences from only one state,
whereas BAR requires energy differences from both directions.  Since
you are likely running simulations at a range of states anyway,
there's no extra cost.

> My question is, for a mutation of a -CH3 group to a -H group, is it
> better to simply run:
> [+ from (Lambda=0 ,  R-CH3, full charges and interactions -STATE A)
> --> (Lambda=1, R-CH, full charges and interactions -STATE B)]
>
> OR
>
> [1) from (Lambda=0 ,  R-CH3, STATE A : Charges and LJ Interactions: ON)
> 2)  (Lambda=?, -CH3 Charges: OFF ,LJ Interactions: ON and unmutated)
> 3)  (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ Interactions: OFF)
> 4)  (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ interactions: ON and Mutated to -H)
> 5)  (Lambda=1, R-CH3, -CH3  STATE B : Charges and LJ Interactions: ON)

Currently, you'd need to run three simulations to do this (will likely
be simplified somewhat in 4.6).

Set 1: turn R-CH3 charges off in a way that preserves the total charge.
Set 2: change CH3 LJ to H LJ
Set 3: Turn R-H charges on in a way that preserves the total charge.

Note that the first and last transformations will need only one two
states -- their energies will be VERY small.

> Reason I'm asking is because when I try the first choice to do it
> STATE A to STATE B in one step, when I reach Lambda=0.85 and above on
> the NVT equilibration right after EM, I receive errors saying that
> bonds are moving way to far off their constraints which leads me to
> believe that the system is moving too far from where it was energy
> minimized.  Errors such as:
>
> Step 188, time 0.376 (ps)
> LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.000017, max 0.000636 (between atoms 9 and 68)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>      9     68   31.2    0.1111   0.1110      0.1111
>
> Step 188, time 0.376 (ps)  LINCS WARNING
> relative constraint deviation after LINCS:
> rms 0.000015, max 0.000531 (between atoms 9 and 68)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>      9     68   31.0    0.1111   0.1110      0.1111

You probably have problems with the charge + LJ terms on the hydrogen.
 Note that you could possibly use soft core for the Coloumb
interactions for hydrogen -- it's sometimes causes sampling problems,
but for H's it probably doesn't matter.  I haven't done this approach
much so, so I can't tell for sure if it will work.