[gmx-users] Gromacs 4.6 & 4.5.3 qualitative differences & 4.6 instability in polarizable force field vacuum/liquid mixture interface simulations
ploetz
ploetz at ksu.edu
Sat Nov 2 18:38:37 CET 2013
Dear Gromacs Users,
I am trying to simulate a system consisting of a vacuum/condensed phase
interface in which a 6x6x12nm condensed phase region is flanked on both ends
(in the z-dimension) by a 6x6x12nm vacuum region to form overall box
dimensions of 6x6x36 nm. The system is a binary liquid mixture of methanol
(0.125 mole fraction methanol) in water using a polarizable (charge on a
spring) force field (COS/M methanol and COS/G2 water) at 300K and 1bar. The
system is stable in Gromacs 4.5.3; however, mdrun gives a segmentation fault
in Gromacs 4.6 when attempting to do dynamics (energy minimization completes
with no apparent problems). If I remove the vacuum region, mdrun works. If I
incrementally add 2 Angstroms to the z-dimension until I reached a vacuum
region of 34 Angstroms total (17 Angstroms on both sides of the condensed
phase region) and try to simulate these systems, mdrun works every time.
When I reach 36 Angstroms, the segmentation fault re-appears. Although not
the system I am actually interested in, I did some simulations using Gromacs
4.6 with the 34 Angstrom vacuum region system and observed an undulating and
very turbulant "vacuum"/condensed phase interface in which a column of
water/methanol mixture came out of the condensed phase region to connect the
two interfaces. Also, the center of mass motion of the system appeared to
not have been removed. In contrast, using Gromacs 4.5.3, the interface is
not undulating, but is "calm" and qualitatively planar, no column forms to
connect the interfaces, and there is no problem with the center of mass
motion removal. Some of the molecules do enter the "vacuum" region when
running with 4.5.3, but this appears to be due to the movement of individual
molecules, not a collective motion of many molecules. This system runs fine
with a non-polarizable force field in Gromacs 4.6.
I have also compared several properties (using g_energy) of the bulk system
(no vacuum region) using Gromacs 4.6 and 4.5.3 and they are not the same for
the polarizable force field, but they are the same for the non-polarizable
force field. Specifically, with the polarizable force field, the LJ(SR)
energy is more positive with 4.6, the LJ(LR) energy is more negative with
4.6, the Coulomb(SR) energy is more negative with 4.6, the Could. recip.
energy is more negative with 4.6, the polarization energy is more positive
with 4.6, the potential energy is more negative in 4.6, the average kinetic
energy (and temperature) is the same but the fluctuations are greater in
4.6, the total energy is more negative in 4.6, the pressure looks fine, and
the volume looks fine. Where I've noted differences, these are all
statistically significant differences.
I would like to know if I can just use 4.5.3 and assume the differences
between the results of 4.5.3 and 4.6 are due to some problem in 4.6. I am
using all the same input files and commands with both versions, only
different executables. I ran the regression tests for 4.6 when I installed
it, and passed them all.
Sincerely,
Elizabeth
-----
ADDITIONAL DETAILS:
Below is where the problem appears for the interface system when z-dimension
of the vacuum region is >= 36 Angstroms total. eq4 is my first attempt at
dynamics, after three successful energy minimizations (1st: charges screened
and no bond constraints, 2nd: charges felt but no bond constraints, 3rd:
charges felt and bond constraints on)
[ploetz at cluster AddRemainingVacuumBack]$ grompp -f eq4.mdp -c em3.pdb -o
eq4.tpr -n index.ndx -p sys.top -nice 0
:-) G R O M A C S (-:
Green Red Orange Magenta Azure Cyan Skyblue
:-) 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.
:-) grompp (-:
Option Filename Type Description
------------------------------------------------------------
-f eq4.mdp Input grompp input file with MD parameters
-po mdout.mdp Output grompp input file with MD parameters
-c em3.pdb Input Structure file: gro g96 pdb tpr etc.
-r conf.gro Input, Opt. Structure file: gro g96 pdb tpr etc.
-rb conf.gro Input, Opt. Structure file: gro g96 pdb tpr etc.
-n index.ndx Input, Opt! Index file
-p sys.top Input Topology file
-pp processed.top Output, Opt. Topology file
-o eq4.tpr Output Run input file: tpr tpb tpa
-t traj.trr Input, Opt. Full precision trajectory: trr trj cpt
-e ener.edr Input, Opt. Energy file
-ref rotref.trr In/Out, Opt. Full precision trajectory: trr trj cpt
Option Type Value Description
------------------------------------------------------
-[no]h bool no Print help info and quit
-[no]version bool no Print version info and quit
-nice int 0 Set the nicelevel
-[no]v bool no Be loud and noisy
-time real -1 Take frame at or first after this time.
-[no]rmvsbds bool yes Remove constant bonded interactions with virtual
sites
-maxwarn int 0 Number of allowed warnings during input
processing. Not for normal use and may generate
unstable systems
-[no]zero bool no Set parameters for bonded interactions without
defaults to zero instead of generating an error
-[no]renum bool yes Renumber atomtypes and minimize number of
atomtypes
Ignoring obsolete mdp entry 'cpp'
Ignoring obsolete mdp entry 'domain-decomposition'
Ignoring obsolete mdp entry 'andersen_seed'
Ignoring obsolete mdp entry 'nstcheckpoint'
Replacing old mdp entry 'unconstrained-start' by 'continuation'
NOTE 1 [file eq4.mdp]:
The Berendsen thermostat does not generate the correct kinetic energy
distribution. You might want to consider using the V-rescale thermostat.
Generated 80 of the 91 non-bonded parameter combinations
Excluding 2 bonded neighbours molecule type 'COSM'
turning all bonds into constraints...
Excluding 2 bonded neighbours molecule type 'COSG2'
turning all bonds into constraints...
Velocities were taken from a Maxwell distribution at 300.15 K
Cleaning up constraints and constant bonded interactions with virtual sites
Number of degrees of freedom in T-Coupling group MOH is 8243.62
Number of degrees of freedom in T-Coupling group SOL is 57717.38
Largest charge group radii for Van der Waals: 0.118, 0.118 nm
Largest charge group radii for Coulomb: 0.118, 0.118 nm
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 48x48x288, spacing 0.120 0.120 0.120
Estimate for the relative computational load of the PME mesh part: 0.37
This run will generate roughly 622 Mb of data
There was 1 note
gcq#228: "Bailed Out Of Edge Synchronization After 10,000 Iterations"
(X/Motif)
[ploetz at cluster AddRemainingVacuumBack]$ more eq4.mdp
; VARIOUS PREPROCESSING OPTIONS
title =
; Preprocessor - specify a full path if necessary.
cpp = /lib/cpp
include =
define =
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 500000
; 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 =
; Selection of energy groups
energygrps = MOH SOL
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 10
; 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 = 1.0
domain-decomposition = no
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = PME
rcoulomb-switch = 0
rcoulomb = 1.0
; Relative dielectric constant for the medium and the reaction field
epsilon-r = 1
epsilon_rf = 1
; Method for doing Van der Waals
vdw-type = Cut-off
; cut-off lengths
rvdw-switch = 0
rvdw = 1.5
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = No
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; Seperate tables between energy group pairs
energygrp_table =
; 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
; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent = No
; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = berendsen
; Groups to couple separately
tc-grps = MOH SOL
; Time constant (ps) and reference temperature (K)
tau-t = 0.1 0.1
ref-t = 300.15 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
; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl =
[ploetz at cluster AddRemainingVacuumBack]$
[ploetz at cluster AddRemainingVacuumBack]$ mdrun -v -s eq4.tpr -c eq4.pdb -g
eq4.log -o eq4.trr -nice 0
:-) G R O M A C S (-:
GRoups of Organic Molecules in ACtion for Science
:-) 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 (-:
Option Filename Type Description
------------------------------------------------------------
-s eq4.tpr Input Run input file: tpr tpb tpa
-o eq4.trr Output Full precision trajectory: trr trj cpt
-x traj.xtc Output, Opt. Compressed trajectory (portable xdr
format)
-cpi state.cpt Input, Opt. Checkpoint file
-cpo state.cpt Output, Opt. Checkpoint file
-c eq4.pdb Output Structure file: gro g96 pdb etc.
-e ener.edr Output Energy file
-g eq4.log Output Log file
-dhdl dhdl.xvg Output, Opt. xvgr/xmgr file
-field field.xvg Output, Opt. xvgr/xmgr file
-table table.xvg Input, Opt. xvgr/xmgr file
-tabletf tabletf.xvg Input, Opt. xvgr/xmgr file
-tablep tablep.xvg Input, Opt. xvgr/xmgr file
-tableb table.xvg Input, Opt. xvgr/xmgr file
-rerun rerun.xtc Input, Opt. Trajectory: xtc trr trj gro g96 pdb cpt
-tpi tpi.xvg Output, Opt. xvgr/xmgr file
-tpid tpidist.xvg Output, Opt. xvgr/xmgr file
-ei sam.edi Input, Opt. ED sampling input
-eo edsam.xvg Output, Opt. xvgr/xmgr file
-j wham.gct Input, Opt. General coupling stuff
-jo bam.gct Output, Opt. General coupling stuff
-ffout gct.xvg Output, Opt. xvgr/xmgr file
-devout deviatie.xvg Output, Opt. xvgr/xmgr file
-runav runaver.xvg Output, Opt. xvgr/xmgr file
-px pullx.xvg Output, Opt. xvgr/xmgr file
-pf pullf.xvg Output, Opt. xvgr/xmgr file
-ro rotation.xvg Output, Opt. xvgr/xmgr file
-ra rotangles.log Output, Opt. Log file
-rs rotslabs.log Output, Opt. Log file
-rt rottorque.log Output, Opt. Log file
-mtx nm.mtx Output, Opt. Hessian matrix
-dn dipole.ndx Output, Opt. Index file
-multidir rundir Input, Opt., Mult. Run directory
-membed membed.dat Input, Opt. Generic data file
-mp membed.top Input, Opt. Topology file
-mn membed.ndx Input, Opt. Index file
Option Type Value Description
------------------------------------------------------
-[no]h bool no Print help info and quit
-[no]version bool no Print version info and quit
-nice int 0 Set the nicelevel
-deffnm string Set the default filename for all file options
-xvg enum xmgrace xvg plot formatting: xmgrace, xmgr or none
-[no]pd bool no Use particle decompostion
-dd vector 0 0 0 Domain decomposition grid, 0 is optimize
-ddorder enum interleave DD node order: interleave, pp_pme or
cartesian
-npme int -1 Number of separate nodes to be used for PME, -1
is guess
-nt int 0 Total number of threads to start (0 is guess)
-ntmpi int 0 Number of thread-MPI threads to start (0 is
guess)
-ntomp int 0 Number of OpenMP threads per MPI process/thread
to start (0 is guess)
-ntomp_pme int 0 Number of OpenMP threads per MPI process/thread
to start (0 is -ntomp)
-pin enum auto Fix threads (or processes) to specific cores:
auto, on or off
-pinoffset int 0 The starting logical core number for pinning to
cores; used to avoid pinning threads from
different mdrun instances to the same core
-pinstride int 0 Pinning distance in logical cores for threads,
use 0 to minimize the number of threads per
physical core
-gpu_id string List of GPU id's to use
-[no]ddcheck bool yes Check for all bonded interactions with DD
-rdd real 0 The maximum distance for bonded interactions
with
DD (nm), 0 is determine from initial coordinates
-rcon real 0 Maximum distance for P-LINCS (nm), 0 is estimate
-dlb enum auto Dynamic load balancing (with DD): auto, no or
yes
-dds real 0.8 Minimum allowed dlb scaling of the DD cell size
-gcom int -1 Global communication frequency
-nb enum auto Calculate non-bonded interactions on: auto, cpu,
gpu or gpu_cpu
-[no]tunepme bool yes Optimize PME load between PP/PME nodes or
GPU/CPU
-[no]testverlet bool no Test the Verlet non-bonded scheme
-[no]v bool yes Be loud and noisy
-[no]compact bool yes Write a compact log file
-[no]seppot bool no Write separate V and dVdl terms for each
interaction type and node to the log file(s)
-pforce real -1 Print all forces larger than this (kJ/mol nm)
-[no]reprod bool no Try to avoid optimizations that affect binary
reproducibility
-cpt real 15 Checkpoint interval (minutes)
-[no]cpnum bool no Keep and number checkpoint files
-[no]append bool yes Append to previous output files when continuing
from checkpoint instead of adding the simulation
part number to all file names
-nsteps int -2 Run this number of steps, overrides .mdp file
option
-maxh real -1 Terminate after 0.99 times this time (hours)
-multi int 0 Do multiple simulations in parallel
-replex int 0 Attempt replica exchange periodically with this
period (steps)
-nex int 0 Number of random exchanges to carry out each
exchange interval (N^3 is one suggestion). -nex
zero or not specified gives neighbor replica
exchange.
-reseed int -1 Seed for replica exchange, -1 is generate a seed
-[no]ionize bool no Do a simulation including the effect of an X-Ray
bombardment on your system
Reading file eq4.tpr, VERSION 4.6 (single precision)
Using 12 MPI threads
starting mdrun 'COSM in water x_COSM 0.125'
500000 steps, 1000.0 ps.
Segmentation fault
[ploetz at cluster AddRemainingVacuumBack]$
[ploetz at cluster AddRemainingVacuumBack]$ more sys.top
#include "polff.itp"
#include "cosm.itp"
#include "cosg2.itp"
[ system ]
; Name
COSM in water x_COSM 0.125
[ molecules ]
; Compound #mols
COSM 1374
COSG2 9620
[ploetz at cluster AddRemainingVacuumBack]$ more polff.itp
;Polarizable force field atom types
;Includes types for COS/G2 (water) COS/M (MeOH 1-site) and CPC (MeOH 2-site)
[ defaults ]
LJ Geometric
[ atomtypes ]
;name mass charge ptype c6 c12
WO 15.99940 0.0 A 0.0 0.0
WH 1.00800 0.0 A 0.0 0.0
WD 0.00000 0.0 D 0.0 0.0
CSMO 15.99940 0.0 A 0.0 0.0
CSMH 1.00800 0.0 A 0.0 0.0
CSMC 15.03500 0.0 A 0.0 0.0
CPCO 15.99940 0.0 A 0.0 0.0
CPCH 1.00800 0.0 A 0.0 0.0
CPCC 15.03500 0.0 A 0.0 0.0
WS 0.00000 0.0 S 0.0 0.0
CSMS 0.00000 0.0 S 0.0 0.0
CPCSC 0.00000 0.0 S 0.0 0.0
CPCSO 0.00000 0.0 S 0.0 0.0
[ nonbond_params ]
;NB: mix of COSM and CPC not included
;NB: includes special COSMC - COS/G2 C12 parameters
WO WO 1 3.24434e-03 3.45765e-06
WO CSMO 1 2.71125e-03 3.00305e-06
WO CSMC 1 5.36612e-03 7.71936e-06
WO CPCO 1 2.68847e-03 2.74273e-06
WO CPCC 1 6.35721e-03 10.57112e-06
CSMO CSMO 1 2.26576e-03 2.60823e-06
CSMO CSMC 1 4.48440e-03 7.10600e-06
CSMC CSMC 1 8.87552e-03 19.36000e-06
CPCO CPCO 1 2.22784e-03 2.17563e-06
CPCO CPCC 1 5.26799e-03 8.38538e-06
CPCC CPCC 1 12.45679e-03 32.31923e-06
[ploetz at cluster AddRemainingVacuumBack]$
[ploetz at cluster AddRemainingVacuumBack]$ more cosm.itp
; Topology for COS/M methanol model
; Yu et al. J., Comput. Chem. 27 (1494-1504), 2006
[ moleculetype ]
; molname nrexcl
COSM 2
[ atoms ]
; id at type res nr residu name at name cg nr charge
1 CSMO 1 MOH OM 1 7.470
2 CSMH 1 MOH HM 1 0.360
3 CSMC 1 MOH CM 1 0.170
4 CSMS 1 MOH SOM 1 -8.000
[ polarization ]
;Center Shell funct alpha (nm^3)
1 4 1 0.001320
[ constraints ]
; ai aj funct length
1 2 2 0.1000
1 3 2 0.1430
2 3 2 0.1988
[ exclusions ]
; iatom excluded from interaction with i
1 2 3 4
2 1 3 4
3 1 2 4
4 1 2 3
#ifdef POSRES
[ position_restraints ]
; iatom type fx fy fz
1 1 100 100 100
3 1 100 100 100
#endif
[ploetz at cluster AddRemainingVacuumBack]$
[ploetz at cluster AddRemainingVacuumBack]$ more cosg2.itp
; Topology for COS/G2 water model
; Yu and Van Gunsteren, J. Chem. Phys. 121 (9549-9564), 2004
[ moleculetype ]
; molname nrexcl
COSG2 2
[ atoms ]
; id at type res nr residu name at name cg nr charge
1 WO 1 SOL OW 1 0.0000
2 WH 1 SOL HW1 1 0.5265
3 WH 1 SOL HW2 1 0.5265
4 WD 1 SOL MW 1 6.9470
5 WS 1 SOL SW 1 -8.0000
[ polarization ]
;Center Shell funct alpha (nm^3)
4 5 1 0.001255
[ settles ]
; i funct dOH dHH
1 1 0.09572 0.15139
[ dummies3 ]
; The position of the dummies is computed as follows:
;
; O
;
; D
;
; H H
;
; 2 * b = distance (OD) / [ cos (angle(DOH)) * distance (OH) ]
; 0.022 nm / [ cos (104.52 / 2 deg) * 0.09572 nm ]
; 0.022 nm / 0.058588228 nm
; Dummy pos x4 = x1 + a*(x2-x1) + b*(x3-X1)
; Dummy from funct a b
4 1 2 3 1 0.187751028 0.187751028
[ exclusions ]
; iatom excluded from interaction with i
1 2 3 4 5
2 1 3 4 5
3 1 2 4 5
4 1 2 3 5
5 1 2 3 4
#ifdef POSRES
[ position_restraints ]
; iatom type fx fy fz
1 1 100 100 100
#endif
[ploetz at cluster AddRemainingVacuumBack]$
--
View this message in context: http://gromacs.5086.x6.nabble.com/Gromacs-4-6-4-5-3-qualitative-differences-4-6-instability-in-polarizable-force-field-vacuum-liquid-ms-tp5012174.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
More information about the gromacs.org_gmx-users
mailing list