[gmx-users] Pulling Mechanics

Justin Lemkul jalemkul at vt.edu
Tue Jan 17 13:46:14 CET 2017



On 1/16/17 3:04 PM, Alexander Yang wrote:
> On 1/13/17 6:37 PM, Alexander Yang wrote:
>>> Hi everyone,
>>>
>>> I am pulling a water molecule into a bilayer using an absolute reference
>>> (mdp file below). I have tried to adapt methodology in Justin Lemkul's
>>> umbrella sampling tutorial, but I have encountered a couple issues:
>>>
>>
>> If you're trying to get water to move into a bilayer, why are you using an
>> absolute restraint?  Intuitively, the groups should be the bilayer and the
>> water.
>>
>
> My goal is to construct an excess free energy profile of water (in the
> z-direction) in these bilayer systems to compare permeability. I'm using a
> z-constraining method that involves pulling a single water molecule to
> various z-coordinates (depths in the bilayer), constraining that molecule,
> and recording the forces on that molecule. My idea has been to use an
> absolute reference to anchor the water molecule to a particular
> z-coordinate while still being free to roam in the xy-plane. I'm not trying
> to introduce extra forces on the bilayer system itself to pull a water
> molecule inside.
>
>
>>> 1) The bilayer (already well-equilibrated) is either being pulled or
>>> drifting. Without the pull code, the bilayer does not move. I've defined
>>> center of mass motion removal for the non-waters, which I would expect to
>>> keep the bilayer still. I've copied the mdp file below (Tracer2117
>>> corresponds to 2117SOL, but I've defined the Tracer2117 group in the
>> index
>>> file as the 3 atoms that correspond to this water)
>>>
>>
>> Applying COM motion removal to one group and not another is highly
>> artificial.
>> Couple that with your use of an absolute restraint (also somewhat
>> artificial)
>> and you're layering bad things upon bad things.
>>
>
> Because I was using an absolute reference, I was hoping to keep the
> equilibrated bilayer still so the profile could be accurately constructed
> for a profile of z-coordinates.
>

But this is not what you will accomplish.  You'll get COM resetting in the 
bilayer, then your water will suffer from "flying ice cube" effects and spurious 
forces due to artificial COM motion.  This undercuts your desired approach.

-Justin

>>
>> -Justin
>>
>>> 2) I'm getting an error that the distance between the pull groups is
>> larger
>>> than 0.49 times the box size (the box is about 6.12 nm x 5.17 nm x 7.42
>>> nm). The origin of my absolute reference is only 0.1 nm away from the
>> water
>>> molecule. The way I've specified my pull-coord-rate, at the timestep when
>>> the pull distance is larger than half the box size, the distance between
>>> the origin of the reference at time 0 and the water molecule at time of
>>> crash exceeds half the box size. Is the reference position actually
>> moving,
>>> or does the pull-coord-rate specify a change to a quantity besides the
>>> reference position such that the water molecule can be pulled while the
>>> reference position stays where it originally was specified?
>>> @ time 0 ns: Water (3.794, 1.447, 1.150), reference origin (3.794, 1.447,
>>> 1.250)
>>> @ crash (time 47.19ns): Water(3.997, 0.032, 4.751). My pull-coord-rate is
>>> 7.7e-5 nm/ps, so would the reference be at (3.794, 1.447, 4.883)? If this
>>> is the case, I'm not sure how this distance violates the box size rule.
>>> However, the distance between reference @ t=0 and water @ t = 47.19ns,
>>> violates the 0.49 box size rule.
>>>
>>>
>>>
>>> ******* Mdp File ***********
>>>
>>>  Run MDP parameters
>>>
>>> integrator                = md
>>>
>>> dt                        = 0.002
>>>
>>> nsteps                    = 25000000
>>>
>>> comm-mode                 = Linear
>>>
>>> nstcomm                   = 1
>>>
>>> comm-grps                 = non-water
>>>
>>>
>>> ; Output parameters
>>>
>>> nstxout                   = 0
>>>
>>> nstvout                   = 0
>>>
>>> nstxtcout                 = 5000
>>>
>>> nstenergy                 = 5000
>>>
>>> nstlog                    = 5000
>>>
>>> nstfout                   = 0
>>>
>>>
>>> ; Bond parameters
>>>
>>> continuation              = yes
>>>
>>> constraint-algorithm      = lincs
>>>
>>> constraints               = all-bonds
>>>
>>> lincs-iter                = 1
>>>
>>> lincs-order               = 4
>>>
>>>
>>> ; Neighbor searching
>>>
>>> cutoff-scheme             = Verlet
>>>
>>> nstlist                   = 10
>>>
>>> rcoulomb                  = 1.4
>>>
>>> rvdw                      = 1.4
>>>
>>>
>>> ; Electrostatics
>>>
>>> coulombtype               = PME
>>>
>>> fourierspacing            = 0.16
>>>
>>> pme_order                 = 4
>>>
>>>
>>> ; Temperature coupling
>>>
>>> tcoupl                    = nose-hoover
>>>
>>> tc_grps                   = non-water   water
>>>
>>> tau_t                     = 0.4         0.4
>>>
>>> ref_t                     = 305         305
>>>
>>>
>>> ; Pressure coupling
>>>
>>> pcoupl                    = Parrinello-Rahman
>>>
>>> pcoupltype                = isotropic
>>>
>>> tau_p                     = 2.0
>>>
>>> ref_p                     = 1.0
>>>
>>> compressibility           = 4.5e-05
>>>
>>> refcoord_scaling          = com
>>>
>>>
>>> ; Misc stuff
>>>
>>> gen_vel                   = no
>>>
>>> pbc                       = xyz
>>>
>>> DispCorr                  = EnerPres
>>>
>>>
>>> ; Pull parameters
>>>
>>> pull                      = yes
>>>
>>> pull-nstxout              = 5000
>>>
>>> pull-nstfout              = 5000
>>>
>>> pull-ngroups              = 1
>>>
>>> pull-ncoords              = 1
>>>
>>> pull-group1-name          = Tracer2117
>>>
>>> pull-coord1-groups        = 0 1
>>>
>>> pull-coord1-type          = umbrella
>>>
>>> pull-coord1-geometry      = direction
>>>
>>> pull-coord1-origin        = 3.794    1.447    1.250
>>>
>>> pull-coord1-dim           = N N Y
>>>
>>> pull-coord1-rate          = 7.7e-05
>>>
>>> pull-coord1-vec = 0 0 1
>>>
>>> pull-coord1-k             = 500.0
>>>
>>> pull-coord1-start         = no
>>>
>>> *************
>>>
>>> I've attached URLS to snapshots below (time0 and time47). Red is water,
>>> grey is the lipid bilayer, yellow is the water molecule (enlarged for
>>> visibility), and I've placed axes in the center. Before the crash, the
>>> water molecule was definitely moving in the correct direction (positive
>> z).
>>> The bilayer drifted in the direction the water molecule was moving.
>>>
>>> (time0) http://imgur.com/a/jWu33
>>>
>>> (time47) http://imgur.com/a/nekmL
>>>
>>> Thanks for the help,
>>>
>>> Alex
>>>
>>
>> --
>> ==================================================
>>
>> 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
>>
>> ==================================================
>>
>>
>> ------------------------------
>>
>> Message: 4
>> Date: Sun, 15 Jan 2017 12:50:52 -0500
>> From: Justin Lemkul <jalemkul at vt.edu>
>> To: gmx-users at gromacs.org
>> Subject: Re: [gmx-users] Protein-Ligand Complex MD simulation ;
>>         Command for Number of Hydrogen bond plo
>> Message-ID: <f01a0916-e7c5-a06a-4699-0a93db6aa1ed at vt.edu>
>> Content-Type: text/plain; charset=windows-1252; format=flowed
>>
>>
>>
>> On 1/13/17 10:51 PM, Adarsh V. K. wrote:
>>> Dear all,
>>>
>>> Protein-ligand complex MD simulation using Gromacs 5.1.4
>>>
>>> Can you please tell me command for Number of Hydrogen bond plot? and
>>> Other interactions between protein and ligand?
>>>
>>
>> gmx hbond calculates hydrogen bonds.  Select the protein and the ligand to
>> get
>> the number of hydrogen bonds between them.
>>
>> -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
>>
>> ==================================================
>>
>>
>> ------------------------------
>>
>> Message: 5
>> Date: Sun, 15 Jan 2017 12:52:00 -0500
>> From: Justin Lemkul <jalemkul at vt.edu>
>> To: gmx-users at gromacs.org
>> Subject: Re: [gmx-users] permeation events- ion channels in bilayer
>> Message-ID: <cff0b0c0-0414-83ff-61f3-2dbf1d671b83 at vt.edu>
>> Content-Type: text/plain; charset=windows-1252; format=flowed
>>
>>
>>
>> On 1/14/17 7:18 AM, Alex Mathew wrote:
>>> Dear GMX users,
>>>
>>> We are looking for water permeation events across the ion-channel, which
>> is
>>> embedded in lipid bilayer. After one micro second MD simulation in
>>> Gromacs2016.1 i used the below script for calculating the  permeation
>>> events,
>>>
>>> http://www.ks.uiuc.edu/Training/Tutorials/science/
>>> nanotubes/files/permeation.tcl
>>>
>>> Unfortunately all the time it gives zero value (I have tried for
>> different
>>> upperEnd  and lowerEnd ), but i can visualize the water movements in VMD
>>> screen.
>>>
>>> Kindly provide your valuable suggestions, since i don't have much
>> expertise
>>> in coding/programming.
>>>
>>
>> If you define permeation as a water moving into the core of the membrane or
>> channel, you can use gmx select to calculate appropriate upper and lower
>> boundaries (based on, e.g. upper and lower leaflet P coordinates or some
>> suitable amino acids that define the boundaries of the channel) and then
>> calculate the number of waters that have an oxygen with a z-coordinate
>> between
>> these boundaries.
>>
>> -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
>>
>> ==================================================
>>
>>
>> ------------------------------
>>
>> Message: 6
>> Date: Sun, 15 Jan 2017 12:54:01 -0500
>> From: Justin Lemkul <jalemkul at vt.edu>
>> To: gmx-users at gromacs.org
>> Subject: Re: [gmx-users] Adding new residue - modified LYS
>> Message-ID: <8676dd45-cd80-5ffb-9fbb-0e40cfc17fcc at vt.edu>
>> Content-Type: text/plain; charset=utf-8; format=flowed
>>
>>
>>
>> On 1/14/17 2:38 PM, Dan Si wrote:
>>> Hi,
>>>
>>>
>>>
>>> I'm struggling with adding a new residue to my MD simulation. I would
>> like
>>> to create a residue (let's call it NML) that would be an original lysine
>>> (LYS) residue with methyl group attached to it's amide bond nitrogen. I'm
>>> using amber99sb-ildn force field for my simulations. Unfortunately, I'm
>> very
>>> new in using gromacs, so I would appreciate really !step-by-step!
>> tutorial :
>>> -)
>>>
>>>
>>>
>>>
>>> Do I need to perform QM calculation to obtain some of the parameters? If
>>> yes, should I use zwitter-ionic structure NH3-CH-(R)-C-(O)-O in the
>>> optimization?
>>>
>>
>> What do you learn from reading the AMBER99sb-ILDN paper (and preceding
>> literature) regarding AMBER parametrization methodology?
>>
>>>
>>>
>>>
>>> I understand, I need to specify the topology of the new residue and
>>> parameters according to amber99sb-ildn force field, but how exactly am I
>>> suppose to do it? Which files I need to modify before running pdb2gmx and
>>> how?
>>>
>>
>> http://www.gromacs.org/Documentation/How-tos/Adding_
>> a_Residue_to_a_Force_Field
>>
>> -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.
>>
>> End of gromacs.org_gmx-users Digest, Vol 153, Issue 72
>> ******************************************************
>>

-- 
==================================================

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

==================================================


More information about the gromacs.org_gmx-users mailing list