Prof Mobley, Sorry to bother you, but do you have any suggestions for my previous post?
Thanks again, Sai On Sun, Mar 31, 2013 at 9:52 PM, Sai Kumar Ramadugu <sramad...@gmail.com>wrote: > Dear Prof Mobley, > > I have some additional questions in understanding the free energy > calculations. > > 1) When doing a mutation from a larger aminoacid to a smaller one, one > necessarily ends up with dummy atoms and atoms that have different mass. In > my case when the glutamate is converted into a glycine or an alanine a > carbon will become a hydrogen (change in mass) and the rest of the excess > atoms will become dummies with the mass they had in Glu. > > The change in mass has an effect in the kinetic energy part of the > Hamiltonian contribution to the free energy and this change is properly > computed by Gromacs as <dKe/dL>. However it is clear from this procedure > that because the dummies have mass, the mass of Ala or Gly will not be the > correct and therefore <dKe/dL> would be quite meaningless. I read papers > that suggest only integrating the potential energy component <dV/dL>. In > fact for example, the Amber11 manual indicates that Amber mutations only > change the potential but keep all masses the same as they were before the > alchemical change. > > Clearly having dummies with no mass would fix this issue but it would be > problematic because F=m*a and with m=0 for a given finite F, the > acceleration on the dummies would be infinite. > > In Monte Carlo all of this is irrelevant because one only computes the > potential contribution to the free energy difference but in MD the story is > different. What is the customary approach to deal with this? Just use > <dV/dL> and discard <dKe/dL>? > > > > 2) The g_bar analysis appears to assume that the interpolation between > states a and b is linear (g_bar issues a warning "using the derivative data > (dH/dlambda) to extrapolate delta H values. This will only work if the > Hamiltonian is linear in lambda." ). However when using softcore potentials > lambda also enters in the expression for "r". It would appear that the > interpolation is not really linear even though the LJ transformation also > looks like > Va *lambda+ (1- lambda)* Vb > Nonetheless I see g_bar being used for such transformations. Could you > shed light on this? > > 3) Also it is common to use small lambda spacings (close to lambda=0 and > lambda=1) but somewhat larger spacing between 0.1 and 0.9. How does g_bar > handle this? Is the integration and error analysis done correctly even when > the lambda spacing may be different at different lambda values? > > Thank you for your time. > > Regards > Sai > > > On Tue, Dec 4, 2012 at 5:27 PM, Sai Kumar Ramadugu <sramad...@gmail.com>wrote: > >> Dear Prof Mobley, >> >> Thank you very much for your time. >> >> While I was waiting for your opinion on my topology files, I did several >> tests to see what happens. >> The tests I did are: >> >> 1. Mutate only the amino acid (Glutamate -1 charge) and no mutation of K+ >> ion, but have positional restraints on the ion. >> 2. Mutate only the amino acid (Glutamate -1 charge) and no mutation of K+ >> ion and no positional restraints or constraints on the system >> 3. The actual simulation that caused me trouble, ie with position >> restraints, constraints and mutating the ion concurrently. >> >> Finally I decided to use the simulation data from no restraints or >> constraints and no mutation of K+ ion. >> With this approach now I am getting a positive deldelG value for the >> mutation going from Glutamate to Alanine. >> >> On Mon, Dec 3, 2012 at 4:07 PM, David Mobley <dmob...@gmail.com> wrote: >> >>> Hi, >>> >>> I had a look at your topologies and don't see anything obviously wrong. >>> You may however want to check that none of the bonded parameters involving >>> the transformed atoms are being unexpectedly perturbed. (In our setups we >>> typically repeat the A state parameters in the B state column for the >>> bonds, angles, and torsions, as some GROMACS versions would otherwise >>> assume these should be turned to zero). >>> >> >> I read in the gromacs manual that with opls_aa force field, the state A >> parameters are repeated if you dont have any parameters for state B. I saw >> those warnings issued by grompp and in order to avoid any problems I did >> repeat the state A parameters in state B columns as you do it. >> >> >>> >>> Sorry I can't be more help. I'm pretty swamped on this end. >>> >>> >>> >>> >> I once again, I thank you very much for the time you took for my problem. >> >> Regards >> Sai >> >> >> >> >>> >>> On Tue, Oct 30, 2012 at 10:35 AM, Sai Kumar Ramadugu < >>> sramad...@gmail.com> wrote: >>> >>>> Prof Mobley, >>>> >>>> I just want to add that in the topology file named Ala_charge_off.top, >>>> I am *turning off* the charges on atoms of alanine 40 residue. The >>>> total integral of <dV/dl> from 0 to 1 gave me a value of 274.73 kJ/mol >>>> and I took the negative of this as this should be same as *turning on*the >>>> charges. To justify my approach I did add the partial charges on Ala40 >>>> residue in a separate set of simulations (which finished today morning) and >>>> the total integral of <dV/dl> from 0 to 1 gave me a value of -274.705 >>>> kJ/mol. >>>> >>>> >>>> Regards >>>> Sai >>>> >>>> >>>> On Sat, Oct 27, 2012 at 11:00 AM, Sai Kumar Ramadugu < >>>> sramad...@gmail.com> wrote: >>>> >>>>> Dear Prof Mobley, >>>>> >>>>> I have attached the topologies for all three steps in presence of >>>>> ligand. The case for the absence of ligand is same. >>>>> The Glutamate happens to be residue 40 in the topology file (around >>>>> line 615). >>>>> In the case of LJ transformation, for the bonds with dummies, I am >>>>> using the same bond length, angles and dihedral angles as it is for their >>>>> respective original atoms. >>>>> I am using the two papers by Martin Karplus to justify this approach. >>>>> "The role of bonded terms in free energy simulations 1: Theoretical >>>>> Analysis", JPCA, 1999, 103, 103-118. >>>>> >>>>> >>>>> On Fri, Oct 26, 2012 at 1:32 PM, David Mobley <dmob...@gmail.com>wrote: >>>>> >>>>>> Hi, >>>>>> >>>>>> A few responses below. I'm afraid I don't have any brilliant insights >>>>>> right now but there are some things to look at... >>>>>> >>>>>> On Thu, Oct 25, 2012 at 11:05 AM, Sai Kumar Ramadugu < >>>>>> sramad...@gmail.com> wrote: >>>>>> >>>>>>> Dear Prof Mobley, >>>>>>> I am sorry for the delay in answering your questions to my original >>>>>>> post. >>>>>>> The reason being, I was doing some calculations to answer some of >>>>>>> your questions. >>>>>>> >>>>>>> *My question*: I am getting negative values for vdw transformation >>>>>>> in gromacs when I am mutating sidechain of Glu to Ala. The overall free >>>>>>> energy difference should be positive but I am getting negative and this >>>>>>> is >>>>>>> mostly effected due to the negative free energy values from LJ >>>>>>> calculation. >>>>>>> Is my LJ transformation correct? >>>>>>> >>>>>> >>>>>> I don't know, honestly. If you're getting a sign that's different >>>>>> from what you expect, it certainly suggests a problem, though it's always >>>>>> possible that the force field disagrees with experiment. At this point, >>>>>> though, I'd be looking for problems with your topologies if I were you (I >>>>>> can't do this since you haven't sent them). >>>>>> >>>>>> >>>>>>> The values themselves: >>>>>>> delG1 (in presence of ligand) Step1 = 1041.31 kJ/mol; Step2: = >>>>>>> -27.08 kJ/mol; Step3 = -274.73 kJ/mol Total = 739.5 >>>>>>> delG2 (in absence of ligand) Step1 = 1037.38 kJ/mol; Step2: = >>>>>>> -20.93 kJ/mol; Step3 = -271.86 kJ/mol Total = 744.59 >>>>>>> >>>>>>> The relative free energy difference: delG1 - delG2 = 739.5 - 744.59 >>>>>>> = -5.09 kJ/mol. For the transformation of Glu to Ala this value has to >>>>>>> be >>>>>>> positive. >>>>>>> >>>>>>> I'm skeptical that the differences in step 1 and step 3 are >>>>>> necessarily meaningful, though they do seem to cancel out. I think you're >>>>>> right to focus on step 2. >>>>>> >>>>>>> >>>>>>> My thermodynamic cycle: >>>>>>> delG1 >>>>>>> Protein-Glutamate-Ligand--------------> Protein-Alanine_ligand >>>>>>> | | >>>>>>> |delG3 |delG4 >>>>>>> | | >>>>>>> Protein-Glutamate ----------------> Protein-Ala >>>>>>> delG2 >>>>>>> I am using delG1 - delG2 for the relative binding free energy of >>>>>>> ligand bound to protein with glutamate to protein with alanine. >>>>>>> >>>>>>> Initially I did a two step approach i.e., I have been doing the >>>>>>> uncharging and charging mutation of Glu--->Ala in one step. >>>>>>> Later I read a post by Prof Shirts in response to Fabian in this >>>>>>> users-list and followed a three step approach. ( >>>>>>> http://gromacs.5086.n6.nabble.com/FEP-td4931063.html). These >>>>>>> calculations just started when you responded. My apologies! >>>>>>> >>>>>>> Step1: Uncharging of Glu >>>>>>> Step2: Vdw mutation of side chain of Glu to Ala >>>>>>> Step3: Charging of Ala >>>>>>> >>>>>>> If you want me to troubleshoot, it'd be nice to see topologies for >>>>>> these, especially step 2, but actually all of them. Most problems result >>>>>> from not setting up the right transformation. >>>>>> >>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>>> >>>>>>>> >>>>>>>>> The topology section of mutating the potassium ion K+ to K0 atom >>>>>>>>> is pasted below: >>>>>>>>> >>>>>>>>> [ moleculetype ] >>>>>>>>> ; Name nrexcl >>>>>>>>> KM 1 >>>>>>>>> >>>>>>>>> [ atoms ] >>>>>>>>> ; nr type resnr residue atom cgnr charge >>>>>>>>> mass typeB chargeB massB >>>>>>>>> 1 opls-408 1 KM KM 1 1.0 >>>>>>>>> 39.0983 DUM_408 0.0 39.0983 >>>>>>>>> >>>>>>>>> >>>>>>>>> I am doing this calculation in presence and absence of the ligand. >>>>>>>>> >>>>>>>>> I take it your position restraints are imposed via posre.itp or >>>>>>>> something similar? Also, I hope you are holding the ion and protein >>>>>>>> apart >>>>>>>> from one another, rather than just holding the ion fixed and allowing >>>>>>>> the >>>>>>>> protein to drift? How exactly are you handlign restraints? >>>>>>>> >>>>>>> >>>>>>> For the position restraints, I am using the restraints below the >>>>>>> [molecule] and [atoms] directive with ifdef command and giving a value >>>>>>> of >>>>>>> 1000 for fx, fy and fz. In the production mdp input file, I am using >>>>>>> -DPOSRES to implement the same. I am pasting the relevant part of the >>>>>>> topology for position restraints. >>>>>>> >>>>>>> #ifdef POSRES_ION >>>>>>> ; Position restraint for each Potassium ion >>>>>>> [ position_restraints ] >>>>>>> ; i funct fcx fcy fcz >>>>>>> 1 1 1000 1000 1000 >>>>>>> #endif >>>>>>> >>>>>>> >>>>>>> Further I am not restraining the protein, but in all the cases, the >>>>>>> protein is not drifting during the simulation time I have which is 2ns. >>>>>>> I >>>>>>> have checked this with distance between the atoms of mutated amino acid >>>>>>> and >>>>>>> the mutating ion (and also other ions in the system). >>>>>>> >>>>>> >>>>>> Well, that's fine as far as it goes I suppose, but you've set up a >>>>>> case where the ion will drift with respect to the protein (as the protein >>>>>> moves) which will likely cause convergence problems in longer >>>>>> simulations. >>>>>> I don't think this is the problem, but it's certainly A (probably minor) >>>>>> problem. >>>>>> >>>>> >>>>> I am thinking to restart my simulations with restraints on the C-alpha >>>>> atoms of the protein. >>>>> >>>>>> >>>>>> >>>>>>> >>>>>>>> >>>>>>>>> Step 2: LJ transformation calculation >>>>>>>>> >>>>>>>>> The relevant part of the topology for LJ transformation is pasted >>>>>>>>> below: >>>>>>>>> >>>>>>>>> ; residue 40 GLU rtp GLU q -1.0 >>>>>>>>> 552 opls_238 40 GLU N 191 0.0 >>>>>>>>> 14.0067 ; >>>>>>>>> 553 opls_241 40 GLU H 191 0.0 >>>>>>>>> 1.008 ; >>>>>>>>> 554 opls_224B 40 GLU CA 191 0.0 12.011 >>>>>>>>> ; >>>>>>>>> 555 opls_140 40 GLU HA 191 0.0 >>>>>>>>> 1.008 ; >>>>>>>>> 556 opls_136 40 GLU CB 192 0.0 >>>>>>>>> 12.011 ; >>>>>>>>> 557 opls_140 40 GLU HB1 192 0.0 1.008 >>>>>>>>> ; >>>>>>>>> 558 opls_140 40 GLU HB2 192 0.0 1.008 >>>>>>>>> ; >>>>>>>>> 559 opls_274 40 GLU CG 193 0.0 >>>>>>>>> 12.011 opls_140 0.0 1.008 ; qtot 0.78 >>>>>>>>> 560 opls_140 40 GLU HG1 193 0.0 1.008 >>>>>>>>> DUM_140 0.0 1.008 ; qtot 0.84 >>>>>>>>> 561 opls_140 40 GLU HG2 193 0.0 1.008 >>>>>>>>> DUM_140 0.0 1.008 ; qtot 0.9 >>>>>>>>> 562 opls_271 40 GLU CD 194 0.0 >>>>>>>>> 12.011 DUM_271 0.0 12.011 ; qtot 1.6 >>>>>>>>> 563 opls_272 40 GLU OE1 194 0.0 15.9994 >>>>>>>>> DUM_272 0.0 15.9994 ; qtot 0.8 >>>>>>>>> 564 opls_272 40 GLU OE2 194 0.0 15.9994 >>>>>>>>> DUM_272 0.0 15.9994 ; qtot 0 >>>>>>>>> 565 opls_235 40 GLU C 195 0.0 >>>>>>>>> 12.011 ; >>>>>>>>> 566 opls_236 40 GLU O 195 0.0 >>>>>>>>> 15.9994 ; >>>>>>>>> >>>>>>>>> Since I am only changing the side chain of Glu to side chain of >>>>>>>>> Ala, I have mutated only the part of the side chain that is different >>>>>>>>> between the two amino acids. >>>>>>>>> >>>>>>>> >>>>>>>> Have you now turned off the LJ interactions on your dummy atoms? >>>>>>>> Or, when do you finish turning off your atoms? >>>>>>>> >>>>>>> >>>>>>> I have turned off the sigma and epsilon at this stage in the >>>>>>> ffnonbonded.itp file. For all the dummy atoms in my system, sigma and >>>>>>> epsilon are zero. >>>>>>> >>>>>> >>>>>> Didn't you say in your LJ transformation you are turning Glu into >>>>>> Ala? Just to confirm, you are turning a carbon into a hydrogen at this >>>>>> stage as well as changing atoms into dummies, right? >>>>>> >>>>> >>>>> Yes, I am mutating the C-gamma to Hydrogen and the rest of the side >>>>> chain atoms of Glutamate to dummies. >>>>> >>>>>> >>>>>>> >>>>>>> >>>>>>>> >>>>>>>> >>>>>>>>> At the same time I am also mutating the K0 atom to a dummy atom. >>>>>>>>> The topology section of mutating the potassium atom to a dummy atom is >>>>>>>>> pasted below: >>>>>>>>> >>>>>>>>> [ moleculetype ] >>>>>>>>> ; Name nrexcl >>>>>>>>> KM 1 >>>>>>>>> >>>>>>>>> [ atoms ] >>>>>>>>> ; nr type resnr residue atom cgnr charge >>>>>>>>> mass typeB chargeB massB >>>>>>>>> 1 opls-408 1 KM KM 1 0.0 >>>>>>>>> 39.0983 DUM_408 0.0 39.0983 >>>>>>>>> >>>>>>>>> >>>>>>>>> I am doing this transformation in presence and absence of the >>>>>>>>> ligand. >>>>>>>>> >>>>>>>>> After these two steps: >>>>>>>>> For the analysis I am just using the values of dV/dl printed for >>>>>>>>> every 10 steps from the simulation from 0 to 1 in lambda and >>>>>>>>> integrating >>>>>>>>> the dV/dl w.r.t. lambda. >>>>>>>>> >>>>>>>>> Step 1 charge in presence of Ligand = 790.109 kJ/mol >>>>>>>>> Step 2 vdw in presence of Ligand = -29.324 kJ/mol >>>>>>>>> >>>>>>>>> The sum of two steps in presence of ligand = 760.785 kJ/mol >>>>>>>>> >>>>>>>>> Step1 charge in absence of Ligand = 787.33 kJ/mol >>>>>>>>> Step 2 vdw in absence of Ligand = -21.8127 kJ/mol >>>>>>>>> >>>>>>>>> The sum of two steps in absence of ligand = 765.517 kJ/mol >>>>>>>>> >>>>>>>>> I'm confused. It looks to me like you should be taking three steps >>>>>>>> Glu -> Ala: Turn off charges on atoms in Glu being deleted while >>>>>>>> turning >>>>>>>> off charges on ion; second, turn off LJ on any atoms in Glu being >>>>>>>> deleted; >>>>>>>> third, change any charges on atoms retained in Ala from what they were >>>>>>>> in >>>>>>>> Glu to what they should be in Ala. I think you are missing some of that >>>>>>>> unless there's something you haven't explained. >>>>>>>> >>>>>>> >>>>>>> Here is the part where I did only two steps and later read a post by >>>>>>> Prof Shirts and started doing the three step approach as you wrote. The >>>>>>> delay for the reply is because of doing new set of calculations. For >>>>>>> step 1 >>>>>>> I have topology and coordinates corresponding to Glu residue and for >>>>>>> step 3 >>>>>>> I have topology and coordinates for Ala residue and in both cases I am >>>>>>> changing only the partial charges on the atoms. Did I misunderstand >>>>>>> anything here? For step 2, I have topology corresponding to Glu and I am >>>>>>> deleting the additional atoms in glutamate to make it alanine. >>>>>>> >>>>>>> >>>>>>>> >>>>>>>> >>>>>>>>> The relative free energy of the mutation Glu-->Ala = -4.732 kJ/mol >>>>>>>>> >>>>>>>>> My main concern is with respect to the LJ transformation. Is my >>>>>>>>> approach correct in terms of modifying the side chain of glutamate to >>>>>>>>> alanine? The doubt arises because the relative free energy difference >>>>>>>>> is >>>>>>>>> negative where as the experimental value is close to 20 kJ/mol. I am >>>>>>>>> way >>>>>>>>> under-predicting the value and with a negative sign. >>>>>>>>> >>>>>>>> >>>>>>>> You haven't discussed how you get the total free energy of this. I >>>>>>>> assume you're doing some solvent calculation as well? >>>>>>>> >>>>>>> >>>>>>> I am using the <dV/dl> vs lambda for lambda 0 to 1 at various lambda >>>>>>> values ( I have 38 lambdas for vdw calculations and 21 lambda values for >>>>>>> charging/uncharging calculations) and integrating the curve from 0 to >>>>>>> 1. I >>>>>>> have tried the gbar approach too. But my integrated values and gbar >>>>>>> values >>>>>>> for free energy were almost same. But I can use gbar if needed. >>>>>>> >>>>>>> I am using the thermodynamic cycle as shown above (similar to >>>>>>> gromacs manual figure 3.10 A). What additional solvent calculation do I >>>>>>> need to do when I am taking the difference delG1 - delG2? >>>>>>> >>>>>> >>>>>> Sorry, now that I see your thermodynamic cycle I see you don't need a >>>>>> "solvent" calculation (or rather you have a "solvent" calculation where >>>>>> it's just the protein in solution with no ligand). >>>>>> >>>>>>> >>>>>>> >>>>>>>> Thanks. >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>>> When I did these same calculations using AMBER 11, I am getting >>>>>>>>> 6.945 kcal/mol which is still less but it does not have a negative >>>>>>>>> sign. >>>>>>>>> The other thing I observed is the coulomb and the vdw dV/dl vs >>>>>>>>> lambda curves for OPLS-AA (gromacs) and AMBER 99SB (AMBER11) have a >>>>>>>>> very >>>>>>>>> similar trend only shifted in the values of dV/dl in the y-axis. >>>>>>>>> >>>>>>>>> I can attach the graphs and include more details, if needed. >>>>>>>>> Let me know. >>>>>>>>> >>>>>>>>> >>>>>>>>> Thanks for your time, >>>>>>>>> >>>>>>>>> Regards >>>>>>>>> Sai >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> -- >>>>>>>> David Mobley >>>>>>>> dmob...@gmail.com >>>>>>>> 949-385-2436 >>>>>>>> >>>>>>>> >>>>>>> >>>>>> >>>>>> >>>>> >>>>> Thank you for your time. >>>>> >>>>> Regards >>>>> Sai >>>>> >>>>>> >>>>>> -- >>>>>> David Mobley >>>>>> dmob...@gmail.com >>>>>> 949-385-2436 >>>>>> >>>>>> >>>>> >>>> >>> >>> >>> -- >>> David Mobley >>> dmob...@gmail.com >>> 949-385-2436 >>> >>> >> > -- 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