[gmx-users] Re: Glutamate to Alanine Mutation
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.comwrote: 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.comwrote: 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
[gmx-users] Re: Glutamate to Alanine Mutation
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.comwrote: 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
[gmx-users] DelG of LJ Transformation Negative Glu---Ala Mutation
Hi Gmx Community, I am doing Glutamate to Alanine mutation in presence and absence of a ligand. The Coulomb transformation is yielding the following results: Glu---Ala + ligand = 790.109 kJ/mol Glu---Ala + no ligand = 787.33 kJ/mol During the LJ transformation Glu---Ala + ligand = - 24.87 kJ/mol Glu---Ala + no ligand = - 12.4403 The reason I am asking this is because in experiments, the Glu---Ala mutation yields a positive deldelG but in TI approach, I am getting a negative deldelG. The major contribution is from the LJ transformation. With the above values, the deldelG from TI approach is -9.65 kJ/mol Is it possible to have a negative deldelG for a mutation like Glutamate (which is charged) to Alanine (which is neutral)? If it is system dependent and you need more information than what provided in this email, I will add in happily. But I am just curious to know your opinions. Regards Sai -- gmx-users mailing listgmx-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
Re: [gmx-users] FEP
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 residatom cgnr chargemasstotal_charge 1CG331 1 SIM C1 1-0.27 12.0110 CG331 0.00 12.0110 2 HGA3 1 SIM H11 2 0.09 1.00800 HGA30.00 1.00800 3 HGA3 1 SIM H12 3 0.09 1.00800 HGA30.00 1.00800 4 HGA3 1 SIM H13 4 0.09 1.00800 HGA30.00 1.00800 5CG331 1 SIM C2 5-0.27 12.0110 CG331 -0.1812.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 HGA30.00 1.00800 9 DUM 1 SIM H31 9 0.00 1.00800 10 DUM 1 SIM H3210 0.00 1.00800 11 DUM 1 SIM H3311 0.00 1.00800 STEP 2. (mutation of groups. CH3-CH2-HD3 becomes D3H-CH2-CH3) [ atoms ] ; nr type resnr residatom cgnr chargemasstotal_charge 1CG331 1 SIM C1 1 0.00 12.0110 HGA30.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.001.00800 5CG331 1 SIM C2 5 -0.18 6 HGA31 SIM H21 6 0.09 7 HGA31 SIM H22 7 0.09 8 HGA31 SIM H23 8 0.00 1.00800 CG331 0.00 12.0110 9 DUM1 SIM H31 9 0.00 1.00800 HGA30.00 1.00800 10 DUM1 SIM H3210 0.00 1.00800 HGA30.00 1.00800 11 DUM1 SIM H3311 0.00 1.00800 HGA30.00 1.00800 STEP3. (opposite of step1) [ atoms ] ; nr type resnr residatom cgnr chargemasstotal_charge 1 HGA3 1 SIM C1 1 0.00 12.0110 HGA3 0.091.00800 2 DUM 1 SIMH11 2 0.00 1.00800 DUM0.00 1.00800 3 DUM 1 SIMH12 3 0.00 1.00800 DUM0.00 1.00800 4 DUM 1 SIMH13 4 0.00 1.00800 DUM0.00 1.00800 5CG331 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 8CG331 1 SIM H23 8 0.00 1.00800 CG331 -0.27 12.0110 9 HGA3 1 SIM H31 9 0.00 1.00800 HGA30.09 1.00800 10 HGA3 1 SIM H3210 0.00 1.00800 HGA30.09 1.00800 11 HGA3 1 SIM H3311 0.00 1.00800 HGA30.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,
[gmx-users] Re: Thermodynamic Integration Glu - Ala Mutation
Any suggestions on this topic? Thanks again for your time. Dear Gromacs Users, I am trying to find out the relative free energy difference of binding of a ligand with wild type protein (Glutamate residue) and mutant protein (Alanine residue). For charge part of the mutation, this is what I'm planning to do: ; residue 40 GLU rtp GLU q -1.0 552 opls_238 40GLU N 191 -0.514.0067 opls_238-0.5 14.0067 ; qtot 0.5 553 opls_241 40GLU H 1910.3 1.008 opls_241 0.3 1.008 ; qtot 0.8 554 opls_224B 40GLU CA191 0.14 12.011 opls_224B0.14 12.011 ; qtot 0.94 555 opls_140 40GLU HA191 0.06 1.008 opls_140 0.061.008 ; qtot 1 556 opls_136 40GLU CB192 -0.12 12.011 opls_135 -0.18 12.011 ; qtot 0.88 557 opls_140 40GLUHB1192 0.06 1.008 opls_140 0.061.008 ; qtot 0.94 558 opls_140 40GLUHB2192 0.06 1.008 opls_140 0.061.008 ; qtot 1 559 opls_274 40GLU CG193 -0.22 12.011 opls_140 0.061.008 ; qtot 0.78 560 opls_140 40GLUHG1193 0.06 1.008 DUM_1400.0 1.008 ; qtot 0.84 561 opls_140 40GLUHG2193 0.06 1.008 DUM_1400.0 1.008 ; qtot 0.9 562 opls_271 40GLU CD 1940.7 12.011 DUM_2710.0 12.011 ; qtot 1.6 563 opls_272 40GLUOE1194 -0.815.9994 DUM_272 0.0 15.9994 ; qtot 0.8 564 opls_272 40GLUOE2194 -0.815.9994 DUM_2720.0 15.9994 ; qtot 0 565 opls_235 40GLU C 1950.5 12.011 opls_235 0.5 12.011 ; qtot 0.5 566 opls_236 40GLU O 195 -0.515.9994 opls_236-0.515.9994 ; qtot 0 I added the DUM_140, DUM_271,DUM_272 atomtypes in ffnonbonded.itp. Further I added the bondtypes, angletypes and dihedraltypes in ffbonded.itp for state B. In order to maintain the electroneutrality, I am mutating a K+ ion to dummy as well. For the K+ going to dummy, I added a DUM_408 atomtype as well. *The questions I have are as follows:* *During the charge mutation will my dummy atoms have sigma and epsilon as 0.0? Since Ala residue does not have Cgamma, Cdelta and the oxygen as well as the Hgamma, all the dummy atoms should have sigma and epsilon as zero. Am I correct for this assumption?* *Further as I am mutating one ion, I dont want the mutating ion to come close to other ions (I have three K+ ions as my protein has -3 charge) as well as the protein atoms. Hence I'm having position restraints on the three ions. * *For the ion that is mutating to dummy, should I have have position restraints on dummy atom as well for the B state topology?* The section of topology for the K+ going to dummy is as follows: [ moleculetype ] ; Namenrexcl KM 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls-408 1 KMKM1 139.0983 DUM_408 0.0 39.0983 #ifdef POSRES_ION ; Position restraint for each Potassium ion [ position_restraints ] ; i funct fcxfcyfcz 11 1000 1000 1000 0 0 0 #endif Here I defined the KM as K+ that will go to dummy and other two K+ ions are represented as regular K+ ions. I'm pasting the regular K+ ion part of the topology below. [ moleculetype ] ; Namenrexcl K 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_408 1 K K1 139.0983 #ifdef POSRES_ION ; Position restraint for each Potassium ion [ position_restraints ] ; i funct fcxfcyfcz 11 1000 1000 1000 #endif Let me know if you need more information. Thanks for your time, Regards Sai -- gmx-users mailing listgmx-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] Thermodynamic Integration Glu - Ala Mutation
Dear Gromacs Users, I am trying to find out the relative free energy difference of binding of a ligand with wild type protein (Glutamate residue) and mutant protein (Alanine residue). For charge part of the mutation, this is what I'm planning to do: ; residue 40 GLU rtp GLU q -1.0 552 opls_238 40GLU N 191 -0.514.0067 opls_238-0.5 14.0067 ; qtot 0.5 553 opls_241 40GLU H 1910.3 1.008 opls_241 0.3 1.008 ; qtot 0.8 554 opls_224B 40GLU CA191 0.14 12.011 opls_224B0.14 12.011 ; qtot 0.94 555 opls_140 40GLU HA191 0.06 1.008 opls_140 0.061.008 ; qtot 1 556 opls_136 40GLU CB192 -0.12 12.011 opls_135 -0.18 12.011 ; qtot 0.88 557 opls_140 40GLUHB1192 0.06 1.008 opls_140 0.061.008 ; qtot 0.94 558 opls_140 40GLUHB2192 0.06 1.008 opls_140 0.061.008 ; qtot 1 559 opls_274 40GLU CG193 -0.22 12.011 opls_140 0.061.008 ; qtot 0.78 560 opls_140 40GLUHG1193 0.06 1.008 DUM_1400.0 1.008 ; qtot 0.84 561 opls_140 40GLUHG2193 0.06 1.008 DUM_1400.0 1.008 ; qtot 0.9 562 opls_271 40GLU CD 1940.7 12.011 DUM_2710.0 12.011 ; qtot 1.6 563 opls_272 40GLUOE1194 -0.815.9994 DUM_272 0.0 15.9994 ; qtot 0.8 564 opls_272 40GLUOE2194 -0.815.9994 DUM_2720.0 15.9994 ; qtot 0 565 opls_235 40GLU C 1950.5 12.011 opls_235 0.5 12.011 ; qtot 0.5 566 opls_236 40GLU O 195 -0.515.9994 opls_236-0.515.9994 ; qtot 0 I added the DUM_140, DUM_271,DUM_272 atomtypes in ffnonbonded.itp. Further I added the bondtypes, angletypes and dihedraltypes in ffbonded.itp for state B. In order to maintain the electroneutrality, I am mutating a K+ ion to dummy as well. For the K+ going to dummy, I added a DUM_408 atomtype as well. *The questions I have are as follows:* *During the charge mutation will my dummy atoms have sigma and epsilon as 0.0? Since Ala residue does not have Cgamma, Cdelta and the oxygen as well as the Hgamma, all the dummy atoms should have sigma and epsilon as zero. Am I correct for this assumption?* *Further as I am mutating one ion, I dont want the mutating ion to come close to other ions (I have three K+ ions as my protein has -3 charge) as well as the protein atoms. Hence I'm having position restraints on the three ions. * *For the ion that is mutating to dummy, should I have have position restraints on dummy atom as well for the B state topology?* The section of topology for the K+ going to dummy is as follows: [ moleculetype ] ; Namenrexcl KM 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls-408 1 KMKM1 139.0983 DUM_408 0.0 39.0983 #ifdef POSRES_ION ; Position restraint for each Potassium ion [ position_restraints ] ; i funct fcxfcyfcz 11 1000 1000 1000 0 0 0 #endif Here I defined the KM as K+ that will go to dummy and other two K+ ions are represented as regular K+ ions. I'm pasting the regular K+ ion part of the topology below. [ moleculetype ] ; Namenrexcl K 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeB chargeB massB 1 opls_408 1 K K1 139.0983 #ifdef POSRES_ION ; Position restraint for each Potassium ion [ position_restraints ] ; i funct fcxfcyfcz 11 1000 1000 1000 #endif Let me know if you need more information. Thanks for your time, Regards Sai -- gmx-users mailing listgmx-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] Glutamate Chi1 Dihedral Not Symmetrical?
Dear Gromacs Users, I am trying to plot the Ryckaert-Bellemans energy for rotating the chi1 dihedral of glutamate in protein. I tried to change the ch1 dihedral from 0-360 degrees at increments of 1 degree and saved the pdbs. Then I used the pdb's to obtain the corresponding gro files and a single topology file. I used opls force field and no water and a 0 step minimization. *My em.mdp file:* * * cpp = /usr/bin/cpp define = -DFLEXIBLE integrator = steep nsteps = 0 constraints = none emtol= 1000.0 emstep = 0.01 nstcomm = 1 ns_type = grid nstlist = 1 rlist= 0 rcoulomb = 0 rvdw = 0 Tcoupl = no Pcoupl = no gen_vel = no nstxout = 1 pbc = no ;energy groups energygrps = Protein After the minimization, when I plot the dihedral energy vs dihedral angle, I do not get a symmetrical curve. After careful observation of the gro files that I obtained using pdb2gmx, the 0-180 values are not exactly same as -180-0 degrees or 180-360 degrees. Is it because of the slightly different values of the dihedral 0-180 and 180-360 that I'm not getting a symmetrical curve or am I doing something wrong? I have attached the dihedral energy vs dihedral angle curve. Regards Sai attachment: chi1_gromacs.png-- gmx-users mailing listgmx-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] Free Energy calcualtions
Hi Gromacs Users, I want to mutate a glutamate in my protein to alanine in presence of a ligand. With glutamate, the protein charge is -3. To neutralize the system, I added 3K+ ions. Now when I mutate GLU to ALA, the charge in state_B will be +1 (protein -2 + 3K+). Right now I'm in the charge part of the mutation. Once this is successful, I will include the VDW mutation too. I have added the mutation details of Glu to Ala residue in the topology in [atoms], [bonds], [angles] and [dihedrals] sections. After I run the grompp command, the result says that my State_B topology has +1 charge since I am not including mutation of one K+ ion to neutral K+ ion. How can I mutate 1 particular K+ to K atom and subsequently to a dummy atom? Since I'm using OPLS force field, we have ions.itp to be used in the topology file, changing one K+ is making it troublesome for me to implement in the topology file. Any suggestions are helpful. Thanks for your time. Regards Sai -- gmx-users mailing listgmx-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
Re: [gmx-users] Half double pair list method in GROMACS [update]
Hi Stephane Chris, I followed all the threads posted by you two. I have a protein using ff99SB and GLYCAM for sugars. I have a disaccharide bound to protein. In xleap of AMBERTOOLS, I use the GLYCAM and ff99SB to generate the topology and coordinate files. I did the tests as Chris' original posts and by Stephane. The 1-4 interaction terms match for sugar alone and protein alone with HEDP method. When I generate topology and coordinate files with xleap of AMBER, there are three atom types that are common to protein and sugar. The atom types for my case are H1, OH and HO. Since for protein pairtypes using ff99SB, the epsilon has to be divided by 10 and pairs section be replicated five times and for sugar the epsilon in the pairtypes be divided by six and replicated six times, I am a little concerned about the three atomtypes that are common. So what I did was to change the atomtypes of sugar as H1S, OHS HOS for sugar and H1, OH and HO for protein and I made the changes accordingly in the pairtypes section for protein and sugar. Is this a valid approach? Any suggestion will be helpful. To make things little more clear: H1H10. 0. A 2.47135e-01 6.56888e-02 ;originally obtained using amb2gmx.pl H1S H1S 0. 0. A 2.47135e-01 6.56888e-02 ;Glycam Hydrogen of Sugar ( I changed this so that the common atom types be separated) In the ffnonbonded_complex_mod.itp: ;;using combination rule of 2 [ pairtypes ] ;;for protein H1 H1 1 0.247135000 0.006568880 ;the epsilon is divided by 10 ;;for sugar H1S H1S 1 0.24713500 0.010948133 ;the epsilon is divided by 6 Thanks for your time, Regards Sai On Mon, Sep 5, 2011 at 11:33 AM, ABEL Stephane 175950 stephane.a...@cea.frwrote: Dear All, Below a little update and results about the application of half double pair list method to scale properly the Coulombic 1-4 interactions in case of a system where the AMBER99SB (fudgeLJ=0.5 and fudgeLJ=0.8333) and GLYCAM06 (fudgeLJ=1.0 and fudgeLJ=1.0) force fields are combined. I have followed the 4 steps described in [1] and used the following values in my forcefield.itp file [ defaults ] ; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ 1 2 yes 1.0 0.16 #include ffnonbonded_mod.itp ;#include ffnonbonded.itp #include ffbonded.itp I used two different topology files for the glycolipid (bDM) and the peptide. One with (*_mod.itp) with pair list parameters duplicated 6 times (bDM) and 5 times (peptide) and with single pair list (*_no_mod.itp) as decribed in [1]. TESTING: Three 3 different systems were examined: A. A first system containing 1 glycolipid (bDM) in water cubic box B. A second system with 1 peptide in TIP3P water C. And a third system with 1 peptide and 1 glycolipid in water cubic box To obtain the glycolipid and peptide energy pairs, I did one step of MD in NVT ensemble with the *.mdp file given in [2] with different energygrps and tc_grps. For 1. energygrps and tc_grps = bDM SOL For 2. energygrps and tc_grps = Protein SOL For 3. energygrps and = Protein bDM SOL bDM/water system Test_A1 ## Control with GLYCAM force field fudgeLJ fudgeQQ parameters and the *_no_mod.itp file : ##; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ ## 1 2 yes 1.0 1.0 Epot (kJ/mol)Coul-SR LJ-SRCoul-14 LJ-14 bDM-bDM -4.00855e+02 -3.71401e+012.03406e+032.79234e+02 Test_A2 with the topology *_mod.itp file and the directive ##; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ ## 1 2 yes 1.0 0.16 Epot (kJ/mol)Coul-SR LJ-SRCoul-14 LJ-14 bDM-bDM -4.00855e+02 -3.71401e+012.03406e+032.79234e+02 Test_A3 with the topology *_no_mod.itp file and the directive ##; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ ## 1 2 yes 1.0 0.16 Epot (kJ/mol)Coul-SR LJ-SRCoul-14 LJ-14 bDM-bDM -4.00855e+02 -3.71401e+013.39010e+022.79234e+02 Coul-14 energy for bDM-bDM in Test_A1 = Coul-14 energy in Test_A2 -- OK ! Coul-14 energy for bDM-bDM Test_A3 is 6 smaller than Coul-14 energy in Test_A1 and Test 2 --- OK ! peptide/water system Test_B1 Control with AMBER force field fudgeLJ fudgeQQ parameters, the *_no_mod.itp file ##; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ ## 1 2 yes 0.5 0.8333 Epot (kJ/mol)Coul-SR LJ-SRCoul-14 LJ-14 Protein-Protein -1.49026e+03 -4.12114e+023.80551e+036.55321e+02 Test_B2 Control with the peptide *_no_mod.itp file and the directive ##; nbfunc
Re: [gmx-users] pulling force vs free energy
I agree with Justin. I have tried myself several SMD simulations for ligand binding studies. I tried 500 simulations and not sure if they are enough. Further, the path dependence is very important part. For different paths that you can try, look at McCammon group's paper in JACS 128, 3019-3026. There are several other papers too, but this discusses about various methods that you can try to get the pmf using Jarzynski On Wed, Nov 16, 2011 at 9:25 AM, Justin A. Lemkul jalem...@vt.edu wrote: Vijayaraj wrote: Message: 4 Date: Wed, 16 Nov 2011 09:34:02 -0500 From: Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu Subject: Re: [gmx-users] pulling force vs free energy To: Discussion list for GROMACS users gmx-users@gromacs.org mailto:gmx-users@gromacs.org** Message-ID: 4ec3c9da.7040...@vt.edu mailto:4EC3C9DA.7040400@vt.**edu4ec3c9da.7040...@vt.edu Content-Type: text/plain; charset=ISO-8859-1; format=flowed Vijayaraj wrote: Hi, What is the relation between pulling force and free energy of binding. can we relate the maximum pulling force with the free energy. for example, 2 systems has the maximum pulling force and free energy as below from umbrella sampling and g_wham analysis, max. forcefree energy system 1 147042 system 2 164732 system 2 has higher pulling force than system 1 and the free energy result is different from this trend. How did you obtain the maximum force, just a single SMD trajectory? If so, I wouldn't put a lot of faith in it necessarily. Umbrella sampling is a more robust method than a single pull. You can use large numbers of pulling simulations and apply Jarzynski's equation to calculate free energy, but there are distinct caveats (although I suppose there are caveats with any method). -Justin yes. the max force is obtained from single SMD trajectory. So in this case we dont have to worry about the correction between max. force and free energy. I found one of my free energy result is 2 times larger than the previous result, where they have applied Jarzynski's equation. SMD is path-dependent, while a true DeltaG is a path-independent quantity. Hence why you cannot easily connect the two. Convergence in sampling and limitations in each technique make it sometimes hard to compare the results that others have obtained with other methods. Proper data collection for Jarzynski's method requires exhaustive sampling, which is often hard to obtain (not a blind criticism of others' work, just a fact). -Justin -- ==**== Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justinhttp://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ==**== -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/**mailman/listinfo/gmx-usershttp://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/** Support/Mailing_Lists/Searchhttp://www.gromacs.org/Support/Mailing_Lists/Searchbefore 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_Listshttp://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-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
Re: [gmx-users] How to calculate work
Hi Sanku, I was using the pullf.xvg and multiplying it with pulling rate and time. f*v*dt = W and getting the total work for each SMD simulation. I'm not sure if this is the best/correct way to do it. But from original Jarzynski's article (PRL (2007) 78(14), 2690-2693) this is what I deduced. I asked this question but I think very few users of Gromacs do Jarzynski Equality (JE) to get free energy differences. If anyone can comment on this, it will be useful for the community. PS: If you use AMBER, it gives the work profile as the output. You can use the total work which is printed at the last line of the output and get the exponential average of beta*W. PPS: May be you are aware of this issue, but the JE suffers from the fact that you do the exponential average and the smaller work values determine everything. Regards Sai On Fri, Nov 4, 2011 at 7:46 PM, Sanku M msank...@yahoo.com wrote: Hi, I am performing steered MD simulation using gromacs. I was wondering how one can get the time profile of the irreversible work from the gromacs pull-code out put . From constant pulling-rate SMD, we get time profile of force pullf.xvg and pullx.xvg. I wonder does multiplying the the value from pullx.xvg and value from pullf.xvg will provide the work . Or, will it be force ( obtained from pullf.xvg) multiplied by pulling rate multiplied by time ? Sanku -- gmx-users mailing listgmx-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 listgmx-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] Jarzynski and PMF
Dear All, I have tried to calculate the free energy of binding using Jarzynski equality. I employed the following procedure. --I did force pulling simulations along Z-direction as exemplified in Justin's umbrella sampling simulations. I did 50, 250 and 500 pulling simulations to test the convergence and also two different pulling rates (2nm/ns and 4nm/ns), two different force constants (2000 kJ mol^-1 nm^-2 and 4000 kJ mol^-1 nm^-2). --From the pulling simulation, I get the pullf.xvg and pullx.xvg. --I use the force from pullf.xvg and calculate the work at each time step by multiplying with v*dt and calculate average work for each simulation. --Once I have the work for each simulation, I calculate the exp(-beta*W) for each simulation and calculate the average of exp(-beta*W). --Once I have the average of exp(-beta*W), I calculate the free energy delF = -kB*T ln(exp(-beta*W)). The question I have is whether my approach for calculating the work is correct/feasible or not. Does the pullf.xvg written by mdrun during a pulling simulation contain the pull force only? Do I need to add any corrections? I'm doing NPT simulations, so can I add the pdV correction as suggested in one of the paper (Macromolecules, 2008, 41 (6), pp 2283–2289) by Berk Hess? Or did I understand the paper incorrectly? For a publication level of work, I'll use NVT ensemble. But for now my biggest concern is calculation of work. Any suggestions will be greatly helpful. Thanks for you time, Sai Ramadugu Dept of Chemistry, University of Iowa. -- gmx-users mailing listgmx-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] PMF for Carbohydrate-Protein
Dear All, I'm running set of umbrella sampling simulations to get the PMF of a disaccharide binding to a protein. I followed the tutorial provided by Justin and changed the values of necessary parameters according to my system. The mdp file for umbrella sampling simulations is as follows: --- ; Pull code pull= umbrella pull_geometry = distance ; can't get PMF with direction pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = active_site pull_pbcatom0 = 575 pull_group1 = new pull_pbcatom1 = 1500 pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 750 ; kJ mol^-1 nm^-2 pull_nstxout= 1000 ; every 2 ps pull_nstfout = 1000 ; every 2 ps Active site in the pull_group0 being the residues in the binding site and new in the pull_group1 being the disaccharide. Initially I did an SMD run where I separated the disaccharide from the binding upto 1.5nm in distance. When I started running the umbrella sampling simulations, I used the conformations with approximately 0.1nm difference starting from an initial separation of 0.3nm between the binding site residues and disaccharide. My question is: when I run g_wham and get the profile, the pmf reaches a plateau, but its not completely horizontal. It drops down a little. Is it acceptable to have the pmf reach a plateau but not completely horizontal? In the histogram, the first two windows have good overlap and are narrow distributions where as the later windows are spread out. I have tried to decrease the distance from 0.1 nm to 0.05 nm and still the sampling around 0.5-0.75 of Xi region does not sample well. Is there anything that you could suggest? I have uploaded the plots of profile and histograms. PMF Plot link: http://img801.imageshack.us/i/pmf.png/ Histogram link: http://img217.imageshack.us/i/histoy.png/ Thanks for your time Sai Ramadugu -- gmx-users mailing listgmx-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
Re: [gmx-users] soft-core
Hi Well the tutorial in the following website has mentioned about the importance of sc_alpha and sc-power. Also the manual gives you more information. http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial. Regards Sai On Sun, Jun 13, 2010 at 6:56 PM, Justin A. Lemkul jalem...@vt.edu wrote: fancy2012 wrote: Dear GMX users, when I generate the input file of MD simulation using grompp, I get this error: ERROR: The soft-core power is 0 and can only be 1 or 2. I don't know how to figure it out? Could somebody give me a hand? Thanks very much in advance! You haven't specified a value for sc-power, so the default of zero is taken. Manual section 7.3.23. -Justin Here is the mdp file: cpp = /lib/cpp constraints = all-bonds integrator = md tinit = 0.0 dt = 0.002 nsteps = 25000 nstcomm = 1 nstxout = 1 nstvout = 0 nstfout = 0 nstlog = 500 nstenergy = 250 nstxtcout = 1000 xtc_precision = 1000 xtc_grps= Protein nstlist = 5 energygrps = Protein SOL ns_type = grid rlist = 0.9 coulombtype = Reaction-Field epsilon_r= 78.0 rcoulomb= 1.4 rvdw= 1.4 Tcoupl = Berendsen tc-grps = Protein SOL ref_t = 300 300 tau_t = 0.1 0.1 Pcoupl = Berendsen tau_p = 1.0 compressibility = 4.6e-5 ref_p;= 1.0 free_energy = yes init_lambda = 0.00 sc-alpha= 1.51 All the best, fancy -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] tutorial of free energy calculations
HI Its not working for me either. But for free energy tutorials using GROMACS, I think the one written by Prof Mobley is like Bible. The link is as follows: http://www.dillgroup.ucsf.edu/group/wiki/index.php/Free_Energy:_Tutorial But the one written by Prof Alan Mark (who is now in Queensland Univ, Australia) is also very useful. Try searching for it. Regards Sai 2010/6/7 fancy2012 fancy2...@yeah.net Dear GMX users, I don't know why I can't open the tutorial of free energy calculations on this website http://www.gromacs.org/Documentation/Tutorials. Can you open it? Thanks very much! All the best, fancy -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] Thr - Met in TI calculation
Hi Fancy, The following link shows the change from p-cresol to Toluene. So follow the steps in the same by converting threonine to methionine. http://compbio.chemistry.uq.edu.au/education/Free-Energy_Course/0.introduction.html Regards Sai 2010/6/5 fancy2012 fancy2...@yeah.net Dear GMX users, I want to calculate the relative binding free energy between a small molecule binding to a protein of Wild Type and Thr - Met using TI, so how should I prepare the topology files of the protein using for TI calculation? Any suggestions will be highly appreciated? All the best, fancy -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: [gmx-users] PBC
Hi Morteza, I had this problem when I was running a trimeric protein attached to an oligosaccharide. *I have used the following command:* * * *trjconv -s .tpr -f .xtc -o -boxcenter tric -pbc mol* * * *but it is not working.* Do the following. It worked for me. In the first round, run trjconv -f *.xtc -s *.tpr -center -boxcenter tric -pbc mol -o 1.xtc and choose backbone (which is usually 4 from the index) for centring. In the second round, use 1.xtc obtained from the first step and run the command as follows trjconv -f 1.xtc -s *.tpr -center -boxcenter tric -pbc a -o 2.xtc and this time also choose the backbone for centring. This should bring back the dimeric protein system in the box. Sometimes the water molecules around the edge are broken. Hope this is not a big worry. By the way is your box cubic or truncated octahedron? My case was a cubic box. Hope this helps. Regards Sai On Tue, Jun 8, 2010 at 2:23 AM, shahab shariati shahab.shari...@gmail.comwrote: Morteza Khabiri wrote: Dear users I have a dimer protein in the water box. It was run for 30ns. during the simulation dimer split to two monomer. This things happen bc of PBC. ( I checked it by vmd pbc option ) to have a two monomer together during trajectories (for visualization) I have used the following command: trjconv -s .tpr -f .xtc -o -boxcenter tric -pbc mol but it is not working. Is there any other method or command which I could implement pbc in trajectory. Thanks in advance Morteza Shahab Shariati wrote: You can use other flags of trjconv command as follows: Trjconv –f *.xtc –s **.tpr –o ***.xtc –pbc nojump –ur compact -center -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Wrapping molecules for octahedron box
Dear Gromacs Users, I have a simulation of a protein in Truncated Octahedron box. I want to calculate the water residence times and mean square displacements (of water) around the active site residues. We developed our own algorithm for doing the same. In earlier simulations I had cubic box and anint() function was working for the periodic boundary conditions. Now that I have truncated octahedron box, is there any function to my rescue? I searched the gromacs forum and could not find any function. Is there any way to circumvent this problem? Sorry if this is a repost. Thanks for your time. Regards Ying-Hua -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Large dVpot/dlamda values at lambda = 1
Hi All, I'm trying to run free energy simulations for a linkage change in an oligosaccharide when bound to its protein. The linkage changes from beta 1-4 to beta 1-3. I followed the tutorial by Prof Mobley. For the lambda values I took 0 to 1 at intervals of 0.05. The values of dVpot/dlambda from 0.00 to 0.95 seem to be in the range of E3-4. But when I check the values of dVpot/dlambda for lambda = 1, it goes from E2 to E24 at the end. This is only for the LBFGS minimization. Does any one suspect a mistake in my topology or my approach? Thanks for your time. Sample .log file of lambda = 1 is pasted below. Step Time Lambda 100 100.01.0 Energies (kJ/mol) Bond AngleProper Dih. Ryckaert-Bell. LJ-14 1.06593e+048.20816e+032.93142e+024.26304e+038.30759e+03 Coulomb-14LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. 4.13439e+041.92653e+05 -1.08257e+04 -1.46637e+06 -2.72128e+05 Potential Pressure (bar) dVpot/dlambda dEkin/dlambda dG/dl constr. -1.48359e+06 -4.15979e+03 -1.01327e+150.0e+000.0e+00 Step Time Lambda 101 101.01.0 Energies (kJ/mol) Bond AngleProper Dih. Ryckaert-Bell. LJ-14 1.06477e+048.21153e+032.93330e+024.26328e+038.31182e+03 Coulomb-14LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. 4.13437e+041.92648e+05 -1.08257e+04 -1.49776e+06 -2.72128e+05 Potential Pressure (bar) dVpot/dlambda dEkin/dlambda dG/dl constr. -1.51500e+06 -4.33201e+03 -2.78144e+180.0e+000.0e+00 Step Time Lambda 103 103.01.0 Energies (kJ/mol) Bond AngleProper Dih. Ryckaert-Bell. LJ-14 1.06776e+048.21126e+032.93576e+024.26328e+038.31168e+03 Coulomb-14LJ (SR) Disper. corr. Coulomb (SR) Coul. recip. 4.13437e+041.92656e+05 -1.08257e+04 -1.76397e+06 -2.72128e+05 Potential Pressure (bar) dVpot/dlambda dEkin/dlambda dG/dl constr. -1.78117e+06 -5.80136e+03 -1.76071e+240.0e+000.0e+00 Regards Sai -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Free Energy Calculations for linkage change
Dear All, I have an disaccharide bound to a C-type lectin. The sequence of disaccharide is beta-galactose1-4beta-GlcNAc. The 3'OH and 4'OH of beta-galactose coordinate to the Ca2+ at the active site. FRET experiments have shown that if the disaccharide linkage is changed to beta-galactose1-3beta-GlcNAc then the binding affinity drops to 20 fold. I have attached both disaccharides in the figure. The left side residue is galactose and right side residue is N-Acetylgalactoseamine. . In order to do a free energy perturbation and find the relative free energy of binding for these two ligands, I have to slowly change the linkage from beta1-4 to beta1-3. If you look at the models of both disaccharides from the figure, to convert the beta 1-4 linked disaccharide to beta 1-3 linked disaccharide, I have to make the the C1 atom of Galactose attached to O4 of N-Acetylgalactoseamine as H atom and rest of the galactose residue to dummy atoms. On the other hand I have to grow the residue from H attached to O3 of N-acetylgalactoseamine to a galactose. The number of atoms involved is 21 atoms. Further the constraints I have is that the sugar should be bound in right conformation at the active site. So may be I need distance restraints. Does my approach make any sense? I have seen some papers by others where they add a fragment or delete a fragment between two ligands. But in my case the difference is the linkage between two residues. Thanks for your time in advance. Regards Sai Ramadugu attachment: disaccharide.gif-- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] Free Energy Calculations for linkage change in disaccharide
Dear All, In my previous email the text in the message was not shown. Thats why I'm writing this again. I request to see the attachment of earlier message. I'm sorry for this!! I have an disaccharide bound to a C-type lectin. The sequence of disaccharide is beta-galactose1-4beta-GlcNAc. The 3'OH and 4'OH of beta-galactose coordinate to the Ca2+ at the active site. FRET experiments have shown that if the disaccharide linkage is changed to beta-galactose1-3beta-GlcNAc then the binding affinity drops to 20 fold. I have attached both disaccharides in the figure. The left side residue is galactose and right side residue is N-Acetylgalactoseamine. . In order to do a free energy perturbation and find the relative free energy of binding for these two ligands, I have to slowly change the linkage from beta1-4 to beta1-3. If you look at the models of both disaccharides from the figure, to convert the beta 1-4 linked disaccharide to beta 1-3 linked disaccharide, I have to make the the C1 atom of Galactose attached to O4 of N-Acetylgalactoseamine as H atom and rest of the galactose residue to dummy atoms. On the other hand I have to grow the residue from H attached to O3 of N-acetylgalactoseamine to a galactose. The number of atoms involved is 21 atoms. Further the constraints I have is that the sugar should be bound in right conformation at the active site. So may be I need distance restraints. Does my approach make any sense? I have seen some papers by others where they add a fragment or delete a fragment between two ligands. But in my case the difference is the linkage between two residues. Thanks for your time in advance. Regards Sai -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
[gmx-users] free energy calculations for two ligands
Dear All, I have a trimeric protein system which binds to two different ligands. The ligands are oligosaccharides. The difference between the two is first ligand has beta 1-4 linkage and the second ligand has beta 1-3 linkage between Galactose (Gal) and N-Acetyl Glucosamine (GlcNAc). Since the linkage changes I have changes in atoms, bonds, pairs, angles and dihedrals. In the pair and dihedrals, some of the pairs are present in the first ligand and not in the second ligand. How do I make the topology for the change from ligand A to ligand B in the free energy simulations? The same is the case with dihedrals. Angles, bonds and atoms have no problem. Can someone please help me in this regard? Thanks for your time. Regards Sai -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
