[gmx-users] MM dihedral scanning
Mohsen Ramezanpour
ramezanpour.mohsen at gmail.com
Wed Jan 18 18:56:57 CET 2017
On Wed, Jan 18, 2017 at 10:35 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> On 1/17/17 10:37 PM, Mohsen Ramezanpour 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
>>
>
> What are you restraining here?
The specific dihedral of interest for scan. explained more at the end.
>
>
> integrator = steep
>> emtol = 1000.0
>> nsteps = 5000
>> vdwtype = Cut-off
>> vdw-modifier = Force-switch
>>
>
> There is no force switching because the cutoffs are infinite. You'll
> probably get an error.
>
Thanks. I belive I should use "None" then. Based on manual: "None: Use an
unmodified Van der Waals potential. With the *group* scheme this means *no
exact cut-off* is used, energies and forces are calculated for all pairs in
the neighborlist.
"
And I have used the cutoff-scheme = group belove.
>
> 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
>>
>
> See above.
>
> 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
>>
>
> You don't need -nsteps.
Sure.
Regarding the restraints (I forward my last email here):
If I understood correctly, NONE of the following options are necessary (not
present) in Gromacs version 5.1.3:
dihre = yes
dihre_fc = 100 ; kJ/(mol rad^2)
dihre_tau = 0.0
nstdihreout = 250
And there is not any similar parameter in this version too.
However, the kfac in the following section should be the "real force" for
restraining.
So, for dihedral restraining during em part, including following lines in
the topol.top (or at the end of molecule.itp) should be okay with above
em.mdp parameters.
#ifdef POSRES
[ dihedral_restraints ]
; ai aj ak al type label phi dphi kfac
A B C D 1 1 120 0 10000
#endif
Which means that I want to apply a force of 10000 kj/mol/rad^2 to a
specific dihedral of ABCD to an angle of 120 for the purpose of scanning.
I chose dphi as "0" and I did not include "power" (because of version 5)
Cheers
>
>
> -Justin
>
>
>
>> 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.
>>>
>>>
>>
>>
>>
> --
> ==================================================
>
> 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