[gmx-users] Group NS scheme parameter for 2016.3 version

Du, Yu duyu at sioc.ac.cn
Tue Sep 26 15:37:41 CEST 2017


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-----------------------
----------------------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.

> 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.



More information about the gromacs.org_gmx-users mailing list