[gmx-users] Steered MD simulation of drug molecules across the plasma membrane in gromacs 5.0.4
Christopher Neale
chris.neale at alum.utoronto.ca
Wed Apr 29 16:30:26 CEST 2015
Your first test after I suggested mdp options was not exactly what I had suggested (you used a negative vector and a negative rate (not that this is inherently disallowed, just that it was not what I suggested). If you're pulling from negative z to positive z then the latest mdp options that you provided seem reasonable. i.e.,
pull = umbrella
pull_geometry = direction
pull_dim = N N Y
pull_coord1_vec =0 0 1
pull_start = yes
pull_ngroups =2
pull-ncoords =1
pull_coord1_groups = 1 2
pull_group1_name = NPROT_ref
pull_group2_name = LIG
pull_coord1_rate = 0.002 ; 0.002 nm per ps = 2 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_nstxout = 50 ; every 1 ps
pull_nstfout = 50 ; every 1 ps
Here, I expect that you should be able to pull entirely across the bilayer in a single continuous steered MD run. I can't comment on the definitions of the groups themselves, because that has changed format since I previously used it.
Now that I understand the intended usage of the new code, I expect that your US runs will also be fine with this format. However, if you want to check (and I would), then create a bizarre test system in which you have a box of water and place a new atom type in the center with that atom having zero charge and zero epsilon for LJ, and put a strong absolute position restraint on this atom. Now pull a single water molecule from negative z to positive z with reference to this dummy atom and use these conformations to start US and from US you should get a PMF that tends to perfect flatness as your simulation time increases and the statistical fluctuations are averaged out. Technically you should also be able to do this without the dummy particle for reference and just use absolute coordinates for pulling, though (a) I am not sure you can do that in gromacs v5, I have not checked, and (b) it's not quite as close a test system to what you want. I have no reason to suspect that you won't get exactly what you want at the bilayer center, just suggesting a possible test if you are concerned.
Chris.
________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of mahender singh <pharmbiochem at live.com>
Sent: 23 April 2015 04:12
To: gromacs user_list
Subject: Re: [gmx-users] Steered MD simulation of drug molecules across the plasma membrane in gromacs 5.0.4
hello
sorry for writing again and again
as from the previous parameters, now I had removed the -ve sign from the pull_coord1_vec =0 0 -1 and pull_coord1_rate = -0.002 as follow
; Pull code
pull = umbrella
pull_geometry = direction
pull_dim = N N Y
pull_coord1_vec =0 0 1
pull_start = yes
pull_ngroups =2
pull-ncoords =1
pull_coord1_groups = 1 2
pull_group1_name = NPROT_ref
pull_group2_name = LIG
pull_coord1_rate = 0.002 ; 0.002 nm per ps = 2 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_nstxout = 50 ; every 1 ps
pull_nstfout = 50 ; every 1 ps
Now I think simulation and the results looking correct and no strange behavior was observed.
So can you kindly comments on the setting which I have mentioned above and their correctness. (there may be chances that I may have understood the parameters wrongly), and from this simulation I want to do the umbrella sampling for free energy calculation.
with regards
Mahender Singh
> Hi
> Christopher Neale and Justin
>
> I used now following pull code (what i understood from the previous emails) to pull the drug across the membrane
> ; Pull code
> pull = umbrella
> pull_geometry = direction
> pull_dim = N N Y
> pull_coord1_vec =0 0 -1
> pull_start = yes
> pull_ngroups =2
> pull-ncoords =1
> pull_coord1_groups = 1 2
> pull_group1_name = NPROT_ref
> pull_group2_name = LIG
> pull_coord1_rate = -0.002 ; 0.002 nm per ps = 2 nm per ns
> pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
> pull_nstxout = 50 ; every 1 ps
> pull_nstfout = 50 ; every 1 ps
>
> there was no error and the molecule was able to cross the membrane from one side to the other side in a single SMD.
> But my concern here is the pullf.xvg output, because in the graph all the force is in negative sign. but if I consider the magnitude of the force, it's giving me the expected graph as there will be increase in the force when drug will enter in the membrane. Can you kindly comment of the statement.
>
> with regards
> Mahender Singh
>
>
> > Thank you dear Justin and Christopher Neale for your help.
> >
> > Solute is at -z relative to the bilayer and wants to pull toward +z from the bulk water on the -Z to the +z bulk water through the membrane.
> > I observed that the + and the - sign in the pull_coord1_rate is associated with increase or decrease in the COM distance between the pull group and the reference group, respectively.
> > I am trying to split simulation into two parts i.e. each leaflet as separate SMD simulation.
> > When solute is moving from the -z to the center of the bilayer, it will stop here as pull_coord1_rate = -0.0020 (decrease in the COM, pull force will be - in sign, am I right that sign in force values is negative because decrease in the distance between COM of two group with time is happening), still in the last frame, COM of the molecule will not moves to the +z side and I am not able to do the second part of the simulation by setting pull_start=yes , as if I will change the pull_coord1_rate = 0.0020 it will pull it back in the -z direction (pull force will be in positive sign) (pull_geometry= distance).
> >
> > @Christopher Neale
> > I will try the following setting and will let you know the results.
> >
> > @ Justin
> > When I am pulling the molecule from the center of bilayer to the bulk water and the same molecule from the bulk water to the center of bilayer (by changing the sign in the pull_coord1_rate = -/+0.0020 ). pull-force v/s time graphs are looking like following (first graph is the pulling of molecule from water to the center of bilayer and second graph is from the center of bilayer to the bulk water). My query was the - and the + sign of the pull force, which I think is due to the increase and decrease in the COM distance between the pull group and the reference group, am I right?
> >
> >
> >
> > with regards
> > Mahender singh
> >
> >
> > > Dear Justin:
> > >
> > > I think you are correct. It may not even really be that much of a hassle though one may need to be very careful in system setup for umbrellas near the center of the bilayer when using pull_start=yes rather than pull_coord1_init (i.e., one may need to make sure that the solute COM is on the side of the bilayer that one specifies with pull_coord1_vec, though having not used it myself perhaps this is not even an issue?). Separately, I presume that pull_coord1_init=0 should be identical with pull_coord1_vec = 0 0 1 or pull_coord1_vec = 0 0 -1, though if I was using the new code that would be something I would check at the outset.
> > >
> > > Dear Mahender:
> > >
> > > Based on Justin's suggestion and your initial mdp file, you should try something like the following. I'm still using gromacs 4.6.7 so I can not verify the accuracy of your other parameters (or even be entirely sure about these ones, but they are worth a try).
> > >
> > > ;; Presuming that you start with the solute at +z relative to the bilayer and want to pull toward -z
> > > pull = umbrella
> > > pull_geometry = direction
> > > pull_dim = N N Y
> > > pull_coord1_vec = 0 0 1
> > > pull_coord1_rate = -0.002
> > > pull_start = yes
> > >
> > > ;; I presume that
> > > ;; pull_coord1_vec = 0 0 1 and pull_coord1_rate = -0.002
> > > ;; is identical to
> > > ;; pull_coord1_vec = 0 0 -1 and pull_coord1_rate = 0.002
> > >
> > > ;; I am not sure if pull_dim = N N Y is necessary, but I don't see how it could hurt.
> > > ;; Presumably the distance and the vector on which the pull_coord1_rate acts are already going to be entirely along z due to pull_coord1_vec , but I'd test that too if not using pull_dim = N N Y
> > >
> > > Please report back whether this works or whether you get strange behaviour when you hit the bilayer center.
> > >
> > > ________________________________________
> > > From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Justin Lemkul <jalemkul at vt.edu>
> > > Sent: 22 April 2015 08:12
> > > To: gmx-users at gromacs.org
> > > Subject: Re: [gmx-users] Steered MD simulation of drug molecules across the plasma membrane in gromacs 5.0.4
> > >
> > > On 4/21/15 6:54 AM, Christopher Neale wrote:
> > > > Dear Mahender:
> > > >
> > > > It's a real pity that the pull_geometry = position has been removed. So now it's impossible to do umbrella sampling of a solute across a lipid bilayer properly? Anyway, I see a note about this removal here ( http://redmine.gromacs.org/issues/1346 ), though no reason is given. I'd suggest that you revert to a more functional version of gromacs ;)
> > > >
> > > > Hopefully I'm missing something and somebody else can point you to a better solution.
> > > >
> > >
> > > I think the behavior can be recovered with "direction" geometry, but that
> > > requires each half of the symmetric sampling windows to be set up separately, I
> > > think, rather than just a continuous vector like "position" used to provide.
> > > That's a shame if true. It's a workaround, but it's much less convenient.
> > >
> > > -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.
>
>
>
> ------------------------------
>
> Message: 4
> Date: Thu, 23 Apr 2015 01:18:01 -0600
> From: Alex <nedomacho at gmail.com>
> To: "gromacs.org_gmx-users at maillist.sys.kth.se"
> <gromacs.org_gmx-users at maillist.sys.kth.se>, <gmx-users at gromacs.org>
> Subject: Re: [gmx-users] Problems in the modelling of interaction
> between peptide and copper
> Message-ID: <151669015.20150423011801 at gmail.com>
> Content-Type: text/plain; charset=us-ascii
>
> Hi,
>
> Your model is not reasonable. First of all, let's forget about
> simulations for a moment, because GMX is a wildly inappropriate choice.
> Even if the energy was not zero (which in your case is probably due to
> some LJ parameters set up improperly, or some really good GMX code
> giving you a hint that this is a bad idea), the results would describe a
> protein next to a collection of LJ entities put together in a crystal
> lattice. Here is why:
>
> When exposed to air or water, copper surface
> would undergo quick oxidation. That aside, let us say we have a magic
> copper surface, which is not oxidated. This is even more problematic
> for MD-type methods, because metal surfaces don't really have
> conventional van der Waals due to presence of Fermi gas in
> metallic crystals, as opposed to the type of polarization in covalent
> dielectrics or even semiconductors.
>
> The only method somewhat reasonable for metals in MD I am aware of is
> the Modified Embedded Atom Method (MEAM), which, I believe, is implemented
> in LAMMPS. At the same time, it won't handle the protein.
>
> Alternatively, you can use DFTMD in CP2K and actually place
> the short region of your protein close to the surface of copper and
> see what happens. This in fact appears to be a half-decent idea and
> I'd try it. Depending on the exchange correlation functions and the
> basis you choose, you could probably do about a 100 atoms for a few
> picoseconds within reasonable time with parallelization over 40-50
> decent cores. One thing to keeo in mind is that
> the last time a checked (a few months back), CP2K had no electron
> energy sampling aside from the Gamma point (VASP can do that, but it
> is expensive). If that doesn't bother you, I would recommend going the DFT route. You could even try to
> build the energy surface for your protein-metal interaction and then
> try to use it to parameterize a lower level method such as MD.
>
> Hope this helps.
>
> Alex
>
>
> MT> Dear all,
>
> MT> I got a problem when modelling the interaction of a peptide and
> MT> copper. I am using the force field of opls-aa. The charge of
> MT> copper atoms are defined as zero and the non-bonded parameters are
> MT> obtained from the CU ions, which can be
> MT> found in opls-aa.
>
> MT> I applied smd and found that the peptide is approaching the
> MT> copper surface and stopped when it is very close to the surface,
> MT> due to some abnormal interaction. The whole process seems quite
> MT> reasonable, however, when I rerun the
> MT> simulation to obtain the interaction energy, for the vdw
> MT> interaction between peptide and copper is zero, absolute zero, in the whole process.
>
> MT> Can anyone give me some suggestions? Is my model reasonable? I am
> MT> not so confident with modelling strategy of using the LJ
> MT> parameters for copper bulk material. Also the zero interaction energy is another big issue.
>
> MT> Cheers,
> MT> Tng
>
>
>
>
>
> ------------------------------
>
> Message: 5
> Date: Thu, 23 Apr 2015 01:18:01 -0600
> From: Alex <nedomacho at gmail.com>
> To: "gromacs.org_gmx-users at maillist.sys.kth.se"
> <gromacs.org_gmx-users at maillist.sys.kth.se>, <gmx-users at gromacs.org>
> Subject: Re: [gmx-users] Problems in the modelling of interaction
> between peptide and copper
> Message-ID: <151669015.20150423011801 at gmail.com>
> Content-Type: text/plain; charset=us-ascii
>
> Hi,
>
> Your model is not reasonable. First of all, let's forget about
> simulations for a moment, because GMX is a wildly inappropriate choice.
> Even if the energy was not zero (which in your case is probably due to
> some LJ parameters set up improperly, or some really good GMX code
> giving you a hint that this is a bad idea), the results would describe a
> protein next to a collection of LJ entities put together in a crystal
> lattice. Here is why:
>
> When exposed to air or water, copper surface
> would undergo quick oxidation. That aside, let us say we have a magic
> copper surface, which is not oxidated. This is even more problematic
> for MD-type methods, because metal surfaces don't really have
> conventional van der Waals due to presence of Fermi gas in
> metallic crystals, as opposed to the type of polarization in covalent
> dielectrics or even semiconductors.
>
> The only method somewhat reasonable for metals in MD I am aware of is
> the Modified Embedded Atom Method (MEAM), which, I believe, is implemented
> in LAMMPS. At the same time, it won't handle the protein.
>
> Alternatively, you can use DFTMD in CP2K and actually place
> the short region of your protein close to the surface of copper and
> see what happens. This in fact appears to be a half-decent idea and
> I'd try it. Depending on the exchange correlation functions and the
> basis you choose, you could probably do about a 100 atoms for a few
> picoseconds within reasonable time with parallelization over 40-50
> decent cores. One thing to keeo in mind is that
> the last time a checked (a few months back), CP2K had no electron
> energy sampling aside from the Gamma point (VASP can do that, but it
> is expensive). If that doesn't bother you, I would recommend going the DFT route. You could even try to
> build the energy surface for your protein-metal interaction and then
> try to use it to parameterize a lower level method such as MD.
>
> Hope this helps.
>
> Alex
>
>
> MT> Dear all,
>
> MT> I got a problem when modelling the interaction of a peptide and
> MT> copper. I am using the force field of opls-aa. The charge of
> MT> copper atoms are defined as zero and the non-bonded parameters are
> MT> obtained from the CU ions, which can be
> MT> found in opls-aa.
>
> MT> I applied smd and found that the peptide is approaching the
> MT> copper surface and stopped when it is very close to the surface,
> MT> due to some abnormal interaction. The whole process seems quite
> MT> reasonable, however, when I rerun the
> MT> simulation to obtain the interaction energy, for the vdw
> MT> interaction between peptide and copper is zero, absolute zero, in the whole process.
>
> MT> Can anyone give me some suggestions? Is my model reasonable? I am
> MT> not so confident with modelling strategy of using the LJ
> MT> parameters for copper bulk material. Also the zero interaction energy is another big issue.
>
> MT> Cheers,
> MT> Tng
>
>
>
>
>
> ------------------------------
>
> --
> 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 132, Issue 106
> *******************************************************
--
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.
More information about the gromacs.org_gmx-users
mailing list