[gmx-users] FEP
Fabian Casteblanco
fabian.casteblanco at gmail.com
Thu Apr 26 17:43:23 CEST 2012
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
More information about the gromacs.org_gmx-users
mailing list