[gmx-users] Group NS scheme parameter for 2016.3 version
Mark Abraham
mark.j.abraham at gmail.com
Tue Sep 26 15:55:09 CEST 2017
Hi,
On Tue, Sep 26, 2017 at 3:38 PM Du, Yu <duyu at sioc.ac.cn> wrote:
> Thank Mark's reply.
>
> > -----Original Messages-----
> > From: "Mark Abraham" <mark.j.abraham at gmail.com>
> > Sent Time: 2017-09-26 19:44:41 (Tuesday)
> > To: gmx-users at gromacs.org
> > Cc:
> > Subject: Re: [gmx-users] Group NS scheme parameter for 2016.3 version
> >
> > Hi,
> >
> > On Tue, Sep 26, 2017 at 4:02 AM Du, Yu <duyu at sioc.ac.cn> wrote:
> >
> > > > -----Original Messages-----
> > > > From: "Mark Abraham" <mark.j.abraham at gmail.com>
> > > > Sent Time: 2017-09-26 06:38:40 (Tuesday)
> > > > To: "Discussion list for GROMACS users" <gmx-users at gromacs.org>,
> > > "Discussion list for GROMACS users" <
> > > gromacs.org_gmx-users at maillist.sys.kth.se>
> > > > Cc:
> > > > Subject: Re: [gmx-users] Group NS scheme parameter for 2016.3 version
> > > >
> > > > Hi,
> > > >
> > > > On Mon, 25 Sep 2017 14:11 Du, Yu <duyu at sioc.ac.cn> wrote:
> > > >
> > > > > Dear GMX Users,
> > > > >
> > > > >
> > > > > Indeed making tabulized potential between groups with Verlet is not
> > > > > trivial as also discussed in OpenMM issue #1765.
> > > > >
> > > > > I'm using gmx 2016.3. If I want use the group NS scheme as precise
> as
> > > > > possible, what parameters in mdp file should I set and what are
> their
> > > > > recommended values for protein and ligand system.
> > > > >
> > > >
> > > > It's straightforward to implement a tabulated interaction kernel, but
> > > > rather less clear how best to let the user describe to mdrun how they
> > > want
> > > > the calculation to work, eg whether interaction types should be
> governed
> > > by
> > > > atom numbers, or types, or names, and how that should be expressed
> in the
> > > > topology. Suggestions for problem types people want to be able to
> handle
> > > > are most welcome.
> > > >
> > >
> > > I want to decompose the interaction between protein and ligand, and
> only
> > > between them, not others.
> > > So I think "vdwtype=user, energygrps = Protein Ligand, energygrp-table
> =
> > > Protein Ligand and table.xvg table_Protein_ligand.xvg" is the only
> solution
> > > in gromacs to achieve it.
> > >
> >
> > Yes. So in future, some method that would select the regions to use the
> > different tabulated interactions based on atom number (ie index groups)
> > would work for you.
>
> Today I experimented on the tabulated potential with group cut-off scheme.
> I used custom potential (no coulomb and no attractive vdw) between protein
> and ligand. I guess it worked. Could you please have a look at the md.log
> and md.mdp file? and help me assure that there is no question.
> ----------------------md.log-----------------------
> Table routines are used for coulomb: true
> Table routines are used for vdw: true
> Cut-off's: NS: 1 Coulomb: 1 LJ: 1
> Long Range LJ corr.: <C6> 3.4946e-04
> System total charge: -0.000
> Read user tables from md_0_3.xvg with 1001 data points.
> Tabscale = 500 points/nm
> Read user tables from md_0_3_Protein_MEL.xvg with 1001 data points.
> Tabscale = 500 points/nm
> Read user tables from md_0_3.xvg with 1001 data points.
> Tabscale = 500 points/nm
> Read user tables from md_0_3.xvg with 1001 data points.
> Tabscale = 500 points/nm
> ...............................
> Statistics over 100001 steps using 1001 frames
>
> Energies (kJ/mol)
> Bond Angle Proper Dih. Improper Dih. LJ-14
> 4.09718e+03 1.11493e+04 1.26265e+04 8.33816e+02 4.52508e+03
> Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) Potential
> 4.81316e+04 4.01305e+04 -2.28798e+03 -4.91463e+05 -3.72257e+05
> Kinetic En. Total Energy Conserved En. Temperature Pres. DC (bar)
> 8.98371e+04 -2.82420e+05 2.61416e+06 3.37891e+02 -2.51191e+02
> Pressure (bar) Constr. rmsd
> 5.71728e+02 0.00000e+00
>
> Total Virial (kJ/mol)
> 2.46727e+04 1.69577e+01 5.39006e+01
> 1.69965e+01 2.48000e+04 -4.34063e+01
> 5.38208e+01 -4.28545e+01 2.47416e+04
>
> Pressure (bar)
> 5.78692e+02 -5.85172e+00 -3.91771e+00
> -5.85597e+00 5.67516e+02 7.54656e+00
> -3.90895e+00 7.48597e+00 5.68976e+02
>
> Epot (kJ/mol) Coul-SR LJ-SR Coul-14 LJ-14
> Protein-Protein -8.49851e+04 -8.95750e+03 4.83260e+04 4.44865e+03
> Protein-MEL 0.00000e+00 1.57386e+01 0.00000e+00 0.00000e+00
> Protein-rest -3.16047e+04 -1.50902e+03 0.00000e+00 0.00000e+00
> MEL-MEL -3.28886e+02 -2.02938e+01 -1.94359e+02 7.64322e+01
> MEL-rest -1.16305e+03 -6.90968e+01 0.00000e+00 0.00000e+00
> rest-rest -3.73381e+05 5.06707e+04 0.00000e+00 0.00000e+00
>
> T-MEL_NA_Protein T-SOL_CL
> 3.41647e+02 3.37007e+02
> ----------------------md.log-----------------------
>
Seems plausible.
> ----------------------md.mdp-----------------------
> ; Run parameters
> integrator = md-vv
> nsteps = 100000 ; 200ps
> dt = 0.002
>
> ; Output parameters
> nstxout = 0
> nstvout = 0
> nstenergy = 500 ; 1ps
> nstlog = 500
> nstxout-compressed = 500
> compressed-x-grps = System
> energygrps = Protein MEL
> energygrp-table = Protein MEL
>
> ; Neighbor searching
> nstlist = 5 ; 10fs
> ns-type = Grid
> pbc = xyz
> rlist = 1.0
> cutoff-scheme = group
> nstcalclr = -1
>
> ; Electrostatics and VdW
> coulombtype = user
> rcoulomb = 1.0
> vdwtype = user
> rvdw = 1.0
>
> ; Temperature coupling
> Tcoupl = nose-hoover
> tc-grps = MEL_NA_Protein SOL_CL
> tau_t = 1.0 1.0
> ref_t = 300 300
>
> ; Pressure coupling
> Pcoupl = no
> ; Dispersion correction
> DispCorr = EnerPres
>
> ; Velocity generation
> gen_vel = no
>
> ; Bond parameters
> continuation = yes
> constraint-algorithm = lincs
> constraints = h-bonds
> lincs-order = 6
> lincs-iter = 2
> ----------------------md.mdp-----------------------
>
> >
> > > How to set up the nonbonded scheme for a given force field varies with
> the
> > > > force field, either based on how it was parameterized or been shown
> to
> > > work
> > > > in practice. This is standard practice. A major defect of the group
> > > scheme
> > > > is that it is inefficient with a buffered list, but that is your
> tradeoff
> > > > to make.
> > > >
> > >
> > > Yes, this is my major confusion. Amber FF don't use charge group not
> like
> > > GROMOS FF. I need to work with Amber FF, so is it proper to use
> > > amber99sb-ildn and gaff with group cut-off scheme?
> > >
> >
> > Using charge groups of size 1 is the same as not using charge groups,
> which
> > is how the AMBER force fields were already implemented, long before the
> > Verlet scheme was implemented. You should also be reading other people's
> > work that used GROMACS and the AMBER force fields...
> >
>
> Yes, I'm really cramming.
>
> >
> > > And now the manual and user guide is dominated by Verlet cut-off
> scheme,
> > > where can I find the guide of using group cut-off scheme? I can't tell
> > > which options is for group scheme in the Molecular dynamics parameters
> > > (.mdp).
> > >
> >
> > Nothing's been removed, and nearly everything works the same way.
> >
> http://manual.gromacs.org/documentation/2016.4/user-guide/mdp-options.html#neighbor-searching
> > and
> > the coulomb and vdw sections are the only bits that matter, and they note
> > the points of difference.
> >
>
> Thanks for the clue, I'll note that part.
>
> >
> > > The following is what I found in
> > >
> regressiontests-2016.3/kernel/nb_kernel_ElecCSTab_VdwCSTab_GeomW4W4/grompp.mdp
> > > ; NEIGHBORSEARCHING PARAMETERS
> > > ; cut-off scheme (group: using charge groups, Verlet: particle based
> > > cut-offs)
> > > cutoff-scheme = Group
> > > ; nblist update frequency
> > > nstlist = 1
> > > ; ns algorithm (simple or grid)
> > > ns-type = Grid
> > > ; Periodic boundary conditions: xyz, no, xy
> > > pbc = xyz
> > > periodic-molecules = no
> > > ; Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,
> > > ; a value of -1 means: use rlist
> > > verlet-buffer-drift = 0.005
> > > ; nblist cut-off
> > > rlist = 0.99
> > > nstcalclr = -1
> > >
> > >
> > > Although it use group cut-off scheme, there is also a
> verlet-buffer-drift,
> > > which doesn't show in the 2016.3 version .mdp option.
> > >
> >
> > Don't use a test input for designing a real simulation. Those are
> designed
> > to do very specific things, including hit corner cases that we think
> nobody
> > should use. And that parameter is an old name that is not a problem when
> > the group scheme is chosen.
>
> Of course, I won't use the test .mdp file for real simulation. I use it
> for learning operational options, which I also found in older version
> Gromacs manual. Thanks for Gromacs' good archive. :)
>
> >
> > You should find some literature that shows the force field and software
> > working together, and test it yourself also. Unfortunately I know of
> nobody
> > who's stepped up and published "this is how you should run *these kinds*
> of
> > simulations" for AMBER ff in GROMACS.
>
> Yeah, this part is definitely my job. Thanks for your points.
>
> >
> > So, lots of thanks for anyone's advice for setting group scheme.
> > >
> > > > Given that you are modifying the force field to add different
> > > interactions,
> > > > you will anyway have to show that the modified form of the force
> field
> > > > works suitably. Even if the Verlet scheme had support for tabulated
> > > > interactions, you would still be limited to having only one set of
> > > > nonbonded parameters per atom type. Is your mix of interaction types
> > > > feasible to use with such a restriction?
> > > >
> > >
> > > I'm really looking forward to the combination of Verlet cut-off scheme
> > > with the tabulated potential. :)
> > >
> >
> > That wasn't the question. Can your ligand atoms have *one* set of vdw
> > parameters that can interact with solvent, ligand and protein atoms
> > according to different functional forms? People sometimes ask for such
>
> Only the protein and ligand interaction is under the effect of tabulated
> potential. If I fully understand the Tabulated potential tutorial, what I
> set doesn't affact the normal interaction between other groups.
>
Yes. But the shape of the interaction is only one of the relevant things.
The parameters, which scale the shape, are another. You have to have the
same *values*. Maybe that's for free with WCA, but it wouldn't be with say
Buckingham. There you would have to do something creative with scaling of
the shape.
Mark
> > things without thinking through whether their model makes sense (and if
> it
> > does, whether they would require schizophrenic atoms).
>
> Turning off the specific interaction between protein and ligand not other
> groups. I think this is really what a crazy me need. :)
>
> Yu
>
> >
> > Mark
> >
> > > Mark
> > > >
> > > > >
> > > > --
> > > > 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.
> > >
> > > --
> > > 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.
> > >
> > --
> > 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.
>
> --
> 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