[gmx-users] Re: Use pull code to restrain the COM
Bin Liu
fdusuperstring at gmail.com
Fri Aug 9 18:48:24 CEST 2013
Hi Justin,
Thank you for your quick reply.
In Setting A, grompp reports
Pull group natoms pbc atom distance at start
reference at t=0
0 0 0
1 146 23185 -0.000 -0.000 -3.268
6.668 6.693 6.492
The values of reference at t=0 are from the .mdp file. They are the COM of
pull group 1 calculated by g_traj -ox -com. But I can not understand the
values of distance at start. Where does -3.268 come from? If pull group 0,
the reference group is set to be 6.668 6.693 6.492, and grompp can
calculate the COM of pull group 1, the distance should be 0.0 0.0 0.0
instead of -0.000 -0.000 -3.268.
In Setting B, grompp reports
Pull group natoms pbc atom distance at start
reference at t=0
0 0 0
1 146 23185 -0.000 -0.000 -3.268
0.000 0.000 -3.268
It seems Setting B doesn't agree with Setting A. Then pull-start=yes
doesn't have the effect as I guessed.
My situation is a bit complicated and challenging. I am simulating the
oligomerization of lipo-peptides inside a lipid bilayer. There are several
lipo-peptides in the system. Initially I placed them in the water, i.e.,
outside the bilayer. Then I found it takes too long (at the scale of
microseconds even seconds) for them to spontaneously insert into the
bilayer and form an oligomer. Then I artificially placed the lipopeptides
into the bilayer, and put them together to create an artificial oligomer. I
want to see if they can stay there. But obviously the artificial insertion
and oligomerization is unfavourable from the viewpoint of energy. So the
lipopeptides have a strong tendency to leave the bilayer. Now I want to
restrain the COM of each lipopeptide (at least in z direction) for some
time, but give them conformational and orientational degrees of freedom as
much as possible, to give rise to the possibility that the artificial
oligomer can find a favourable conformation and stay there after the
removal of restraints. Unfortunately GROMACS can not restrain the COMs of
several groups simultaneously. Then in one short simulation (say 100ps), I
can only restrain the COM of one lipopeptide and freeze the others (in z
direction). The whole process proceeds as:
MD (restrain COM of peptide A, freeze BCD) -> CG EM -> MD (restrain COM of
B, freeze ACD) -> CG EM
-> MD (restrain COM of C, freeze ABD) -> CG EM -> MD (restrain COM of D,
freeze ABC) -> Next cyle -> ....
>From my current experience, the EM step is necessary as in MD steps the
freezed lipopeptides developed a large amount of strain. If the strain can
not get released, the MD runs will explode in several hundreds of
picoseconds (GROMACS reports failures of LINCS and eventually segmentation
fault, an obvious indication that the system has exploded.) And pull-k1= 50
is indeed very small. The restrained lipopeptides can displace in z
direction by a quite large amount in only hundreds of picoseconds.
Could you suggest a different approach, or some suggestions to refine this
process? Thank you very much.
Best Regards
Bin
2013/8/9 <gmx-users-request at gromacs.org>
> Send gmx-users mailing list submissions to
> gmx-users at gromacs.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> or, via email, send a message with subject or body 'help' to
> gmx-users-request at gromacs.org
>
> You can reach the person managing the list at
> gmx-users-owner at gromacs.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of gmx-users digest..."
>
>
> Today's Topics:
>
> 1. Aw: [gmx-users] Umbrella sampling - position restraints
> (lloyd riggs)
> 2. Use pull code to restrain the COM (Bin Liu)
> 3. Re: Use pull code to restrain the COM (Justin Lemkul)
> 4. Re: Trying to replicate Aqvist's results (solvation free
> energy). (Andr? Farias de Moura)
> 5. Re: Trying to replicate Aqvist's results (solvation free
> energy). (Heymman)
> 6. g_wham -sym (Shima Arasteh)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Thu, 8 Aug 2013 23:19:59 +0200 (CEST)
> From: "lloyd riggs" <lloyd.riggs at gmx.ch>
> Subject: Aw: [gmx-users] Umbrella sampling - position restraints
> To: "Discussion list for GROMACS users" <gmx-users at gromacs.org>
> Message-ID:
>
> <trinity-790b128a-a41d-4a6e-87ac-908747401fc1-1375996798996 at 3capp-gmx-bs33
> >
>
> Content-Type: text/plain; charset="us-ascii"
>
> An HTML attachment was scrubbed...
> URL:
> http://lists.gromacs.org/pipermail/gmx-users/attachments/20130808/fe9e3133/attachment-0001.html
>
> ------------------------------
>
> Message: 2
> Date: Thu, 8 Aug 2013 18:43:09 -0400
> From: Bin Liu <fdusuperstring at gmail.com>
> Subject: [gmx-users] Use pull code to restrain the COM
> To: gmx-users at gromacs.org
> Message-ID:
> <
> CAHA8nLjJPYdupTamu50TzxkhtWd+gkPUR87QkCNzWJ3v37Siew at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1
>
> Hi All,
>
> I would like to restrain the COM of a molecule, say a protein, in my
> simulation. I found the pull code can do the job. However, I am not sure
> about several parameters in the .mdp file.
>
> For example, if I want to restrain the COM of the protein in only z
> direction, but not rigidly, I can specify the following parameters:
>
> *; Setting A*
> *pull = umbrella*
> *pull-geometry = position*
> *pull-dim = N N Y ; only applied to Z direction*
> *pull-ngroups = 1
> *
> *pull-group1 = Protein*
> *pull-start = no*
> *pull-init1 = 3.0 3.0 3.0 ; suppose the initial position of the COM
> is 3.0 3.0 3.0*
> *pull-k1 = 50*
>
> I am not confident about my understanding of the settings. Here I didn't
> specify *pull-group0*. Therefore the position of the reference point is 0.0
> 0.0 0.0. But *pull-init1* moves the reference point to *3.0 3.0 3.0. *I
> don't know the reasonable range for *pull-k1 *too. I want the COM can only
> have a small deviation from its initial position. Is 50 too large or too
> small?
>
> And I figured out another approach.
>
> *; Setting B*
> *pull = umbrella*
> *pull-geometry = position*
> *pull-dim = N N Y ; only applied to Z direction*
> *pull-ngroups = 1
> *
> *pull-group1 = Protein*
> *pull-start = yes*
> *pull-k1 = 50*
> *
> *
> The GROMACS manual says *pull-start=yes *adds the COM distance of the
> starting conformation to *pull-init. *Does it mean Setting B is equivalent
> to Setting A? If I set *pull-geometry= position, *does *pull-start=yes *add
> the initial position of the COM to *pull-init* as a vector? If it's true,
> then I don't need to calculate the initial position of the COM.
>
> Thanks a lot.
>
> Best
>
> Bin
>
>
> ------------------------------
>
> Message: 3
> Date: Thu, 08 Aug 2013 18:53:54 -0400
> From: Justin Lemkul <jalemkul at vt.edu>
> Subject: Re: [gmx-users] Use pull code to restrain the COM
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <52042182.1020702 at vt.edu>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
>
>
> On 8/8/13 6:43 PM, Bin Liu wrote:
> > Hi All,
> >
> > I would like to restrain the COM of a molecule, say a protein, in my
> > simulation. I found the pull code can do the job. However, I am not sure
> > about several parameters in the .mdp file.
> >
> > For example, if I want to restrain the COM of the protein in only z
> > direction, but not rigidly, I can specify the following parameters:
> >
> > *; Setting A*
> > *pull = umbrella*
> > *pull-geometry = position*
> > *pull-dim = N N Y ; only applied to Z direction*
> > *pull-ngroups = 1
> > *
> > *pull-group1 = Protein*
> > *pull-start = no*
> > *pull-init1 = 3.0 3.0 3.0 ; suppose the initial position of the COM
> > is 3.0 3.0 3.0*
> > *pull-k1 = 50*
> >
> > I am not confident about my understanding of the settings. Here I didn't
> > specify *pull-group0*. Therefore the position of the reference point is
> 0.0
> > 0.0 0.0. But *pull-init1* moves the reference point to *3.0 3.0 3.0. *I
>
> What information makes you say this? What is the grompp output?
>
> > don't know the reasonable range for *pull-k1 *too. I want the COM can
> only
> > have a small deviation from its initial position. Is 50 too large or too
> > small?
> >
>
> That depends on your definition of "small." I would suspect that a pull_k1
> value of only 50 is rather weak. There are no hard and fast rules here,
> unfortunately.
>
> > And I figured out another approach.
> >
> > *; Setting B*
> > *pull = umbrella*
> > *pull-geometry = position*
> > *pull-dim = N N Y ; only applied to Z direction*
> > *pull-ngroups = 1
> > *
> > *pull-group1 = Protein*
> > *pull-start = yes*
> > *pull-k1 = 50*
> > *
> > *
> > The GROMACS manual says *pull-start=yes *adds the COM distance of the
> > starting conformation to *pull-init. *Does it mean Setting B is
> equivalent
> > to Setting A? If I set *pull-geometry= position, *does *pull-start=yes
> *add
> > the initial position of the COM to *pull-init* as a vector? If it's true,
> > then I don't need to calculate the initial position of the COM.
> >
>
> With pull_geometry = position, restraints are determined based on
> pull_vec1,
> which defaults to (0,0,0) so the default may be appropriate.
>
> Out of curiosity, why is it necessary to restrain the protein in a given
> position? If there's a homogeneous solvent around it, the restraint is
> pointless and does cause problems with COM motion (as noted in the
> manual). If
> there are other molecules in the system, you can probably define a better
> reference point.
>
> -Justin
>
> --
> ==================================================
>
> Justin A. Lemkul, Ph.D.
> Postdoctoral Fellow
>
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 601
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
>
> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>
> ==================================================
>
>
> ------------------------------
>
> Message: 4
> Date: Thu, 8 Aug 2013 21:58:59 -0300
> From: Andr? Farias de Moura <moura at ufscar.br>
> Subject: Re: [gmx-users] Trying to replicate Aqvist's results
> (solvation free energy).
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID:
> <CAKTwKgq--YgMCvXC7rpTZCjGDGv7jxeBVzvOnmC=
> ekU8eH91cA at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1
>
> first thing that I notice is that you are using PBC whereas Aqvist employed
> a spherical boundary model. I also think that PME is not applicable for a
> spherical boundary, so the reciprocal space contributions may be considered
> as spurious if you want to reproduce the original results. Maybe you should
> read Aqvist's paper methodology section and compare carefully the
> description of the original model and yours.
>
> cheers
>
> Andre
>
> On Thu, Aug 8, 2013 at 5:37 PM, Francis de Lasalle <
> francis.de.lasalle at umontreal.ca> wrote:
>
> > Hello,
> >
> > As we're trying to get familiar with Gromacs, in order to be able to
> > perform solvation free energy calculations of ions, we're trying to
> > replicate the results obtained by Aqvist (the author of the article from
> > which were obtained Lennard-Jones parameters for alkaline and alkaline
> > earth ions, for the OPLS-AA force field). The solvation free energies
> we're
> > computing are all about 50 kJ/mol different from the values reported by
> > Aqvist. We verified that the Lennard-Jones parameters for our ions and
> for
> > the water are identical to that of Aqvist and we're also using SPC water.
> > We're wondering why we can't reproduce these results and any input would
> be
> > greatly appreciated. Here's the article :
> >
> > http://pubs.acs.org/doi/pdf/**10.1021/j100384a009<
> http://pubs.acs.org/doi/pdf/10.1021/j100384a009>
> >
> > We're using 200 lambda states with 2000000 steps incremented by 0.002 ps
> > (4000 ps). Here's our run file for Li+ :
> >
> > integrator = sd
> > nsteps = 2000000
> > dt = 0.002
> > nstenergy = 1
> > nstlog = 1
> > nstxout = 0
> > nstvout = 0
> > cutoff-scheme = group
> > rlist = 1.0
> > dispcorr = EnerPres
> > vdw-type = cut-off
> > rvdw = 1.0
> > coulombtype = pme
> > rcoulomb = 1.0
> > fourierspacing = 0.12
> > constraint-algorithm = SHAKE
> > continuation = no
> > shake-tol = 0.00001
> > constraints = all-angles
> > tcoupl = v-rescale
> > tc-grps = system
> > tau-t = 0.2
> > ref-t = 298
> > ref-p = 1
> > compressibility = 4.5e-5
> > tau-p = 5
> > free-energy = yes
> > couple-moltype = lithium
> > sc-power = 1
> > sc-sigma = 0.3
> > sc-alpha = 1.0
> > couple-intramol = no
> > couple-lambda1 = vdwq
> > couple-lambda0 = none
> > init-lambda-state = $LAMBDA
> > fep-lambdas = $LAMBDAS
> >
> > The box has been equilibrated for 1 ns at 298 K. The box is a
> dodecahedron
> > with 240 water molecules.
> >
> > Thanks a lot.
> >
> > Francis
> >
> >
> > --
> > gmx-users mailing list gmx-users at gromacs.org
> > http://lists.gromacs.org/**mailman/listinfo/gmx-users<
> http://lists.gromacs.org/mailman/listinfo/gmx-users>
> > * Please search the archive at http://www.gromacs.org/**
> > Support/Mailing_Lists/Search<
> http://www.gromacs.org/Support/Mailing_Lists/Search>before posting!
> > * Please don't post (un)subscribe requests to the list. Use the www
> > interface or send it to gmx-users-request at gromacs.org.
> > * Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<
> http://www.gromacs.org/Support/Mailing_Lists>
> >
> >
>
>
> --
> _____________
>
> Prof. Dr. André Farias de Moura
> Department of Chemistry
> Federal University of São Carlos
> São Carlos - Brazil
> phone: +55-16-3351-8090
>
>
> ------------------------------
>
> Message: 5
> Date: Fri, 9 Aug 2013 00:02:23 -0700 (PDT)
> From: Heymman <francis.de.lasalle at umontreal.ca>
> Subject: [gmx-users] Re: Trying to replicate Aqvist's results
> (solvation free energy).
> To: gmx-users at gromacs.org
> Message-ID: <1376031743678-5010412.post at n6.nabble.com>
> Content-Type: text/plain; charset=us-ascii
>
> Thank you for posting. I also noticed they were using a spherical boundary
> model; you think using PBC would justify a 50 kJ/mol difference ? As for
> PME, we wanted to use the classical Ewald summation method, but it turns
> out
> that Gromacs hasn't implemented it, yet, for free energy of solvation
> calculations. Other than that, the only noticeable difference I could find
> concerns the choice of cutoff radius for ion-water and water-water
> interactions.
>
>
>
> --
> View this message in context:
> http://gromacs.5086.x6.nabble.com/Trying-to-replicate-Aqvist-s-results-solvation-free-energy-tp5010407p5010412.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
>
>
> ------------------------------
>
> Message: 6
> Date: Fri, 9 Aug 2013 02:48:34 -0700 (PDT)
> From: Shima Arasteh <shima_arasteh2001 at yahoo.com>
> Subject: [gmx-users] g_wham -sym
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID:
> <1376041714.42130.YahooMailNeo at web162406.mail.bf1.yahoo.com>
> Content-Type: text/plain; charset=iso-8859-1
>
> Hi,
>
> I use the
> g_wham -it tpr-files.dat -if pullf-files.dat -o -hist -unit kca -sym
>
>
> I' d like to know if it is possible to symmetrize the profile around a
> non-zero point? forexample z=60?
>
> Thanks for your suggestions in advance.
>
>
> Sincerely,
> Shima
>
>
> ------------------------------
>
> --
> gmx-users mailing list
> gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>
> End of gmx-users Digest, Vol 112, Issue 28
> ******************************************
>
More information about the gromacs.org_gmx-users
mailing list