And oh, none of the Kenno's web pages for tutorial and software is working at the moment.
On Tue, Jan 17, 2017 at 8:37 PM, Mohsen Ramezanpour < ramezanpour.moh...@gmail.com> wrote: > Hi Gromacs users, > > I want to follow recommended procedure (by Justin) for MM scanning of > rotatable dihedrals and do it manually by Gromacs. > > I have made the following em.mdp (for restraining desired rotatable bond > and relaxing the rest of molecule) and md-zero.mdp for the single-point > energy calculation. > I was wondering if these parameters seems reasonable. > > > *em.mdp:* > > define = -DPOSRES > integrator = steep > emtol = 1000.0 > nsteps = 5000 > vdwtype = Cut-off > vdw-modifier = Force-switch > coulombtype = pme > ; > constraints = h-bonds > constraint_algorithm = LINCS > ; to do a MM scan > nstlist = 0 > ns-type = simple > pbc = no > cutoff-scheme = group > rcoulomb = 0 > rvdw = 0 > rvdw_switch = 0 > rlist = 0 > > > > > *md-zero.mdp:* > > define = > integrator = md > dt = 0.001 > nsteps = 0 > nstlog = 0 > nstxout = 0 > nstvout = 0 > nstfout = 0 > nstcalcenergy = 0 > nstenergy = 0 > ; > nstxout-compressed = 100000 > compressed-x-precision = 1000 > coulombtype = pme > vdwtype = Cut-off > vdw-modifier = Force-switch > constraints = h-bonds > constraint_algorithm = LINCS > nstcomm = 100 > comm_mode = linear > comm_grps = system > refcoord_scaling = com > ; > nstlist = 0 > ns-type = simple > pbc = no > cutoff-scheme = group > rcoulomb = 0 > rvdw = 0 > rvdw_switch = 0 > rlist = 0 > > *I will use this command to calculate it:* > > mdrun -s input.tpr -rerun configuration.pdb -nsteps 0 > > Thanks in advance for your comments, > Mohsen > > On Fri, Jan 13, 2017 at 2:04 PM, Justin Lemkul <jalem...@vt.edu> wrote: > >> >> >> On 1/13/17 3:26 PM, Mohsen Ramezanpour wrote: >> >>> Thanks. interesting! >>> So, for now, lets focus on one of those dihedrals which has another >>> equivalent. >>> If I do PES on a rotatable bond in one ethyl, the other equivalent >>> rotatable bond in other ethyl will be free to move. If I understood >>> correctly, you meant it does not matter and I can do PES with one of the >>> two. >>> >>> Looking at QM and MM profiles by GAAMP, I figured out something >>> interesting: >>> If I shift the MM graph for 120 in phi direction, MM and QM will be >>> similar >>> in their general behaviour and it is much much better that before. >>> It works for shift of 240 too, but 120 results in better agreement with >>> QM >>> PES. >>> >> >> This sounds like the phase and/or multiplicity are off. Perhaps a bad >> initial guess of the parameters threw off the fitting. >> >> So, there might be two possible reasons for this: >>> 1) parameters are not good and this is how it is >>> 2) GAAMP has started rotation from different starting points (with a >>> difference of 120) >>> 3) GAAMP made a mistake for taking one dihedral for QM and took the >>> equivalent one for MM PES. (i.e. taking one of oxygen atoms in QM and >>> other >>> Oxygen in MM for the same rotatable bond) >>> >>> I should read the deatils of how GAAMP did these exactly to rule out the >>> options 2 and 3 >>> >> >> That's the beauty of GAAMP - it gives you full input and output for all >> the QM and MM processes. It should be quite apparent what the programs >> did. Perhaps it did not deal elegantly with the symmetry of the molecule. >> Its programs are quite robust but no black box is perfect (automated >> parametrization is a long-standing goal but in truth there's a lot of work >> left to be done). >> >> -Justin >> >> >> Cheers >>> >>> >>> >>> On Fri, Jan 13, 2017 at 5:45 AM, Justin Lemkul <jalem...@vt.edu> wrote: >>> >>> >>>> >>>> On 1/12/17 11:20 PM, Mohsen Ramezanpour wrote: >>>> >>>> On Thu, Jan 12, 2017 at 8:36 PM, Justin Lemkul <jalem...@vt.edu> wrote: >>>>> >>>>> >>>>> >>>>>> On 1/12/17 10:21 PM, Mohsen Ramezanpour wrote: >>>>>> >>>>>> On Thu, Jan 12, 2017 at 6:31 PM, Justin Lemkul <jalem...@vt.edu> >>>>>> wrote: >>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> On 1/12/17 7:47 PM, Mohsen Ramezanpour wrote: >>>>>>>> >>>>>>>> On Thu, Jan 12, 2017 at 4:24 PM, Justin Lemkul <jalem...@vt.edu> >>>>>>>> wrote: >>>>>>>> >>>>>>>> >>>>>>>>> >>>>>>>>> >>>>>>>>> On 1/12/17 3:18 PM, Mohsen Ramezanpour wrote: >>>>>>>>> >>>>>>>>>> >>>>>>>>>> Dear Gromacs users, >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> For parameterization of a molecule in Charmm36, I have got the QM >>>>>>>>>>> scanning >>>>>>>>>>> and partial charges from GAMMP server. However, the fitted >>>>>>>>>>> parameters >>>>>>>>>>> are >>>>>>>>>>> not good enough. >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> That's very surprising. What's wrong with what GAAMP gave you? >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>> The dihedral has two local minima in both QM and fitted ones both >>>>>>>>>> from >>>>>>>>>> >>>>>>>>>> GAAMP. The angle for the minima are okay but the corresponding >>>>>>>>>> depths >>>>>>>>>> >>>>>>>>> are >>>>>>>>> not. >>>>>>>>> In fact, the depth for the first local minimum is larger than >>>>>>>>> second >>>>>>>>> one >>>>>>>>> in >>>>>>>>> QM, while the situation is reverse in MM profile with fitted >>>>>>>>> parameters. >>>>>>>>> This makes the dihedral to be more (statistically) in wrong angle >>>>>>>>> (in >>>>>>>>> local >>>>>>>>> minimum which is not the most favourable one). >>>>>>>>> GAAMP, unfortunately, did not work well with my case (some critical >>>>>>>>> partial >>>>>>>>> charges and critical dihedrals). >>>>>>>>> >>>>>>>>> >>>>>>>>> I decided to do the MM scanning and try to get better parameters >>>>>>>>> for >>>>>>>>> the >>>>>>>>> >>>>>>>>> >>>>>>>>> dihedral. >>>>>>>>>> >>>>>>>>>> >>>>>>>>>>> Unfortunately, I do not have any experience with this part, and I >>>>>>>>>>> could >>>>>>>>>>> not >>>>>>>>>>> find any tutorial for how to do this in Gromacs. >>>>>>>>>>> I was wondering if you are aware of any tutorial which could >>>>>>>>>>> help me >>>>>>>>>>> to >>>>>>>>>>> overcome this challenging step. >>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> Tutorial (CGenFF theory is the same as CHARMM, by design): >>>>>>>>>>> >>>>>>>>>>> http://mackerell.umaryland.edu/~kenno/cgenff/download.php#tutor >>>>>>>>>>> >>>>>>>>>> >>>>>>>>>> Fitting program and other resources: >>>>>>>>>> http://mackerell.umaryland.edu/~kenno/lsfitpar/ >>>>>>>>>> http://mackerell.umaryland.edu/ff_dev.shtml >>>>>>>>>> >>>>>>>>>> Obviously, these are all CHARMM-centric approaches and frankly the >>>>>>>>>> modules >>>>>>>>>> within CHARMM make parametrization rather straightforward (not >>>>>>>>>> "easy," >>>>>>>>>> mind >>>>>>>>>> you, but straightforward). Since I began working with CHARMM, it >>>>>>>>>> has >>>>>>>>>> become indispensable in my daily routine. >>>>>>>>>> >>>>>>>>>> If you want to do things in GROMACS, the main issue is that you >>>>>>>>>> will >>>>>>>>>> have >>>>>>>>>> to do MM scans in a more manual fashion, by restraining the target >>>>>>>>>> dihedrals (very strongly) in a series of configurations >>>>>>>>>> (typically at >>>>>>>>>> intervals of 15 degrees over a full rotation) while allowing the >>>>>>>>>> rest >>>>>>>>>> of >>>>>>>>>> the molecule to relax to match the QM. >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> If I got it right, I must do EM on the molecule while only the >>>>>>>>>> >>>>>>>>> desired >>>>>>>>> dihedral is fixed in a specific angle. Which aspect should match >>>>>>>>> with >>>>>>>>> QM? >>>>>>>>> If you mean structurally, what is the criteria for matching >>>>>>>>> (RMSD?.)? >>>>>>>>> >>>>>>>>> >>>>>>>>> "Match" in this case means "treat equivalently," therefore only one >>>>>>>>> >>>>>>>>> constraint (restraint in the MM) should be applied while allowing >>>>>>>> the >>>>>>>> rest >>>>>>>> of the molecule to relax freely. You do have a difficult case >>>>>>>> because >>>>>>>> each >>>>>>>> of your molecules is symmetric; this means the same dihedral term is >>>>>>>> affecting multiple torsions. >>>>>>>> >>>>>>>> So, probably I should 4 dihedrals and try to optimize all at once?!( >>>>>>>> >>>>>>>> because all 4 dihedrals of O-C-CH2-CH3 seems equivalent to me). Am I >>>>>>> right? >>>>>>> >>>>>>> >>>>>>> No, there are 2 O-C-CH2-CH3 dihedrals, not 4. >>>>>>> >>>>>> >>>>>> >>>>>> How come? :-) there are two ethyls and two oxygens in the ring. Each >>>>> ethyl >>>>> forms two dihedrals (one with each oxygen). >>>>> Thus, there would be 4 dihedrals. And all of these 4 should be >>>>> equivalent. >>>>> >>>>> >>>>> Sorry, mentally jumping ahead. By atom typing, yes, there are 4 that >>>> are >>>> the same, but there are still only 2 equivalent rotatable bonds that >>>> need >>>> to be addressed. So you can define the rotation in the PES as either of >>>> the two for a given rotatable bond, they're redundant. So choose one >>>> and >>>> fit to it. You are, in any case, doing a PES of the entire molecule, so >>>> the fitting accounts for the redundancy. >>>> >>>> -Justin >>>> >>>> >>>> >>>> >>>>>> >>>>>> >>>>>>>> Deactivate the restraint, obtain the potential energy of the >>>>>>>> molecule >>>>>>>> >>>>>>>>> via >>>>>>>>> >>>>>>>>> mdrun -rerun and plot as a function of the dihedral. >>>>>>>>> >>>>>>>>> >>>>>>>>>> >>>>>>>>>> This should be a zero step EM, right? The molecule should not be >>>>>>>>>> >>>>>>>>> allowed >>>>>>>>> to >>>>>>>>> change its conformation. >>>>>>>>> >>>>>>>>> >>>>>>>>> No, a zero-step MD. EM actually changes the coordinates before >>>>>>>>> step >>>>>>>>> >>>>>>>>> zero. >>>>>>>> >>>>>>>> A bit of shell scripting and careful topology modification and >>>>>>>> this >>>>>>>> can >>>>>>>> >>>>>>>> >>>>>>>> be done. >>>>>>>>> >>>>>>>>> >>>>>>>>>> >>>>>>>>>> Two more questions on this part: >>>>>>>>>> >>>>>>>>>> 1) I am using the QM scanning data and partial charges from GAAMP. >>>>>>>>> When >>>>>>>>> I >>>>>>>>> do this MM scanning, do I need to exclude any 1-4 interactions or I >>>>>>>>> can >>>>>>>>> behave this dihedral as other dihedrals? >>>>>>>>> >>>>>>>>> >>>>>>>>> 1-4 interactions are always at full strength in CHARMM. >>>>>>>>> >>>>>>>> >>>>>>>> 2) this dihedral is part of a lipid. >>>>>>>> >>>>>>>> Do I need to do these on only "one Lipid in vacuum" or >>>>>>>> >>>>>>>>> OR >>>>>>>>> on all lipids in "a bilayer in vacuum" or in a "bilayer in solvent" >>>>>>>>> >>>>>>>>> I think it should be "one Lipid in vacuum". >>>>>>>>> >>>>>>>>> >>>>>>>>> One model compound in vacuum, from which you will construct the >>>>>>>>> lipid. >>>>>>>>> >>>>>>>>> >>>>>>>> How if this compound (which is small part of the lipid) is charged? >>>>>>>> >>>>>>>> Should >>>>>>> it be still in vacuum? >>>>>>> >>>>>>> >>>>>>> Yes. The fact that a molecule is irrelevant except for when >>>>>> considering >>>>>> the dipole moment. >>>>>> >>>>>> Is there any specific consideration to be made in zero-step MD? e.g. a >>>>>> big >>>>>> >>>>>> simulation box or specific parameter in .mdp file? >>>>>>> >>>>>>> >>>>>>> Infinite cutoffs, no PBC. >>>>>>> >>>>>> >>>>>> -Justin >>>>>> >>>>>> >>>>>> -- >>>>>> ================================================== >>>>>> >>>>>> Justin A. Lemkul, Ph.D. >>>>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow >>>>>> >>>>>> Department of Pharmaceutical Sciences >>>>>> School of Pharmacy >>>>>> Health Sciences Facility II, Room 629 >>>>>> University of Maryland, Baltimore >>>>>> 20 Penn St. >>>>>> Baltimore, MD 21201 >>>>>> >>>>>> jalem...@outerbanks.umaryland.edu | (410) 706-7441 >>>>>> http://mackerell.umaryland.edu/~jalemkul >>>>>> >>>>>> ================================================== >>>>>> -- >>>>>> Gromacs Users mailing list >>>>>> >>>>>> * Please search the archive at http://www.gromacs.org/Support >>>>>> /Mailing_Lists/GMX-Users_List before posting! >>>>>> >>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >>>>>> >>>>>> * For (un)subscribe requests visit >>>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >>>>>> send a mail to gmx-users-requ...@gromacs.org. >>>>>> >>>>>> >>>>>> >>>>> >>>>> >>>>> -- >>>> ================================================== >>>> >>>> Justin A. Lemkul, Ph.D. >>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow >>>> >>>> Department of Pharmaceutical Sciences >>>> School of Pharmacy >>>> Health Sciences Facility II, Room 629 >>>> University of Maryland, Baltimore >>>> 20 Penn St. >>>> Baltimore, MD 21201 >>>> >>>> jalem...@outerbanks.umaryland.edu | (410) 706-7441 >>>> http://mackerell.umaryland.edu/~jalemkul >>>> >>>> ================================================== >>>> -- >>>> Gromacs Users mailing list >>>> >>>> * Please search the archive at http://www.gromacs.org/Support >>>> /Mailing_Lists/GMX-Users_List before posting! >>>> >>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >>>> >>>> * For (un)subscribe requests visit >>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >>>> send a mail to gmx-users-requ...@gromacs.org. >>>> >>>> >>> >>> >>> >> -- >> ================================================== >> >> Justin A. Lemkul, Ph.D. >> Ruth L. Kirschstein NRSA Postdoctoral Fellow >> >> Department of Pharmaceutical Sciences >> School of Pharmacy >> Health Sciences Facility II, Room 629 >> University of Maryland, Baltimore >> 20 Penn St. >> Baltimore, MD 21201 >> >> jalem...@outerbanks.umaryland.edu | (410) 706-7441 >> http://mackerell.umaryland.edu/~jalemkul >> >> ================================================== >> -- >> Gromacs Users mailing list >> >> * Please search the archive at http://www.gromacs.org/Support >> /Mailing_Lists/GMX-Users_List before posting! >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> * For (un)subscribe requests visit >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >> send a mail to gmx-users-requ...@gromacs.org. >> > > > > -- > *Rewards work better than punishment ...* > -- *Rewards work better than punishment ...* -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.