[gmx-developers] Group cutoff-scheme simulation temperature affected by nh-chain-length

Du, Yu duyu at sioc.ac.cn
Fri Nov 3 02:59:56 CET 2017


Hi Mark and Dear GMX-Developers, 


Mark, thank you very much. I got your point. The difference is about the combination rule (5.3 and 5.4).
Geometric mean of c6 (comb-rule=1) is different from the arithmetic mean of sigma (comb-rule=2)


amber99sb-ildn is designed under comb-rule=2, althought the average difference is little.


Is there any means to adapt the source code of src/gromacs/tables/forcetable.cpp or others to get tabulized potential work with amber.ff whose comb-rule=2?


Yu


-----Original Messages-----
From:"Mark Abraham" <mark.j.abraham at gmail.com>
Sent Time:2017-11-03 09:24:19 (Friday)
To: gmx-developers at gromacs.org
Cc: GMX-developer <gromacs.org_gmx-developers at maillist.sys.kth.se>
Subject: Re: [gmx-developers] Group cutoff-scheme simulation temperature affected by nh-chain-length


Hi,


On Fri, Nov 3, 2017 at 2:12 AM Du, Yu <duyu at sioc.ac.cn> wrote:

Hi Mark and Dear GMX-Developers, 


I tried the gmx check -s1 -s2, it gives comparison each recording steps.


That sounds like your comparing the trajectories. I intended you to compare the .tpr files, because the combinations are computed in that. Yours may differ.
 
But I need the reason of not changing comb-rule 2to1 and sigma,epsilon to c6,c12 in amber ff.



Could you explain why the combined parameters are different and whose combined parameters are different from? 


See the equations in section 5.3.2 of the reference manual. Only in some cases with the parameters for a pair of atom types be the same in the two rules. Whether that applies is up to you, but if so then the gmx compare will confirm your deduction.
 
I have adapted the sigma and epsilon of ffnonbonded.itp to c6 and c12, comb-rule 2 to 1, I think it's reasonable and these places (forcefield.itp and ffnonbonded.itp) are where one only needs to adapt.  


In so doing, you surely assumed that you were treating an A-A interaction. But the A-B interactions are the relevant ones for *combining* parameters.
 
But Prof. Hess said I should never change the comb-rule of ff. Thanks for any detailed information.


Because it does not work for e.g. proteins.


Mark
 
Best, 
Yu


-----Original Messages-----
From:"Mark Abraham" <mark.j.abraham at gmail.com>
Sent Time:2017-11-02 22:52:54 (Thursday)
To:gmx-developers at gromacs.org
Cc: GMX-developer <gromacs.org_gmx-developers at maillist.sys.kth.se>

Subject: Re: [gmx-developers] Group cutoff-scheme simulation temperature affected by nh-chain-length


Hi,


The general reason is that the combined parameters are different. Whether that matters depends on your system composition. You can see if that affects you by doing gmx compare -s1 -s2


Mark


On Thu, Nov 2, 2017 at 3:42 PM Du, Yu <duyu at sioc.ac.cn> wrote:

Dear Prof. Hess, 


I compared the MD results of comb-rule=1 and 2 amber99sb-ildn ff and found that there is little difference between them. The MD results of original amber ff and adapted are pasted in the end of this email.


Could you give the reason that change should never be made to the comb-rule and the coorsponding ffnonbonded.itp file?


Thanks a lot.


Regards, 
Yu


##########################original amber ff############################
        Statistics over 25001 steps using 251 frames


   Energies (kJ/mol)
          Angle    Proper Dih.  Improper Dih.          LJ-14     Coulomb-14
    4.02507e+03    5.15818e+03    2.75373e+02    2.11916e+03    1.69171e+04
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    6.20228e+04   -5.36969e+05    3.15506e+03   -4.43296e+05    8.45579e+04
   Total Energy  Conserved En.    Temperature Pressure (bar)   Constr. rmsd
   -3.58738e+05   -4.39556e+05    3.00296e+02   -1.81660e+02    0.00000e+00


   Total Virial (kJ/mol)
    3.01440e+04    1.75821e+01   -1.58821e+01
    1.84730e+01    3.00339e+04   -9.97169e+01
   -1.46410e+01   -1.00240e+02    3.00330e+04


   Pressure (bar)
   -1.84467e+02   -3.37456e+00    2.10386e+00
   -3.46046e+00   -1.75964e+02    1.04734e+01
    1.98421e+00    1.05238e+01   -1.84549e+02


      T-Protein  T-non-Protein
    2.99844e+02    3.00324e+02




        M E G A - F L O P S   A C C O U N T I N G


 NB=Group-cutoff nonbonded kernels    NxN=N-by-N cluster Verlet kernels
 RF=Reaction-Field  VdW=Van der Waals  QSTab=quadratic-spline table
 W3=SPC/TIP3p  W4=TIP4p (single or pairs)
 V&F=Potential and force  V=Potential only  F=Force only


 Computing:                               M-Number         M-Flops  % Flops
-----------------------------------------------------------------------------
 Pair Search distance check           14291.569978      128624.130     0.5
 NxN Ewald Elec. + LJ [F]            184675.631168    12188591.657    51.1
 NxN Ewald Elec. + LJ [V&F]            2250.700664      240824.971     1.0
 NxN Ewald Elec. [F]                 160299.489552     9778268.863    41.0
 NxN Ewald Elec. [V&F]                 1953.462840      164090.879     0.7
 1,4 nonbonded interactions             127.655106       11488.960     0.0
 Calc Weights                          2540.801628       91468.859     0.4
 Spread Q Bspline                     54203.768064      108407.536     0.5
 Gather F Bspline                     54203.768064      325222.608     1.4
 3D-FFT                               69761.190336      558089.523     2.3
 Solve PME                               48.401936        3097.724     0.0
 Reset In Box                            42.311124         126.933     0.0
 CG-CoM                                  42.412752         127.238     0.0
 Angles                                  88.678547       14897.996     0.1
 Propers                                139.355574       31912.426     0.1
 Impropers                               10.650426        2215.289     0.0
 Virial                                  10.413396         187.441     0.0
 Stop-CM                                  8.536752          85.368     0.0
 Calc-Ekin                              169.413876        4574.175     0.0
 Lincs                                  165.281870        9916.912     0.0
 Lincs-Mat                             3504.083832       14016.335     0.1
 Constraint-V                          1040.480273        8323.842     0.0
 Constraint-Vir                          12.528823         300.692     0.0
 Settle                                 583.493470      188468.391     0.8
-----------------------------------------------------------------------------
 Total                                                23873328.747   100.0
##########################original amber ff############################


##########################comb-rule=1 amber ff############################
        Statistics over 25001 steps using 251 frames


   Energies (kJ/mol)
          Angle    Proper Dih.  Improper Dih.          LJ-14     Coulomb-14
    4.01117e+03    5.20573e+03    2.78324e+02    2.02413e+03    1.68747e+04
        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    6.22595e+04   -5.37413e+05    3.13120e+03   -4.43628e+05    8.45492e+04
   Total Energy  Conserved En.    Temperature Pressure (bar)   Constr. rmsd
   -3.59079e+05   -4.40129e+05    3.00265e+02   -1.82026e+02    0.00000e+00


   Total Virial (kJ/mol)
    3.01455e+04   -9.34378e+01   -1.54004e+02
   -9.20285e+01    3.00302e+04   -4.38835e+01
   -1.53975e+02   -4.42383e+01    3.00377e+04


   Pressure (bar)
   -1.89288e+02    8.10551e+00    1.64126e+01
    7.96965e+00   -1.76137e+02    7.49214e-01
    1.64098e+01    7.83424e-01   -1.80654e+02


      T-Protein  T-non-Protein
    3.00317e+02    3.00262e+02




        M E G A - F L O P S   A C C O U N T I N G


 NB=Group-cutoff nonbonded kernels    NxN=N-by-N cluster Verlet kernels
 RF=Reaction-Field  VdW=Van der Waals  QSTab=quadratic-spline table
 W3=SPC/TIP3p  W4=TIP4p (single or pairs)
 V&F=Potential and force  V=Potential only  F=Force only


 Computing:                               M-Number         M-Flops  % Flops
-----------------------------------------------------------------------------
 Pair Search distance check           14608.461222      131476.151     0.6
 NxN Ewald Elec. + LJ [F]            183977.095776    12142488.321    51.1
 NxN Ewald Elec. + LJ [V&F]            2242.460760      239943.301     1.0
 NxN Ewald Elec. [F]                 159476.727584     9728080.383    40.9
 NxN Ewald Elec. [V&F]                 1943.414296      163246.801     0.7
 1,4 nonbonded interactions             127.655106       11488.960     0.0
 Calc Weights                          2540.801628       91468.859     0.4
 Spread Q Bspline                     54203.768064      108407.536     0.5
 Gather F Bspline                     54203.768064      325222.608     1.4
 3D-FFT                               69761.190336      558089.523     2.3
 Solve PME                               48.401936        3097.724     0.0
 Reset In Box                            42.345000         127.035     0.0
 CG-CoM                                  42.412752         127.238     0.0
 Angles                                  88.678547       14897.996     0.1
 Propers                                139.355574       31912.426     0.1
 Impropers                               10.650426        2215.289     0.0
 Virial                                  10.413396         187.441     0.0
 Stop-CM                                  8.536752          85.368     0.0
 Calc-Ekin                              169.413876        4574.175     0.0
 Lincs                                  168.515866       10110.952     0.0
 Lincs-Mat                             3570.772848       14283.091     0.1
 Constraint-V                          1043.984536        8351.876     0.0
 Constraint-Vir                          12.566057         301.585     0.0
 Settle                                 583.673648      188526.588     0.8
-----------------------------------------------------------------------------
 Total                                                23778711.227   100.0
##########################comb-rule=1 amber ff############################


-----Original Messages-----
From:"Du, Yu" <duyu at sioc.ac.cn>

Sent Time:2017-11-02 19:45:40 (Thursday)

To: GMX-developer <gromacs.org_gmx-developers at maillist.sys.kth.se>
Cc:
Subject: Re: [gmx-developers] Group cutoff-scheme simulation temperature affected by nh-chain-length

Thanks for Prof. Berk's prompt fix.


You advised that "You should never change the combination rule of a force field" in Bug #2286 and because tabulized potential needs comb-rule=1, what should I do to use the amber99sb-ildn ff? (ff needed for protein-ligand simulation system) Or in Gromacs, I can only use the group cutoff-scheme and tabulized potential with gromos ff?

Could you please explain the consequence I may get when I use the comb-func=1 amber99sb-ildn ff which is changed by myself?


Best, 
Yu

-----Original Messages-----
From:"Du, Yu" <duyu at sioc.ac.cn>
Sent Time:2017-11-02 14:00:40 (Thursday)
To: GMX-developer <gromacs.org_gmx-developers at maillist.sys.kth.se>
Cc:
Subject: Re: [gmx-developers] Group cutoff-scheme simulation temperature affected by nh-chain-length

Dear GMX-Developers and Shirts, 


To reproduce the the problem, I use the material in GROMACS Tutorials Lysozyme in Water.



Yes, Prof. Shirts. If I use standard vdwtype=cut-off, coulombtype=pme and cutoff-scheme=group, the simulation temperature accords with the ref_t.
If I use gromos ff and tabulized  potential, there is no temperature problem.


Now the problem is pinned on the relation between tabulized potential and thermostat under amber99sb-ildn ff.


I filed Bug #2286 in redmine.
https://redmine.gromacs.org/issues/2286


Hope you could give some advice about this issue.


Yu

-----Original Messages-----
From:"Michael R Shirts" <Michael.Shirts at Colorado.EDU>
Sent Time:2017-10-23 19:43:38 (Monday)
To: "gmx-developers at gromacs.org" <gmx-developers at gromacs.org>, GMX-developer <gromacs.org_gmx-developers at maillist.sys.kth.se>
Cc:
Subject: Re: [gmx-developers] Group cutoff-scheme simulation temperature affected by nh-chain-length



If you can file a bug report in redmine (http://redmine.gromacs.org/) with example files, that would be useful. Before you file it can you check if problem happens if you use standard cutoffs, instead of tabulated cutoffs? I wonder if there is some sort of code path that it’s going down.  I will bet that combination of inputs has only very rarely been used.  Also, you said the temperature jumps. Does it come back down to 310, or does it stay there?

 

Generally, I would recommend using v-rescale instead of Nose-Hoover.  Nose-Hoover has some known flaws; increasing the chain length is supposed to help things, but doesn’t really help that much for realistic systems, just with very simple ones.

 

However, any Nose-Hoover change length should in principle give the same average, and did when the code was written.  If it doesn’t, the code was broken somewhere, possibly because of lack of coverage in the regression tests and something getting out of sync.

 

Best,

~~~~~~~~~~~~~~~~

Michael Shirts

Associate Professor

michael.shirts at colorado.edu

http://www.colorado.edu/lab/shirtsgroup/

Phone: (303) 735-7860

Office: JSCBB C123

Department of Chemical and Biological Engineering

University of Colorado Boulder

 

From: <gromacs.org_gmx-developers-bounces at maillist.sys.kth.se> on behalf of "Du, Yu" <duyu at sioc.ac.cn>
Reply-To: "gmx-developers at gromacs.org" <gmx-developers at gromacs.org>
Date: Monday, October 23, 2017 at 12:48 AM
To: GMX-developer <gromacs.org_gmx-developers at maillist.sys.kth.se>
Subject: [gmx-developers] Group cutoff-scheme simulation temperature affected by nh-chain-length

 

Dear GMX Developer and Prof. Hess, 

 

I think this question is not common to GMX users, so I post it here.

 

Prof. Hess is very familiar with the tabulized potential and group cutoff scheme, so I also CC this email to you. Any advice will be appreciated. 

 

The Backgroud:

I simulated protein-ligand system (ref_t=300) with amber99sb-ildn ff. I used the normal table6-12.xvg tabulized potential to reproduce this problem.

 

The Problem:

In NVT simulation, integrator=md-vv, Tcoupl=nose-hoover, coulombtype=pme-user, vdwtype=user. I found that at the start of the simulation, there is a temperature jump to around 310K (average 310K in 100ps) higher than the ref_t (300K) under nh-chain-length=10.

 

1)I tuned and spotted the nh-chain-length parameter. If I used nh-chain-length=1 there is no temperature jump and temperature can converge very well to 300K.

 

2)I replicated the same experiment in verlet scheme, the temperature converged to ref_t 300K in 100ps and no jump in nh-chain-length=2, which can 50K temperature jump (to 350K).

 

The Question:

1)How will the nh-chain-length parameter affect the simulation temperature under the group cut-off scheme?

2)Is it valid to use NH (nh-chain-length=1) not NH chain (nh-chain-length>1) to simulate the protein-ligand system?

 

Thanks for any advice.

 

--

Du, Yu
PhD Student,
Shanghai Institute of Organic Chemistry

345 Ling Ling Rd., Shanghai, China. 

Zip: 200032, Tel: (86) 021 5492 5275

--
Gromacs Developers mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_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-developers or send a mail to gmx-developers-request at gromacs.org.

--
Du, Yu
PhD Student,
Shanghai Institute of Organic Chemistry
345 Ling Ling Rd., Shanghai, China. 
Zip: 200032, Tel: (86) 021 5492 5275
--
Gromacs Developers mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_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-developers or send a mail to gmx-developers-request at gromacs.org.

--
Du, Yu
PhD Student,
Shanghai Institute of Organic Chemistry
345 Ling Ling Rd., Shanghai, China. 
Zip: 200032, Tel: (86) 021 5492 5275
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20171103/d86ea97e/attachment-0001.html>


More information about the gromacs.org_gmx-developers mailing list