[gmx-users] Simualtion blows up with Verlet cutoff scheme

Bu Wang bw2 at alfred.edu
Thu Mar 14 03:02:55 CET 2013


Hello,

We are having a problem with the non-bond energy calculation in MD
simulations using Gromacs 6.4.1. Our system is a simple alumina crystal and
the non-bond interactions are defined with our own buckingham potentials.
Everything works fine when the group cutoff scheme is used. But with
similar settings in Verlet scheme, the simulation would blow up right after
step 0 showing absurdly high buckingham energy and kinetic energy.

I've tried modifying parameters related to the Verlet scheme, such as
verlet-buffer-drift, vdw-modifier and rvdw-switch, but with no luck. I also
compared the outputs from group and Verlet schemes but couldn't see
anything weird. I suspect that we made some mistakes in setting up the
topology, but couldn't figure out what might be wrong. I attached the
topology and pasted below the parameters copied from the output. I would
greatly appreciate any suggestions. Thanks in advance!

Input Parameters:
   integrator           = md
   nsteps               = 1000000
   init-step            = 0
   cutoff-scheme        = Verlet
   ns_type              = Grid
   nstlist              = 10
   ndelta               = 2
   nstcomm              = 500
   comm-mode            = Linear
   nstlog               = 500
   nstxout              = 0
   nstvout              = 1000
   nstfout              = 0
   nstcalcenergy        = 100
   nstenergy            = 500
   nstxtcout            = 1000
   init-t               = 0
   delta-t              = 0.0005
   xtcprec              = 1e+08
   fourierspacing       = 0.12
   nkx                  = 32
   nky                  = 32
   nkz                  = 32
   pme-order            = 6
   ewald-rtol           = 1e-05
   ewald-geometry       = 0
   epsilon-surface      = 0
   optimize-fft         = TRUE
   ePBC                 = xyz
   bPeriodicMols        = FALSE
   bContinuation        = FALSE
   bShakeSOR            = FALSE
   etc                  = V-rescale
   bPrintNHChains       = FALSE
   nsttcouple           = 1
   epc                  = Berendsen
   epctype              = Isotropic
   nstpcouple           = 1
   tau-p                = 0.01
   ref-p (3x3):
      ref-p[    0]={ 1.00000e+00,  0.00000e+00,  0.00000e+00}
      ref-p[    1]={ 0.00000e+00,  1.00000e+00,  0.00000e+00}
      ref-p[    2]={ 0.00000e+00,  0.00000e+00,  1.00000e+00}
   compress (3x3):
      compress[    0]={ 5.00000e-06,  0.00000e+00,  0.00000e+00}
      compress[    1]={ 0.00000e+00,  5.00000e-06,  0.00000e+00}
      compress[    2]={ 0.00000e+00,  0.00000e+00,  5.00000e-06}
   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                = 1.219
   rlistlong            = 1.219
   nstcalclr            = 1
   rtpi                 = 0.05
   coulombtype          = PME
   coulomb-modifier     = None
   rcoulomb-switch      = 0
   rcoulomb             = 1.2
   vdwtype              = Cut-off
   vdw-modifier         = Potential-shift
   rvdw-switch          = 1.1
   rvdw                 = 1.2
   epsilon-r            = 1
   epsilon-rf           = inf
   tabext               = 1
   implicit-solvent     = No
   gb-algorithm         = Still
   gb-epsilon-solvent   = 80
   nstgbradii           = 1
   rgbradii             = 1
   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           = 1
   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:       10959
   ref-t:         700
   tau-t:         0.1
anneal:          No
ann-npoints:           0
   acc:           0           0           0
   nfreeze:           N           N           N
   energygrp-flags[  0]: 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
Using 1 OpenMP thread

Detecting CPU-specific acceleration.
Present hardware specification:
Vendor: GenuineIntel
Brand:  Intel(R) Xeon(R) CPU           E5645  @ 2.40GHz
Family:  6  Model: 44  Stepping:  2
Features: aes apic clfsh cmov cx8 cx16 htt lahf_lm mmx msr nonstop_tsc pcid
pclmuldq pdcm pdpe1gb popcnt pse rdtscp sse2 sse3 sse4.1 sse4.2 ssse3
Acceleration most likely to fit this hardware: SSE4.1
Acceleration selected at GROMACS compile time: SSE4.1

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.384195 nm for Ewald
Cut-off's:   NS: 1.219   Coulomb: 1.2   BHAM: 1.2
Determining largest Buckingham b parameter for table
Buckingham b parameters, min: 0, max: 45.6204
System total charge: 0.000
Generated table with 4438 data points for Ewald.
Tabscale = 2000 points/nm
Generated table with 4438 data points for LJ6.
Tabscale = 2000 points/nm
Generated table with 4438 data points for EXPMIN.
Tabscale = 43.84 points/nm

Using SSE4.1 4x2 non-bonded kernels

Using full Lennard-Jones parameter combination matrix

Potential shift: LJ r^-12: 0.112 r^-6 0.335, Ewald 0.000e+00
Initialized non-bonded Ewald correction tables, spacing: 7.23e-04 size: 3072

Removing pbc first time
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 ++++
G. Bussi, D. Donadio and M. Parrinello
Canonical sampling through velocity rescaling
J. Chem. Phys. 126 (2007) pp. 014101
-------- -------- --- Thank You --- -------- --------

There are: 3654 Atoms
Initial temperature: 699.927 K

Started mdrun on node 0 Wed Mar 13 18:32:37 2013

           Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
  Buck.ham (SR)   Coulomb (SR)   Coul. recip.      Potential    Kinetic En.
    5.03375e+11   -1.27701e+07    2.74229e+03    5.03363e+11    3.10221e+15
   Total Energy    Temperature Pressure (bar)
    3.10272e+15    6.80916e+13    8.55442e+14

--
Bu


More information about the gromacs.org_gmx-users mailing list