[gmx-users] MM dihedral scanning
Mohsen Ramezanpour
ramezanpour.mohsen at gmail.com
Wed Jan 18 04:39:22 CET 2017
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.mohsen at 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 <jalemkul at 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 <jalemkul at vt.edu> wrote:
>>>
>>>
>>>>
>>>> On 1/12/17 11:20 PM, Mohsen Ramezanpour wrote:
>>>>
>>>> On Thu, Jan 12, 2017 at 8:36 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>>>
>>>>>
>>>>>
>>>>>> On 1/12/17 10:21 PM, Mohsen Ramezanpour wrote:
>>>>>>
>>>>>> On Thu, Jan 12, 2017 at 6:31 PM, Justin Lemkul <jalemkul at vt.edu>
>>>>>> wrote:
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On 1/12/17 7:47 PM, Mohsen Ramezanpour wrote:
>>>>>>>>
>>>>>>>> On Thu, Jan 12, 2017 at 4:24 PM, Justin Lemkul <jalemkul at 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
>>>>>>
>>>>>> jalemkul at 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-request at 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
>>>>
>>>> jalemkul at 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-request at 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
>>
>> jalemkul at 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-request at gromacs.org.
>>
>
>
>
> --
> *Rewards work better than punishment ...*
>
--
*Rewards work better than punishment ...*
More information about the gromacs.org_gmx-users
mailing list