[gmx-users] Group NS scheme parameter for 2016.3 version
Du, Yu
duyu at sioc.ac.cn
Tue Sep 26 16:40:18 CEST 2017
Hi,
> -----Original Messages-----
> From: "Mark Abraham" <mark.j.abraham at gmail.com>
> Sent Time: 2017-09-26 21:54:54 (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 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.
Thanks for checking.
>
>
> > ----------------------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.
>
Absolutely, I guess WCA is the best solution for vdw decomposition. But even Gromacs can't make WCA in tabulated potential because the h(r) or g(r) has to contain the sigma from the result of combination rule. As far as I know this is not supported in current xvg file format.
So here I propose adding topology variables and operator in tabulated xvg file, e.g. $sigma, $epsilon, $c6, $c12 and +, -, *, /, ^..., which will add more flexibility to Gromacs' custom function. Excited! :)
Yu
> 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.
> >
> --
> 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