set(OUTPUT_SIMD "IBM_QPX")
elseif(CPU_DETECTION_FEATURES MATCHES " neon_asimd ")
set(OUTPUT_SIMD "ARM_NEON_ASIMD")
- elseif(CPU_DETECTION_FEATURES MATCHES " neon ")
+ elseif(CPU_DETECTION_FEATURES MATCHES " neon " AND NOT GMX_DOUBLE)
set(OUTPUT_SIMD "ARM_NEON")
endif()
endif()
# => compile sm_20, sm_30, sm_35, sm_37, sm_50, & sm_52 SASS, and compute_52 PTX
# - with CUDA >=8.0 CC 6.0-6.2 is supported (but we know nothing about CC 6.2, so we won't generate code or it)
# => compile sm_20, sm_30, sm_35, sm_37, sm_50, sm_52, sm_60, sm_61 SASS, and compute_60 and compute_61 PTX
- #
+ # - with CUDA >=9.0 CC 7.0 is supported and CC 2.0 is no longer supported
+ # => compile sm_30, sm_35, sm_37, sm_50, sm_52, sm_60, sm_61, sm_70 SASS, and compute_70 PTX
#
# Note that CUDA 6.5.19 second patch release supports cc 5.2 too, but
# CUDA_VERSION does not contain patch version and having PTX 5.0 JIT-ed is
# equally fast as compiling with sm_5.2 anyway.
# First add flags that trigger SASS (binary) code generation for physical arch
- list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_20,code=sm_20")
+ if(CUDA_VERSION VERSION_LESS "9.00") # < 9.0
+ list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_20,code=sm_20")
+ endif()
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_30,code=sm_30")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_35,code=sm_35")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_60,code=sm_60")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_61,code=sm_61")
endif()
+ if(NOT CUDA_VERSION VERSION_LESS "9.0") # >= 9.0
+ list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_70,code=sm_70")
+ endif()
# Next add flags that trigger PTX code generation for the newest supported virtual arch
# that's useful to JIT to future architectures
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_50,code=compute_50")
elseif(CUDA_VERSION VERSION_LESS "8.0")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_52,code=compute_52")
- else() # version >= 8.0
+ elseif(CUDA_VERSION VERSION_LESS "9.0")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_60,code=compute_60")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_61,code=compute_61")
+ else() # version >= 9.0
+ list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_70,code=compute_70")
endif()
endif()
elseif(GMX_SIMD_ACTIVE STREQUAL "ARM_NEON")
+ if (GMX_DOUBLE)
+ message(FATAL_ERROR "ARM_NEON SIMD support is not available for a double precision build because the architecture lacks double-precision support")
+ endif()
+
gmx_find_flags(
"#include<arm_neon.h>
int main(){float32x4_t x=vdupq_n_f32(0.5);x=vmlaq_f32(x,x,x);return vgetq_lane_f32(x,0)>0;}"
user-guide/mdrun-features.rst
user-guide/mdrun-performance.rst
user-guide/mdp-options.rst
+ user-guide/run-time-errors.rst
user-guide/file-formats.rst
user-guide/cmdline.rst
user-guide/environment-variables.rst
+ user-guide/terminology.rst
user-guide/plotje.gif
user-guide/xvgr.gif
conf.py
the run input file, one should use the minimized conformation from
the full precision trajectory file, as the structure file is not
accurate enough.
-{\tt \normindex{g_nmeig}} does the diagonalization and
+{\tt \normindex{gmx nmeig}} does the diagonalization and
the sorting of the normal modes according to their frequencies.
-Both {\tt mdrun} and {\tt g_nmeig} should be run in double precision.
-The normal modes can be analyzed with the program {\tt g_anaeig}.
+Both {\tt mdrun} and {\tt gmx nmeig} should be run in double precision.
+The normal modes can be analyzed with the program {\tt gmx anaeig}.
Ensembles of structures at any temperature and for any subset of
-normal modes can be generated with {\tt \normindex{g_nmens}}.
+normal modes can be generated with {\tt \normindex{gmx nmens}}.
An overview of normal-mode analysis and the related principal component
analysis (see \secref{covanal}) can be found in~\cite{Hayward95b}.
% } % Brace matches ifthenelse test for gmxlite
{\gromacs} now also supports the use of Bennett's Acceptance Ratio~\cite{Bennett1976}
for calculating values of $\Delta$G for transformations from state A to state B using
-the program {\tt \normindex{g_bar}}. The same data can also be used to calculate free
+the program {\tt \normindex{gmx bar}}. The same data can also be used to calculate free
energies using MBAR~\cite{Shirts2008}, though the analysis currently requires external tools from
the external {\tt pymbar} package, at https://SimTK.org/home/pymbar.
The procedure for essential dynamics sampling or flooding is as follows.
First, the eigenvectors and eigenvalues need to be determined
-using covariance analysis ({\tt g_covar})
-or normal-mode analysis ({\tt g_nmeig}).
+using covariance analysis ({\tt gmx covar})
+or normal-mode analysis ({\tt gmx nmeig}).
Then, this information is fed into {\tt make_edi},
which has many options for selecting vectors and setting parameters,
see {\tt gmx make_edi -h}.
%
% This file is part of the GROMACS molecular simulation package.
%
-% Copyright (c) 2013,2014, by the GROMACS development team, led by
+% Copyright (c) 2013,2014,2017, by the GROMACS development team, led by
% Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
% and including many others, as listed in the AUTHORS file in the
% top-level source directory and at http://www.gromacs.org.
if we check this by inserting $m=1$ we get back \eqnref{simplevar0}
\subsection{Summing energy terms}
-The {\tt \normindex{g_energy}} program can also sum energy terms into one, {\eg}
+The {\tt \normindex{gmx energy}} program can also sum energy terms into one, {\eg}
potential + kinetic = total. For the partial averages this is again easy
if we have $S$ energy components $s$:
\beq
%
% This file is part of the GROMACS molecular simulation package.
%
-% Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
+% Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
% Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
% and including many others, as listed in the AUTHORS file in the
% top-level source directory and at http://www.gromacs.org.
\left( \int_{t_0}^{{t_0}+t} P_{xz}(t') \mbox{d} t' \right)^2
\right\rangle_{t_0}
\eeq
-This can be done with {\tt g_energy}.
+This can be done with {\tt gmx energy}.
This method converges very slowly~\cite{Hess2002a}, and as such
a nanosecond simulation might not
be long enough for an accurate determination of the viscosity.
mdrun-performance
managing-simulations
mdp-options
+ run-time-errors
file-formats
cmdline
environment-variables
+ terminology
cases.
``-tunepme``
- Defaults to "on." If "on," will optimize various aspects of the
- PME and DD algorithms, shifting load between ranks and/or GPUs to
- maximize throughput
+ Defaults to "on." If "on," a Verlet-scheme simulation will
+ optimize various aspects of the PME and DD algorithms, shifting
+ load between ranks and/or GPUs to maximize throughput. Some
+ mdrun features are not compatible with this, and these ignore
+ this option.
``-dlb``
Can be set to "auto," "no," or "yes."
--- /dev/null
+Run-time errors
+===============
+
+The vast majority of error messages generated by |Gromacs| are descriptive,
+informing the user where the exact error lies. Some errors that arise are noted
+below, along with more details on what the issue is and how to solve it.
+
+General
+-------
+
+Cannot allocate memory
+^^^^^^^^^^^^^^^^^^^^^^
+
+The executed program has attempted to assign memory to be used in the
+calculation, but could not get a sufficient amount of memory from the operating
+system.
+
+Possible solutions are the following.
+
+* Install more memory in the computer.
+* Use a computer with more memory.
+* Reduce the scope of the number of atoms selected for analysis.
+* Reduce the length of trajectory file being processed.
+* In some cases confusion between Ångström and nm may lead to users wanting
+ to generate a water box that is 10\ :sup:`3` times larger than what they
+ think it is (e.g. :ref:`gmx solvate`).
+
+The user should bear in mind that the cost in time and/or memory for various
+activities will scale with the number of atoms/groups/residues N or the
+simulation length T as order N, N log N, or N\ :sup:`2` (or maybe worse!) and
+the same for T, depending on the type of activity. If it takes a long time,
+think about what you are doing, and the underlying algorithm (see the
+`reference manual`_, man page, or use the ``-h`` flag for the utility), and see
+if there's something sensible you can do that has better scaling properties.
+
+``gmx pdb2gmx``
+---------------
+
+Residue XXX not found in residue topology database
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This means that the force field you have selected while running
+:ref:`gmx pdb2gmx` does not have an entry in the residue topology database
+(:ref:`rtp` file) for XXX. The residue database entry is necessary both for
+stand-alone molecules (e.g. formaldehyde) or a peptide (standard or
+non-standard). This entry defines the atom types, connectivity, bonded and
+non-bonded interaction types for the residue and is necessary to use
+:ref:`gmx pdb2gmx` to build a :ref:`top` file. A residue database entry may be
+missing simply because the database does not contain the residue at all, or
+because the name is different.
+
+For new users, this error appears because they are running :ref:`gmx pdb2gmx`
+blindly on a ref:`pdb` file they have without consideration of the contents of
+the file. A force field is not something that is magical, it can only deal with
+molecules or residues (building blocks) that are provided in the residue
+database or included otherwise.
+
+If you want to use :ref:`gmx pdb2gmx` to automatically generate your topology,
+you have to ensure that the appropriate :ref:`rtp` entry is present within the
+desired force field and has the same name as the building block you are trying
+to use. If you call your molecule "HIS," then :ref:`gmx pdb2gmx` will not
+magically build a random molecule; it will try to build histidine, based on the
+``[ HIS ]`` entry in the :ref:`rtp` file, so it will look for the exact atomic
+entries for histidine, no more no less.
+
+If you want a topology for an arbitrary molecule, you cannot use
+:ref:`gmx pdb2gmx` (unless you build the :ref:`rtp` entry yourself). You will
+have to build it by hand, or use another program (such as :ref:`gmx x2top` or
+one of the scripts contributed by users) to build the :ref:`top` file.
+
+If there is not an entry for this residue in the database, then the options for
+obtaining the force field parameters are:
+
+* see if there is a different name being used for the residue in the residue
+ database and rename as appropriate,
+* parameterize the residue / molecule yourself (lots of work, even for an
+ expert),
+* find a topology file for the molecule, convert it to an :ref:`itp` file and
+ include it in your :ref:`top` file,
+* use another force field which has parameters available for this,
+* search the primary literature for publications for parameters for the residue
+ that are consistent with the force field that is being used.
+
+Once you have determined the parameters and topology for your residue, see
+Adding a residue to a force field for instructions on how to proceed.
+
+Long bonds and/or missing atoms
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+There are probably atoms missing earlier in the :ref:`pdb` file which makes
+:ref:`gmx pdb2gmx` go crazy. Check the screen output of :ref:`gmx pdb2gmx`, as
+as it will tell you which one is missing. Then add the atoms in your :ref:`pdb`
+file, energy minimization will put them in the right place, or fix the side
+chain with e.g. the `WhatIF program <http://swift.cmbi.ru.nl/whatif/>`__.
+
+Chain identifier 'X' was used in two non-sequential blocks
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This means that within the coordinate file fed to :ref:`gmx pdb2gmx`, the X
+chain has been split, possibly by the incorrect insertion of one molecule within
+another. The solution is simple: move the inserted molecule to a location within
+the file so that it is not splitting another molecule.
+
+This message may also mean that the same chain identifier has been used for two
+separate chains. In that case, rename the second chain to a unique identifier.
+
+Atom X is missing in residue XXX Y in the pdb file
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Related to the long bonds/missing atoms error above, this error is usually
+quite obvious in its meaning. That is, :ref:`gmx pdb2gmx` expects certain atoms
+within the given residue, based on the entries in the force field :ref:`rtp`
+file. There are several cases to which this error applies:
+
+* Missing hydrogen atoms; the error message may be suggesting that an entry in
+ the :ref:`hdb` file is missing. More likely, the nomenclature of your hydrogen
+ atoms simply does not match what is expected by the :ref:`rtp` entry. In this
+ case, use ``-ignh`` to allow :ref:`gmx pdb2gmx` to add the correct hydrogens
+ for you, or re-name the problematic atoms.
+* A terminal residue (usually the N-terminus) is missing H atoms; this usually
+ suggests that the proper ``-ter`` option has not been supplied or chosen
+ properly. In the case of the AMBER force fields, nomenclature is typically the
+ problem. N-terminal and C-terminal residues must be prefixed by N and C,
+ respectively. For example, an N-terminal alanine should not be listed in the
+ :ref:`pdb` file as ``ALA``, but rather ``NALA``, as specified in the
+ `ffamber instructions <http://ffamber.cnsm.csulb.edu/ffamber.php>`__.
+* Atoms are simply missing in the structure file provided to :ref:`gmx pdb2gmx`;
+ look for ``REMARK 465`` and ``REMARK 470`` entries in the :ref:`pdb` file.
+ These atoms will have to be modeled in using external software. There is no
+ |Gromacs| tool to re-construct incomplete models at present.
+
+Contrary to what the error message says, the use of the option ``-missing`` is
+almost always inappropriate. The ``-missing`` option should only be used to
+generate specialized topologies for amino acid-like molecules to take advantage
+of :ref:`rtp` entries. If you find yourself using ``-missing`` in order to
+generate a topology for a protein or nucleic acid, don't; the topology produced
+is likely physically unrealistic.
+
+Atom X in residue YYY not found in rtp entry
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+If you are attempting to assemble a topology using :ref:`gmx pdb2gmx`, the atom
+names are expected to match those found in the :ref:`rtp` file that define the
+building block(s) in your structure. In most cases, the problem arises from a
+naming mismatch, so simply re-name the atoms in your coordinate file
+appropriately. In other cases, you may be supplying a structure that has
+residues that do not conform to the expectations of the force field, in which
+case you should investigate why such a difference is occurring and make a
+decision based on what you find: use a different force field, manually edit the
+structure, etc.
+
+No force fields found (files with name 'forcefield.itp' in subdirectories ending on '.ff')
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This means your environment is not configured to use |Gromacs| properly, because
+:ref:`gmx pdb2gmx` cannot find its databases of forcefield information. This
+could happen because a |Gromacs| installation was moved from one location to
+another. Either follow the instructions about getting access to |Gromacs| after
+installation or re-install |Gromacs| before doing so.
+
+``gmx grompp``
+--------------
+
+Found a second defaults directive file
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This is caused by the ``[ defaults ]`` directive appearing more than once in the
+topology or force field files for the system -- it can only appear once. A
+typical cause of this is a second defaults being set in an included topology
+file (:ref:`itp`), that has been sourced from somewhere else. For
+specifications on how the topology files work, see `Reference Manual`_.::
+
+ [ defaults ]
+ ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
+ 1 1 no 1.0 1.0
+
+One solution is to simply comment out (or delete) the lines of code out in the
+file where it is included for the second time i.e.,::
+
+ ;[ defaults ]
+ ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
+ ;1 1 no 1.0 1.0
+
+A better approach to finding a solution is to re-think what you are doing. The
+``[ defaults ]`` directive should only be appearing at the top of your
+:ref:`top` file where you choose the force field. If you are trying to mix two
+force fields, then you are asking for trouble. If a molecule :ref:`itp` file
+tries to choose a force field, then whoever produced it is asking for trouble.
+
+Invalid order for directive xxx
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The directives in the :ref:`top` and :ref:`itp` files have rules about the order
+in which they can appear, and this error is seen when the order is violated.
+Consider the examples and discussion in chapter 5 of the |Gromacs| manual,
+and/or from tutorial material. The include file mechanism cannot be used to
+``#include`` a file in just any old location, because they contain directives
+and these have to be properly placed.
+
+In particular, "Invalid order for directive defaults" is a result of defaults
+being set in the topology or force field files in the inappropriate location;
+the ``[ defaults ]`` section can only appear once and must be the first
+directive in the topology. The ``[ defaults ]`` directive is typically present
+in the force field file (``forcefield.itp``), and is added to the topology when
+you ``#include`` this file in the system topology.
+
+If the directive in question is atomtypes (which is the most common source of
+this error) or any other bonded or nonbonded ``[ *types ]`` directive, typically
+the user is adding some non-standard species (ligand, solvent, etc.) that
+introduces new atom types or parameters into the system. As indicated above,
+these new types and parameters must appear before any ``[ moleculetype ]``
+directive. The force field has to be fully constructed before any molecules can
+be defined.
+
+Atom index n in position_restraints out of bounds
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+A common problem is placing position restraint files for multiple molecules out
+of order. Recall that a position restraint :ref:`itp` file containing a
+``[ position_restraints ]`` block can only belong to the ``[ moleculetype ]``
+block that contains it. For example,
+
+WRONG:::
+
+ #include "topol_A.itp"
+ #include "topol_B.itp"
+ #include "ligand.itp"
+
+ #ifdef POSRES
+ #include "posre_A.itp"
+ #include "posre_B.itp"
+ #include "ligand_posre.itp"
+ #endif
+
+RIGHT:::
+
+ #include "topol_A.itp"
+ #ifdef POSRES
+ #include "posre_A.itp"
+ #endif
+
+ #include "topol_B.itp"
+ #ifdef POSRES
+ #include "posre_B.itp"
+ #endif
+
+ #include "ligand.itp"
+ #ifdef POSRES
+ #include "ligand_posre.itp"
+ #endif
+
+Further, the atom index of each ``[ position_restraint ]`` must be relative to
+the ``[ moleculetype ]``, not relative to the system (because the parsing has
+not reached ``[ molecules ]`` yet, there is no such concept as "system"). So
+you cannot use the output of a tool like :ref:`gmx genrestr` blindly (as ``gmx
+genrestr -h`` warns).
+
+System has non-zero total charge
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Notifies you that counter-ions may be required for the system to neutralize the
+charge or there may be problems with the topology.
+
+If the charge is a non-integer, then this indicates that there is a problem
+with the topology. If :ref:`gmx pdb2gmx` has been used, then look at the right hand
+comment column of the atom listing, which lists the cumulative charge. This
+should be an integer after every residue (and/or charge group where
+applicable). This will assist in finding the residue where things start
+departing from integer values. Also check the capping groups that have been
+used.
+
+If the charge is already close to an integer, then the difference is caused by
+rounding errors and not a major problem.
+
+Note for PME users: It is possible to use a uniform neutralizing background
+charge in PME to compensate for a system with a net background charge. This may
+however, especially for non-homogeneous systems, lead to unwanted artifacts, as
+shown in Hub, J. S., de Groot, B. L., Grubmüller, H. & Groenhof, G. Quantifying
+artifacts in Ewald simulations of inhomogeneous systems with a net charge.
+*J. Chem. Theory Comput.* **10**, 381–390 (2014). Nevertheless, it is a standard
+practice to actually add counter-ions to make the system net neutral.
+
+Incorrect number of parameters
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Look at the topology file for the system. You've not given enough parameters
+for one of the bonded definitions. Sometimes this also occurs if you've mangled
+the include file mechanism or the topology file format (see `Reference Manual`_)
+when you edited the file.
+
+Number of coordinates in coordinate file does not match topology
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This is pointing out that, based on the information provided in the topology
+file, :ref:`top`, the total number of atoms or particles within the system does
+not match exactly with what is provided within the coordinate file, often a
+:ref:`gro` or a :ref:`pdb`.
+
+The most common reason for this is simply that the user has failed to update
+the topology file after solvating or adding additional molecules to the system,
+or made a typographical error in the number of one of the molecules within the
+system. Ensure that the end of the topology file being used contains something
+like the following, that matches exactly with what is within the coordinate
+file being used, in terms of both numbers and order of the molecules:::
+
+ [ molecules ]
+ ; Compound #mol
+ Protein 1
+ SOL 10189
+ NA+ 10
+
+Fatal error: No such moleculetype XXX
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Each type of molecule in your ``[ molecules ]`` section of your :ref:`top` file
+must have a corresponding ``[ moleculetype ]`` section defined previously,
+either in the :ref:`top` file or an included :ref:`itp` file. See
+`Reference Manual`_ for the syntax description. Your :ref:`top` file doesn't
+have such a definition for the indicated molecule. Check the contents of the
+relevant files, how you have named your molecules, and how you have tried to
+refer to them later. Pay attention to the status of ``#ifdef`` and/or
+``#include`` statements.
+
+T-Coupling group XXX has fewer than 10% of the atoms
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+It is possible to specify separate thermostats (temperature coupling groups)
+for every molecule type within a simulation. This is a particularly bad practice
+employed by many new users to molecular dynamics simulations. Doing so is a bad
+idea, as you can introduce errors and artifacts that are hard to predict. In
+some cases it is best to have all molecules within a single group, using system.
+If separate coupling groups are required to avoid the "hot solvent cold solute"
+problem, then ensure that they are of "sufficient size" and combine molecule
+types that appear together within the simulation. For example, for a protein in
+water with counter-ions, one would likely want to use ``Protein`` and
+``Non-Protein``.
+
+The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease rlist
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This error is generated in the cases as noted within the message. The
+dimensions of the box are such that an atom will interact with itself (when
+using periodic boundary conditions), thus violating the minimum image
+convention. Such an event is totally unrealistic and will introduce some
+serious artefacts. The solution is again what is noted within the message,
+either increase the size of the simulation box so that it is at an absolute
+minimum twice the cut-off length in all three dimensions (take care here if are
+using pressure coupling, as the box dimensions will change over time and if
+they decrease even slightly, you will still be violating the minimum image
+convention) or decrease the cut-off length (depending on the force field
+utilised, this may not be an option).
+
+Unknown left-hand XXXX in parameter file
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+:ref:`gmx grompp` has found an unknown term in the :ref:`mdp` file fed to it.
+You should check the spelling of XXXX and look for typographical errors. Be
+aware that quite a few run parameters changed between some |Gromacs| versions
+and the output from grompp will sometimes offer helpful commentary about these
+situations.
+
+Atom index (1) in bonds out of bounds
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This kind of error looks like:::
+
+ Fatal error:
+ [ file spc.itp, line 32 ]
+ Atom index (1) in bonds out of bounds (1-0).
+ This probably means that you have inserted topology
+ section "settles" in a part belonging to a different
+ molecule than you intended to. in that case move the
+ "settles" section to the right molecule.
+
+This error is fairly self-explanatory. You should look at your :ref:`top` file
+and check that all of the ``[ molecules ]`` sections contain all of the data
+pertaining to that molecule, and no other data. That is, you cannot
+``#include`` another molecule type (an :ref:`itp` file) before the previous
+``[ moleculetype ]`` has ended. Consult the examples in chapter 5 of the manual
+for information on the required ordering of the different ``[ sections ]``. Pay
+attention to the contents of any files you have included with ``#include``
+directives.
+
+This error can also arise if you are using a water model that is not enabled
+for use with your chosen force field by default. For example, if you are
+attempting to use the SPC water model with an AMBER force field, you will see
+this error. The reason is that, in spc.itp, there is no ``#ifdef`` statement
+defining atom types for any of the AMBER force fields. You can either add this
+section yourself, or use a different water model.
+
+XXX non-matching atom names
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This error usually indicates that the order of the topology file does not match
+that of the coordinate file. When running :ref:`gmx grompp`, the program reads
+through the topology, mapping the supplied parameters to the atoms in the
+coordinate file. If there is a mismatch, this error is generated. To remedy the
+problem, make sure that the contents of your ``[ molecules ]`` directive
+matches the exact order of the atoms in the coordinate file.
+
+In some cases, the error is harmless. For example, when running simulations
+with the MARTINI force field, the workflow relies on :ref:`gmx grompp` to apply
+the correct names, which are not previously assigned. Also, perhaps you are
+using a coordinate file that has the old (pre-version 4.5) ion nomenclature. In
+this particular case, allowing :ref:`gmx grompp` to re-assign names is harmless.
+For just about any other situation, when this error comes up, it should not be
+ignored. Just because the ``-maxwarn`` option is available does not mean you
+should use it in the blind hope of your simulation working. It will undoubtedly
+blow up.
+
+The sum of the two largest charge group radii (X) is larger than rlist - rvdw/rcoulomb
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This error warns that some combination of settings will result in poor energy
+conservation at the longest cutoff, which occurs when charge groups move in or
+out of neighborlist range. The error can have the following two sources.
+
+* Your charge groups encompass too many atoms. Most charge groups should be
+ less than 4 atoms or less.
+* Your :ref:`mdp` settings are incompatible with the chosen algorithms. For
+ switch or shift functions, rlist must be larger than the longest cutoff
+ (``rvdw`` or ``rcoulomb``) to provide buffer space for charge groups that move
+ beyond the neighbor searching radius. If set incorrectly, you may miss
+ interactions, contributing to poor energy conservation.
+
+A similar error ("The sum of the two largest charge group radii (X) is larger
+than rlist") can arise under two following circumstances.
+
+* The charge groups are inappropriately large or rlist is set too low.
+* Molecules are broken across periodic boundaries, which is not a problem in a
+ periodic system. In this case, the sum of the two largest charge groups will
+ correspond to a value of twice the box vector along which the molecule is
+ broken.
+
+Invalid line in coordinate file for atom X
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This error arises if the format of the :ref:`gro` file is broken in some way.
+The most common explanation is that the second line in the :ref:`gro` file
+specifies an incorrect number of atoms, causing :ref:`gmx grompp` to continue
+searching for atoms but finding box vectors.
+
+``gmx mdrun``
+-------------
+
+Stepsize too small, or no change in energy. Converged to machine precision, but not to the requested precision
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This is not an error as such. It is simply informing you that during the energy
+minimization process it reached the limit possible to minimize the structure
+with your current parameters. It does not mean that the system has not been
+minimized fully, but in some situations that may be the case. If the system has
+a significant amount of water present, then an Epot of the order of -105 to
+-106 (in conjunction with an Fmax between 10 and 1000 kJ mol-1 nm-1) is
+typically a reasonable value for starting most MD simulations from the
+resulting structure. The most important result is likely the value of Fmax, as
+it describes the slope of the potential energy surface, i.e. how far from an
+energy minimum your structure lies. Only for special purposes, such as normal
+mode analysis type of calculations, it may be necessary to minimize further.
+
+Further minimization may be achieved by using a different energy minimization
+method or by making use of double precision-enabled |Gromacs|.
+
+LINCS/SETTLE/SHAKE warnings
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Sometimes, when running dynamics, mdrun may suddenly stop (perhaps after
+writing several :ref:`pdb` files) after a series of warnings about the
+constraint algorithms (e.g. LINCS, SETTLE or SHAKE) are written to the
+:ref:`log` file. These algorithms often used to constrain bond lengths and/or
+angles. When a system is blowing up (i.e. exploding due to diverging forces),
+the constraints are usually the first thing to fail. This doesn't necessarily
+mean you need to troubleshoot the constraint algorithm. Usually it is a sign
+of something more fundamentally wrong (physically unrealistic) with your system.
+See also the advice here about diagnosing unstable systems.
+
+1-4 interaction not within cut-off
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Some of your atoms have moved so two atoms separated by three bonds are
+separated by more than the cut-off distance. This is very bad. Most importantly,
+do not increase your cut-off! This error actually indicates that the atoms have
+very large velocities, which usually means that (part of) your molecule(s) is
+(are) blowing up. If you are using LINCS for constraints, you probably also
+already got a number of LINCS warnings. When using SHAKE this will give rise to
+a SHAKE error, which halts your simulation before the "1-4 not within cutoff"
+error can appear.
+
+There can be a number of reasons for the large velocities in your system. If it
+happens at the beginning of the simulation, your system might be not
+equilibrated well enough (e.g. it contains some bad contacts). Try a(nother)
+round of energy minimization to fix this. Otherwise you might have a very high
+temperature, and/or a timestep that is too large. Experiment with these
+parameters until the error stops occurring. If this doesn't help, check the
+validity of the parameters in your topology!
+
+Simulation running but no output
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Not an error as such, but mdrun appears to be chewing up CPU time but nothing
+is being written to the output files. There are a number of reasons why this
+may occur.
+
+* Your simulation might simply be (very) slow, and since output is buffered, it
+ can take quite some time for output to appear in the respective files. If you
+ are trying to fix some problems and you want to get output as fast as
+ possible, you can set the environment variable ``GMX_LOG_BUFFER`` to 0.
+* Something might be going wrong in your simulation, causing e.g. not-a-numbers
+ (NAN) to be generated (these are the result of e.g. division by zero).
+ Subsequent calculations with NAN's will generate floating point exceptions
+ which slow everything down by orders of magnitude.
+* You might have all ``nst*`` parameters (see your :ref:`mdp` file) set to 0,
+ this will suppress most output.
+* Your disk might be full. Eventually this will lead to mdrun crashing, but
+ since output is buffered, it might take a while for mdrun to realize it can't
+ write.
+
+Can not do Conjugate Gradients with constraints
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This error means you can't do energy minimization with the conjugate gradient
+algorithm if your topology has constraints defined (see `Reference Manual`_
+for details).
+
+Pressure scaling more than 1%
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This error tends to be generated when the simulation box begins to oscillate
+(due to large pressures and/or small coupling constants), the system starts
+to resonate and then crashes. This can mean that the system isn't equilibrated
+sufficiently before using pressure coupling. Therefore, better/more
+equilibration may fix the issue.
+
+It is recommended to observe the system trajectory prior and during the crash.
+This may indicate if a particular part of the system/structure is the problem.
+
+In some cases, if the system has been equilibrated sufficiently, this error can
+mean that the pressure coupling constant, ``tau_p``, is too small (particularly
+when using the Berendsen weak coupling method). Increasing that value will slow
+down the response to pressure changes and may stop the resonance from occuring.
+
+This error can also appear when using a timestep that is too large, e.g. 5 fs,
+in the absence of constraints and/or virtual sites.
+
+Range Checking error
+^^^^^^^^^^^^^^^^^^^^
+
+This usually means your simulation is blowing up. Probably you need to do
+better energy minimization and/or equilibration and/or topology design.
+
+X particles communicated to PME node Y are more than a cell length out of the domain decomposition cell of their charge group
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This is another way that mdrun tells you your system is blowing up. In
+|Gromacs| version 4.0, domain decomposition was introduced to divide the system
+into regions containing nearby atoms (see `Reference Manual`_ or Hess, B.,
+Kutzner, C., Van Der Spoel, D. & Lindahl, E. GROMACS 4: algorithms for highly
+efficient, load-balanced, and scalable molecular simulation.
+*J. Chem. Theory Comput.* **4**, 435–447 (2008) for more details). If you have
+particles that are flying across the system, you will get this fatal error.
+The message indicates that some piece of your system is tearing apart (hence
+out of the "cell of their charge group"). Refer to the Blowing Up page for
+advice on how to fix this issue.
+
+A charge group moved too far between two domain decomposition steps
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+As immediately above.
+
+Software inconsistency error: Some interactions seem to be assigned multiple times
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+As immediately above.
+
+There is no domain decomposition for n nodes that is compatible with the given box and a minimum cell size of x nm
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This means you tried to run a parallel calculation, and when mdrun tried to
+partition your simulation cell into chunks for each processor, it couldn't. The
+minimum cell size is controlled by the size of the largest charge group or
+bonded interaction and the largest of rvdw, rlist and rcoulomb, some other
+effects of bond constraints, and a safety margin. Thus it is not possible to
+run a small simulation with large numbers of processors. So, if :ref:`gmx grompp` warned
+you about a large charge group, pay attention and reconsider its size. mdrun
+prints a breakdown of how it computed this minimum size in the :ref:`log` file, so
+you can perhaps find a cause there.
+
+If you didn't think you were running a parallel calculation, be aware that,
+from version 4.5 onwards, |Gromacs| uses thread-based parallelism by default.
+To prevent this, you can either give mdrun the ``-nt 1`` command line option, or
+build |Gromacs| so that it will not use threads. Otherwise, you might be using
+an MPI-enabled |Gromacs| and not be aware of the fact.
+
+.. _reference manual: gmx-manual-parent-dir_
--- /dev/null
+Terminology
+===========
+
+Pressure
+--------
+
+The pressure in molecular dynamics can be computed from the kinetic energy and
+the virial.
+
+Fluctuation
+^^^^^^^^^^^
+
+Whether or not pressure coupling is used within a simulation, the pressure
+value for the simulation box will oscillate significantly. Instantaneous
+pressure is meaningless, and not well-defined. Over a picosecond time scale it
+usually will not be a good indicator of the true pressure. This variation is
+entirely normal due to the fact that pressure is a macroscopic property and can
+only be measured properly as time average, while it is being measured and/or
+adjusted with pressure coupling on the microscopic scale. How much it varies
+and the speed at which it does depends on the number of atoms in the system,
+the type of pressure coupling used and the value of the coupling constants.
+Fluctuations of the order of hundreds of bar are typical. For a box of 216
+waters, fluctuations of 500-600 bar are standard. Since the fluctuations go
+down with the square root of the number of particles, a system of 21600 water
+molecules (100 times larger) will still have pressure fluctuations of 50-60 bar.
+
+Periodic boundary conditions
+----------------------------
+
+Periodic boundary conditions (PBC) are used in molecular dynamics simulations
+to avoid problems with boundary effects caused by finite size, and make the
+system more like an infinite one, at the cost of possible periodicity effects.
+
+Beginners visualizing a trajectory sometimes think they are observing a problem
+when
+
+* the molecule(s) does not stay in the centre of the box, or
+* it appears that (parts of) the molecule(s) diffuse out of the box, or
+* holes are created, or
+* broken molecules appear, or
+* their unit cell was a rhombic dodecahedron or cubic octahedron but it looks
+ like a slanted cube after the simulation, or
+* crazy bonds all across the simulation cell appear.
+
+This is not a problem or error that is occuring, it is what you should expect.
+
+The existence of PBC means that any atom that leaves a simulation box by, say,
+the right-hand face, then enters the simulation box by the left-hand face. In
+the example of a large protein, if you look at the face of the simulation box
+that is opposite to the one from which the protein is protruding, then a hole
+in the solvent will be visible. The reason that the molecule(s) move from where
+they were initially located within the box is (for the vast majority of
+simulations) they are free to diffuse around. And so they do. They are not held
+in a magic location of the box. The box is not centered around anything while
+performing the simulation. Molecules are not made whole as a matter of course.
+Moreover, any periodic cell shape can be expressed as a parallelepiped (a.k.a.
+triclinic cell), and |Gromacs| does so internally regardless of the initial
+shape of the box.
+
+These visual issues can be fixed after the conclusion of the simulation by
+judicious use of the optional inputs to :ref:`gmx trjconv` to process the
+trajectory files. Similarly, analyses such as RMSD of atomic positions can be
+flawed when a reference structure is compared with a structure that needs
+adjusting for periodicity effects, and the solution with :ref:`gmx trjconv`
+follows the same lines. Some complex cases needing more than one operation will
+require more than one invocation of :ref:`gmx trjconv` in order to work.
+
+For further information, see `Reference Manual`_.
+
+Suggested workflow
+^^^^^^^^^^^^^^^^^^
+
+Fixing periodicity effects with :ref:`gmx trjconv` to suit visualization or
+analysis can be tricky. Multiple invocations can be necessary. You may need to
+create custom index groups (e.g. to keep your ligand with your protein)
+Following the steps below in order (omitting those not required) should help
+get a pleasant result. You will need to consult ``gmx trjconv -h`` to find out
+the details for each step. That's deliberate -- there is no magic "do what I
+want" recipe. You have to decide what you want, first. :-)
+
+#. First make your molecules whole if you want them whole.
+#. Cluster your molecules/particles if you want them clustered.
+#. If you want jumps removed, extract the first frame from the trajectory to
+ use as the reference, and then use ``-pbc nojump`` with that first
+ frame as reference.
+#. Center your system using some criterion. Doing so shifts the system, so
+ don't use ``-pbc nojump`` after this step.
+#. Perhaps put everything in some box with the other ``-pbc`` or ``-ur``
+ options.
+#. Fit the resulting trajectory to some (other) reference structure (if
+ desired), and don't use any PBC related option afterwards.
+
+With point three, the issue is that :ref:`gmx trjconv` removes the jumps from
+the first frame using the reference structure provided with -s. If the reference
+structure (run input file) is not clustered/whole, using ``-pbc nojump``
+will undo steps 1 and 2.
+
+Thermostats
+-----------
+
+Thermostats are designed to help a simulation sample from the correct ensemble
+(i.e. NVT or NPT) by modulating the temperature of the system in some fashion.
+First, we need to establish what we mean by temperature. In simulations, the
+"instantaneous (kinetic) temperature" is usually computed from the kinetic
+energy of the system using the equipartition theorem. In other words, the
+temperature is computed from the system's total kinetic energy.
+
+So, what's the goal of a thermostat? Actually, it turns out the goal is not to
+keep the temperature constant, as that would mean fixing the total kinetic
+energy, which would be silly and not the aim of NVT or NPT. Rather, it's to
+ensure that the average temperature of a system be correct.
+
+To see why this is the case, imagine a glass of water sitting in a room.
+Suppose you can look very closely at a few molecules in some small region of
+the glass, and measure their kinetic energies. You would not expect the kinetic
+energy of this small number of particles to remain precisely constant; rather,
+you'd expect fluctuations in the kinetic energy due to the small number of
+particles. As you average over larger and larger numbers of particles, the
+fluctuations in the average get smaller and smaller, so finally by the time you
+look at the whole glass, you say it has "constant temperature".
+
+Molecular dynamics simulations are often fairly small compared to a glass of
+water, so we have bigger fluctuations. So it's really more appropriate here to
+think of the role of a thermostat as ensuring that we have
+
+(a) the correct average temperature, and
+(b) the fluctuations of the correct size.
+
+See `Reference Manual`_ for details on how temperature coupling is applied and
+the types currently available.
+
+What to do
+^^^^^^^^^^
+
+Some hints on practices that generally are a good idea:
+
+* Preferably, use a thermostat that samples the correct distribution of
+ temperatures (for examples, see the corresponding manual section), in addition
+ to giving you the correct average temperature.
+* At least: use a thermostat that gives you the correct average temperature,
+ and apply it to components of your system for which they are justified (see
+ the first bullet in `What not to do`_). In some cases, using
+ ``tc-grps = System`` may lead to the "hot solvent/cold solute" problem
+ described in the 3rd reference in `Further reading`_.
+
+What not to do
+^^^^^^^^^^^^^^
+
+Some hints on practices that generally not a good idea to use:
+
+* Do not use separate thermostats for every component of your system. Some
+ molecular dynamics thermostats only work well in the thermodynamic limit. A
+ group must be of sufficient size to justify its own thermostat. If you use one
+ thermostat for, say, a small molecule, another for protein, and another for
+ water, you are likely introducing errors and artifacts that are hard to
+ predict. In particular, do not couple ions in aqueous solvent in a separate
+ group from that solvent. For a protein simulation, using ``tc-grps = Protein
+ Non-Protein`` is usually best.
+* Do not use thermostats that work well only in the limit of a large number of
+ degrees of freedom for systems with few degrees of freedom. For example, do
+ not use Nosé-Hoover or Berendsen thermostats for types of free energy
+ calculations where you will have a component of the system with very few
+ degrees of freedom in an end state (i.e. a noninteracting small molecule).
+
+Further reading
+^^^^^^^^^^^^^^^
+
+#. Cheng, A. & Merz, K. M. Application of the nosé- hoover chain algorithm to
+ the study of protein dynamics. *J. Phys. Chem.* **100** (5), 1927–1937
+ (1996).
+#. Mor, A., Ziv, G. & Levy, Y. Simulations of proteins with inhomogeneous
+ degrees of freedom: the effect of thermostats. *J. Comput. Chem.* **29**
+ (12), 1992–1998 (2008).
+#. Lingenheil, M., Denschlag, R., Reichold, R. & Tavan, P. The
+ "hot-solvent/cold-solute" problem revisited. *J. Chem. Theory Comput.* **4**
+ (8), 1293–1306 (2008).
+
+Energy conservation
+-------------------
+
+In principle, a molecular dynamics simulation should conserve the total energy,
+the total momentum and (in a non-periodic system) the total angular momentum. A
+number of algorithmic and numerical issues make that this is not always the
+case:
+
+* Cut-off treatment and/or long-range electrostatics treatment (see Van Der
+ Spoel, D. & van Maaren, P. J. The origin of layer structure artifacts in
+ simulations of liquid water. *J. Chem. Theor. Comp.* **2**, 1–11 (2006).)
+* Treatment of neighborlists,
+* Constraint algorithms (see e.g. Hess, B. P-LINCS: A parallel linear constraint
+ solver for molecular simulation. *J. Chem. Theor. Comp.* **4**, 116–122
+ (2008).)
+* The integration timestep,
+* Temperature coupling and pressure coupling,
+* Round-off error (in particular in single precision), for example subtracting
+ large numbers (Lippert, R. A. et al. A common, avoidable source of error in
+ molecular dynamics integrators. *J. Chem. Phys.* **126**, 046101 (2007).)
+* The choice of the integration algorithm (in |Gromacs| this is normally
+ leap-frog),
+* Removal of center of mass motion: when doing this in more than one group the
+ conservation of energy will be violated.
+
+Average structure
+-----------------
+
+Various |Gromacs| utilities can compute average structures. Presumably the idea
+for this comes from something like an ensemble-average NMR structure. In some
+cases, it makes sense to calculate an average structure (as a step on the way
+to calculating root-mean-squared fluctuations (RMSF), for example, one needs
+the average position of all of the atoms).
+
+However, it's important to remember that an average structure isn't necessarily
+meaningful. By way of analogy, suppose I alternate holding a ball in my left
+hand, then in my right hand. What's the average position of the ball? Halfway
+in between -- even though I always have it either in my left hand or my right
+hand. Similarly, for structures, averages will tend to be meaningless anytime
+there are separate metastable conformational states. This can happen on a
+sidechain level, or for some regions of backbone, or even whole helices or
+components of the secondary structure.
+
+Thus, if you derive an average structure from a molecular dynamics simulation,
+and find artifacts like unphysical bond lengths, weird structures, etc., this
+doesn't necessarily mean something is wrong. It just shows the above: an
+average structure from a simulation is not necessarily a physically meaningful
+structure.
+
+Blowing up
+----------
+
+*Blowing up* is a highly technical term used to describe a common sort of
+simulation failure. In brief, it describes a failure typically due to an
+unacceptably large force that ends up resulting in a failure of the integrator.
+
+To give a bit more background, it's important to remember that molecular
+dynamics numerically integrates Newton's equations of motion by taking small,
+discrete timesteps, and using these timesteps to determine new velocities and
+positions from velocities, positions, and forces at the previous timestep. If
+forces become too large at one timestep, this can result in extremely large
+changes in velocity/position when going to the next timestep. Typically, this
+will result in a cascade of errors: one atom experiences a very large force one
+timestep, and thus goes shooting across the system in an uncontrolled way in
+the next timestep, overshooting its preferred location or landing on top of
+another atom or something similar. This then results in even larger forces the
+next timestep, more uncontrolled motions, and so on. Ultimately, this will
+cause the simulation package to crash in some way, since it can't cope with
+such situations. In simulations with constraints, the first symptom of this
+will usually be some LINCS or SHAKE warning or error -- not because the
+constraints are the source of the problem, but just because they're the first
+thing to crash. Similarly, in simulations with domain decomposition, you may
+see messages about particles being more than a cell length out of the domain
+decomposition cell of their charge group, which are symptomatic of your
+underlying problem, and not the domain decomposition algorithm itself. Likewise
+for warnings about tabulated or 1-4 interactions being outside the distance
+supported by the table. This can happen on one computer system while another
+resulted in a stable simulation because of the impossibility of numerical
+reproducibility of these calculations on different computer systems.
+
+Possible causes include:
+
+* you didn't minimize well enough,
+* you have a bad starting structure, perhaps with steric clashes,
+* you are using too large a timestep (particularly given your choice of
+ constraints),
+* you are doing particle insertion in free energy calculations without using
+ soft core,
+* you are using inappropriate pressure coupling (e.g. when you are not in
+ equilibrium, Berendsen can be best while relaxing the volume, but you will
+ need to switch to a more accurate pressure-coupling algorithm later),
+* you are using inappropriate temperature coupling, perhaps on inappropriate
+ groups, or
+* your position restraints are to coordinates too different from those present
+ in the system, or
+* you have a single water molecule somewhere within the system that is
+ isolated from the other water molecules, or
+* you are experiencing a bug in :ref:`gmx mdrun`.
+
+Because blowing up is due, typically, to forces that are too large for a
+particular timestep size, there are a couple of basic solutions:
+
+* make sure the forces don't get that large, or
+* use a smaller timestep.
+
+Better system preparation is a way to make sure that forces don't get large, if
+the problems are occurring near the beginning of a simulation.
+
+Diagnosing an unstable system
+-----------------------------
+
+Troubleshooting a system that is blowing up can be challenging, especially for
+an inexperienced user. Here are a few general tips that one may find useful
+when addressing such a scenario:
+
+#. If the crash is happening relatively early (within a few steps), set
+ ``nstxout`` (or ``nstxout-compressed``) to 1, capturing all possible frames.
+ Watch the resulting trajectory to see which atoms/residues/molecules become
+ unstable first.
+#. Simplify the problem to try to establish a cause:
+
+ * If you have a new box of solvent, try minimizing and simulating a single
+ molecule to see if the instability is due to some inherent problem with
+ the molecule's topology or if instead there are clashes in your starting
+ configuration.
+ * If you have a protein-ligand system, try simulating the protein alone in
+ the desired solvent. If it is stable, simulate the ligand in vacuo to see
+ if its topology gives stable configurations, energies, etc.
+ * Remove the use of fancy algorithms, particularly if you haven't
+ equilibrated thoroughly first
+
+#. Monitor various components of the system's energy using :ref:`gmx energy`.
+ If an intramolecular term is spiking, that may indicate improper bonded
+ parameters, for example.
+#. Make sure you haven't been ignoring error messages (missing atoms when
+ running :ref:`gmx pdb2gmx`, mismatching names when running :ref:`gmx grompp`,
+ etc.) or using work-arounds (like using ``gmx grompp -maxwarn`` when you
+ shouldn't be) to make sure your topology is intact and being interpreted
+ correctly.
+#. Make sure you are using appropriate settings in your :ref:`mdp` file for the
+ force field you have chosen and the type of system you have. Particularly
+ important settings are treatment of cutoffs, proper neighbor searching
+ interval (``nstlist``), and temperature coupling. Improper settings can lead
+ to a breakdown in the model physics, even if the starting configuration of
+ the system is reasonable.
+
+If using implicit solvation, starting your equilibration with a smaller time
+step than your production run can help energy equipartition more stably.
+
+There are several common situations in which instability frequently arises,
+usually in the introduction of new species (ligands or other molecules) into
+the system. To determine the source of the problem, simplify the system (e.g.
+the case of a protein-ligand complex) in the following way.
+
+#. Does the protein (in water) minimize adequately by itself? This is a test of
+ the integrity of the coordinates and system preparation. If this fails,
+ something probably went wrong when running :ref:`gmx pdb2gmx` (see below), or
+ maybe :ref:`gmx genion` placed an ion very close to the protein (it is
+ random, after all).
+#. Does the ligand minimize in vacuo? This is a test of the topology. If it
+ does not, check your parameterization of the ligand and any implementation of
+ new parameters in force field files.
+#. (If previous item is successful) Does the ligand minimize in water, and/or
+ does a short simulation of the ligand in water succeed?
+
+Other sources of possible problems are in the biomolecule topology itself.
+
+#. Did you use ``-missing`` when running :ref:`gmx pdb2gmx`? If so, don't.
+ Reconstruct missing coordinates rather than ignoring them.
+#. Did you override long/short bond warnings by changing the lengths? If so,
+ don't. You probably have missing atoms or some terrible input geometry.
+
+.. _reference manual: gmx-manual-parent-dir_
gmx_bool bUseGPU,
gmx_bool *bPrinting)
{
+ GMX_RELEASE_ASSERT(ir->cutoff_scheme != ecutsGROUP, "PME tuning is not supported with cutoff-scheme=group (because it contains bugs)");
+
pme_load_balancing_t *pme_lb;
real spm, sp;
int d;
}
else
{
+ /* TODO Remove these lines and pme_lb->cutoff_scheme */
tmpr_coulomb = set->rcut_coulomb + pme_lb->rbufOuter_coulomb;
tmpr_vdw = pme_lb->rcut_vdw + pme_lb->rbufOuter_vdw;
+ /* Two (known) bugs with cutoff-scheme=group here:
+ * - This modification of rlist results in incorrect DD comunication.
+ * - We should set fr->bTwinRange = (fr->rlistlong > fr->rlist).
+ */
set->rlistOuter = std::min(tmpr_coulomb, tmpr_vdw);
set->rlistInner = set->rlistOuter;
}
/* TODO: centralize the code that sets the potentials shifts */
if (ic->coulomb_modifier == eintmodPOTSHIFT)
{
- ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+ GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero");
+ ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb) / ic->rcoulomb;
}
if (EVDW_PME(ic->vdwtype))
{
"The diffusion coefficient is determined by linear regression of the MSD,",
"where, unlike for the normal output of D, the times are weighted",
"according to the number of reference points, i.e. short times have",
- "a higher weight. Also when [TT]-beginfit[tt]=-1,fitting starts at 10%",
- "and when [TT]-endfit[tt]=-1, fitting goes to 90%.",
+ "a higher weight. Also when [TT]-beginfit[tt] is -1, fitting starts at 10%",
+ "and when [TT]-endfit[tt] is -1, fitting goes to 90%.",
"Using this option one also gets an accurate error estimate",
"based on the statistics between individual molecules.",
"Note that this diffusion coefficient and error estimate are only",
#ifndef CUDA_ARCH_UTILS_CUH_
#define CUDA_ARCH_UTILS_CUH_
+#include "config.h"
+
/*! \file
* \brief CUDA arch dependent definitions.
*
*/
static const int warp_size = 32;
static const int warp_size_log2 = 5;
+/*! \brief Bitmask corresponding to all threads active in a warp.
+ * NOTE that here too we assume 32-wide warps.
+ */
+static const unsigned int c_fullWarpMask = 0xffffffff;
+
+/* Below are backward-compatibility wrappers for CUDA 9 warp-wide intrinsics. */
+
+/*! \brief Compatibility wrapper around the CUDA __syncwarp() instrinsic. */
+static __forceinline__ __device__
+void gmx_syncwarp(const unsigned int activeMask = c_fullWarpMask)
+{
+#if GMX_CUDA_VERSION < 9000
+ /* no sync needed on pre-Volta. */
+ GMX_UNUSED_VALUE(activeMask);
+#else
+ __syncwarp(activeMask);
+#endif
+}
+
+/*! \brief Compatibility wrapper around the CUDA __ballot()/__ballot_sync() instrinsic. */
+static __forceinline__ __device__
+unsigned int gmx_ballot_sync(const unsigned int activeMask,
+ const int pred)
+{
+#if GMX_CUDA_VERSION < 9000
+ GMX_UNUSED_VALUE(activeMask);
+ return __ballot(pred);
+#else
+ return __ballot_sync(activeMask, pred);
+#endif
+}
+
+/*! \brief Compatibility wrapper around the CUDA __any()/__any_sync() instrinsic. */
+static __forceinline__ __device__
+int gmx_any_sync(const unsigned int activeMask,
+ const int pred)
+{
+#if GMX_CUDA_VERSION < 9000
+ GMX_UNUSED_VALUE(activeMask);
+ return __any(pred);
+#else
+ return __any_sync(activeMask, pred);
+#endif
+}
+
+/*! \brief Compatibility wrapper around the CUDA __shfl_up()/__shfl_up_sync() instrinsic. */
+template <typename T>
+static __forceinline__ __device__
+T gmx_shfl_up_sync(const unsigned int activeMask,
+ const T var,
+ unsigned int offset)
+{
+#if GMX_CUDA_VERSION < 9000
+ GMX_UNUSED_VALUE(activeMask);
+ return __shfl_up(var, offset);
+#else
+ return __shfl_up_sync(activeMask, var, offset);
+#endif
+}
+
+/*! \brief Compatibility wrapper around the CUDA __shfl_down()/__shfl_down_sync() instrinsic. */
+template <typename T>
+static __forceinline__ __device__
+T gmx_shfl_down_sync(const unsigned int activeMask,
+ const T var,
+ unsigned int offset)
+{
+#if GMX_CUDA_VERSION < 9000
+ GMX_UNUSED_VALUE(activeMask);
+ return __shfl_down(var, offset);
+#else
+ return __shfl_down_sync(activeMask, var, offset);
+#endif
+}
/*! \brief Allow disabling CUDA textures using the GMX_DISABLE_CUDA_TEXTURES macro.
*
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2012,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2012,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
/*! Launches asynchronous host to device memory copy in stream s. */
int cu_copy_D2H_async(void * /*h_dest*/, void * /*d_src*/, size_t /*bytes*/, cudaStream_t /*s = 0*/);
-
/*! Launches synchronous host to device memory copy. */
int cu_copy_H2D(void * /*d_dest*/, void * /*h_src*/, size_t /*bytes*/);
ic->epsfac = fr->epsfac;
ic->ewaldcoeff_q = fr->ewaldcoeff_q;
- if (fr->coulomb_modifier == eintmodPOTSHIFT)
+ if (EEL_PME_EWALD(ic->eeltype) && ic->coulomb_modifier == eintmodPOTSHIFT)
{
- ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+ GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero");
+ ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb) / ic->rcoulomb;
}
else
{
*
* Note that the current kernel implementation only supports NTHREAD_Z > 1 with
* shuffle-based reduction, hence CC >= 3.0.
+ *
+ *
+ * NOTEs / TODO on Volta / CUDA 9 support extensions:
+ * - the current way of computing active mask using ballot_sync() should be
+ * reconsidered: we can compute all masks with bitwise ops iso ballot and
+ * secondly, all conditionals are warp-uniform, so the sync is not needed;
+ * - reconsider the use of __syncwarp(): its only role is currently to prevent
+ * WAR hazard due to the cj preload; we should try to replace it with direct
+ * loads (which may be faster given the improved L1 on Volta).
*/
/* Kernel launch bounds for different compute capabilities. The value of NTHREAD_Z
#endif /* CALC_ENERGIES */
#ifdef EXCLUSION_FORCES
- const int nonSelfInteraction = !(nb_sci.shift == CENTRAL & tidxj <= tidxi);
+ const int nonSelfInteraction = !(nb_sci.shift == CENTRAL & tidxj <= tidxi);
#endif
+ int j4LoopStart = cij4_start + tidxz;
+ unsigned int j4LoopThreadMask = gmx_ballot_sync(c_fullWarpMask, j4LoopStart < cij4_end);
/* loop over the j clusters = seen by any of the atoms in the current super-cluster */
- for (j4 = cij4_start + tidxz; j4 < cij4_end; j4 += NTHREAD_Z)
+ for (j4 = j4LoopStart; j4 < cij4_end; j4 += NTHREAD_Z)
{
wexcl_idx = pl_cj4[j4].imei[widx].excl_ind;
imask = pl_cj4[j4].imei[widx].imask;
wexcl = excl[wexcl_idx].pair[(tidx) & (warp_size - 1)];
+ unsigned int imaskSkipConditionThreadMask = j4LoopThreadMask;
#ifndef PRUNE_NBL
+ imaskSkipConditionThreadMask = gmx_ballot_sync(j4LoopThreadMask, imask);
if (imask)
#endif
{
{
cjs[tidxi + tidxj * c_nbnxnGpuJgroupSize/c_splitClSize] = pl_cj4[j4].cj[tidxi];
}
+ gmx_syncwarp(imaskSkipConditionThreadMask);
/* Unrolling this loop
- with pruning leads to register spilling;
Tested with up to nvcc 7.5 */
for (jm = 0; jm < c_nbnxnGpuJgroupSize; jm++)
{
- if (imask & (superClInteractionMask << (jm * c_numClPerSupercl)))
+ const unsigned int jmSkipCondition = imask & (superClInteractionMask << (jm * c_numClPerSupercl));
+ const unsigned int jmSkipConditionThreadMask = gmx_ballot_sync(imaskSkipConditionThreadMask, jmSkipCondition);
+ if (jmSkipCondition)
{
mask_ji = (1U << (jm * c_numClPerSupercl));
#endif
for (i = 0; i < c_numClPerSupercl; i++)
{
- if (imask & mask_ji)
+ const unsigned int iInnerSkipCondition = imask & mask_ji;
+ const unsigned int iInnerSkipConditionThreadMask = gmx_ballot_sync(jmSkipConditionThreadMask, iInnerSkipCondition);
+ if (iInnerSkipCondition)
{
ci = sci * c_numClPerSupercl + i; /* i cluster index */
/* If _none_ of the atoms pairs are in cutoff range,
the bit corresponding to the current
cluster-pair in imask gets set to 0. */
- if (!__any(r2 < rlist_sq))
+ if (!gmx_any_sync(iInnerSkipConditionThreadMask, r2 < rlist_sq))
{
imask &= ~mask_ji;
}
}
/* reduce j forces */
- reduce_force_j_warp_shfl(fcj_buf, f, tidxi, aj);
+ reduce_force_j_warp_shfl(fcj_buf, f, tidxi, aj, jmSkipConditionThreadMask);
}
}
#ifdef PRUNE_NBL
pl_cj4[j4].imei[widx].imask = imask;
#endif
}
+ // avoid shared memory WAR hazards between loop iterations
+ gmx_syncwarp(j4LoopThreadMask);
+ // update thread mask for next loop iteration
+ j4LoopThreadMask = gmx_ballot_sync(j4LoopThreadMask, (j4 + NTHREAD_Z) < cij4_end);
}
/* skip central shifts when summing shift forces */
ai = (sci * c_numClPerSupercl + i) * c_clSize + tidxi;
reduce_force_i_warp_shfl(fci_buf[i], f,
&fshift_buf, bCalcFshift,
- tidxj, ai);
+ tidxj, ai, c_fullWarpMask);
}
/* add up local shift forces into global mem, tidxj indexes x,y,z */
#ifdef CALC_ENERGIES
/* reduce the energies over warps and store into global memory */
- reduce_energy_warp_shfl(E_lj, E_el, e_lj, e_el, tidx);
+ reduce_energy_warp_shfl(E_lj, E_el, e_lj, e_el, tidx, c_fullWarpMask);
#endif
}
#endif /* FUNCTION_DECLARATION_ONLY */
/*! \brief Log of the i and j cluster size.
* change this together with c_clSize !*/
-static const int c_clSizeLog2 = 3;
+static const int c_clSizeLog2 = 3;
/*! \brief Square of cluster size. */
-static const int c_clSizeSq = c_clSize*c_clSize;
+static const int c_clSizeSq = c_clSize*c_clSize;
/*! \brief j-cluster size after split (4 in the current implementation). */
-static const int c_splitClSize = c_clSize/c_nbnxnGpuClusterpairSplit;
+static const int c_splitClSize = c_clSize/c_nbnxnGpuClusterpairSplit;
/*! \brief Stride in the force accumualation buffer */
-static const int c_fbufStride = c_clSizeSq;
+static const int c_fbufStride = c_clSizeSq;
/*! \brief i-cluster interaction mask for a super-cluster with all c_numClPerSupercl=8 bits set */
-static const unsigned superClInteractionMask = ((1U << c_numClPerSupercl) - 1U);
+static const unsigned superClInteractionMask = ((1U << c_numClPerSupercl) - 1U);
-static const float c_oneSixth = 0.16666667f;
-static const float c_oneTwelveth = 0.08333333f;
+static const float c_oneSixth = 0.16666667f;
+static const float c_oneTwelveth = 0.08333333f;
/*! Convert LJ sigma,epsilon parameters to C6,C12. */
#if GMX_PTX_ARCH >= 300 || GMX_PTX_ARCH == 0
static __forceinline__ __device__
void reduce_force_j_warp_shfl(float3 f, float3 *fout,
- int tidxi, int aidx)
+ int tidxi, int aidx,
+ const unsigned int activemask)
{
- f.x += __shfl_down(f.x, 1);
- f.y += __shfl_up (f.y, 1);
- f.z += __shfl_down(f.z, 1);
+ f.x += gmx_shfl_down_sync(activemask, f.x, 1);
+ f.y += gmx_shfl_up_sync (activemask, f.y, 1);
+ f.z += gmx_shfl_down_sync(activemask, f.z, 1);
if (tidxi & 1)
{
f.x = f.y;
}
- f.x += __shfl_down(f.x, 2);
- f.z += __shfl_up (f.z, 2);
+ f.x += gmx_shfl_down_sync(activemask, f.x, 2);
+ f.z += gmx_shfl_up_sync (activemask, f.z, 2);
if (tidxi & 2)
{
f.x = f.z;
}
- f.x += __shfl_down(f.x, 4);
+ f.x += gmx_shfl_down_sync(activemask, f.x, 4);
if (tidxi < 3)
{
static __forceinline__ __device__
void reduce_force_i_warp_shfl(float3 fin, float3 *fout,
float *fshift_buf, bool bCalcFshift,
- int tidxj, int aidx)
+ int tidxj, int aidx,
+ const unsigned int activemask)
{
- fin.x += __shfl_down(fin.x, c_clSize);
- fin.y += __shfl_up (fin.y, c_clSize);
- fin.z += __shfl_down(fin.z, c_clSize);
+ fin.x += gmx_shfl_down_sync(activemask, fin.x, c_clSize);
+ fin.y += gmx_shfl_up_sync (activemask, fin.y, c_clSize);
+ fin.z += gmx_shfl_down_sync(activemask, fin.z, c_clSize);
if (tidxj & 1)
{
fin.x = fin.y;
}
- fin.x += __shfl_down(fin.x, 2*c_clSize);
- fin.z += __shfl_up (fin.z, 2*c_clSize);
+ fin.x += gmx_shfl_down_sync(activemask, fin.x, 2*c_clSize);
+ fin.z += gmx_shfl_up_sync (activemask, fin.z, 2*c_clSize);
if (tidxj & 2)
{
static __forceinline__ __device__
void reduce_energy_warp_shfl(float E_lj, float E_el,
float *e_lj, float *e_el,
- int tidx)
+ int tidx,
+ const unsigned int activemask)
{
int i, sh;
#pragma unroll 5
for (i = 0; i < 5; i++)
{
- E_lj += __shfl_down(E_lj, sh);
- E_el += __shfl_down(E_el, sh);
+ E_lj += gmx_shfl_down_sync(activemask, E_lj, sh);
+ E_el += gmx_shfl_down_sync(activemask, E_el, sh);
sh += sh;
}
"Checkpointing is merely intended for plain continuation of runs. "
"For safety reasons you must specify all file names (e.g. with -deffnm), "
"and all these files must match the names used in the run prior to checkpointing "
- "since we will append to them by default. If the files are not available, you "
+ "since we will append to them by default. If you used -deffnm and the files listed above as not "
+ "present are in fact present, try explicitly specifying them in respective mdrun options. "
+ "If the files are not available, you "
"can add the -noappend flag to mdrun and write separate new parts. "
"For mere concatenation of files, you should use the gmx trjcat tool instead.",
nfiles-nexist, nfiles);
static inline SimdFloat gmx_simdcall
maskzRsqrt(SimdFloat x, SimdFBool m)
{
+ // The result will always be correct since we mask the result with m, but
+ // for debug builds we also want to make sure not to generate FP exceptions
#ifndef NDEBUG
x.simdInternal_ = vbslq_f32(m.simdInternal_, x.simdInternal_, vdupq_n_f32(1.0f));
#endif
static inline SimdFloat gmx_simdcall
maskzRcp(SimdFloat x, SimdFBool m)
{
+ // The result will always be correct since we mask the result with m, but
+ // for debug builds we also want to make sure not to generate FP exceptions
#ifndef NDEBUG
x.simdInternal_ = vbslq_f32(m.simdInternal_, x.simdInternal_, vdupq_n_f32(1.0f));
#endif
operator<<(SimdFInt32 a, int n)
{
return {
- vshlq_n_s32(a.simdInternal_, n)
+ vshlq_s32(a.simdInternal_, vdupq_n_s32(n >= 32 ? 32 : n))
};
}
operator>>(SimdFInt32 a, int n)
{
return {
- vshrq_n_s32(a.simdInternal_, n)
+ vshlq_s32(a.simdInternal_, vdupq_n_s32(n >= 32 ? -32 : -n))
};
}
operator<<(SimdDInt32 a, int n)
{
return {
- vshl_n_s32(a.simdInternal_, n)
+ vshl_s32(a.simdInternal_, vdup_n_s32(n >= 32 ? 32 : n))
};
}
operator>>(SimdDInt32 a, int n)
{
return {
- vshr_n_s32(a.simdInternal_, n)
+ vshl_s32(a.simdInternal_, vdup_n_s32(n >= 32 ? -32 : -n))
};
}
replExParams);
}
- /* PME tuning is only supported with PME for Coulomb. Is is not supported
- * with only LJ PME, or for reruns.
- */
+ /* PME tuning is only supported in the Verlet scheme, with PME for
+ * Coulomb. It is not supported with only LJ PME, or for
+ * reruns. */
bPMETune = ((Flags & MD_TUNEPME) && EEL_PME(fr->eeltype) && !bRerunMD &&
- !(Flags & MD_REPRODUCIBLE));
+ !(Flags & MD_REPRODUCIBLE) && ir->cutoff_scheme != ecutsGROUP);
if (bPMETune)
{
pme_loadbal_init(&pme_loadbal, cr, mdlog, ir, state->box,
{ "-nstlist", FALSE, etINT, {&nstlist_cmdline},
"Set nstlist when using a Verlet buffer tolerance (0 is guess)" },
{ "-tunepme", FALSE, etBOOL, {&bTunePME},
- "Optimize PME load between PP/PME ranks or GPU/CPU" },
+ "Optimize PME load between PP/PME ranks or GPU/CPU (only with the Verlet cut-off scheme)" },
{ "-v", FALSE, etBOOL, {&bVerbose},
"Be loud and noisy" },
{ "-pforce", FALSE, etREAL, {&pforce},
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff