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.castebla...@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@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-requ...@gromacs.org. > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >
-- gmx-users mailing list gmx-users@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-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
