[gmx-users] table-extension distance appears to be over-ridden by .mdp file cut-off values

Elizabeth Ploetz ploetz at ksu.edu
Wed Dec 18 21:39:44 CET 2013


Dear Gromacs Users,

I'm trying to use tabulated VdW potentials that are non-zero until image-distance/2.0, but they seem to be over-ridden by the cut-offs in the .mdp file. I'm hoping someone with experience with the use of tabulated potentials plus the table-extension .mdp file option can tell me if I have not made all of the necessary changes so that the tables are correctly implemented. Here are the details of what I've done.

I have five systems each consisting of a single Lennard-Jones (6-12) sphere in TIP3P water (rhombic dodecahedron, 15nm image-distance for each system). The Lennard-Jones parameters for the infinitely-dilute LJ-sphere are as follows:
System Sigma (nm) Epsilon (kJ/mol)
1             0.3             0.1
2             1.5             0.1
3             3.0             0.1
4             4.5             0.1
5             6.0             0.1
For each system, rlist=rvdw=rcoulomb=0.8 (nm). vdw-type=User, coulombtype=PME, table-extension=6.7 (nm). 6.7 nm was chosen because 6.7nm+0.8nm=7.5nm=half the image-distance. A tabulated potential is only used so that the LJ-sphere - water interactions can have a different (larger) cut-off than the water-water interactions since the sigma/2.0 values for the LJ-sphere are greater than the 0.8 nm cut-off for Systems 3-5. 

When I run these systems (version 4.6), the LJ-sphere-to-water center-of-mass RDFs look correct for Systems 1 and 2 (the two LJ-spheres that have a radius less than the rlist=rvdw=rcoulomb cut-off). The RDFs look "wrong" for all larger spheres, and worse and worse as the Lennard-Jones sigma values increase. By "wrong," I mean that the RDFs take on positive values at too short of radial distances (i.e., water is being "sucked in" too close to the Lennard-Jones sphere). This makes sense if my tabulated VdW-potentials are not really being used (or, equivalently, not being used beyond 0.8 nm), and consequently the waters do not "see" the LJ-sphere until they are <0.8 nm away from the sphere. The simulations for the larger LJ-spheres are also unstable, which would also make sense, because waters that could get so close to the sphere would be well within the repulsive region of the sphere-to-water potential once they finally suddenly "saw" the sphere at <0.8 nm. The types of errors I get when the systems crash (usually after around 20 ps) are domain decomposition errors, specifically:

Fatal error:
A charge group moved too far between two domain decomposition steps
This usually means that your system is not well equilibrated
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

Or (when I change to particle-decomposition or use only one processor), the simulations eventually crash due to LINCS. 

The details of my .mdp file and tabulated potentials are below my signature, as well as an example of the .log file output.

Have I missed any extra flags that need to be defined so that the tabulated potentials are really used [the .log file output suggests they are really being used and so does mdrun.debug (not attached, but I can provide)]? Or, does anyone know of a way to check whether or not the tabulated potentials are really being used up to the cut-off + table-extension and not just up to the cut-off?

Please note that, as a test, I did a series of energy minimizations on a simpler system, the Lennard-Jones sphere (sigma = 6 nm) and a single water molecule, when varying the coordinates of the water molecule systematically between 7.5 nm and 0.3 nm away from the LJ-sphere. The LJ(SR) sphere-water potential energy was zero for all distances except when the water molecule was < 0.8 nm from the LJ-sphere. This, again, indicates that the rlist=rcoul=rvdw=0.8 nm cut-off distance is being used for the sphere-water interactions instead of the extended distance.

Sincerely,

Elizabeth Ploetz
-----
EXTRA INFORMATION:

(1) Relevant portion of my .mdp file:

define                   = -DFLEXIBLE
; RUN CONTROL PARAMETERS
integrator               = md
; Start time and timestep in ps
tinit                    = 0
dt                       = 0.002
nsteps                   = 50000
; For exact run continuation or redoing part of a run
init_step                = 0
; mode for center of mass motion removal
comm-mode                = Linear
; number of steps for center of mass motion removal
nstcomm                  = 500
; group(s) for center of mass motion removal
comm-grps                = 

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                  = 500
nstvout                  = 0
nstfout                  = 0
; Checkpointing helps you continue after crashes
nstcheckpoint            = 5000
; Output frequency for energies to log file and energy file
nstlog                   = 500
nstenergy                = 500 
; Output frequency and precision for xtc file
nstxtcout                = 0
xtc-precision            = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps                 = 

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist                  = 1
nstcalclr                = 1
; ns algorithm (simple or grid)
ns-type                  = Grid
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
pbc                      = xyz
; nblist cut-off        
rlist                    = 0.8

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype              = PME
rcoulomb-switch          = 0
rcoulomb                 = 0.8
; Relative dielectric constant for the medium and the reaction field
epsilon-r                = 1
epsilon_rf               = 1
; Method for doing Van der Waals
vdw-type                 = User
; cut-off lengths       
rvdw-switch              = 0
rvdw                     = 0.8
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = No
; Extension of the potential lookup tables beyond the cut-off
; Seperate tables between energy group pairs
energygrps               = LJ SOL
energygrp_table          = LJ SOL
table-extension          = 6.7
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx               = 0
fourier_ny               = 0
fourier_nz               = 0
; EWALD/PME/PPPM parameters
pme_order                = 4
ewald_rtol               = 1e-05
ewald_geometry           = 3d
epsilon_surface          = 0
optimize_fft             = no

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling  
tcoupl                   = berendsen  
; Groups to couple separately
tc-grps                  = System  
; Time constant (ps) and reference temperature (K)
tau-t                    = 0.1      
ref-t                    = 300.15          
; Pressure coupling     
Pcoupl                   = no 
Pcoupltype               = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p                    = 0.5
compressibility          = 4.5e-5
ref-p                    = 1
; Random seed for Andersen thermostat
andersen_seed            = 815131

; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel                  = yes
gen-temp                 = 300.15
gen-seed                 = 173529

; OPTIONS FOR BONDS    
constraints              = all-bonds
; Type of constraint algorithm
constraint-algorithm     = Lincs
; Do not constrain the start configuration
unconstrained-start      = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR                = no
; Relative tolerance of shake
shake-tol                = 0.0001
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter               = 4
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle          = 30
; Convert harmonic bonds to morse potentials
morse                    = no


(2) I have table.xvg and table_LJ_SOL.xvg in my local directory when using mdrun.
table_LJ_SOL.xvg was made with this fortran code:

      program gen_table
      implicit none
      real,parameter :: delr=0.002,rcut=15.00
      real :: r
      integer :: nbins,j
      nbins=int((rcut)/delr)+1
      do j=0,nbins
        r=delr*j
c       if(r.le.0.04.or.r.gt.0.8)then
        if(r.le.0.04)then 
          write(6,'(e16.5,6(2x,e16.10))')r,0.0000000000e+00,
     *                                    0.0000000000e+00,
     *                                    0.0000000000e+00,
     *                                    0.0000000000e+00,
     *                                    0.0000000000e+00,
     *                                    0.0000000000e+00
        else
          write(6,'(e16.5,6(2x,e16.10))')r,1/r,1/(r*r),
     *      -1/(r**6),-6/(r**7),
     *      1/(r**12),12/(r**13)
        end if
      end do
      stop
      end

table.xvg was also made with that code, except that "if(r.le.0.04)then" was commented out instead of "if(r.le.0.04.or.r.gt.0.8)then"


(3) Beginning of an example log file:
Log file opened on Wed Dec 18 10:51:26 2013
Host: compute-0-34.local  pid: 26844  nodeid: 0  nnodes:  1
Gromacs version:    VERSION 4.6
Precision:          single
Memory model:       64 bit
MPI library:        thread_mpi
OpenMP support:     enabled
GPU support:        disabled
invsqrt routine:    gmx_software_invsqrt(x)
CPU acceleration:   SSE2
FFT library:        fftw-3.3.3-sse2
Large file support: enabled
RDTSCP usage:       enabled
Built on:           Mon Mar 11 23:28:46 CDT 2013
Built by:           root at cluster.chem.ksu.edu [CMAKE]
Build OS/arch:      Linux 2.6.18-164.6.1.el5 x86_64
Build CPU vendor:   AuthenticAMD
Build CPU brand:    Six-Core AMD Opteron(tm) Processor 2427
Build CPU family:   16   Model: 8   Stepping: 0
Build CPU features: apic clfsh cmov cx8 cx16 htt lahf_lm misalignsse mmx msr nonstop_tsc pdpe1gb
 popcnt pse rdtscp sse2 sse3 sse4a
C compiler:         /raid/gromacs/bin/gcc GNU gcc (GCC) 4.7.2
C compiler flags:   -msse2   -Wextra -Wno-missing-field-initializers -Wno-sign-compare -Wall -Wn
o-unused -Wunused-value   -fomit-frame-pointer -funroll-all-loops -fexcess-precision=fast  -O3 -
DNDEBUG


                         :-)  G  R  O  M  A  C  S  (-:

           Glycine aRginine prOline Methionine Alanine Cystine Serine

                             :-)  VERSION 4.6  (-:

        Contributions from Mark Abraham, Emile Apol, Rossen Apostolov, 
           Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar,  
     Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans, 
        Peter Kasson, Carsten Kutzner, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
         Copyright (c) 2001-2012,2013, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
       modify it under the terms of the GNU Lesser General Public License
        as published by the Free Software Foundation; either version 2.1
             of the License, or (at your option) any later version.

                                :-)  mdrun  (-:


++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
B. Hess and C. Kutzner and D. van der Spoel and E. Lindahl
GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable
molecular simulation
J. Chem. Theory Comput. 4 (2008) pp. 435-447
-------- -------- --- Thank You --- -------- --------


++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C.
Berendsen
GROMACS: Fast, Flexible and Free
J. Comp. Chem. 26 (2005) pp. 1701-1719
-------- -------- --- Thank You --- -------- --------


++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
E. Lindahl and B. Hess and D. van der Spoel
GROMACS 3.0: A package for molecular simulation and trajectory analysis
J. Mol. Mod. 7 (2001) pp. 306-317
-------- -------- --- Thank You --- -------- --------


++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
H. J. C. Berendsen, D. van der Spoel and R. van Drunen
GROMACS: A message-passing parallel molecular dynamics implementation
Comp. Phys. Comm. 91 (1995) pp. 43-56
-------- -------- --- Thank You --- -------- --------

Input Parameters:
   integrator           = md
   nsteps               = 50000
   init-step            = 0
   cutoff-scheme        = Group
   ns_type              = Grid
   nstlist              = 1
   ndelta               = 2
   nstcomm              = 500
   comm-mode            = Linear
   nstlog               = 500
   nstxout              = 500
   nstvout              = 0
   nstfout              = 0
   nstcalcenergy        = 100
   nstenergy            = 500
   nstxtcout            = 0
   init-t               = 0
   delta-t              = 0.002
   xtcprec              = 1000
   fourierspacing       = 0.12
   nkx                  = 128
   nky                  = 128
   nkz                  = 128
   pme-order            = 4
   ewald-rtol           = 1e-05
   ewald-geometry       = 0
   epsilon-surface      = 0
   optimize-fft         = FALSE
   ePBC                 = xyz
   bPeriodicMols        = FALSE
   bContinuation        = FALSE
   bShakeSOR            = FALSE
   etc                  = Berendsen
   bPrintNHChains       = FALSE
   nsttcouple           = 1
   epc                  = No
   epctype              = Isotropic
   nstpcouple           = -1
   tau-p                = 0.5
   ref-p (3x3):
      ref-p[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      ref-p[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      ref-p[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
   compress (3x3):
      compress[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      compress[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      compress[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
   refcoord-scaling     = No
   posres-com (3):
      posres-com[0]= 0.00000e+00
      posres-com[1]= 0.00000e+00
      posres-com[2]= 0.00000e+00
   posres-comB (3):
      posres-comB[0]= 0.00000e+00
      posres-comB[1]= 0.00000e+00
      posres-comB[2]= 0.00000e+00
   verlet-buffer-drift  = 0.005
   rlist                = 0.8
   rlistlong            = 0.8
   nstcalclr            = 0
   rtpi                 = 0.05
   coulombtype          = PME
   coulomb-modifier     = None
   rcoulomb-switch      = 0
   rcoulomb             = 0.8
   vdwtype              = User
   vdw-modifier         = None
   rvdw-switch          = 0
   rvdw                 = 0.8
   epsilon-r            = 1
   epsilon-rf           = 1
   tabext               = 6.7
   implicit-solvent     = No
   gb-algorithm         = Still
   gb-epsilon-solvent   = 80
   nstgbradii           = 1
   rgbradii             = 2
   gb-saltconc          = 0
   gb-obc-alpha         = 1
   gb-obc-beta          = 0.8
   gb-obc-gamma         = 4.85
   gb-dielectric-offset = 0.009
   sa-algorithm         = Ace-approximation
   sa-surface-tension   = 2.05016
   DispCorr             = No
   bSimTemp             = FALSE
   free-energy          = no
   nwall                = 0
   wall-type            = 9-3
   wall-atomtype[0]     = -1
   wall-atomtype[1]     = -1
   wall-density[0]      = 0
   wall-density[1]      = 0
   wall-ewald-zfac      = 3
   pull                 = no
   rotation             = FALSE
   disre                = No
   disre-weighting      = Conservative
   disre-mixed          = FALSE
   dr-fc                = 1000
   dr-tau               = 0
   nstdisreout          = 100
   orires-fc            = 0
   orires-tau           = 0
   nstorireout          = 100
   dihre-fc             = 0
   em-stepsize          = 0.01
   em-tol               = 10
   niter                = 20
   fc-stepsize          = 0
   nstcgsteep           = 1000
   nbfgscorr            = 10
   ConstAlg             = Lincs
   shake-tol            = 0.0001
   lincs-order          = 4
   lincs-warnangle      = 30
   lincs-iter           = 4
   bd-fric              = 0
   ld-seed              = 1993
   cos-accel            = 0
   deform (3x3):
      deform[    0]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      deform[    1]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
      deform[    2]={ 0.00000e+00,  0.00000e+00,  0.00000e+00}
   adress               = FALSE
   userint1             = 0
   userint2             = 0
   userint3             = 0
   userint4             = 0
   userreal1            = 0
   userreal2            = 0
   userreal3            = 0
   userreal4            = 0
grpopts:
   nrdf:      521885
   ref-t:      300.15
   tau-t:         0.1
anneal:          No
ann-npoints:           0
   acc:	           0           0           0
   nfreeze:           N           N           N
   energygrp-flags[  0]: 0 2
   energygrp-flags[  1]: 2 0
   efield-x:
      n = 0
   efield-xt:
      n = 0
   efield-y:
      n = 0
   efield-yt:
      n = 0
   efield-z:
      n = 0
   efield-zt:
      n = 0
   bQMMM                = FALSE
   QMconstraints        = 0
   QMMMscheme           = 0
   scalefactor          = 1
qm-opts:
   ngQM                 = 0
Using 1 MPI thread

Detecting CPU-specific acceleration.
Present hardware specification:
Vendor: AuthenticAMD
Brand:  Six-Core AMD Opteron(tm) Processor 2427
Family: 16  Model:  8  Stepping:  0
Features: apic clfsh cmov cx8 cx16 htt lahf_lm misalignsse mmx msr nonstop_tsc pdpe1gb popcnt ps
e rdtscp sse2 sse3 sse4a
Acceleration most likely to fit this hardware: SSE2
Acceleration selected at GROMACS compile time: SSE2

Table routines are used for coulomb: FALSE
Table routines are used for vdw:     TRUE
Will do PME sum in reciprocal space.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen 
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- --- Thank You --- -------- --------

Will do ordinary reciprocal space Ewald sum.
Using a Gaussian width (1/beta) of 0.25613 nm for Ewald
Cut-off's:   NS: 0.8   Coulomb: 0.8   LJ: 0.8
System total charge: 0.000
Read user tables from table.xvg with 7501 data points.
Tabscale = 500 points/nm
Generated table with 3750 data points for Ewald.
Tabscale = 500 points/nm
Read user tables from table_LJ_SOL.xvg with 7501 data points.
Tabscale = 500 points/nm
Generated table with 3750 data points for Ewald.
Tabscale = 500 points/nm

Enabling SPC-like water optimization for 74555 molecules.

Potential shift: LJ r^-12: 0.000 r^-6 0.000, Ewald 0.000e+00
Initialized non-bonded Ewald correction tables, spacing: 5.90e-04 size: 25421

Removing pbc first time

Initializing LINear Constraint Solver

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
B. Hess and H. Bekker and H. J. C. Berendsen and J. G. E. M. Fraaije
LINCS: A Linear Constraint Solver for molecular simulations
J. Comp. Chem. 18 (1997) pp. 1463-1472
-------- -------- --- Thank You --- -------- --------

The number of constraints is 149110
Center of mass motion removal mode is Linear
We have the following groups for center of mass motion removal:
  0:  rest

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
H. J. C. Berendsen, J. P. M. Postma, A. DiNola and J. R. Haak
Molecular dynamics with coupling to an external bath
J. Chem. Phys. 81 (1984) pp. 3684-3690
-------- -------- --- Thank You --- -------- --------

There are: 223666 Atoms
Max number of connections per atom is 4
Total number of connections is 596440
Max number of graph edges per atom is 2
Total number of graph edges is 298220

Constraining the starting coordinates (step 0)

Constraining the coordinates at t0-dt (step 0)
RMS relative constraint deviation after constraining: 2.24e-06
Initial temperature: 300.862 K

Started mdrun on node 0 Wed Dec 18 10:51:28 2013

           Step           Time         Lambda
              0        0.00000        0.00000


Grid: 22 x 22 x 16 cells
   Energies (kJ/mol)
          Angle        LJ (SR)   Coulomb (SR)   Coul. recip.      Potential
    6.00247e+04    6.43606e+05   -3.87278e+06   -5.92226e+05   -3.76137e+06
    Kinetic En.   Total Energy    Temperature Pressure (bar)   Constr. rmsd
    6.52742e+05   -3.10863e+06    3.00857e+02   -3.66814e+03    1.82007e-06

...etc...


More information about the gromacs.org_gmx-users mailing list