[gmx-users] MM dihedral scanning
Mohsen Ramezanpour
ramezanpour.mohsen at gmail.com
Wed Jan 18 04:37:21 CET 2017
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 ...*
More information about the gromacs.org_gmx-users
mailing list