[gmx-users] FEP
Sai Kumar Ramadugu
sramadugu at gmail.com
Fri Jun 15 00:52:29 CEST 2012
Hi Fabian,
I am trying something similar with Glutamate to Alanine mutation.
Does your dummy atoms i.e., DUM1 have a value of 0.0 for sigma and epsilon
during all three steps or only in step 2?
Thanks for the time,
Sai
On Thu, Apr 26, 2012 at 10:43 AM, Fabian Casteblanco <
fabian.casteblanco at gmail.com> wrote:
> Hello all,
>
> This is in reply to Michael shirts a while ago on a FEP of a R-CH3 to
> an R-H group. Below is the orignal email. I recently tested out
> mutating a CH3-CH3-(3dummy atoms) molecule on both sides in order to
> test out that a peturbation would give you a total of 0. forcefield
> used was CGenFF
>
> So for procedure was that I turned the CH3 group on the left to an H
> with 3 dummy atoms and I turned an H and 3 dummy atoms on the right to
> a CH3 group, essentially mutating the molecule to another version of
> itself (should give you a total dG of 0 if it works).
>
> STEP 1. (uncharging everything except the -CH2 group inbetween both
> sides .... CH3-CH2-HD3) *D are dummy atoms
> [ atoms ]
> ; nr type resnr resid atom cgnr charge mass total_charge
> 1 CG331 1 SIM C1 1 -0.27
> 12.0110 CG331 0.00 12.0110
> 2 HGA3 1 SIM H11 2 0.09
> 1.00800 HGA3 0.00 1.00800
> 3 HGA3 1 SIM H12 3 0.09
> 1.00800 HGA3 0.00 1.00800
> 4 HGA3 1 SIM H13 4 0.09
> 1.00800 HGA3 0.00 1.00800
> 5 CG331 1 SIM C2 5 -0.27 12.0110
> CG331
> -0.18 12.0110
> 6 HGA3 1 SIM H21 6 0.09
> 7 HGA3 1 SIM H22 7 0.09
> 8 HGA3 1 SIM H23 8 0.09 1.00800
> HGA3 0.00 1.00800
> 9 DUM 1 SIM H31 9 0.00 1.00800
> 10 DUM 1 SIM H32 10 0.00 1.00800
> 11 DUM 1 SIM H33 11 0.00 1.00800
>
> STEP 2. (mutation of groups. CH3-CH2-HD3 becomes D3H-CH2-CH3)
> [ atoms ]
> ; nr type resnr resid atom cgnr charge mass total_charge
> 1 CG331 1 SIM C1 1 0.00
> 12.0110 HGA3 0.00 1.00800
> 2 HGA3 1 SIM H11 2 0.00
> 1.00800 DUM 0.00 1.00800
> 3 HGA3 1 SIM H12 3 0.00
> 1.00800 DUM 0.00 1.00800
> 4 HGA3 1 SIM H13 4 0.00 1.00800
> DUM
> 0.00 1.00800
> 5 CG331 1 SIM C2 5 -0.18
> 6 HGA3 1 SIM H21 6 0.09
> 7 HGA3 1 SIM H22 7 0.09
> 8 HGA3 1 SIM H23 8 0.00
> 1.00800 CG331 0.00 12.0110
> 9 DUM 1 SIM H31 9 0.00
> 1.00800 HGA3 0.00
> 1.00800
> 10 DUM 1 SIM H32 10 0.00 1.00800
> HGA3 0.00
> 1.00800
> 11 DUM 1 SIM H33 11 0.00
> 1.00800 HGA3 0.00 1.00800
>
> STEP3. (opposite of step1)
> [ atoms ]
> ; nr type resnr resid atom cgnr charge mass total_charge
> 1 HGA3 1 SIM C1 1 0.00 12.0110
> HGA3
> 0.09 1.00800
> 2 DUM 1 SIM H11 2 0.00
> 1.00800 DUM 0.00 1.00800
> 3 DUM 1 SIM H12 3 0.00
> 1.00800 DUM 0.00 1.00800
> 4 DUM 1 SIM H13 4 0.00
> 1.00800 DUM 0.00 1.00800
> 5 CG331 1 SIM C2 5 -0.18
> 12.0110 CG331
> -0.27 12.0110
> 6 HGA3 1 SIM H21 6 0.09
> 7 HGA3 1 SIM H22 7 0.09
> 8 CG331 1 SIM H23 8 0.00
> 1.00800 CG331 -0.27 12.0110
> 9 HGA3 1 SIM H31 9 0.00
> 1.00800 HGA3 0.09 1.00800
> 10 HGA3 1 SIM H32 10 0.00
> 1.00800 HGA3 0.09 1.00800
> 11 HGA3 1 SIM H33 11 0.00
> 1.00800 HGA3 0.09 1.00800
>
> So in the end, the procedure worked using the g_bar method. I took 20
> simulations for steps 1 and 3 and 50 simulations for step 2. Here are
> the results:
>
> Step1: -37.62 +/- 0
> Step2: -0.12 +/- 0.24
> Step3: 37.68 +/- 0.13
> TOTAL ~ -0.06 kJ/mol which is expected
>
> My question is that I did something similar for a larger molecule (65
> atoms) and someone said before that step 1 and 3 should be very small.
> I'm getting values in the near -20 kJ/mol for step 1 and 3. I feel
> it should be straightforward since I'm simply uncharging a CH3 group
> and I don't get why I'm getting such a large number. If it is suppose
> to be very small, is it safe to assume that step2 should be a majority
> of the final value.
>
> If anyone has any experience in this area all help is greatly
> appreciated. Thank you very much.
>
> -Fabian
>
>
>
>
> > 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.
>
>
> --
> Best regards,
>
> Fabian F. Casteblanco
> Rutgers University --
> Chemical Engineering
> --
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20120614/41bbc6e6/attachment.html>
More information about the gromacs.org_gmx-users
mailing list