Hi Fabian,<div>I am trying something similar with Glutamate to Alanine mutation.</div><div>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?</div><div><br>
</div><div>Thanks for the time,</div><div><br></div><div>Sai<br><br><div class="gmail_quote">On Thu, Apr 26, 2012 at 10:43 AM, Fabian Casteblanco <span dir="ltr"><<a href="mailto:fabian.casteblanco@gmail.com" target="_blank">fabian.casteblanco@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello all,<br>
<br>
This is in reply to Michael shirts a while ago on a FEP of a R-CH3 to<br>
an R-H group. Below is the orignal email. I recently tested out<br>
mutating a CH3-CH3-(3dummy atoms) molecule on both sides in order to<br>
test out that a peturbation would give you a total of 0. forcefield<br>
used was CGenFF<br>
<br>
So for procedure was that I turned the CH3 group on the left to an H<br>
with 3 dummy atoms and I turned an H and 3 dummy atoms on the right to<br>
a CH3 group, essentially mutating the molecule to another version of<br>
itself (should give you a total dG of 0 if it works).<br>
<br>
STEP 1. (uncharging everything except the -CH2 group inbetween both<br>
sides .... CH3-CH2-HD3) *D are dummy atoms<br>
[ atoms ]<br>
; nr type resnr resid atom cgnr charge mass total_charge<br>
1 CG331 1 SIM C1 1 -0.27<br>
12.0110 CG331 0.00 12.0110<br>
2 HGA3 1 SIM H11 2 0.09<br>
1.00800 HGA3 0.00 1.00800<br>
3 HGA3 1 SIM H12 3 0.09<br>
1.00800 HGA3 0.00 1.00800<br>
4 HGA3 1 SIM H13 4 0.09<br>
1.00800 HGA3 0.00 1.00800<br>
5 CG331 1 SIM C2 5 -0.27 12.0110 CG331<br>
-0.18 12.0110<br>
6 HGA3 1 SIM H21 6 0.09<br>
7 HGA3 1 SIM H22 7 0.09<br>
8 HGA3 1 SIM H23 8 0.09 1.00800 HGA3 0.00 1.00800<br>
9 DUM 1 SIM H31 9 0.00 1.00800<br>
10 DUM 1 SIM H32 10 0.00 1.00800<br>
11 DUM 1 SIM H33 11 0.00 1.00800<br>
<br>
STEP 2. (mutation of groups. CH3-CH2-HD3 becomes D3H-CH2-CH3)<br>
[ atoms ]<br>
; nr type resnr resid atom cgnr charge mass total_charge<br>
1 CG331 1 SIM C1 1 0.00<br>
12.0110 HGA3 0.00 1.00800<br>
2 HGA3 1 SIM H11 2 0.00<br>
1.00800 DUM 0.00 1.00800<br>
3 HGA3 1 SIM H12 3 0.00<br>
1.00800 DUM 0.00 1.00800<br>
4 HGA3 1 SIM H13 4 0.00 1.00800 DUM<br>
0.00 1.00800<br>
5 CG331 1 SIM C2 5 -0.18<br>
6 HGA3 1 SIM H21 6 0.09<br>
7 HGA3 1 SIM H22 7 0.09<br>
8 HGA3 1 SIM H23 8 0.00<br>
1.00800 CG331 0.00 12.0110<br>
9 DUM 1 SIM H31 9 0.00 1.00800 HGA3 0.00<br>
1.00800<br>
10 DUM 1 SIM H32 10 0.00 1.00800 HGA3 0.00<br>
1.00800<br>
11 DUM 1 SIM H33 11 0.00<br>
1.00800 HGA3 0.00 1.00800<br>
<br>
STEP3. (opposite of step1)<br>
[ atoms ]<br>
; nr type resnr resid atom cgnr charge mass total_charge<br>
1 HGA3 1 SIM C1 1 0.00 12.0110 HGA3<br>
0.09 1.00800<br>
2 DUM 1 SIM H11 2 0.00<br>
1.00800 DUM 0.00 1.00800<br>
3 DUM 1 SIM H12 3 0.00<br>
1.00800 DUM 0.00 1.00800<br>
4 DUM 1 SIM H13 4 0.00<br>
1.00800 DUM 0.00 1.00800<br>
5 CG331 1 SIM C2 5 -0.18 12.0110 CG331<br>
-0.27 12.0110<br>
6 HGA3 1 SIM H21 6 0.09<br>
7 HGA3 1 SIM H22 7 0.09<br>
8 CG331 1 SIM H23 8 0.00<br>
1.00800 CG331 -0.27 12.0110<br>
9 HGA3 1 SIM H31 9 0.00<br>
1.00800 HGA3 0.09 1.00800<br>
10 HGA3 1 SIM H32 10 0.00<br>
1.00800 HGA3 0.09 1.00800<br>
11 HGA3 1 SIM H33 11 0.00<br>
1.00800 HGA3 0.09 1.00800<br>
<br>
So in the end, the procedure worked using the g_bar method. I took 20<br>
simulations for steps 1 and 3 and 50 simulations for step 2. Here are<br>
the results:<br>
<br>
Step1: -37.62 +/- 0<br>
Step2: -0.12 +/- 0.24<br>
Step3: 37.68 +/- 0.13<br>
TOTAL ~ -0.06 kJ/mol which is expected<br>
<br>
My question is that I did something similar for a larger molecule (65<br>
atoms) and someone said before that step 1 and 3 should be very small.<br>
I'm getting values in the near -20 kJ/mol for step 1 and 3. I feel<br>
it should be straightforward since I'm simply uncharging a CH3 group<br>
and I don't get why I'm getting such a large number. If it is suppose<br>
to be very small, is it safe to assume that step2 should be a majority<br>
of the final value.<br>
<br>
If anyone has any experience in this area all help is greatly<br>
appreciated. Thank you very much.<br>
<br>
-Fabian<br>
<br>
<br>
<br>
<br>
> So BAR is only<br>
> referring to the mathematical code used to calculate the overall free<br>
> energy for the FEP, correct?<br>
<br>
Yes. The information required is the same, with the exception that<br>
exponential averaging requires energy differences from only one state,<br>
whereas BAR requires energy differences from both directions. Since<br>
you are likely running simulations at a range of states anyway,<br>
there's no extra cost.<br>
<br>
> My question is, for a mutation of a -CH3 group to a -H group, is it<br>
> better to simply run:<br>
> [+ from (Lambda=0 , R-CH3, full charges and interactions -STATE A)<br>
> --> (Lambda=1, R-CH, full charges and interactions -STATE B)]<br>
><br>
> OR<br>
><br>
> [1) from (Lambda=0 , R-CH3, STATE A : Charges and LJ Interactions: ON)<br>
> 2) (Lambda=?, -CH3 Charges: OFF ,LJ Interactions: ON and unmutated)<br>
> 3) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ Interactions: OFF)<br>
> 4) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ interactions: ON and Mutated to -H)<br>
> 5) (Lambda=1, R-CH3, -CH3 STATE B : Charges and LJ Interactions: ON)<br>
<br>
Currently, you'd need to run three simulations to do this (will likely<br>
be simplified somewhat in 4.6).<br>
<br>
Set 1: turn R-CH3 charges off in a way that preserves the total charge.<br>
Set 2: change CH3 LJ to H LJ<br>
Set 3: Turn R-H charges on in a way that preserves the total charge.<br>
<br>
Note that the first and last transformations will need only one two<br>
states -- their energies will be VERY small.<br>
<br>
> Reason I'm asking is because when I try the first choice to do it<br>
> STATE A to STATE B in one step, when I reach Lambda=0.85 and above on<br>
> the NVT equilibration right after EM, I receive errors saying that<br>
> bonds are moving way to far off their constraints which leads me to<br>
> believe that the system is moving too far from where it was energy<br>
> minimized. Errors such as:<br>
><br>
> Step 188, time 0.376 (ps)<br>
> LINCS WARNING<br>
> relative constraint deviation after LINCS:<br>
> rms 0.000017, max 0.000636 (between atoms 9 and 68)<br>
> bonds that rotated more than 30 degrees:<br>
> atom 1 atom 2 angle previous, current, constraint length<br>
> 9 68 31.2 0.1111 0.1110 0.1111<br>
><br>
> Step 188, time 0.376 (ps) LINCS WARNING<br>
> relative constraint deviation after LINCS:<br>
> rms 0.000015, max 0.000531 (between atoms 9 and 68)<br>
> bonds that rotated more than 30 degrees:<br>
> atom 1 atom 2 angle previous, current, constraint length<br>
> 9 68 31.0 0.1111 0.1110 0.1111<br>
<br>
You probably have problems with the charge + LJ terms on the hydrogen.<br>
Note that you could possibly use soft core for the Coloumb<br>
interactions for hydrogen -- it's sometimes causes sampling problems,<br>
but for H's it probably doesn't matter. I haven't done this approach<br>
much so, so I can't tell for sure if it will work.<br>
<br>
<br>
--<br>
Best regards,<br>
<br>
Fabian F. Casteblanco<br>
Rutgers University --<br>
Chemical Engineering<br>
<span class="HOEnZb"><font color="#888888">--<br>
gmx-users mailing list <a href="mailto:gmx-users@gromacs.org">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/Search" target="_blank">http://www.gromacs.org/Support/Mailing_Lists/Search</a> before posting!<br>
Please don't post (un)subscribe requests to the list. Use the<br>
www interface or send it to <a href="mailto:gmx-users-request@gromacs.org">gmx-users-request@gromacs.org</a>.<br>
Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" target="_blank">http://www.gromacs.org/Support/Mailing_Lists</a><br>
</font></span></blockquote></div><br></div>