Project

General

Profile

Bug #1376

Gromacs 4.6 & 4.5.3 qualitative differences & 4.6 instability in polarizable force field vacuum/liquid mixture interface simulations

Added by Elizabeth Ploetz almost 4 years ago. Updated over 3 years ago.

Status:
Closed
Priority:
Normal
Category:
mdrun
Target version:
-
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

Hi Professor van der Spoel,

I'll try to think of ways to simplify the system. In the mean time, here is what I posted to the gmx-user group.

I've added the following files:
cosg2.itp
cosm.itp
em3.pdb
eq4-4.5.3.edr
eq4-4.5.3.log
eq4-4.5.3.tpr
eq4-4.5.3.trr
eq4-4.6.edr
eq4-4.6.log
eq4-4.6.tpr
eq4-4.6.trr
eq4.mdp
index.ndx
polff.itp
sys.top
in Redmine.tar

which were generated from
grompp -f eq4.mdp -c em3.pdb -o eq4-4.6.tpr -n index.ndx -p sys.top -nice 0
mdrun -v -s eq4-4.6.tpr -c eq4-4.6.pdb -g eq4-4.6.log -o eq4-4.6.trr -e eq4-4.6.edr -nice 0
or
/raid/ploetz/gromacs-4.5.3/bin/grompp -f eq4.mdp -c em3.pdb -o eq4-4.5.3.tpr -n index.ndx -p sys.top -nice 0
/raid/ploetz/gromacs-4.5.3/bin/mdrun -v -s eq4-4.5.3.tpr -c eq4-4.5.3.pdb -g eq4-4.5.3.log -o eq4-4.5.3.trr -e eq4-4.5.3.edr -nice 0

Elizabeth

On 2013-11-02 18:38, ploetz wrote:

Dear Gromacs Users,

Please start a redmine.gromacs.org issue and assign it to me, but try to
simplify the system as much as possible. You can cut and paste all the
information to the redmine issue.

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
energy is more positive with 4.6, the LJ 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@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
------------------------------------------------------
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@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@cluster AddRemainingVacuumBack]$

[ploetz@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
------------------------------------------------------
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
testverlet bool no Test the Verlet non-bonded scheme
[no]v bool yes Be loud and noisy
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
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@cluster AddRemainingVacuumBack]$

[ploetz@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@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@cluster AddRemainingVacuumBack]$

[ploetz@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@cluster AddRemainingVacuumBack]$
[ploetz@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@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.

--
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
[hidden email] http://folding.bmc.uu.se
--
gmx-users mailing list [hidden email]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Redmine.tar (8.91 MB) Redmine.tar Elizabeth Ploetz, 11/02/2013 10:33 PM

Associated revisions

Revision 34eddd6e (diff)
Added by David van der Spoel over 3 years ago

Implemented fatal error for nstcalcenergy!=1 with shells.

The combination of using shell particles with a value for
nstcalcenergy that is different from 1 leads to energies
not being communicated around on a parallel machine and hence
a meaningless shell relaxation procedure. This patch will just
bail out with a fatal error if the combination is detected.
Fixes #1376.

Change-Id: Ia54f1480b2358355438acb5de31dfa56b7f8d603

Revision 236439e4 (diff)
Added by David van der Spoel over 3 years ago

Added check in grompp for shells and inputrec issues.

Some combinations of inputrec settings do not work with shells,
in particular nstcalcenergy > 1 or use of a twin range cutoff.
This is now checked for in grompp. Fixes #1376.

Change-Id: I4382bcf5231920b22c725a7cf8e7b4e17c2526d9

Revision 69837f7c (diff)
Added by David van der Spoel over 3 years ago

Added check for domain decomposition and shells.

Fixes #1376 in as much as that it prevents incorrect calculations.
Can be left out when merging to 5.0 that does not support
particle decomposition. In 5.0 we can write a message about OpenMP
if that works with shells.
Change-Id: Ic07ae70cd3cd5298b6b2d579719da4d267b2d8fb

History

#1 Updated by Erik Lindahl over 3 years ago

David, feedback?

#2 Updated by David van der Spoel over 3 years ago

Both 4.5 and 4.6 give me unstable simulations, in 4.5:

MDStep= 0/ 0 EPot: -3.73548875e+05, rmsF: 1.19e+04
MDStep= 0/ 1 EPot: -5.07857375e+05, rmsF: 2.74e+03
MDStep= 0/ 2 EPot: -5.15182875e+05, rmsF: 5.69e+02
MDStep= 0/ 3 EPot: -5.15478812e+05, rmsF: 1.47e+02
MDStep= 0/ 4 EPot: -5.15493375e+05, rmsF: 5.74e+01
MDStep= 0/ 5 EPot: -5.15493281e+05, rmsF: 2.31e+01
MDStep= 0/ 6 EPot: -5.15489750e+05, rmsF: 9.04e+00

^Mstep 0MDStep= 1/ 0 EPot: -4.40416465e+21, rmsF: 4.47e+02
MDStep= 1/ 1 EPot: -4.40416465e+21, rmsF: 1.04e+02
MDStep= 1/ 2 EPot: -4.40416465e+21, rmsF: 2.80e+01
MDStep= 1/ 3 EPot: -4.40416465e+21, rmsF: 8.17e+00
MDStep= 2/ 0 EPot: -4.40416465e+21, rmsF: 4.23e+02
MDStep= 2/ 1 EPot: -4.40416465e+21, rmsF: 9.63e+01
MDStep= 2/ 2 EPot: -4.40416465e+21, rmsF: 2.59e+01
MDStep= 2/ 3 EPot: -4.40416465e+21, rmsF: 7.88e+00

4.6:

MDStep= 0/ 0 EPot: -3.73551031e+05, rmsF: 1.19e+04
MDStep= 0/ 1 EPot: -5.09971531e+05, rmsF: 2.74e+03
MDStep= 0/ 2 EPot: -5.19412812e+05, rmsF: 5.69e+02
MDStep= 0/ 3 EPot: -5.21814094e+05, rmsF: 1.47e+02
MDStep= 0/ 4 EPot: -5.23941594e+05, rmsF: 5.75e+01
MDStep= 0/ 5 EPot: -5.26050562e+05, rmsF: 2.32e+01
MDStep= 0/ 6 EPot: -5.28164188e+05, rmsF: 9.08e+00

step 0MDStep= 1/ 0 EPot: 1.12479906e+05, rmsF: 4.47e+02
MDStep= 1/ 1 EPot: 1.12160562e+05, rmsF: 1.04e+02
MDStep= 1/ 2 EPot: 1.12043656e+05, rmsF: 2.79e+01
MDStep= 1/ 3 EPot: 1.12031992e+05, rmsF: 8.18e+00
MDStep= 2/ 0 EPot: 1.12095203e+05, rmsF: 4.23e+02

etc.

#3 Updated by David van der Spoel over 3 years ago

It also seems that the output that was uploaded does not correspond to the input, or that there was an issue in gromacs 4.5.3. The system was set up with a twin-range cut off (1.0/1.5) but no long-range LJ is reported in the 4.5.3.log file.

The culprit seems to be nstcalcenergy. This is an evil option. It should be default 1.

Using 4.6 and nstcalcenergy = nstlist = 1 I get a stable simulation. The nstlist = 1 is needed because of the twin range cutoff. However, using rvdw = 1.0 (turning off twin range) we can put back nstlist to 10. Nstlist = 1 is bad, but not as bad as in unpolarizable MD as the neighborlist remains the same between iterations.

In summary this is not a real bug, but we should detect the combination of shells and nstcalcenergy != 1.

#4 Updated by Erik Lindahl over 3 years ago

David, did you add the detection you wrote about?

#5 Updated by Gerrit Code Review Bot over 3 years ago

Gerrit received a related patchset '1' for Issue #1376.
Uploader: David van der Spoel ()
Change-Id: Ia54f1480b2358355438acb5de31dfa56b7f8d603
Gerrit URL: https://gerrit.gromacs.org/3654

#6 Updated by Justin Lemkul over 3 years ago

The fatal error at run-time is the only solution at the moment, though it would be better if grompp could detect this instead. My Drude branch can do this because of the new .mdp options I have added. I am OK with the fix for now, but if we leave this issue open, I can make changes in grompp once the Drude stuff is ready to merge (reasonably soon).

#7 Updated by Erik Lindahl over 3 years ago

  • Status changed from New to Resolved

#8 Updated by Gerrit Code Review Bot over 3 years ago

Gerrit received a related patchset '1' for Issue #1376.
Uploader: David van der Spoel ()
Change-Id: I4382bcf5231920b22c725a7cf8e7b4e17c2526d9
Gerrit URL: https://gerrit.gromacs.org/3658

#9 Updated by David van der Spoel over 3 years ago

In relation to the original issue: both the 4.5 and the 4.6 tpr files had nstcalcenergy > 1 and hence both will have yielded incorrect results, albeit differently.

#10 Updated by Elizabeth Ploetz over 3 years ago

Hi all,

Thanks so much for your attention to my observation/mistake. Do you have any idea why 4.5.3 would be stable and 4.6 would crash since they both had nstcalcenergy > 1 (see previous comment) and yet that was identified as the "culprit" (see earlier comment)? I'm just trying to understand more fully the comment "and hence both will have yielded incorrect results, albeit differently."

Thanks again,

Elizabeth

#11 Updated by David van der Spoel over 3 years ago

You had nstcalcenergy = 10 in 4.5 and 100 in 4.6. The effect when running in parallel could be that the energies were only summed every 10 respectively 100 steps, and hence that the minimization was not global but only local. This is a serious error that we should have prevented earlier, but as you noted it was difficult to detect. I'm afraid I have already published simulations run with these settings.

#12 Updated by Elizabeth Ploetz over 3 years ago

Okay, thanks for helping me understand! I'm really sorry about your simulations (and mine ;) ).

#13 Updated by Elizabeth Ploetz over 3 years ago

Dear Prof. van der Spoel,

I have three questions:

First, if I change rlist=rcoulomb=rvdw=1.5 (so no twin-range) and change nstlist=1/nstcalcenergy=1, then I get a segmentation fault if I use more than three cores [tested on my local machine running Gromacs 4.6 and tested on DiaGrid, a distributed research computing grid provided by IT at Purdue running Gromacs 4.6.1 (to make sure this was not a problem with my Gromacs installation)]. Was your simulation stable in parallel with more than 3 cores? (I also tried nstlist=10/nstcalcenergy=1 and nstlist=5/nstcalcenergy=5, but this gave the same pattern of segmentation faults with >3 cores.)

Second, when you said that "Nstlist = 1 is bad," you mean computationally inefficient (as opposed to incorrect), am I right? I'm new to polarizable simulations, but when I run nonpolarizable simulations, I always have nstlist=10 and a twin-range cut-off. You are not suggesting that this combination is incorrect in nonpolarizable simulations, are you? (I don't see the point of a list if it must be updated every step!)

Third, what specific energy is it that is only calculated with a frequency of nstcalcenergy? I don't quite understand this parameter.

Thanks,

Elizabeth

#14 Updated by David van der Spoel over 3 years ago

For gromacs 4.6 you still need the -pd flag to mdrun when running in parallel. I have reproduced the crash and I will put in another warning to mdrun. -pd triggers the old way of parallellization in gromacs. However with the development gromacs version the calculation crashes in the same manner as with 4.6 without -pd.

Indeed nstlist = 1 is ineffective. The nstcalcenergy option is an optimization where the energy is only computed and summed over all processors every so many steps. This is under the presumption that only the forces are needed for the dynamics. However global communication is needed when the energy is optimized with respect to the shell positions, every step of the dynamics.

#15 Updated by Gerrit Code Review Bot over 3 years ago

Gerrit received a related patchset '1' for Issue #1376.
Uploader: David van der Spoel ()
Change-Id: Ic07ae70cd3cd5298b6b2d579719da4d267b2d8fb
Gerrit URL: https://gerrit.gromacs.org/3688

#16 Updated by Erik Lindahl over 3 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF