set(EXPECTED_DOXYGEN_VERSION 1.8.5)
-set(EXPECTED_SPHINX_VERSION 1.4.1)
+set(EXPECTED_SPHINX_VERSION 1.6.1)
if (DEFINED PYTHON_EXECUTABLE)
# Keep quiet on subsequent runs of cmake
set(HTML_OUTPUT_DIR "${CMAKE_CURRENT_BINARY_DIR}/html")
file(MAKE_DIRECTORY ${HTML_OUTPUT_DIR})
+# Prepare directories for pdf/tex output
+set(TEX_OUTPUT_DIR "${CMAKE_CURRENT_BINARY_DIR}/manual")
+set(SPHINX_LATEX_FILE "${TEX_OUTPUT_DIR}/gromacs.tex")
+file(MAKE_DIRECTORY ${TEX_OUTPUT_DIR})
+
# The directory from which man pages will be installed; if it remains
# empty, they will be silently skipped.
set(MAN_PAGE_DIR)
endif()
add_subdirectory(doxygen)
-add_subdirectory(manual)
-
if (SPHINX_FOUND)
# We need to have all the Sphinx input files in a single directory, and
# since some of them are generated, we copy everything into the build tree,
# to this directory.
set(SPHINX_INPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}/sphinx-input)
set(SPHINX_EXTENSION_PATH ${CMAKE_CURRENT_SOURCE_DIR})
+ # As the manual build now depends also on Sphinx, the inclusion path needs
+ # to be set after we know the basic information for Sphinx.
+ add_subdirectory(manual)
if (SOURCE_MD5SUM STREQUAL "unknown")
# But for testing the webpage build (e.g. from the repo) we
# need a default value.
# this should be set to the matching MD5 sum.
set(REGRESSIONTEST_MD5SUM_STRING "${REGRESSIONTEST_MD5SUM}")
endif()
+ # The reference manual still contains the individual sections from the
+ # LaTeX document, with the files below grouped and annotated by chapter.
+ set(REFERENCEMANUAL_SPHINX_FILES_GENERAL
+ # Main index file, preface and introduction.
+ reference-manual/index.rst
+ reference-manual/preface.rst
+ reference-manual/introduction.rst
+ # Definitions and Units chapter
+ reference-manual/definitions.rst
+ # Topologies chapter
+ reference-manual/topologies.rst
+ reference-manual/topologies/particle-type.rst
+ reference-manual/topologies/parameter-files.rst
+ reference-manual/topologies/molecule-definition.rst
+ reference-manual/topologies/constraint-algorithm-section.rst
+ reference-manual/topologies/pdb2gmx-input-files.rst
+ reference-manual/topologies/topology-file-formats.rst
+ reference-manual/topologies/force-field-organization.rst
+ # File formats chapter
+ reference-manual/file-formats.rst
+ # Run parameters and programs chapter
+ reference-manual/run-parameters.rst
+ # Implementation details chapter
+ reference-manual/details.rst
+ # Averages and fluctations chapter
+ reference-manual/averages.rst
+ # References
+ reference-manual/references.rst
+ # PNG formated plot files that don't need to be converted into PNG
+ # for the web page.
+ reference-manual/plots/peregrine.png
+ reference-manual/plots/adress.png
+ reference-manual/plots/plotje.png
+ reference-manual/plots/xvgr.png
+ )
+ set(REFERENCEMANUAL_SPHINX_FILES_WITH_IMAGES
+ # Algorithms chapter
+ reference-manual/algorithms.rst
+ reference-manual/algorithms/periodic-boundary-conditions.rst
+ reference-manual/algorithms/group-concept.rst
+ reference-manual/algorithms/molecular-dynamics.rst
+ reference-manual/algorithms/shell-molecular-dynamics.rst
+ reference-manual/algorithms/constraint-algorithms.rst
+ reference-manual/algorithms/simulated-annealing.rst
+ reference-manual/algorithms/stochastic-dynamics.rst
+ reference-manual/algorithms/brownian-dynamics.rst
+ reference-manual/algorithms/energy-minimization.rst
+ reference-manual/algorithms/normal-mode-analysis.rst
+ reference-manual/algorithms/free-energy-calculations.rst
+ reference-manual/algorithms/replica-exchange.rst
+ reference-manual/algorithms/essential-dynamics.rst
+ reference-manual/algorithms/expanded-ensemble.rst
+ reference-manual/algorithms/parallelization-domain-decomp.rst
+ # Interaction functions and force fields chapter
+ reference-manual/functions.rst
+ reference-manual/functions/bonded-interactions.rst
+ reference-manual/functions/force-field.rst
+ reference-manual/functions/free-energy-interactions.rst
+ reference-manual/functions/interaction-methods.rst
+ reference-manual/functions/long-range-electrostatics.rst
+ reference-manual/functions/long-range-vdw.rst
+ reference-manual/functions/nonbonded-interactions.rst
+ reference-manual/functions/polarization.rst
+ reference-manual/functions/restraints.rst
+ # Special topics chapter
+ reference-manual/special.rst
+ reference-manual/special/free-energy-implementation.rst
+ reference-manual/special/pulling.rst
+ reference-manual/special/awh.rst
+ reference-manual/special/enforced-rotation.rst
+ reference-manual/special/electric-fields.rst
+ reference-manual/special/comp-electrophys.rst
+ reference-manual/special/free-energy-pmf.rst
+ reference-manual/special/remove-fast-dgf.rst
+ reference-manual/special/viscosity-calculation.rst
+ reference-manual/special/tabulated-interaction-functions.rst
+ reference-manual/special/qmmm.rst
+ reference-manual/special/vmd-imd.rst
+ reference-manual/special/membrane-embedding.rst
+ # Analysis chapter
+ reference-manual/analysis.rst
+ reference-manual/analysis/using-groups.rst
+ reference-manual/analysis/looking-at-trajectory.rst
+ reference-manual/analysis/general-properties.rst
+ reference-manual/analysis/radial-distribution-function.rst
+ reference-manual/analysis/correlation-function.rst
+ reference-manual/analysis/curve-fitting.rst
+ reference-manual/analysis/mean-square-displacement.rst
+ reference-manual/analysis/bond-angle-dihedral.rst
+ reference-manual/analysis/radius-of-gyration.rst
+ reference-manual/analysis/rmsd.rst
+ reference-manual/analysis/covariance-analysis.rst
+ reference-manual/analysis/dihedral-pca.rst
+ reference-manual/analysis/protein-related.rst
+ reference-manual/analysis/interface-related.rst)
+ # The image files have also been ordered by the respective
+ # chapter they are included in in the reference manual
+ set(REFERENCEMANUAL_IMAGE_FILES
+ # General folder
+ reference-manual/plots/decomp.pdf
+ reference-manual/plots/dih.pdf
+ reference-manual/plots/drift-all.pdf
+ reference-manual/plots/f-angle.pdf
+ reference-manual/plots/f-bond.pdf
+ reference-manual/plots/fp-highres.pdf
+ reference-manual/plots/int-mat.pdf
+ reference-manual/plots/mdpar.pdf
+ reference-manual/plots/parsort.pdf
+ reference-manual/plots/ring.pdf
+ reference-manual/plots/shiftf.pdf
+ # Algorithms chapter
+ reference-manual/algorithms/plots/dd-cells.pdf
+ reference-manual/algorithms/plots/dd-tric.pdf
+ reference-manual/algorithms/plots/flowchart.pdf
+ reference-manual/algorithms/plots/free1.pdf
+ reference-manual/algorithms/plots/free2.pdf
+ reference-manual/algorithms/plots/leapfrog.pdf
+ reference-manual/algorithms/plots/lincs.pdf
+ reference-manual/algorithms/plots/maxwell.pdf
+ reference-manual/algorithms/plots/mpmd-pme.pdf
+ reference-manual/algorithms/plots/nstric.pdf
+ reference-manual/algorithms/plots/par-lincs2.pdf
+ reference-manual/algorithms/plots/pbctric.pdf
+ reference-manual/algorithms/plots/rhododec.pdf
+ reference-manual/algorithms/plots/truncoct.pdf
+ reference-manual/algorithms/plots/verlet-drift.pdf
+ # Interaction functions chapter
+ reference-manual/functions/plots/angle.pdf
+ reference-manual/functions/plots/bstretch.pdf
+ reference-manual/functions/plots/chain.pdf
+ reference-manual/functions/plots/dummies.pdf
+ reference-manual/functions/plots/f-bham.pdf
+ reference-manual/functions/plots/fbposres.pdf
+ reference-manual/functions/plots/f-dih.pdf
+ reference-manual/functions/plots/f-dr.pdf
+ reference-manual/functions/plots/fig-02.pdf
+ reference-manual/functions/plots/fig-04.pdf
+ reference-manual/functions/plots/f-imps.pdf
+ reference-manual/functions/plots/f-lj.pdf
+ reference-manual/functions/plots/f-morse.pdf
+ reference-manual/functions/plots/f-pr.pdf
+ reference-manual/functions/plots/f-rbs.pdf
+ reference-manual/functions/plots/ring-imp.pdf
+ reference-manual/functions/plots/softcore.pdf
+ reference-manual/functions/plots/subst-im.pdf
+ reference-manual/functions/plots/tetra-im.pdf
+ reference-manual/functions/plots/vcrf.pdf
+ reference-manual/functions/plots/vsite-4fdn.pdf
+ # Special topics chapter
+ reference-manual/special/plots/awh-invN.pdf
+ reference-manual/special/plots/awh-pmfs.pdf
+ reference-manual/special/plots/awh-sampleweights.pdf
+ reference-manual/special/plots/awh-traj.pdf
+ reference-manual/special/plots/compelsetup.pdf
+ reference-manual/special/plots/dumaro.pdf
+ reference-manual/special/plots/dumtypes.pdf
+ reference-manual/special/plots/equipotential.pdf
+ reference-manual/special/plots/field.pdf
+ reference-manual/special/plots/gaussians.pdf
+ reference-manual/special/plots/pulldirrel.pdf
+ reference-manual/special/plots/pull.pdf
+ reference-manual/special/plots/pullref.pdf
+ reference-manual/special/plots/rotation.pdf
+ # Analysis chapter
+ reference-manual/analysis/plots/dih-def.pdf
+ reference-manual/analysis/plots/distm.pdf
+ reference-manual/analysis/plots/dssp.pdf
+ reference-manual/analysis/plots/hbond-insert.pdf
+ reference-manual/analysis/plots/hbond.pdf
+ reference-manual/analysis/plots/hpr-wheel.pdf
+ reference-manual/analysis/plots/msdwater.pdf
+ reference-manual/analysis/plots/ngmxdump.pdf
+ reference-manual/analysis/plots/phipsi.pdf
+ reference-manual/analysis/plots/rama.pdf
+ reference-manual/analysis/plots/rdfO-O.pdf
+ reference-manual/analysis/plots/rdf.pdf
+ reference-manual/analysis/plots/sgangle.pdf
+ )
set(SPHINX_SOURCE_FILES
index.rst
download.rst
- dev-manual/index.rst
+ conf.py
+ links.dat
dev-manual/build-system.rst
+ dev-manual/change-management.rst
dev-manual/commitstyle.rst
dev-manual/documentation-generation.rst
+ dev-manual/contribute.rst
dev-manual/doxygen.rst
dev-manual/error-handling.rst
dev-manual/formatting.rst
dev-manual/gmxtree.rst
dev-manual/includestyle.rst
+ dev-manual/index.rst
dev-manual/jenkins.rst
dev-manual/language-features.rst
dev-manual/naming.rst
dev-manual/overview.rst
+ dev-manual/physical_validation.rst
dev-manual/redmine-states.png
dev-manual/relocatable-binaries.rst
dev-manual/reportstyle.rst
dev-manual/tools.rst
dev-manual/uncrustify.rst
fragments/doxygen-links.rst
+ how-to/index.rst
+ how-to/beginners.rst
+ how-to/topology.rst
+ how-to/special.rst
+ how-to/visualize.rst
install-guide/index.rst
release-notes/index.rst
+ release-notes/2019/major/highlights.rst
+ release-notes/2019/major/features.rst
+ release-notes/2019/major/performance.rst
+ release-notes/2019/major/tools.rst
+ release-notes/2019/major/bugs-fixed.rst
+ release-notes/2019/major/removed-functionality.rst
+ release-notes/2019/major/deprecated-functionality.rst
+ release-notes/2019/major/portability.rst
+ release-notes/2019/major/miscellaneous.rst
release-notes/2018/2018.4.rst
release-notes/2018/2018.3.rst
release-notes/2018/2018.2.rst
release-notes/2016/major/removed-features.rst
release-notes/2016/major/miscellaneous.rst
release-notes/older/index.rst
- user-guide/index.rst
+ # the entry for user-guide/index.rst should not appear here,
+ # as it will be included conditionally further down depending on
+ # if the documentation will be build with the full reference
+ # manual or without.
+ user-guide/cmdline.rst
user-guide/cutoff-schemes.rst
- user-guide/getting-started.rst
- user-guide/force-fields.rst
+ user-guide/deprecation-policy.rst
+ user-guide/environment-variables.rst
user-guide/faq.rst
- user-guide/flow.rst
user-guide/floating-point.rst
- user-guide/system-preparation.rst
+ user-guide/flow.rst
+ user-guide/force-fields.rst
+ user-guide/getting-started.rst
+ user-guide/index.rst
user-guide/managing-simulations.rst
+ user-guide/mdp-options.rst
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/security.rst
+ user-guide/system-preparation.rst
user-guide/terminology.rst
- user-guide/plotje.gif
- user-guide/xvgr.gif
- conf.py
- links.dat
)
include(SphinxMacros.cmake)
gmx_init_sphinx_setup(${SPHINX_INPUT_DIR})
+ # set temporary variables for doi inclusion
+ # into the manual, plain string + some wrapping
+ # for release builds, and dummy string for non-release
+ # builds
+ if("${GMX_MANUAL_DOI}" STREQUAL "")
+ # empty string means no doi, set dummy text
+ set(GMX_MANUAL_DOI_STRING "This is not a release build of GROMACS, so please reference")
+ set(GMX_MANUAL_DOI_STRING "${GMX_MANUAL_DOI_STRING} one of the GROMACS papers and the base release of the manual.")
+ else()
+ # release version, set plain old boring string
+ set(GMX_MANUAL_DOI_STRING "Please reference this documentation as https://doi.org/${GMX_MANUAL_DOI}.")
+ endif()
+ # same for source doi, but modify the text
+ if("${GMX_SOURCE_DOI}" STREQUAL "")
+ # empty string means no release build
+ set(GMX_SOURCE_DOI_STRING "This is not a release build of GROMACS. Please reference one of the")
+ set(GMX_SOURCE_DOI_STRING "${GMX_SOURCE_DOI_STRING} GROMACS papers, as well as the base release that this version is built from.")
+ set(GMX_SOURCE_DOI_STRING "${GMX_SOURCE_DOI_STRING} Also, please state what modifcations have been performed or where the version")
+ set(GMX_SOURCE_DOI_STRING "${GMX_SOURCE_DOI_STRING} was sourced from.")
+ else()
+ # release version, give them a doi url string
+ set(GMX_SOURCE_DOI_STRING "To cite the source code for this release, please cite")
+ set(GMX_SOURCE_DOI_STRING "${GMX_SOURCE_DOI_STRING} https://doi.org/${GMX_SOURCE_DOI}.")
+ endif()
+
+ if(IMAGE_CONVERT_POSSIBLE)
+ set(IMAGE_CONVERT_STRING "possible")
+ else()
+ set(IMAGE_CONVERT_STRING "impossible")
+ endif()
+
set(SPHINX_CONFIG_VARS_FILE ${SPHINX_INPUT_DIR}/conf-vars.py)
gmx_configure_version_file(conf-vars.py.cmakein ${SPHINX_CONFIG_VARS_FILE}
EXTRA_VARS
SPHINX_EXTENSION_PATH RELENG_PATH
+ IMAGE_CONVERT_STRING
EXPECTED_DOXYGEN_VERSION
EXPECTED_SPHINX_VERSION
CMAKE_MINIMUM_REQUIRED_VERSION REQUIRED_CUDA_VERSION
REQUIRED_CUDA_COMPUTE_CAPABILITY REGRESSIONTEST_VERSION
SOURCE_MD5SUM REGRESSIONTEST_MD5SUM_STRING
GMX_TNG_MINIMUM_REQUIRED_VERSION
- GMX_LMFIT_MINIMUM_REQUIRED_VERSION
+ GMX_LMFIT_REQUIRED_VERSION
+ GMX_MANUAL_DOI_STRING
+ GMX_SOURCE_DOI_STRING
COMMENT "Configuring Sphinx configuration file")
gmx_add_sphinx_input_file(${SPHINX_CONFIG_VARS_FILE})
gmx_add_sphinx_source_files(FILES ${SPHINX_SOURCE_FILES})
dev-manual/releng/jenkins-ui.rst
)
endif()
- gmx_add_sphinx_input_target(sphinx-input)
+ gmx_add_sphinx_source_files(
+ FILES
+ ${REFERENCEMANUAL_SPHINX_FILES_GENERAL})
+ if (IMAGE_CONVERT_POSSIBLE)
+ gmx_add_sphinx_source_files(
+ FILES
+ ${REFERENCEMANUAL_SPHINX_FILES_WITH_IMAGES}
+ ${REFERENCEMANUAL_IMAGE_FILES})
+ gmx_add_sphinx_image_conversion_files(
+ FILES
+ ${REFERENCEMANUAL_IMAGE_FILES})
+ endif()
+ gmx_add_sphinx_input_target(sphinx-input-rst)
+ gmx_add_sphinx_image_conversion_target(sphinx-image-conversion)
+ add_custom_target(sphinx-input)
+ add_dependencies(sphinx-input sphinx-input-rst sphinx-image-conversion)
# Remove other rst files from the build tree, since they confuse Sphinx.
# Skip generated files in onlinehelp/, and fragments.
# The latter do not cause issues with obsolete files, as they
add_custom_target(install-guide
COMMAND
${SPHINX_EXECUTABLE}
- -q -b text
+ -q -E -b text
-w sphinx-install.log
-d ${CMAKE_CURRENT_BINARY_DIR}/install-guide/_doctrees
-c ${SPHINX_INPUT_DIR}
# Sphinx cache with pickled ReST documents
set(SPHINX_CACHE_DIR "${CMAKE_CURRENT_BINARY_DIR}/_doctrees")
-
add_custom_target(webpage-sphinx
+ DEPENDS sphinx-programs
+ DEPENDS sphinx-input
+ DEPENDS sphinx-image-conversion
COMMAND
${CMAKE_COMMAND} -E make_directory ${SPHINX_INPUT_DIR}/_static
COMMAND
${SPHINX_EXECUTABLE}
- -q -b html
+ -q -E -b html
-w sphinx-html.log
-d "${SPHINX_CACHE_DIR}"
"${SPHINX_INPUT_DIR}"
COMMENT "Building HTML documentation with Sphinx"
VERBATIM
)
- add_dependencies(webpage-sphinx sphinx-input sphinx-programs)
add_custom_target(man
COMMAND
${SPHINX_EXECUTABLE}
- -q -b man
+ -q -E -b man
-w sphinx-man.log
-d ${SPHINX_CACHE_DIR}
-t do_man
# manually build the target.
set(MAN_PAGE_DIR ${CMAKE_CURRENT_BINARY_DIR})
endif()
+
else()
+ set(MANUAL_BUILD_IS_POSSIBLE OFF)
+ set(MANUAL_BUILD_NOT_POSSIBLE_REASON "Sphinx version ${EXPECTED_SPHINX_VERSION} is not available")
+
add_custom_target(webpage-sphinx
COMMAND ${CMAKE_COMMAND} -E echo
"HTML pages cannot be built because Sphinx version ${EXPECTED_SPHINX_VERSION} is not available"
COMMAND ${CMAKE_COMMAND} -E echo
"man pages cannot be built because Sphinx version ${EXPECTED_SPHINX_VERSION} is not available"
VERBATIM)
+ add_custom_target(sphinx-create-texman
+ COMMAND ${CMAKE_COMMAND} -E echo
+ "Cannot prepare LaTeX input files because Sphinx version ${EXPECTED_SPHINX_VERSION} is not available"
+ VERBATIM)
endif()
if (MAN_PAGE_DIR)
+.. _user guide:
+
**********
User guide
**********
For getting, building and installing |Gromacs|, see the
:doc:`/install-guide/index`.
For background on algorithms and implementations, see the
-`reference manual`_.
+:ref:`reference manual part <gmx-reference-manual-rst>` of the documentation.
+
+|GMX_MANUAL_DOI_STRING|
-.. _reference manual: gmx-manual-parent-dir_
+|GMX_SOURCE_DOI_STRING|
.. TODO This is going to require more organization now that
we are getting more content available.
mdrun-features
mdrun-performance
run-time-errors
- file-formats
+ cmdline
terminology
environment-variables
floating-point
+ security
+ deprecation-policy
:mdp-value:`integrator=tpic` gives identical results to
single-rank :mdp-value:`integrator=tpic`.
+ .. mdp-value:: mimic
+
+ Enable MiMiC QM/MM coupling to run hybrid molecular dynamics.
+ Keey in mind that its required to launch CPMD compiled with MiMiC as well.
+ In this mode all options regarding integration (T-coupling, P-coupling,
+ timestep and number of steps) are ignored as CPMD will do the integration
+ instead. Options related to forces computation (cutoffs, PME parameters,
+ etc.) are working as usual. Atom selection to define QM atoms is read
+ from :mdp:`QMMM-grps`
+
.. mdp:: tinit
- (0) \[ps\]
+ (0) [ps]
starting time for your run (only makes sense for time-based
integrators)
.. mdp:: dt
- (0.001) \[ps\]
+ (0.001) [ps]
time step for integration (only makes sense for time-based
integrators)
.. mdp:: init-step
(0)
- The starting step. The time at an step i in a run is
+ The starting step. The time at step i in a run is
calculated as: t = :mdp:`tinit` + :mdp:`dt` *
(:mdp:`init-step` + i). The free-energy lambda is calculated
as: lambda = :mdp:`init-lambda` + :mdp:`delta-lambda` *
.. mdp-value:: Angular
- Remove center of mass translational and rotational velocity around
- the center of mass
+ Remove center of mass translational and rotational velocity
.. mdp-value:: Linear-acceleration-correction
Remove center of mass translational velocity. Correct the center of
mass position assuming linear acceleration over :mdp:`nstcomm` steps.
This is useful for cases where an acceleration is expected on the
- center of mass which is nearly constant over mdp:`nstcomm` steps.
+ center of mass which is nearly constant over :mdp:`nstcomm` steps.
This can occur for example when pulling on a group using an absolute
reference.
.. mdp:: nstcomm
- (100) \[steps\]
+ (100) [steps]
frequency for center of mass motion removal
.. mdp:: comm-grps
.. mdp:: bd-fric
- (0) \[amu ps-1\]
+ (0) [amu ps\ :sup:`-1`]
Brownian dynamics friction coefficient. When :mdp:`bd-fric` is 0,
the friction coefficient for each particle is calculated as mass/
:mdp:`tau-t`.
.. mdp:: ld-seed
- (-1) \[integer\]
+ (-1) [integer]
used to initialize random generator for thermal noise for
stochastic and Brownian dynamics. When :mdp:`ld-seed` is set to -1,
a pseudo random seed is used. When running BD or SD on multiple
.. mdp:: emtol
- (10.0) \[kJ mol-1 nm-1\]
+ (10.0) [kJ mol\ :sup:`-1` nm\ :sup:`-1`]
the minimization is converged when the maximum force is smaller
than this value
.. mdp:: emstep
- (0.01) \[nm\]
+ (0.01) [nm]
initial step-size
.. mdp:: nstcgsteep
- (1000) \[steps\]
+ (1000) [steps]
frequency of performing 1 steepest descent step while doing
conjugate gradient energy minimization.
.. mdp:: fcstep
- (0) \[ps^2\]
+ (0) [ps\ :sup:`2`]
the step size for optimizing the flexible constraints. Should be
chosen as mu/(d2V/dq2) where mu is the reduced mass of two
particles in a flexible constraint and d2V/dq2 is the second
.. mdp:: rtpi
- (0.05) \[nm\]
+ (0.05) [nm]
the test particle insertion radius, see integrators
:mdp-value:`integrator=tpi` and :mdp-value:`integrator=tpic`
.. mdp:: nstxout
- (0) \[steps\]
- number of steps that elapse between writing coordinates to output
- trajectory file, the last coordinates are always written
+ (0) [steps]
+ number of steps that elapse between writing coordinates to the output
+ trajectory file (:ref:`trr`), the last coordinates are always written
.. mdp:: nstvout
- (0) \[steps\]
- number of steps that elapse between writing velocities to output
- trajectory, the last velocities are always written
+ (0) [steps]
+ number of steps that elapse between writing velocities to the output
+ trajectory file (:ref:`trr`), the last velocities are always written
.. mdp:: nstfout
- (0) \[steps\]
- number of steps that elapse between writing forces to output
- trajectory.
+ (0) [steps]
+ number of steps that elapse between writing forces to the output
+ trajectory file (:ref:`trr`), the last forces are always written.
.. mdp:: nstlog
- (1000) \[steps\]
+ (1000) [steps]
number of steps that elapse between writing energies to the log
file, the last energies are always written
.. mdp:: nstenergy
- (1000) \[steps\]
- number of steps that else between writing energies to energy file,
+ (1000) [steps]
+ number of steps that elapse between writing energies to energy file,
the last energies are always written, should be a multiple of
:mdp:`nstcalcenergy`. Note that the exact sums and fluctuations
over all MD steps modulo :mdp:`nstcalcenergy` are stored in the
.. mdp:: nstxout-compressed
- (0) \[steps\]
+ (0) [steps]
number of steps that elapse between writing position coordinates
- using lossy compression
+ using lossy compression (:ref:`xtc` file)
.. mdp:: compressed-x-precision
- (1000) \[real\]
+ (1000) [real]
precision with which to write to the compressed trajectory file
.. mdp:: compressed-x-grps
.. mdp:: nstlist
- \(10) \[steps\]
+ (10) [steps]
.. mdp-value:: >0
Use no periodic boundary conditions, ignore the box. To simulate
without cut-offs, set all cut-offs and :mdp:`nstlist` to 0. For
best performance without cut-offs on a single MPI rank, set
- :mdp:`nstlist` to zero and :mdp:`ns-type` =simple.
+ :mdp:`nstlist` to zero and :mdp-value:`ns-type=simple`.
.. mdp-value:: xy
Use periodic boundary conditions in x and y directions
- only. This works only with :mdp:`ns-type` =grid and can be used
+ only. This works only with :mdp-value:`ns-type=grid` and can be used
in combination with walls_. Without walls or with only one wall
the system size is infinite in the z direction. Therefore
pressure coupling or Ewald summation methods can not be
.. mdp:: verlet-buffer-tolerance
- (0.005) \[kJ/mol/ps\]
+ (0.005) [kJ mol\ :sup:`-1` ps\ :sup:`-1`]
Useful only with the :mdp-value:`cutoff-scheme=Verlet` :mdp:`cutoff-scheme`. This sets
the maximum allowed error for pair interactions per particle caused
.. mdp:: rlist
- (1) \[nm\]
+ (1) [nm]
Cut-off distance for the short-range neighbor list. With the
:mdp-value:`cutoff-scheme=Verlet` :mdp:`cutoff-scheme`, this is by default set by the
:mdp:`verlet-buffer-tolerance` option and the value of
:mdp:`fourierspacing`. The relative accuracy of
direct/reciprocal space is controlled by :mdp:`ewald-rtol`.
- NOTE: Ewald scales as O(N^3/2) and is thus extremely slow for
+ NOTE: Ewald scales as O(N\ :sup:`3/2`) and is thus extremely slow for
large systems. It is included mainly for reference - in most
cases PME will perform much better.
:mdp:`fourierspacing` and the interpolation order with
:mdp:`pme-order`. With a grid spacing of 0.1 nm and cubic
interpolation the electrostatic forces have an accuracy of
- 2-3*10^-4. Since the error from the vdw-cutoff is larger than
+ 2-3*10\ :sup:`-4`. Since the error from the vdw-cutoff is larger than
this you might try 0.15 nm. When running in parallel the
interpolation parallelizes better than the FFT, so try
decreasing grid dimensions while increasing interpolation.
.. mdp-value:: Reaction-Field-zero
In |Gromacs|, normal reaction-field electrostatics with
- :mdp:`cutoff-scheme` = :mdp-value:`cutoff-scheme=group` leads to bad energy
+ :mdp-value:`cutoff-scheme=group` leads to bad energy
conservation. :mdp-value:`coulombtype=Reaction-Field-zero` solves this by making
the potential zero beyond the cut-off. It can only be used with
an infinite dielectric constant (:mdp:`epsilon-rf` =0), because
only for that value the force vanishes at the
cut-off. :mdp:`rlist` should be 0.1 to 0.3 nm larger than
- :mdp:`rcoulomb` to accommodate for the size of charge groups
+ :mdp:`rcoulomb` to accommodate the size of charge groups
and diffusion between neighbor list updates. This, and the fact
that table lookups are used instead of analytical functions make
- :mdp-value:`coulombtype=Reaction-Field-zero` computationally more expensive than
+ reaction-field-zero computationally more expensive than
normal reaction-field.
.. mdp-value:: Shift
A combination of PME and a switch function for the direct-space
part (see above). :mdp:`rcoulomb` is allowed to be smaller than
:mdp:`rlist`. This is mainly useful constant energy simulations
- (note that using PME with :mdp:`cutoff-scheme` = :mdp-value:`cutoff-scheme=Verlet`
+ (note that using PME with :mdp-value:`cutoff-scheme=Verlet`
will be more efficient).
.. mdp-value:: PME-User
.. mdp:: rcoulomb-switch
- (0) \[nm\]
+ (0) [nm]
where to start switching the Coulomb potential, only relevant
when force or potential switching is used
.. mdp:: rcoulomb
- (1) \[nm\]
+ (1) [nm]
distance for the Coulomb cut-off
.. mdp:: epsilon-r
.. mdp-value:: Shift
- This functionality is deprecated and replaced by
- :mdp:`vdw-modifier` = Force-switch. The LJ (not Buckingham)
- potential is decreased over the whole range and the forces decay
- smoothly to zero between :mdp:`rvdw-switch` and
+ This functionality is deprecated and replaced by using
+ :mdp-value:`vdwtype=Cut-off` with :mdp-value:`vdw-modifier=Force-switch`.
+ The LJ (not Buckingham) potential is decreased over the whole range and
+ the forces decay smoothly to zero between :mdp:`rvdw-switch` and
:mdp:`rvdw`. The neighbor search cut-off :mdp:`rlist` should
- be 0.1 to 0.3 nm larger than :mdp:`rvdw` to accommodate for the
+ be 0.1 to 0.3 nm larger than :mdp:`rvdw` to accommodate the
size of charge groups and diffusion between neighbor list
updates.
.. mdp-value:: Switch
- This functionality is deprecated and replaced by
- :mdp:`vdw-modifier` = Potential-switch. The LJ (not Buckingham)
- potential is normal out to :mdp:`rvdw-switch`, after which it
- is switched off to reach zero at :mdp:`rvdw`. Both the
+ This functionality is deprecated and replaced by using
+ :mdp-value:`vdwtype=Cut-off` with :mdp-value:`vdw-modifier=Potential-switch`.
+ The LJ (not Buckingham) potential is normal out to :mdp:`rvdw-switch`, after
+ which it is switched off to reach zero at :mdp:`rvdw`. Both the
potential and force functions are continuously smooth, but be
aware that all switch functions will give rise to a bulge
(increase) in the force (since we are switching the
potential). The neighbor search cut-off :mdp:`rlist` should be
- 0.1 to 0.3 nm larger than :mdp:`rvdw` to accommodate for the
+ 0.1 to 0.3 nm larger than :mdp:`rvdw` to accommodate the
size of charge groups and diffusion between neighbor list
updates.
.. mdp:: rvdw-switch
- (0) \[nm\]
-
+ (0) [nm]
where to start switching the LJ force and possibly the potential,
only relevant when force or potential switching is used
.. mdp:: rvdw
- (1) \[nm\]
+ (1) [nm]
distance for the LJ or Buckingham cut-off
.. mdp:: DispCorr
.. mdp:: table-extension
- (1) \[nm\]
+ (1) [nm]
Extension of the non-bonded potential lookup tables beyond the
largest cut-off distance. The value should be large enough to
account for charge group sizes and the diffusion between
.. mdp:: fourierspacing
- (0.12) \[nm\]
+ (0.12) [nm]
For ordinary Ewald, the ratio of the box dimensions and the spacing
determines a lower bound for the number of wave vectors to use in
each (signed) direction. For PME and P3M, that ratio determines a
.. mdp:: ewald-rtol
- (1e-5)
+ (10\ :sup:`-5`)
The relative strength of the Ewald-shifted direct potential at
:mdp:`rcoulomb` is given by :mdp:`ewald-rtol`. Decreasing this
will give a more accurate direct sum, but then you need more wave
.. mdp:: ewald-rtol-lj
- (1e-3)
+ (10\ :sup:`-3`)
When doing PME for VdW-interactions, :mdp:`ewald-rtol-lj` is used
to control the relative strength of the dispersion potential at
:mdp:`rvdw` in the same way as :mdp:`ewald-rtol` controls the
.. mdp-value:: berendsen
- Temperature coupling with a Berendsen-thermostat to a bath with
+ Temperature coupling with a Berendsen thermostat to a bath with
temperature :mdp:`ref-t`, with time constant
:mdp:`tau-t`. Several groups can be coupled separately, these
are specified in the :mdp:`tc-grps` field separated by spaces.
but in this case :mdp:`tau-t` controls the period of the
temperature fluctuations at equilibrium, which is slightly
different from a relaxation time. For NVT simulations the
- conserved energy quantity is written to energy and log file.
+ conserved energy quantity is written to the energy and log files.
.. mdp-value:: andersen
- Temperature coupling by randomizing a fraction of the particles
+ Temperature coupling by randomizing a fraction of the particle velocities
at each timestep. Reference temperature and coupling groups are
selected as above. :mdp:`tau-t` is the average time between
randomization of each molecule. Inhibits particle dynamics
.. mdp-value:: andersen-massive
- Temperature coupling by randomizing all particles at infrequent
- timesteps. Reference temperature and coupling groups are
+ Temperature coupling by randomizing velocities of all particles at
+ infrequent timesteps. Reference temperature and coupling groups are
selected as above. :mdp:`tau-t` is the time between
randomization of all molecules. Inhibits particle dynamics
somewhat, but little or no ergodicity issues. Currently only
.. mdp:: tau-t
- \[ps\]
+ [ps]
time constant for coupling (one for each group in
:mdp:`tc-grps`), -1 means no temperature coupling
.. mdp:: ref-t
- \[K\]
+ [K]
reference temperature for coupling (one for each group in
:mdp:`tc-grps`)
.. mdp-value:: Berendsen
Exponential relaxation pressure coupling with time constant
- :mdp:`tau-p`. The box is scaled every timestep. It has been
+ :mdp:`tau-p`. The box is scaled every :mdp:`nstpcouple` steps. It has been
argued that this does not yield a correct thermodynamic
ensemble, but it is the most efficient way to scale a box at the
beginning of a run.
equilibrium. This is probably a better method when you want to
apply pressure scaling during data collection, but beware that
you can get very large oscillations if you are starting from a
- different pressure. For simulations where the exact fluctation
+ different pressure. For simulations where the exact fluctations
of the NPT ensemble are important, or if the pressure coupling
time is very short it may not be appropriate, as the previous
time step pressure is used in some steps of the |Gromacs|
.. mdp:: tau-p
- (1) \[ps\]
+ (1) [ps]
The time constant for pressure coupling (one value for all
directions).
.. mdp:: compressibility
- \[bar^-1\]
- The compressibility (NOTE: this is now really in bar^-1) For water at 1
- atm and 300 K the compressibility is 4.5e-5 bar^-1. The number of
+ [bar\ :sup:`-1`]
+ The compressibility (NOTE: this is now really in bar\ :sup:`-1`) For water at 1
+ atm and 300 K the compressibility is 4.5e-5 bar\ :sup:`-1`. The number of
required values is implied by :mdp:`pcoupltype`.
.. mdp:: ref-p
- \[bar\]
+ [bar]
The reference pressure for coupling. The number of required values
is implied by :mdp:`pcoupltype`.
Generate velocities in :ref:`gmx grompp` according to a
Maxwell distribution at temperature :mdp:`gen-temp`, with
random seed :mdp:`gen-seed`. This is only meaningful with
- integrator :mdp-value:`integrator=md`.
+ :mdp-value:`integrator=md`.
.. mdp:: gen-temp
- (300) \[K\]
+ (300) [K]
temperature for Maxwell distribution
.. mdp:: gen-seed
- (-1) \[integer\]
+ (-1) [integer]
used to initialize random generator for random velocities,
when :mdp:`gen-seed` is set to -1, a pseudo random seed is
used.
.. mdp:: lincs-warnangle
- (30) \[deg\]
+ (30) [deg]
maximum angle that a bond can rotate before LINCS will complain
.. mdp:: morse
(0)
When set to 1 there is a wall at ``z=0``, when set to 2 there is
also a wall at ``z=z-box``. Walls can only be used with :mdp:`pbc`
- ``=xy``. When set to 2 pressure coupling and Ewald summation can be
+ ``=xy``. When set to 2, pressure coupling and Ewald summation can be
used (it is usually best to use semiisotropic pressure coupling
with the ``x/y`` compressibility set to 0, as otherwise the surface
area will change). Walls interact wit the rest of the system
.. mdp:: wall-r-linpot
- (-1) \[nm\]
+ (-1) [nm]
Below this distance from the wall the potential is continued
linearly and thus the force is constant. Setting this option to a
postive value is especially useful for equilibration when some
.. mdp:: wall-density
- \[nm^-3/nm^-2\]
+ [nm\ :sup:`-3`] / [nm\ :sup:`-2`]
the number density of the atoms for each wall for wall types 9-3
and 10-4
COM pulling
^^^^^^^^^^^
-Note that where pulling coordinate are applicable, there can be more
+Note that where pulling coordinates are applicable, there can be more
than one (set with :mdp:`pull-ncoords`) and multiple related :ref:`mdp`
variables will exist accordingly. Documentation references to things
like :mdp:`pull-coord1-vec` should be understood to apply to to the
-applicable pulling coordinate.
+applicable pulling coordinate, eg. the second pull coordinate is described by
+pull-coord2-vec, pull-coord2-k, and so on.
.. mdp:: pull
.. mdp:: pull-cylinder-r
- (1.5) \[nm\]
- the radius of the cylinder for
- :mdp:`pull-coord1-geometry` = :mdp-value:`pull-coord1-geometry=cylinder`
+ (1.5) [nm]
+ the radius of the cylinder for :mdp-value:`pull-coord1-geometry=cylinder`
.. mdp:: pull-constr-tol
- (1e-6)
+ (10\ :sup:`-6`)
the relative constraint tolerance for constraint pulling
.. mdp:: pull-print-com
.. mdp-value:: no
only print the distance for each pull coordinate
-
+
.. mdp-value:: yes
print the distance and Cartesian components selected in
frequency for writing out the force of all the pulled group
(0 is never)
+.. mdp:: pull-pbc-ref-prev-step-com
+
+ .. mdp-value:: no
+
+ Use the reference atom (:mdp:`pull-group1-pbcatom`) for the
+ treatment of periodic boundary conditions.
+
+ .. mdp-value:: yes
+
+ Use the COM of the previous step as reference for the treatment
+ of periodic boundary conditions. The reference is initialized
+ using the reference atom (:mdp:`pull-group1-pbcatom`), which should
+ be located centrally in the group. Using the COM from the
+ previous step can be useful if one or more pull groups are large.
+
+.. mdp:: pull-xout-average
+
+ .. mdp-value:: no
+
+ Write the instantaneous coordinates for all the pulled groups.
+
+ .. mdp-value:: yes
+
+ Write the average coordinates (since last output) for all the
+ pulled groups. N.b., some analysis tools might expect instantaneous
+ pull output.
+
+.. mdp:: pull-fout-average
+
+ .. mdp-value:: no
+
+ Write the instantaneous force for all the pulled groups.
+
+ .. mdp-value:: yes
+
+ Write the average force (since last output) for all the
+ pulled groups. N.b., some analysis tools might expect instantaneous
+ pull output.
.. mdp:: pull-ngroups
vector. For determining the COM, all atoms in the group are put at
their periodic image which is closest to
:mdp:`pull-group1-pbcatom`. A value of 0 means that the middle
- atom (number wise) is used. This parameter is not used with
+ atom (number wise) is used, which is only safe for small groups.
+ :ref:`gmx grompp` checks that the maximum distance from the reference
+ atom (specifically chosen, or not) to the other atoms in the group
+ is not too large. This parameter is not used with
:mdp:`pull-coord1-geometry` cylinder. A value of -1 turns on cosine
weighting, which is useful for a group of molecules in a periodic
system, *e.g.* a water slab (see Engin et al. J. Chem. Phys. B
.. mdp:: pull-coord1-init
- (0.0) \[nm\] / \[deg\]
- The reference distance at t=0.
+ (0.0) [nm] or [deg]
+ The reference distance or reference angle at t=0.
.. mdp:: pull-coord1-rate
- (0) \[nm/ps\] / \[deg/ps\]
- The rate of change of the reference position.
+ (0) [nm/ps] or [deg/ps]
+ The rate of change of the reference position or reference angle.
.. mdp:: pull-coord1-k
- (0) \[kJ mol-1 nm-2\] / \[kJ mol-1 nm-1\] / \[kJ mol-1 rad-2\] / \[kJ mol-1 rad-1\]
+ (0) [kJ mol\ :sup:`-1` nm\ :sup:`-2`] or [kJ mol\ :sup:`-1` nm\ :sup:`-1`] or
+ [kJ mol\ :sup:`-1` rad\ :sup:`-2`] or [kJ mol\ :sup:`-1` rad\ :sup:`-1`]
The force constant. For umbrella pulling this is the harmonic force
- constant in kJ mol-1 nm-2 (or kJ mol-1 rad-2 for angles). For constant force pulling this is the
+ constant in kJ mol\ :sup:`-1` nm\ :sup:`-2` (or kJ mol\ :sup:`-1` rad\ :sup:`-2`
+ for angles). For constant force pulling this is the
force constant of the linear potential, and thus the negative (!)
- of the constant force in kJ mol-1 nm-1 (or kJ mol-1 rad-1 for angles).
+ of the constant force in kJ mol\ :sup:`-1` nm\ :sup:`-1`
+ (or kJ mol\ :sup:`-1` rad\ :sup:`-1` for angles).
Note that for angles the force constant is expressed in terms of radians
(while :mdp:`pull-coord1-init` and :mdp:`pull-coord1-rate` are expressed in degrees).
.. mdp:: pull-coord1-kB
- (pull-k1) \[kJ mol-1 nm-2\] / \[kJ mol-1 nm-1\] / \[kJ mol-1 rad-2\] / \[kJ mol-1 rad-1\]
+ (pull-k1) [kJ mol\ :sup:`-1` nm\ :sup:`-2`] or [kJ mol\ :sup:`-1` nm\ :sup:`-1`]
+ or [kJ mol\ :sup:`-1` rad\ :sup:`-2`] or [kJ mol\ :sup:`-1` rad\ :sup:`-1`]
As :mdp:`pull-coord1-k`, but for state B. This is only used when
:mdp:`free-energy` is turned on. The force constant is then (1 -
lambda) * :mdp:`pull-coord1-k` + lambda * :mdp:`pull-coord1-kB`.
.. mdp:: awh1-error-init
- (10.0) \[kJ mol-1\]
+ (10.0) [kJ mol\ :sup:`-1`]
Estimated initial average error of the PMF for this bias. This value together with the
given diffusion constant(s) :mdp:`awh1-dim1-diffusion` determine the initial biasing rate.
The error is obviously not known *a priori*. Only a rough estimate of :mdp:`awh1-error-init`
.. mdp-value:: constant
The bias is tuned towards a constant (uniform) coordinate distribution
- in the defined sampling interval (defined by \[:mdp:`awh1-dim1-start`, :mdp:`awh1-dim1-end`\]).
+ in the defined sampling interval (defined by [:mdp:`awh1-dim1-start`, :mdp:`awh1-dim1-end`]).
.. mdp-value:: cutoff
.. mdp:: awh1-target-beta-scaling
- [0] \[\]
+ (0)
For :mdp-value:`awh1-target=boltzmann` and :mdp-value:`awh1-target=local-boltzmann`
it is the unitless beta scaling factor taking values in (0,1).
.. mdp:: awh1-target-cutoff
- [0] \[kJ mol-1\]
+ (0) [kJ mol\ :sup:`-1`]
For :mdp-value:`awh1-target=cutoff` this is the cutoff, should be > 0.
.. mdp:: awh1-user-data
.. mdp:: awh1-ndim
- (1) \[integer\]
+ (1) [integer]
Number of dimensions of the coordinate, each dimension maps to 1 pull coordinate.
The following options should be specified for each such dimension. Below only
the options for dimension number 1 is shown. Options for other dimension indices are
.. mdp:: awh1-dim1-force-constant
- (0) \[kJ/mol/nm^2\] or \[kJ/mol/rad^2\]
+ (0) [kJ mol\ :sup:`-1` nm\ :sup:`-2`] or [kJ mol\ :sup:`-1` rad\ :sup:`-2`]
Force constant for the (convolved) umbrella potential(s) along this
coordinate dimension.
.. mdp:: awh1-dim1-start
- (0.0) \[nm\]/\[rad\]
+ (0.0) [nm] or [rad]
Start value of the sampling interval along this dimension. The range of allowed
values depends on the relevant pull geometry (see :mdp:`pull-coord1-geometry`).
For periodic geometries :mdp:`awh1-dim1-start` greater than :mdp:`awh1-dim1-end`
.. mdp:: awh1-dim1-end
- (0.0) \[nm\]/\[rad\]
+ (0.0) [nm] or [rad]
End value defining the sampling interval together with :mdp:`awh1-dim1-start`.
.. mdp:: awh1-dim1-period
- (0.0) \[nm\]/\[rad\]
+ (0.0) [nm] or [rad]
The period of this reaction coordinate, use 0 when the coordinate is not periodic.
.. mdp:: awh1-dim1-diffusion
- (1e-5) \[nm^2/ps\]/\[rad^2/ps\]
+ (10\ :sup:`-5`) [nm\ :sup:`2`/ps] or [rad\ :sup:`2`/ps]
Estimated diffusion constant for this coordinate dimension determining the initial
biasing rate. This needs only be a rough estimate and should not critically
affect the results unless it is set to something very low, leading to slow convergence,
.. mdp:: awh1-dim1-cover-diameter
- (0.0)) \[nm\]/\[rad\]
+ (0.0) [nm] or [rad]
Diameter that needs to be sampled by a single simulation around a coordinate value
before the point is considered covered in the initial stage (see :mdp-value:`awh1-growth=exp-linear`).
A value > 0 ensures that for each covering there is a continuous transition of this diameter
.. mdp:: rot-pivot0
- (0.0 0.0 0.0)
- Pivot point (nm) for the potentials ``iso``, ``pm``, ``rm``, and ``rm2``.
+ (0.0 0.0 0.0) [nm]
+ Pivot point for the potentials ``iso``, ``pm``, ``rm``, and ``rm2``.
.. mdp:: rot-rate0
- (0)
- Reference rotation rate (degree/ps) of group 0.
+ (0) [degree ps\ :sup:`-1`]
+ Reference rotation rate of group 0.
.. mdp:: rot-k0
- (0)
- Force constant (kJ/(mol*nm^2)) for group 0.
+ (0) [kJ mol\ :sup:`-1` nm\ :sup:`-2`]
+ Force constant for group 0.
.. mdp:: rot-slab-dist0
- (1.5)
- Slab distance (nm), if a flexible axis rotation type was chosen.
+ (1.5) [nm]
+ Slab distance, if a flexible axis rotation type was chosen.
.. mdp:: rot-min-gauss0
.. mdp:: rot-eps0
- (0.0001)
- Value of additive constant epsilon' (nm^2) for ``rm2*`` and ``flex2*`` potentials.
+ (0.0001) [nm\ :sup:`2`]
+ Value of additive constant epsilon for ``rm2*`` and ``flex2*`` potentials.
.. mdp:: rot-fit-method0
distance restraints over an ensemble of molecules in one
simulation box. Normally, one would perform ensemble averaging
- over multiple subsystems, each in a separate box, using ``mdrun
- -multi``. Supply ``topol0.tpr``, ``topol1.tpr``, ... with
- different coordinates and/or velocities. The environment
+ over multiple simulations, using ``mdrun
+ -multidir``. The environment
variable ``GMX_DISRE_ENSEMBLE_SIZE`` sets the number of systems
- within each ensemble (usually equal to the ``mdrun -multi``
- value).
+ within each ensemble (usually equal to the number of directories
+ supplied to ``mdrun -multidir``).
.. mdp:: disre-weighting
.. mdp:: disre-fc
- (1000) \[kJ mol-1 nm-2\]
+ (1000) [kJ mol\ :sup:`-1` nm\ :sup:`-2`]
force constant for distance restraints, which is multiplied by a
(possibly) different factor for each restraint given in the `fac`
column of the interaction in the topology file.
.. mdp:: disre-tau
- (0) \[ps\]
+ (0) [ps]
time constant for distance restraints running average. A value of
zero turns off time averaging.
.. mdp:: nstdisreout
- (100) \[steps\]
+ (100) [steps]
period between steps when the running time-averaged and
instantaneous distances of all atom pairs involved in restraints
are written to the energy file (can make the energy file very
.. mdp-value:: yes
use orientation restraints, ensemble averaging can be performed
- with `mdrun -multi`
+ with ``mdrun -multidir``
.. mdp:: orire-fc
- (0) \[kJ mol\]
+ (0) [kJ mol\ :sup:`-1`]
force constant for orientation restraints, which is multiplied by a
(possibly) different weight factor for each restraint, can be set
to zero to obtain the orientations from a free simulation
.. mdp:: orire-tau
- (0) \[ps\]
+ (0) [ps]
time constant for orientation restraints running average. A value
of zero turns off time averaging.
.. mdp:: nstorireout
- (100) \[steps\]
+ (100) [steps]
period between steps when the running time-averaged and
instantaneous orientations for all restraints, and the molecular
order tensor are written to the energy file (can make the energy
.. mdp:: fep-lambdas
- \[array\]
+ [array]
Zero, one or more lambda values for which Delta H values will be
determined and written to dhdl.xvg every :mdp:`nstdhdl`
steps. Values must be between 0 and 1. Free energy differences
.. mdp:: coul-lambdas
- \[array\]
+ [array]
Zero, one or more lambda values for which Delta H values will be
determined and written to dhdl.xvg every :mdp:`nstdhdl`
steps. Values must be between 0 and 1. Only the electrostatic
.. mdp:: vdw-lambdas
- \[array\]
+ [array]
Zero, one or more lambda values for which Delta H values will be
determined and written to dhdl.xvg every :mdp:`nstdhdl`
steps. Values must be between 0 and 1. Only the van der Waals
.. mdp:: bonded-lambdas
- \[array\]
+ [array]
Zero, one or more lambda values for which Delta H values will be
determined and written to dhdl.xvg every :mdp:`nstdhdl`
steps. Values must be between 0 and 1. Only the bonded interactions
.. mdp:: restraint-lambdas
- \[array\]
+ [array]
Zero, one or more lambda values for which Delta H values will be
determined and written to dhdl.xvg every :mdp:`nstdhdl`
steps. Values must be between 0 and 1. Only the restraint
.. mdp:: mass-lambdas
- \[array\]
+ [array]
Zero, one or more lambda values for which Delta H values will be
determined and written to dhdl.xvg every :mdp:`nstdhdl`
steps. Values must be between 0 and 1. Only the particle masses are
.. mdp:: temperature-lambdas
- \[array\]
+ [array]
Zero, one or more lambda values for which Delta H values will be
determined and written to dhdl.xvg every :mdp:`nstdhdl`
steps. Values must be between 0 and 1. Only the temperatures
.. mdp:: sc-sigma
- (0.3) \[nm\]
+ (0.3) [nm]
the soft-core sigma for particles which have a C6 or C12 parameter
equal to zero or a sigma smaller than :mdp:`sc-sigma`
.. mdp:: sim-temp-low
- (300) \[K\]
+ (300) [K]
Low temperature for simulated tempering.
.. mdp:: sim-temp-high
- (300) \[K\]
+ (300) [K]
High temperature for simulated tempering.
.. mdp:: simulated-tempering-scaling
.. mdp:: accelerate
- (0) \[nm ps^-2\]
+ (0) [nm ps\ :sup:`-2`]
acceleration for :mdp:`acc-grps`; x, y and z for each group
(*e.g.* ``0.1 0.0 0.0 -0.1 0.0 0.0`` means that first group has
- constant acceleration of 0.1 nm ps-2 in X direction, second group
+ constant acceleration of 0.1 nm ps\ :sup:`-2` in X direction, second group
the opposite).
.. mdp:: freezegrps
Groups that are to be frozen (*i.e.* their X, Y, and/or Z position
will not be updated; *e.g.* ``Lipid SOL``). :mdp:`freezedim`
- specifies for which dimension the freezing applies. To avoid
- spurious contibrutions to the virial and pressure due to large
+ specifies for which dimension(s) the freezing applies. To avoid
+ spurious contributions to the virial and pressure due to large
forces between completely frozen atoms you need to use energy group
exclusions, this also saves computing time. Note that coordinates
of frozen atoms are not scaled by pressure-coupling algorithms.
.. mdp:: cos-acceleration
- (0) \[nm ps^-2\]
+ (0) [nm ps\ :sup:`-2`]
the amplitude of the acceleration profile for calculating the
viscosity. The acceleration is in the X-direction and the magnitude
is :mdp:`cos-acceleration` cos(2 pi z/boxheight). Two terms are
.. mdp:: deform
- (0 0 0 0 0 0) \[nm ps-1\]
+ (0 0 0 0 0 0) [nm ps\ :sup:`-1`]
The velocities of deformation for the box elements: a(x) b(y) c(z)
b(x) c(x) c(y). Each step the box elements for which :mdp:`deform`
is non-zero are calculated as: box(ts)+(t-ts)*deform, off-diagonal
alternating and pulsed. The general expression for the field
has the form of a gaussian laser pulse:
- E(t) = E0 exp ( -(t-t0)^2/(2 sigma^2) ) cos(omega (t-t0))
+ E(t) = E0 exp ( -(t-t0)\ :sup:`2`/(2 sigma\ :sup:`2`) ) cos(omega (t-t0))
For example, the four parameters for direction x are set in the
- three fields of :mdp:`electric-field-x` (and similar for y and z)
+ three fields of :mdp:`electric-field-x` (and similar for y and z)
like
electric-field-x = E0 omega t0 sigma
electric field is applied.
More details in Carl Caleman and David van der Spoel: Picosecond
- Melting of Ice by an Infrared Laser Pulse - A Simulation Study
+ Melting of Ice by an Infrared Laser Pulse - A Simulation Study.
Angew. Chem. Intl. Ed. 47 pp. 14 17-1420 (2008)
.. mdp:: QMMM-grps
- groups to be descibed at the QM level
+ groups to be descibed at the QM level (works also in case of MiMiC QM/MM)
.. mdp:: QMMMscheme
.. mdp:: QMcharge
- (0) \[integer\]
+ (0) [integer]
The total charge in `e` of the :mdp:`QMMM-grps`. In case there are
more than one :mdp:`QMMM-grps`, the total charge of each ONIOM
layer needs to be specified separately.
.. mdp:: QMmult
- (1) \[integer\]
+ (1) [integer]
The multiplicity of the :mdp:`QMMM-grps`. In case there are more
than one :mdp:`QMMM-grps`, the multiplicity of each ONIOM layer
needs to be specified separately.
.. mdp:: CASorbitals
- (0) \[integer\]
+ (0) [integer]
The number of orbitals to be included in the active space when
doing a CASSCF computation.
.. mdp:: CASelectrons
- (0) \[integer\]
+ (0) [integer]
The number of electrons to be included in the active space when
doing a CASSCF computation.
CASSCF method.
-Implicit solvent
-^^^^^^^^^^^^^^^^
-
-.. mdp:: implicit-solvent
-
- .. mdp-value:: no
-
- No implicit solvent
-
- .. mdp-value:: GBSA
-
- Do a simulation with implicit solvent using the Generalized Born
- formalism. Three different methods for calculating the Born
- radii are available, Still, HCT and OBC. These are specified
- with the :mdp:`gb-algorithm` field. The non-polar solvation is
- specified with the :mdp:`sa-algorithm` field.
-
-.. mdp:: gb-algorithm
-
- .. mdp-value:: Still
-
- Use the Still method to calculate the Born radii
-
- .. mdp-value:: HCT
-
- Use the Hawkins-Cramer-Truhlar method to calculate the Born
- radii
-
- .. mdp-value:: OBC
-
- Use the Onufriev-Bashford-Case method to calculate the Born
- radii
-
-.. mdp:: nstgbradii
-
- (1) \[steps\]
- Frequency to (re)-calculate the Born radii. For most practial
- purposes, setting a value larger than 1 violates energy
- conservation and leads to unstable trajectories.
-
-.. mdp:: rgbradii
-
- (1.0) \[nm\]
- Cut-off for the calculation of the Born radii. Currently must be
- equal to rlist
-
-.. mdp:: gb-epsilon-solvent
-
- (80)
- Dielectric constant for the implicit solvent
-
-.. mdp:: gb-saltconc
-
- (0) \[M\]
- Salt concentration for implicit solvent models, currently not used
-
-.. mdp:: gb-obc-alpha
-.. mdp:: gb-obc-beta
-.. mdp:: gb-obc-gamma
-
- Scale factors for the OBC model. Default values of 1, 0.78 and 4.85
- respectively are for OBC(II). Values for OBC(I) are 0.8, 0 and 2.91
- respectively
-
-.. mdp:: gb-dielectric-offset
-
- (0.009) \[nm\]
- Distance for the di-electric offset when calculating the Born
- radii. This is the offset between the center of each atom the
- center of the polarization energy for the corresponding atom
-
-.. mdp:: sa-algorithm
-
- .. mdp-value:: Ace-approximation
-
- Use an Ace-type approximation
-
- .. mdp-value:: None
-
- No non-polar solvation calculation done. For GBSA only the polar
- part gets calculated
-
-.. mdp:: sa-surface-tension
-
- \[kJ mol-1 nm-2\]
- Default value for surface tension with SA algorithms. The default
- value is -1; Note that if this default value is not changed it will
- be overridden by :ref:`gmx grompp` using values that are specific
- for the choice of radii algorithm (0.0049 kcal/mol/Angstrom^2 for
- Still, 0.0054 kcal/mol/Angstrom2 for HCT/OBC) Setting it to 0 will
- while using an sa-algorithm other than None means no non-polar
- calculations are done.
-
-
Computational Electrophysiology
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Use these options to switch on and control ion/water position exchanges in "Computational
.. mdp:: coupl-steps
- (\10) Average the number of ions per compartment over these many swap attempt steps.
+ (10) Average the number of ions per compartment over these many swap attempt steps.
This can be used to prevent that ions near a compartment boundary
(diffusing through a channel, e.g.) lead to unwanted back and forth swaps.
.. mdp:: cyl0-r
- (2.0) \[nm\] Radius of the split cylinder #0.
+ (2.0) [nm] Radius of the split cylinder #0.
Two split cylinders (mimicking the channel pores) can optionally be defined
relative to the center of the split group. With the help of these cylinders
it can be counted which ions have passed which channel. The split cylinder
.. mdp:: cyl0-up
- (1.0) \[nm\] Upper extension of the split cylinder #0.
+ (1.0) [nm] Upper extension of the split cylinder #0.
.. mdp:: cyl0-down
- (1.0) \[nm\] Lower extension of the split cylinder #0.
+ (1.0) [nm] Lower extension of the split cylinder #0.
.. mdp:: cyl1-r
- (2.0) \[nm\] Radius of the split cylinder #1.
+ (2.0) [nm] Radius of the split cylinder #1.
.. mdp:: cyl1-up
- (1.0) \[nm\] Upper extension of the split cylinder #1.
+ (1.0) [nm] Upper extension of the split cylinder #1.
.. mdp:: cyl1-down
- (1.0) \[nm\] Lower extension of the split cylinder #1.
+ (1.0) [nm] Lower extension of the split cylinder #1.
User defined thingies
Removed features
^^^^^^^^^^^^^^^^
-This feature has been removed from |Gromacs|, but so that old
+These features have been removed from |Gromacs|, but so that old
:ref:`mdp` and :ref:`tpr` files cannot be mistakenly misused, we still
parse this option. :ref:`gmx grompp` and :ref:`gmx mdrun` will issue a
fatal error if this is set.
(no)
+.. mdp:: implicit-solvent
+
+ (no)
+
.. _reference manual: gmx-manual-parent-dir_
#include "config.h"
-#include <assert.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
+#include <cassert>
#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
#include <algorithm>
+#include <list>
+#include "gromacs/domdec/domdec.h"
#include "gromacs/ewald/ewald-utils.h"
#include "gromacs/fft/parallel_3dfft.h"
#include "gromacs/fileio/pdbio.h"
#include "gromacs/timing/cyclecounter.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/topology/topology.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "pme-spline-work.h"
#include "pme-spread.h"
+/*! \brief Help build a descriptive message in \c error if there are
+ * \c errorReasons why PME on GPU is not supported.
+ *
+ * \returns Whether the lack of errorReasons indicate there is support. */
+static bool
+addMessageIfNotSupported(const std::list<std::string> &errorReasons,
+ std::string *error)
+{
+ bool foundErrorReasons = errorReasons.empty();
+ if (!foundErrorReasons && error)
+ {
+ std::string regressionTestMarker = "PME GPU does not support";
+ // this prefix is tested for in the regression tests script gmxtest.pl
+ *error = regressionTestMarker + ": " + gmx::joinStrings(errorReasons, "; ") + ".";
+ }
+ return foundErrorReasons;
+}
+
+bool pme_gpu_supports_build(std::string *error)
+{
+ std::list<std::string> errorReasons;
+ if (GMX_DOUBLE)
+ {
+ errorReasons.emplace_back("double precision");
+ }
+ if (GMX_GPU == GMX_GPU_NONE)
+ {
+ errorReasons.emplace_back("non-GPU build of GROMACS");
+ }
+ return addMessageIfNotSupported(errorReasons, error);
+}
+
+bool pme_gpu_supports_input(const t_inputrec &ir, const gmx_mtop_t &mtop, std::string *error)
+{
+ std::list<std::string> errorReasons;
+ if (!EEL_PME(ir.coulombtype))
+ {
+ errorReasons.emplace_back("systems that do not use PME for electrostatics");
+ }
+ if (ir.pme_order != 4)
+ {
+ errorReasons.emplace_back("interpolation orders other than 4");
+ }
+ if (ir.efep != efepNO)
+ {
+ if (gmx_mtop_has_perturbed_charges(mtop))
+ {
+ errorReasons.emplace_back("free energy calculations with perturbed charges (multiple grids)");
+ }
+ }
+ if (EVDW_PME(ir.vdwtype))
+ {
+ errorReasons.emplace_back("Lennard-Jones PME");
+ }
+ if (ir.cutoff_scheme == ecutsGROUP)
+ {
+ errorReasons.emplace_back("group cutoff scheme");
+ }
+ if (!EI_DYNAMICS(ir.eI))
+ {
+ errorReasons.emplace_back("not a dynamical integrator");
+ }
+ return addMessageIfNotSupported(errorReasons, error);
+}
+
+/*! \brief \libinternal
+ * Finds out if PME with given inputs is possible to run on GPU.
+ * This function is an internal final check, validating the whole PME structure on creation,
+ * but it still duplicates the preliminary checks from the above (externally exposed) pme_gpu_supports_input() - just in case.
+ *
+ * \param[in] pme The PME structure.
+ * \param[out] error The error message if the input is not supported on GPU.
+ * \returns True if this PME input is possible to run on GPU, false otherwise.
+ */
+static bool pme_gpu_check_restrictions(const gmx_pme_t *pme, std::string *error)
+{
+ std::list<std::string> errorReasons;
+ if (pme->nnodes != 1)
+ {
+ errorReasons.emplace_back("PME decomposition");
+ }
+ if (pme->pme_order != 4)
+ {
+ errorReasons.emplace_back("interpolation orders other than 4");
+ }
+ if (pme->bFEP)
+ {
+ errorReasons.emplace_back("free energy calculations (multiple grids)");
+ }
+ if (pme->doLJ)
+ {
+ errorReasons.emplace_back("Lennard-Jones PME");
+ }
+ if (GMX_DOUBLE)
+ {
+ errorReasons.emplace_back("double precision");
+ }
+ if (GMX_GPU == GMX_GPU_NONE)
+ {
+ errorReasons.emplace_back("non-GPU build of GROMACS");
+ }
+
+ return addMessageIfNotSupported(errorReasons, error);
+}
+
+PmeRunMode pme_run_mode(const gmx_pme_t *pme)
+{
+ GMX_ASSERT(pme != nullptr, "Expecting valid PME data pointer");
+ return pme->runMode;
+}
+
+gmx::PinningPolicy pme_get_pinning_policy()
+{
+ return gmx::PinningPolicy::PinnedIfSupported;
+}
+
/*! \brief Number of bytes in a cache line.
*
* Must also be a multiple of the SIMD and SIMD4 register size, to
/* pme_solve is roughly double the cost of an fft */
- return (n1 + n2 + 3*n3)/(double)(6*pme->nkx*pme->nky*pme->nkz);
+ return (n1 + n2 + 3*n3)/static_cast<double>(6*pme->nkx*pme->nky*pme->nkz);
}
/*! \brief Initialize atom communication data structure */
int ndata,
int commplainsize)
{
- int b, i;
- pme_grid_comm_t *pgc;
gmx_bool bCont;
- int fft_start, fft_end, send_index1, recv_index1;
#if GMX_MPI
MPI_Status stat;
* that belong to higher nodes (modulo nnodes)
*/
- snew(ol->s2g0, ol->nnodes+1);
- snew(ol->s2g1, ol->nnodes);
+ ol->s2g0.resize(ol->nnodes + 1);
+ ol->s2g1.resize(ol->nnodes);
if (debug)
{
fprintf(debug, "PME slab boundaries:");
}
- for (i = 0; i < nnodes; i++)
+ for (int i = 0; i < nnodes; i++)
{
/* s2g0 the local interpolation grid start.
* s2g1 the local interpolation grid end.
* spatially uniform along dimension x or y, we need to round
* s2g0 down and s2g1 up.
*/
- ol->s2g0[i] = ( i *ndata + 0 )/nnodes;
- ol->s2g1[i] = ((i+1)*ndata + nnodes-1)/nnodes + norder - 1;
+ ol->s2g0[i] = (i * ndata + 0) / nnodes;
+ ol->s2g1[i] = ((i + 1) * ndata + nnodes - 1) / nnodes + norder - 1;
if (debug)
{
}
/* Determine with how many nodes we need to communicate the grid overlap */
- b = 0;
+ int testRankCount = 0;
do
{
- b++;
+ testRankCount++;
bCont = FALSE;
- for (i = 0; i < nnodes; i++)
+ for (int i = 0; i < nnodes; i++)
{
- if ((i+b < nnodes && ol->s2g1[i] > ol->s2g0[i+b]) ||
- (i+b >= nnodes && ol->s2g1[i] > ol->s2g0[i+b-nnodes] + ndata))
+ if ((i + testRankCount < nnodes && ol->s2g1[i] > ol->s2g0[i + testRankCount]) ||
+ (i + testRankCount >= nnodes && ol->s2g1[i] > ol->s2g0[i + testRankCount - nnodes] + ndata))
{
bCont = TRUE;
}
}
}
- while (bCont && b < nnodes);
- ol->noverlap_nodes = b - 1;
-
- snew(ol->send_id, ol->noverlap_nodes);
- snew(ol->recv_id, ol->noverlap_nodes);
- for (b = 0; b < ol->noverlap_nodes; b++)
- {
- ol->send_id[b] = (ol->nodeid + (b + 1)) % ol->nnodes;
- ol->recv_id[b] = (ol->nodeid - (b + 1) + ol->nnodes) % ol->nnodes;
- }
- snew(ol->comm_data, ol->noverlap_nodes);
+ while (bCont && testRankCount < nnodes);
+ ol->comm_data.resize(testRankCount - 1);
ol->send_size = 0;
- for (b = 0; b < ol->noverlap_nodes; b++)
+
+ for (size_t b = 0; b < ol->comm_data.size(); b++)
{
- pgc = &ol->comm_data[b];
+ pme_grid_comm_t *pgc = &ol->comm_data[b];
+
/* Send */
- fft_start = ol->s2g0[ol->send_id[b]];
- fft_end = ol->s2g0[ol->send_id[b]+1];
- if (ol->send_id[b] < nodeid)
+ pgc->send_id = (ol->nodeid + (b + 1)) % ol->nnodes;
+ int fft_start = ol->s2g0[pgc->send_id];
+ int fft_end = ol->s2g0[pgc->send_id + 1];
+ if (pgc->send_id < nodeid)
{
fft_start += ndata;
fft_end += ndata;
}
- send_index1 = ol->s2g1[nodeid];
- send_index1 = std::min(send_index1, fft_end);
- pgc->send_index0 = fft_start;
- pgc->send_nindex = std::max(0, send_index1 - pgc->send_index0);
- ol->send_size += pgc->send_nindex;
+ int send_index1 = ol->s2g1[nodeid];
+ send_index1 = std::min(send_index1, fft_end);
+ pgc->send_index0 = fft_start;
+ pgc->send_nindex = std::max(0, send_index1 - pgc->send_index0);
+ ol->send_size += pgc->send_nindex;
/* We always start receiving to the first index of our slab */
+ pgc->recv_id = (ol->nodeid - (b + 1) + ol->nnodes) % ol->nnodes;
fft_start = ol->s2g0[ol->nodeid];
- fft_end = ol->s2g0[ol->nodeid+1];
- recv_index1 = ol->s2g1[ol->recv_id[b]];
- if (ol->recv_id[b] > nodeid)
+ fft_end = ol->s2g0[ol->nodeid + 1];
+ int recv_index1 = ol->s2g1[pgc->recv_id];
+ if (pgc->recv_id > nodeid)
{
recv_index1 -= ndata;
}
#if GMX_MPI
/* Communicate the buffer sizes to receive */
- for (b = 0; b < ol->noverlap_nodes; b++)
+ for (size_t b = 0; b < ol->comm_data.size(); b++)
{
- MPI_Sendrecv(&ol->send_size, 1, MPI_INT, ol->send_id[b], b,
- &ol->comm_data[b].recv_size, 1, MPI_INT, ol->recv_id[b], b,
+ MPI_Sendrecv(&ol->send_size, 1, MPI_INT, ol->comm_data[b].send_id, b,
+ &ol->comm_data[b].recv_size, 1, MPI_INT, ol->comm_data[b].recv_id, b,
ol->mpi_comm, &stat);
}
#endif
/* For non-divisible grid we need pme_order iso pme_order-1 */
- snew(ol->sendbuf, norder*commplainsize);
- snew(ol->recvbuf, norder*commplainsize);
-}
-
-/*! \brief Destroy data structure for communication */
-static void
-destroy_overlap_comm(const pme_overlap_t *ol)
-{
- sfree(ol->s2g0);
- sfree(ol->s2g1);
- sfree(ol->send_id);
- sfree(ol->recv_id);
- sfree(ol->comm_data);
- sfree(ol->sendbuf);
- sfree(ol->recvbuf);
+ ol->sendbuf.resize(norder * commplainsize);
+ ol->recvbuf.resize(norder * commplainsize);
}
int minimalPmeGridSize(int pmeOrder)
bool gmx_pme_check_restrictions(int pme_order,
int nkx, int nky, int nkz,
- int nnodes_major,
+ int numPmeDomainsAlongX,
bool useThreads,
bool errorsAreFatal)
{
std::string message = gmx::formatString(
"pme_order (%d) is larger than the maximum allowed value (%d). Modify and recompile the code if you really need such a high order.",
pme_order, PME_ORDER_MAX);
- GMX_THROW(InconsistentInputError(message));
+ GMX_THROW(gmx::InconsistentInputError(message));
}
const int minGridSize = minimalPmeGridSize(pme_order);
std::string message = gmx::formatString(
"The PME grid sizes need to be >= 2*(pme_order-1) (%d)",
minGridSize);
- GMX_THROW(InconsistentInputError(message));
+ GMX_THROW(gmx::InconsistentInputError(message));
}
/* Check for a limitation of the (current) sum_fftgrid_dd code.
* We only allow multiple communication pulses in dim 1, not in dim 0.
*/
- if (useThreads && (nkx < nnodes_major*pme_order &&
- nkx != nnodes_major*(pme_order - 1)))
+ if (useThreads && (nkx < numPmeDomainsAlongX*pme_order &&
+ nkx != numPmeDomainsAlongX*(pme_order - 1)))
{
if (!errorsAreFatal)
{
return false;
}
gmx_fatal(FARGS, "The number of PME grid lines per rank along x is %g. But when using OpenMP threads, the number of grid lines per rank along x should be >= pme_order (%d) or = pmeorder-1. To resolve this issue, use fewer ranks along x (and possibly more along y and/or z) by specifying -dd manually.",
- nkx/(double)nnodes_major, pme_order);
+ nkx/static_cast<double>(numPmeDomainsAlongX), pme_order);
}
return true;
return (enumerator + denominator - 1)/denominator;
}
-gmx_pme_t *gmx_pme_init(const t_commrec *cr,
- int nnodes_major,
- int nnodes_minor,
- const t_inputrec *ir,
- int homenr,
- gmx_bool bFreeEnergy_q,
- gmx_bool bFreeEnergy_lj,
- gmx_bool bReproducible,
- real ewaldcoeff_q,
- real ewaldcoeff_lj,
- int nthread,
- PmeRunMode runMode,
- PmeGpu *pmeGpu,
- gmx_device_info_t *gpuInfo,
- const gmx::MDLogger & /*mdlog*/)
+gmx_pme_t *gmx_pme_init(const t_commrec *cr,
+ const NumPmeDomains &numPmeDomains,
+ const t_inputrec *ir,
+ int homenr,
+ gmx_bool bFreeEnergy_q,
+ gmx_bool bFreeEnergy_lj,
+ gmx_bool bReproducible,
+ real ewaldcoeff_q,
+ real ewaldcoeff_lj,
+ int nthread,
+ PmeRunMode runMode,
+ PmeGpu *pmeGpu,
+ const gmx_device_info_t *gpuInfo,
+ PmeGpuProgramHandle pmeGpuProgram,
+ const gmx::MDLogger & /*mdlog*/)
{
int use_threads, sum_use_threads, i;
ivec ndata;
fprintf(debug, "Creating PME data structures.\n");
}
- unique_cptr<gmx_pme_t, gmx_pme_destroy> pme(new gmx_pme_t());
+ gmx::unique_cptr<gmx_pme_t, gmx_pme_destroy> pme(new gmx_pme_t());
pme->sum_qgrid_tmp = nullptr;
pme->sum_qgrid_dd_tmp = nullptr;
pme->nnodes = 1;
pme->bPPnode = TRUE;
- pme->nnodes_major = nnodes_major;
- pme->nnodes_minor = nnodes_minor;
+ pme->nnodes_major = numPmeDomains.x;
+ pme->nnodes_minor = numPmeDomains.y;
#if GMX_MPI
- if (nnodes_major*nnodes_minor > 1)
+ if (numPmeDomains.x*numPmeDomains.y > 1)
{
pme->mpi_comm = cr->mpi_comm_mygroup;
MPI_Comm_rank(pme->mpi_comm, &pme->nodeid);
MPI_Comm_size(pme->mpi_comm, &pme->nnodes);
- if (pme->nnodes != nnodes_major*nnodes_minor)
+ if (pme->nnodes != numPmeDomains.x*numPmeDomains.y)
{
gmx_incons("PME rank count mismatch");
}
}
else
{
- if (nnodes_minor == 1)
+ if (numPmeDomains.y == 1)
{
#if GMX_MPI
pme->mpi_comm_d[0] = pme->mpi_comm;
pme->nodeid_minor = 0;
}
- else if (nnodes_major == 1)
+ else if (numPmeDomains.x == 1)
{
#if GMX_MPI
pme->mpi_comm_d[0] = MPI_COMM_NULL;
}
else
{
- if (pme->nnodes % nnodes_major != 0)
+ if (pme->nnodes % numPmeDomains.x != 0)
{
- gmx_incons("For 2D PME decomposition, #PME ranks must be divisible by the number of ranks in the major dimension");
+ gmx_incons("For 2D PME decomposition, #PME ranks must be divisible by the number of domains along x");
}
pme->ndecompdim = 2;
#if GMX_MPI
- MPI_Comm_split(pme->mpi_comm, pme->nodeid % nnodes_minor,
+ MPI_Comm_split(pme->mpi_comm, pme->nodeid % numPmeDomains.y,
pme->nodeid, &pme->mpi_comm_d[0]); /* My communicator along major dimension */
- MPI_Comm_split(pme->mpi_comm, pme->nodeid/nnodes_minor,
+ MPI_Comm_split(pme->mpi_comm, pme->nodeid/numPmeDomains.y,
pme->nodeid, &pme->mpi_comm_d[1]); /* My communicator along minor dimension */
MPI_Comm_rank(pme->mpi_comm_d[0], &pme->nodeid_major);
" PME grid_x (%d) and grid_y (%d) should be divisible by #PME_ranks_x (%d)\n"
" and PME grid_y (%d) and grid_z (%d) should be divisible by #PME_ranks_y (%d)\n"
"\n",
- (int)((imbal-1)*100 + 0.5),
+ gmx::roundToInt((imbal-1)*100),
pme->nkx, pme->nky, pme->nnodes_major,
pme->nky, pme->nkz, pme->nnodes_minor);
}
/* Double-check for a limitation of the (current) sum_fftgrid_dd code.
* Note that gmx_pme_check_restrictions checked for this already.
*/
- if (pme->bUseThreads && pme->overlap[0].noverlap_nodes > 1)
+ if (pme->bUseThreads && (pme->overlap[0].comm_data.size() > 1))
{
gmx_incons("More than one communication pulse required for grid overlap communication along the major dimension while using threads");
}
pme->overlap[0].s2g1[pme->nodeid_major]-pme->overlap[0].s2g0[pme->nodeid_major+1],
pme->overlap[1].s2g1[pme->nodeid_minor]-pme->overlap[1].s2g0[pme->nodeid_minor+1]);
/* This routine will allocate the grid data to fit the FFTs */
- const auto allocateRealGridForGpu = (pme->runMode == PmeRunMode::Mixed) ? gmx::PinningPolicy::CanBePinned : gmx::PinningPolicy::CannotBePinned;
+ const auto allocateRealGridForGpu = (pme->runMode == PmeRunMode::Mixed) ? gmx::PinningPolicy::PinnedIfSupported : gmx::PinningPolicy::CannotBePinned;
gmx_parallel_3dfft_init(&pme->pfft_setup[i], ndata,
&pme->fftgrid[i], &pme->cfftgrid[i],
pme->mpi_comm_d,
}
/* Use atc[0] for spreading */
- init_atomcomm(pme.get(), &pme->atc[0], nnodes_major > 1 ? 0 : 1, TRUE);
+ init_atomcomm(pme.get(), &pme->atc[0], numPmeDomains.x > 1 ? 0 : 1, TRUE);
if (pme->ndecompdim >= 2)
{
init_atomcomm(pme.get(), &pme->atc[1], 1, FALSE);
pme->lb_buf2 = nullptr;
pme->lb_buf_nalloc = 0;
- pme_gpu_reinit(pme.get(), gpuInfo);
+ if (pme_gpu_active(pme.get()))
+ {
+ if (!pme->gpu)
+ {
+ // Initial check of validity of the data
+ std::string errorString;
+ bool canRunOnGpu = pme_gpu_check_restrictions(pme.get(), &errorString);
+ if (!canRunOnGpu)
+ {
+ GMX_THROW(gmx::NotImplementedError(errorString));
+ }
+ }
+
+ pme_gpu_reinit(pme.get(), gpuInfo, pmeGpuProgram);
+ }
pme_init_all_work(&pme->solve_work, pme->nthread, pme->nkx);
}
void gmx_pme_reinit(struct gmx_pme_t **pmedata,
- t_commrec * cr,
+ const t_commrec *cr,
struct gmx_pme_t * pme_src,
const t_inputrec * ir,
const ivec grid_size,
// so we don't expect the actual logging.
// TODO: when PME is an object, it should take reference to mdlog on construction and save it.
GMX_ASSERT(pmedata, "Invalid PME pointer");
- *pmedata = gmx_pme_init(cr, pme_src->nnodes_major, pme_src->nnodes_minor,
+ NumPmeDomains numPmeDomains = { pme_src->nnodes_major, pme_src->nnodes_minor };
+ *pmedata = gmx_pme_init(cr, numPmeDomains,
&irc, homenr, pme_src->bFEP_q, pme_src->bFEP_lj, FALSE, ewaldcoeff_q, ewaldcoeff_lj,
- pme_src->nthread, pme_src->runMode, pme_src->gpu, nullptr, dummyLogger);
+ pme_src->nthread, pme_src->runMode, pme_src->gpu, nullptr, nullptr, dummyLogger);
//TODO this is mostly passing around current values
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
{
gmx_incons("gmx_pme_calc_energy called in parallel");
}
- if (pme->bFEP_q > 1)
+ if (pme->bFEP_q)
{
gmx_incons("gmx_pme_calc_energy with free energy");
}
/*! \brief Calculate initial Lorentz-Berthelot coefficients for LJ-PME */
static void
-calc_initial_lb_coeffs(struct gmx_pme_t *pme, real *local_c6, real *local_sigma)
+calc_initial_lb_coeffs(struct gmx_pme_t *pme, const real *local_c6, const real *local_sigma)
{
int i;
for (i = 0; i < pme->atc[0].n; ++i)
/*! \brief Calculate next Lorentz-Berthelot coefficients for LJ-PME */
static void
-calc_next_lb_coeffs(struct gmx_pme_t *pme, real *local_sigma)
+calc_next_lb_coeffs(struct gmx_pme_t *pme, const real *local_sigma)
{
int i;
real chargeA[], real chargeB[],
real c6A[], real c6B[],
real sigmaA[], real sigmaB[],
- matrix box, t_commrec *cr,
+ matrix box, const t_commrec *cr,
int maxshift_x, int maxshift_y,
- t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ t_nrnb *nrnb, gmx_wallcycle *wcycle,
matrix vir_q, matrix vir_lj,
real *energy_q, real *energy_lj,
real lambda_q, real lambda_lj,
gmx_bool bFirst, bDoSplines;
int fep_state;
int fep_states_lj = pme->bFEP_lj ? 2 : 1;
- const gmx_bool bCalcEnerVir = flags & GMX_PME_CALC_ENER_VIR;
- const gmx_bool bBackFFT = flags & (GMX_PME_CALC_F | GMX_PME_CALC_POT);
- const gmx_bool bCalcF = flags & GMX_PME_CALC_F;
+ const gmx_bool bCalcEnerVir = (flags & GMX_PME_CALC_ENER_VIR) != 0;
+ const gmx_bool bBackFFT = (flags & (GMX_PME_CALC_F | GMX_PME_CALC_POT)) != 0;
+ const gmx_bool bCalcF = (flags & GMX_PME_CALC_F) != 0;
+ /* We could be passing lambda!=1 while no q or LJ is actually perturbed */
+ if (!pme->bFEP_q)
+ {
+ lambda_q = 1;
+ }
+ if (!pme->bFEP_lj)
+ {
+ lambda_lj = 1;
+ }
+
assert(pme->nnodes > 0);
assert(pme->nnodes == 1 || pme->ndecompdim > 0);
{
fprintf(debug, "PME: number of ranks = %d, rank = %d\n",
cr->nnodes, cr->nodeid);
- fprintf(debug, "Grid = %p\n", (void*)grid);
+ fprintf(debug, "Grid = %p\n", static_cast<void*>(grid));
if (grid == nullptr)
{
gmx_fatal(FARGS, "No grid!");
}
}
- where();
if (pme->nnodes == 1)
{
{
wallcycle_start(wcycle, ewcPME_REDISTXF);
do_redist_pos_coeffs(pme, cr, start, homenr, bFirst, x, coefficient);
- where();
wallcycle_stop(wcycle, ewcPME_REDISTXF);
}
if (pme->nnodes > 1)
{
gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_FORWARD);
- where();
}
#endif
{
wallcycle_stop(wcycle, ewcPME_FFT);
}
- where();
/* solve in k-space for our local cells */
if (thread == 0)
if (thread == 0)
{
wallcycle_stop(wcycle, (grid_index < DO_Q ? ewcPME_SOLVE : ewcLJPME));
- where();
inc_nrnb(nrnb, eNR_SOLVEPME, loop_count);
}
}
/* do 3d-invfft */
if (thread == 0)
{
- where();
wallcycle_start(wcycle, ewcPME_FFT);
}
gmx_parallel_3dfft_execute(pfft_setup, GMX_FFT_COMPLEX_TO_REAL,
{
wallcycle_stop(wcycle, ewcPME_FFT);
- where();
if (pme->nodeid == 0)
{
gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_BACKWARD);
}
#endif
- where();
unwrap_periodic_pmegrid(pme, grid);
}
{
/* interpolate forces for our local atoms */
- where();
/* If we are running without parallelization,
* atc->f is the actual force array, not a buffer,
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
}
- where();
inc_nrnb(nrnb, eNR_GATHERFBSP,
pme->pme_order*pme->pme_order*pme->pme_order*pme->atc[0].n);
{
local_c6[i] = atc->coefficient[i];
}
- where();
do_redist_pos_coeffs(pme, cr, start, homenr, FALSE, x, RedistSigma);
local_sigma = pme->lb_buf2;
{
local_sigma[i] = atc->coefficient[i];
}
- where();
wallcycle_stop(wcycle, ewcPME_REDISTXF);
}
pfft_setup = pme->pfft_setup[grid_index];
calc_next_lb_coeffs(pme, local_sigma);
grid = pmegrid->grid.grid;
- where();
if (flags & GMX_PME_SPREAD)
{
if (pme->nnodes > 1)
{
gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_FORWARD);
- where();
}
#endif
copy_pmegrid_to_fftgrid(pme, grid, fftgrid, grid_index);
{
wallcycle_stop(wcycle, ewcPME_FFT);
}
- where();
}
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
if (thread == 0)
{
wallcycle_stop(wcycle, ewcLJPME);
- where();
inc_nrnb(nrnb, eNR_SOLVEPME, loop_count);
}
}
pfft_setup = pme->pfft_setup[grid_index];
grid = pmegrid->grid.grid;
calc_next_lb_coeffs(pme, local_sigma);
- where();
#pragma omp parallel num_threads(pme->nthread) private(thread)
{
try
/* do 3d-invfft */
if (thread == 0)
{
- where();
wallcycle_start(wcycle, ewcPME_FFT);
}
{
wallcycle_stop(wcycle, ewcPME_FFT);
- where();
if (pme->nodeid == 0)
{
gmx_sum_qgrid_dd(pme, grid, GMX_SUM_GRID_BACKWARD);
}
#endif
- where();
unwrap_periodic_pmegrid(pme, grid);
if (bCalcF)
{
/* interpolate forces for our local atoms */
- where();
bClearF = (bFirst && PAR(cr));
scale = pme->bFEP ? (fep_state < 1 ? 1.0-lambda_lj : lambda_lj) : 1.0;
scale *= lb_scale_factor[grid_index-2];
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
}
- where();
inc_nrnb(nrnb, eNR_GATHERFBSP,
pme->pme_order*pme->pme_order*pme->pme_order*pme->atc[0].n);
wallcycle_stop(wcycle, ewcPME_REDISTXF);
}
- where();
if (bCalcEnerVir)
{
sfree(pme->bsp_mod[i]);
}
- destroy_overlap_comm(&pme->overlap[0]);
- destroy_overlap_comm(&pme->overlap[1]);
-
sfree(pme->lb_buf1);
sfree(pme->lb_buf2);
sfree(pme->sum_qgrid_tmp);
sfree(pme->sum_qgrid_dd_tmp);
+ destroy_pme_spline_work(pme->spline_work);
+
if (pme_gpu_active(pme) && pme->gpu)
{
pme_gpu_destroy(pme->gpu);
#include <algorithm>
#include <vector>
+#include "gromacs/compat/make_unique.h"
#include "gromacs/fileio/filetypes.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/gmxfio-xdr.h"
*/
enum tpxv {
tpxv_ComputationalElectrophysiology = 96, /**< support for ion/water position swaps (computational electrophysiology) */
- tpxv_Use64BitRandomSeed, /**< change ld_seed from int to gmx_int64_t */
+ tpxv_Use64BitRandomSeed, /**< change ld_seed from int to int64_t */
tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, /**< potentials for supporting coarse-grained force fields */
tpxv_InteractiveMolecularDynamics, /**< interactive molecular dynamics (IMD) */
tpxv_RemoveObsoleteParameters1, /**< remove optimize_fft, dihre_fc, nstcheckpoint */
tpxv_PullExternalPotential, /**< Added pull type external potential */
tpxv_GenericParamsForElectricField, /**< Introduced KeyValueTree and moved electric field parameters */
tpxv_AcceleratedWeightHistogram, /**< sampling with accelerated weight histogram method (AWH) */
+ tpxv_RemoveImplicitSolvation, /**< removed support for implicit solvation */
+ tpxv_PullPrevStepCOMAsReference, /**< Enabled using the COM of the pull group of the last frame as reference for PBC */
+ tpxv_MimicQMMM, /**< Inroduced support for MiMiC QM/MM interface */
+ tpxv_PullAverage, /**< Added possibility to output average pull force and position */
tpxv_Count /**< the total number of tpxv versions */
};
/* This number should be the most recent backwards incompatible version
* I.e., if this number is 9, we cannot read tpx version 9 with this code.
*/
-static const int tpx_incompatible_version = 30; // GMX3.2 has version 31
+static const int tpx_incompatible_version = 57; // GMX4.0 has version 58
} t_ftupd;
/*
+ * TODO The following three lines make little sense, please clarify if
+ * you've had to work out how ftupd works.
+ *
* The entries should be ordered in:
* 1. ascending function type number
* 2. ascending file version number
+ *
+ * Because we support reading of old .tpr file versions (even when
+ * mdrun can no longer run the simulation), we need to be able to read
+ * obsolete t_interaction_function types. Any data read from such
+ * fields is discarded. Their names have _NOLONGERUSED appended to
+ * them to make things clear.
*/
static const t_ftupd ftupd[] = {
- { 34, F_FENEBONDS },
- { 43, F_TABBONDS },
- { 43, F_TABBONDSNC },
{ 70, F_RESTRBONDS },
{ tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_RESTRANGLES },
{ 76, F_LINEAR_ANGLES },
- { 34, F_QUARTIC_ANGLES },
- { 43, F_TABANGLES },
{ tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_RESTRDIHS },
{ tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_CBTDIHS },
- { 43, F_TABDIHS },
{ 65, F_CMAP },
- { 60, F_GB12 },
- { 61, F_GB13 },
- { 61, F_GB14 },
- { 72, F_GBPOL },
- { 72, F_NPSOLVATION },
- { 41, F_LJC14_Q },
- { 41, F_LJC_PAIRS_NB },
- { 32, F_BHAM_LR_NOLONGERUSED },
- { 32, F_RF_EXCL },
- { 32, F_COUL_RECIP },
+ { 60, F_GB12_NOLONGERUSED },
+ { 61, F_GB13_NOLONGERUSED },
+ { 61, F_GB14_NOLONGERUSED },
+ { 72, F_GBPOL_NOLONGERUSED },
+ { 72, F_NPSOLVATION_NOLONGERUSED },
{ 93, F_LJ_RECIP },
- { 46, F_DPD },
- { 36, F_THOLE_POL },
{ 90, F_FBPOSRES },
- { 49, F_VSITE4FDN },
- { 50, F_VSITEN },
- { 46, F_COM_PULL },
- { 46, F_ECONSERVED },
{ 69, F_VTEMP_NOLONGERUSED},
{ 66, F_PDISPCORR },
- { 54, F_DVDL_CONSTR },
{ 76, F_ANHARM_POL },
{ 79, F_DVDL_COUL },
{ 79, F_DVDL_VDW, },
static void do_pullgrp_tpx_pre95(t_fileio *fio,
t_pull_group *pgrp,
t_pull_coord *pcrd,
- gmx_bool bRead,
- int file_version)
+ gmx_bool bRead)
{
rvec tmp;
pcrd->init = tmp[0];
gmx_fio_do_real(fio, pcrd->rate);
gmx_fio_do_real(fio, pcrd->k);
- if (file_version >= 56)
- {
- gmx_fio_do_real(fio, pcrd->kB);
- }
- else
- {
- pcrd->kB = pcrd->k;
- }
+ gmx_fio_do_real(fio, pcrd->kB);
}
static void do_pull_group(t_fileio *fio, t_pull_group *pgrp, gmx_bool bRead)
gmx_fio_do_int(fio, expand->lmc_forced_nstart);
gmx_fio_do_int(fio, expand->lmc_seed);
gmx_fio_do_real(fio, expand->mc_temp);
- gmx_fio_do_int(fio, expand->bSymmetrizedTMatrix);
+ gmx_fio_do_gmx_bool(fio, expand->bSymmetrizedTMatrix);
gmx_fio_do_int(fio, expand->nstTij);
gmx_fio_do_int(fio, expand->minvarmin);
gmx_fio_do_int(fio, expand->c_range);
snew(fepvals->all_lambda[g], fepvals->n_lambda);
}
gmx_fio_ndo_double(fio, fepvals->all_lambda[g], fepvals->n_lambda);
- gmx_fio_ndo_int(fio, fepvals->separate_dvdl, efptNR);
+ gmx_fio_ndo_gmx_bool(fio, fepvals->separate_dvdl, efptNR);
}
else if (fepvals->init_lambda >= 0)
{
}
}
gmx_fio_do_real(fio, fepvals->sc_alpha);
- if (file_version >= 38)
- {
- gmx_fio_do_int(fio, fepvals->sc_power);
- }
- else
- {
- fepvals->sc_power = 2;
- }
+ gmx_fio_do_int(fio, fepvals->sc_power);
if (file_version >= 79)
{
gmx_fio_do_real(fio, fepvals->sc_r_power);
}
if (file_version >= 79)
{
- gmx_fio_do_int(fio, fepvals->bScCoul);
+ gmx_fio_do_gmx_bool(fio, fepvals->bScCoul);
}
else
{
gmx_fio_do_real(fio, pull->constr_tol);
if (file_version >= 95)
{
- gmx_fio_do_int(fio, pull->bPrintCOM);
+ gmx_fio_do_gmx_bool(fio, pull->bPrintCOM);
/* With file_version < 95 this value is set below */
}
if (file_version >= tpxv_ReplacePullPrintCOM12)
{
- gmx_fio_do_int(fio, pull->bPrintRefValue);
- gmx_fio_do_int(fio, pull->bPrintComp);
+ gmx_fio_do_gmx_bool(fio, pull->bPrintRefValue);
+ gmx_fio_do_gmx_bool(fio, pull->bPrintComp);
}
else if (file_version >= tpxv_PullCoordTypeGeom)
{
int idum;
gmx_fio_do_int(fio, idum); /* used to be bPrintCOM2 */
- gmx_fio_do_int(fio, pull->bPrintRefValue);
- gmx_fio_do_int(fio, pull->bPrintComp);
+ gmx_fio_do_gmx_bool(fio, pull->bPrintRefValue);
+ gmx_fio_do_gmx_bool(fio, pull->bPrintComp);
}
else
{
}
gmx_fio_do_int(fio, pull->nstxout);
gmx_fio_do_int(fio, pull->nstfout);
+ if (file_version >= tpxv_PullPrevStepCOMAsReference)
+ {
+ gmx_fio_do_gmx_bool(fio, pull->bSetPbcRefToPrevStepCOM);
+ }
+ else
+ {
+ pull->bSetPbcRefToPrevStepCOM = FALSE;
+ }
if (bRead)
{
snew(pull->group, pull->ngroup);
{
/* We read and ignore a pull coordinate for group 0 */
do_pullgrp_tpx_pre95(fio, &pull->group[g], &pull->coord[std::max(g-1, 0)],
- bRead, file_version);
+ bRead);
if (g > 0)
{
pull->coord[g-1].group[0] = 0;
bRead, file_version, ePullOld, eGeomOld, dimOld);
}
}
+ if (file_version >= tpxv_PullAverage)
+ {
+ gmx_bool v;
+
+ v = pull->bXOutAverage;
+ gmx_fio_do_gmx_bool(fio, v);
+ pull->bXOutAverage = v;
+ v = pull->bFOutAverage;
+ gmx_fio_do_gmx_bool(fio, v);
+ pull->bFOutAverage = v;
+ }
}
snew(rotg->x_ref, rotg->nat);
}
gmx_fio_ndo_rvec(fio, rotg->x_ref, rotg->nat);
- gmx_fio_do_rvec(fio, rotg->vec);
+ gmx_fio_do_rvec(fio, rotg->inputVec);
gmx_fio_do_rvec(fio, rotg->pivot);
gmx_fio_do_real(fio, rotg->rate);
gmx_fio_do_real(fio, rotg->k);
{
do_swapgroup(fio, &swap->grp[ig], bRead);
}
- gmx_fio_do_int(fio, swap->massw_split[eChannel0]);
- gmx_fio_do_int(fio, swap->massw_split[eChannel1]);
+ gmx_fio_do_gmx_bool(fio, swap->massw_split[eChannel0]);
+ gmx_fio_do_gmx_bool(fio, swap->massw_split[eChannel1]);
gmx_fio_do_int(fio, swap->nstswap);
gmx_fio_do_int(fio, swap->nAverage);
gmx_fio_do_real(fio, swap->threshold);
gmx_fio_do_int(fio, swap->grp[3].nat);
gmx_fio_do_int(fio, swap->grp[eGrpSolvent].nat);
gmx_fio_do_int(fio, swap->grp[eGrpSplit0].nat);
- gmx_fio_do_int(fio, swap->massw_split[eChannel0]);
+ gmx_fio_do_gmx_bool(fio, swap->massw_split[eChannel0]);
gmx_fio_do_int(fio, swap->grp[eGrpSplit1].nat);
- gmx_fio_do_int(fio, swap->massw_split[eChannel1]);
+ gmx_fio_do_gmx_bool(fio, swap->massw_split[eChannel1]);
gmx_fio_do_int(fio, swap->nstswap);
gmx_fio_do_int(fio, swap->nAverage);
gmx_fio_do_real(fio, swap->threshold);
static void do_inputrec(t_fileio *fio, t_inputrec *ir, gmx_bool bRead,
- int file_version, real *fudgeQQ)
+ int file_version)
{
int i, j, k, idum = 0;
- real rdum, bd_temp;
- gmx_bool bdum = 0;
+ real rdum;
+ gmx_bool bdum = false;
if (file_version != tpx_version)
{
ir->init_step = idum;
}
- if (file_version >= 58)
- {
- gmx_fio_do_int(fio, ir->simulation_part);
- }
- else
- {
- ir->simulation_part = 1;
- }
+ gmx_fio_do_int(fio, ir->simulation_part);
if (file_version >= 67)
{
{
ir->nstcalcenergy = 1;
}
- if (file_version < 53)
- {
- /* The pbc info has been moved out of do_inputrec,
- * since we always want it, also without reading the inputrec.
- */
- gmx_fio_do_int(fio, ir->ePBC);
- if (file_version >= 45)
- {
- gmx_fio_do_int(fio, ir->bPeriodicMols);
- }
- else
- {
- if (ir->ePBC == 2)
- {
- ir->ePBC = epbcXYZ;
- ir->bPeriodicMols = TRUE;
- }
- else
- {
- ir->bPeriodicMols = FALSE;
- }
- }
- }
if (file_version >= 81)
{
gmx_fio_do_int(fio, ir->cutoff_scheme);
gmx_fio_do_int(fio, ir->ns_type);
gmx_fio_do_int(fio, ir->nstlist);
gmx_fio_do_int(fio, idum); /* used to be ndelta; not used anymore */
- if (file_version < 41)
- {
- gmx_fio_do_int(fio, idum);
- gmx_fio_do_int(fio, idum);
- }
- if (file_version >= 45)
- {
- gmx_fio_do_real(fio, ir->rtpi);
- }
- else
- {
- ir->rtpi = 0.05;
- }
+
+ gmx_fio_do_real(fio, ir->rtpi);
+
gmx_fio_do_int(fio, ir->nstcomm);
- if (file_version > 34)
- {
- gmx_fio_do_int(fio, ir->comm_mode);
- }
- else if (ir->nstcomm < 0)
- {
- ir->comm_mode = ecmANGULAR;
- }
- else
- {
- ir->comm_mode = ecmLINEAR;
- }
- ir->nstcomm = abs(ir->nstcomm);
+ gmx_fio_do_int(fio, ir->comm_mode);
/* ignore nstcheckpoint */
if (file_version < tpxv_RemoveObsoleteParameters1)
gmx_fio_do_int(fio, dummy_nstcalclr);
}
gmx_fio_do_int(fio, ir->coulombtype);
- if (file_version < 32 && ir->coulombtype == eelRF)
- {
- ir->coulombtype = eelRF_NEC_UNSUPPORTED;
- }
if (file_version >= 81)
{
gmx_fio_do_int(fio, ir->coulomb_modifier);
gmx_fio_do_real(fio, ir->rvdw);
gmx_fio_do_int(fio, ir->eDispCorr);
gmx_fio_do_real(fio, ir->epsilon_r);
- if (file_version >= 37)
- {
- gmx_fio_do_real(fio, ir->epsilon_rf);
- }
- else
- {
- if (EEL_RF(ir->coulombtype))
- {
- ir->epsilon_rf = ir->epsilon_r;
- ir->epsilon_r = 1.0;
- }
- else
- {
- ir->epsilon_rf = 1.0;
- }
- }
+ gmx_fio_do_real(fio, ir->epsilon_rf);
gmx_fio_do_real(fio, ir->tabext);
- gmx_fio_do_int(fio, ir->gb_algorithm);
- gmx_fio_do_int(fio, ir->nstgbradii);
- gmx_fio_do_real(fio, ir->rgbradii);
- gmx_fio_do_real(fio, ir->gb_saltconc);
- gmx_fio_do_int(fio, ir->implicit_solvent);
- if (file_version >= 55)
+ // This permits reading a .tpr file that used implicit solvent,
+ // and later permitting mdrun to refuse to run it.
+ if (bRead)
{
- gmx_fio_do_real(fio, ir->gb_epsilon_solvent);
- gmx_fio_do_real(fio, ir->gb_obc_alpha);
- gmx_fio_do_real(fio, ir->gb_obc_beta);
- gmx_fio_do_real(fio, ir->gb_obc_gamma);
- if (file_version >= 60)
+ if (file_version < tpxv_RemoveImplicitSolvation)
{
- gmx_fio_do_real(fio, ir->gb_dielectric_offset);
- gmx_fio_do_int(fio, ir->sa_algorithm);
+ gmx_fio_do_int(fio, idum);
+ gmx_fio_do_int(fio, idum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_int(fio, idum);
+ ir->implicit_solvent = (idum > 0);
}
else
{
- ir->gb_dielectric_offset = 0.009;
- ir->sa_algorithm = esaAPPROX;
+ ir->implicit_solvent = false;
}
- gmx_fio_do_real(fio, ir->sa_surface_tension);
-
- /* Override sa_surface_tension if it is not changed in the mpd-file */
- if (ir->sa_surface_tension < 0)
+ if (file_version < tpxv_RemoveImplicitSolvation)
{
- if (ir->gb_algorithm == egbSTILL)
- {
- ir->sa_surface_tension = 0.0049 * 100 * CAL2JOULE;
- }
- else if (ir->gb_algorithm == egbHCT || ir->gb_algorithm == egbOBC)
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ if (file_version >= 60)
{
- ir->sa_surface_tension = 0.0054 * 100 * CAL2JOULE;
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_int(fio, idum);
}
+ gmx_fio_do_real(fio, rdum);
}
-
- }
- else
- {
- /* Better use sensible values than insane (0.0) ones... */
- ir->gb_epsilon_solvent = 80;
- ir->gb_obc_alpha = 1.0;
- ir->gb_obc_beta = 0.8;
- ir->gb_obc_gamma = 4.85;
- ir->sa_surface_tension = 2.092;
}
-
if (file_version >= 81)
{
gmx_fio_do_real(fio, ir->fourier_spacing);
gmx_fio_do_rvec(fio, ir->compress[XX]);
gmx_fio_do_rvec(fio, ir->compress[YY]);
gmx_fio_do_rvec(fio, ir->compress[ZZ]);
- if (file_version >= 47)
- {
- gmx_fio_do_int(fio, ir->refcoord_scaling);
- gmx_fio_do_rvec(fio, ir->posres_com);
- gmx_fio_do_rvec(fio, ir->posres_comB);
- }
- else
- {
- ir->refcoord_scaling = erscNO;
- clear_rvec(ir->posres_com);
- clear_rvec(ir->posres_comB);
- }
+ gmx_fio_do_int(fio, ir->refcoord_scaling);
+ gmx_fio_do_rvec(fio, ir->posres_com);
+ gmx_fio_do_rvec(fio, ir->posres_comB);
+
if (file_version < 79)
{
gmx_fio_do_int(fio, ir->andersen_seed);
ir->andersen_seed = 0;
}
- if (file_version < 37)
- {
- gmx_fio_do_real(fio, rdum);
- }
-
gmx_fio_do_real(fio, ir->shake_tol);
- if (file_version < 54)
- {
- // cppcheck-suppress redundantPointerOp
- gmx_fio_do_real(fio, *fudgeQQ);
- }
gmx_fio_do_int(fio, ir->efep);
do_fepvals(fio, ir->fepvals, bRead, file_version);
{
do_expandedvals(fio, ir->expandedvals, ir->fepvals, bRead, file_version);
}
- if (file_version >= 57)
- {
- gmx_fio_do_int(fio, ir->eDisre);
- }
+
+ gmx_fio_do_int(fio, ir->eDisre);
gmx_fio_do_int(fio, ir->eDisreWeighting);
gmx_fio_do_gmx_bool(fio, ir->bDisreMixed);
gmx_fio_do_real(fio, ir->dr_fc);
if (file_version < 79)
{
gmx_fio_do_real(fio, rdum);
- if (file_version < 56)
- {
- gmx_fio_do_real(fio, rdum);
- gmx_fio_do_int(fio, idum);
- }
}
gmx_fio_do_real(fio, ir->em_stepsize);
gmx_fio_do_int(fio, ir->nProjOrder);
gmx_fio_do_real(fio, ir->LincsWarnAngle);
gmx_fio_do_int(fio, ir->nLincsIter);
- if (file_version < 33)
- {
- gmx_fio_do_real(fio, bd_temp);
- }
gmx_fio_do_real(fio, ir->bd_fric);
if (file_version >= tpxv_Use64BitRandomSeed)
{
gmx_fio_do_int(fio, idum);
ir->ld_seed = idum;
}
- if (file_version >= 33)
- {
- for (i = 0; i < DIM; i++)
- {
- gmx_fio_do_rvec(fio, ir->deform[i]);
- }
- }
- else
+
+ for (i = 0; i < DIM; i++)
{
- for (i = 0; i < DIM; i++)
- {
- clear_rvec(ir->deform[i]);
- }
+ gmx_fio_do_rvec(fio, ir->deform[i]);
}
gmx_fio_do_real(fio, ir->cos_accel);
}
/* pull stuff */
- if (file_version >= 48)
{
int ePullOld = 0;
do_pull(fio, ir->pull, bRead, file_version, ePullOld);
}
}
- else
- {
- ir->bPull = FALSE;
- }
if (file_version >= tpxv_AcceleratedWeightHistogram)
{
/* Enforced rotation */
if (file_version >= 74)
{
- gmx_fio_do_int(fio, ir->bRot);
- if (ir->bRot == TRUE)
+ gmx_fio_do_gmx_bool(fio, ir->bRot);
+ if (ir->bRot)
{
if (bRead)
{
/* Interactive molecular dynamics */
if (file_version >= tpxv_InteractiveMolecularDynamics)
{
- gmx_fio_do_int(fio, ir->bIMD);
- if (TRUE == ir->bIMD)
+ gmx_fio_do_gmx_bool(fio, ir->bIMD);
+ if (ir->bIMD)
{
if (bRead)
{
gmx_fio_ndo_real(fio, ir->opts.nrdf, ir->opts.ngtc);
gmx_fio_ndo_real(fio, ir->opts.ref_t, ir->opts.ngtc);
gmx_fio_ndo_real(fio, ir->opts.tau_t, ir->opts.ngtc);
- if (file_version < 33 && ir->eI == eiBD)
- {
- for (i = 0; i < ir->opts.ngtc; i++)
- {
- ir->opts.tau_t[i] = bd_temp;
- }
- }
}
if (ir->opts.ngfrz > 0)
{
gmx_fio_ndo_real(fio, ir->opts.anneal_temp[j], k);
}
/* Walls */
- if (file_version >= 45)
{
gmx_fio_do_int(fio, ir->nwall);
gmx_fio_do_int(fio, ir->wall_type);
- if (file_version >= 50)
- {
- gmx_fio_do_real(fio, ir->wall_r_linpot);
- }
- else
- {
- ir->wall_r_linpot = -1;
- }
+ gmx_fio_do_real(fio, ir->wall_r_linpot);
gmx_fio_do_int(fio, ir->wall_atomtype[0]);
gmx_fio_do_int(fio, ir->wall_atomtype[1]);
gmx_fio_do_real(fio, ir->wall_density[0]);
gmx_fio_do_real(fio, ir->wall_density[1]);
gmx_fio_do_real(fio, ir->wall_ewald_zfac);
}
- else
- {
- ir->nwall = 0;
- ir->wall_type = 0;
- ir->wall_atomtype[0] = -1;
- ir->wall_atomtype[1] = -1;
- ir->wall_density[0] = 0;
- ir->wall_density[1] = 0;
- ir->wall_ewald_zfac = 3;
- }
+
/* Cosine stuff for electric fields */
if (file_version < tpxv_GenericParamsForElectricField)
{
}
/* QMMM stuff */
- if (file_version >= 39)
{
gmx_fio_do_gmx_bool(fio, ir->bQMMM);
gmx_fio_do_int(fio, ir->QMMMscheme);
snew(ir->opts.SAoff, ir->opts.ngQM);
snew(ir->opts.SAsteps, ir->opts.ngQM);
}
- if (ir->opts.ngQM > 0)
+ if (ir->opts.ngQM > 0 && ir->bQMMM)
{
gmx_fio_ndo_int(fio, ir->opts.QMmethod, ir->opts.ngQM);
gmx_fio_ndo_int(fio, ir->opts.QMbasis, ir->opts.ngQM);
* changing the tpr format for every QMMM change.
*/
std::vector<int> dummy(ir->opts.ngQM, 0);
- gmx_fio_ndo_gmx_bool(fio, dummy.data(), ir->opts.ngQM);
- gmx_fio_ndo_gmx_bool(fio, dummy.data(), ir->opts.ngQM);
+ gmx_fio_ndo_int(fio, dummy.data(), ir->opts.ngQM);
+ gmx_fio_ndo_int(fio, dummy.data(), ir->opts.ngQM);
}
/* end of QMMM stuff */
}
case F_ANGRESZ:
gmx_fio_do_real(fio, iparams->pdihs.phiA);
gmx_fio_do_real(fio, iparams->pdihs.cpA);
- if ((ftype == F_ANGRES || ftype == F_ANGRESZ) && file_version < 42)
- {
- /* Read the incorrectly stored multiplicity */
- gmx_fio_do_real(fio, iparams->harmonic.rB);
- gmx_fio_do_real(fio, iparams->harmonic.krB);
- iparams->pdihs.phiB = iparams->pdihs.phiA;
- iparams->pdihs.cpB = iparams->pdihs.cpA;
- }
- else
- {
- gmx_fio_do_real(fio, iparams->pdihs.phiB);
- gmx_fio_do_real(fio, iparams->pdihs.cpB);
- gmx_fio_do_int(fio, iparams->pdihs.mult);
- }
+ gmx_fio_do_real(fio, iparams->pdihs.phiB);
+ gmx_fio_do_real(fio, iparams->pdihs.cpB);
+ gmx_fio_do_int(fio, iparams->pdihs.mult);
break;
case F_RESTRDIHS:
gmx_fio_do_real(fio, iparams->pdihs.phiA);
gmx_fio_do_int(fio, iparams->vsiten.n);
gmx_fio_do_real(fio, iparams->vsiten.a);
break;
- case F_GB12:
- case F_GB13:
- case F_GB14:
- /* We got rid of some parameters in version 68 */
- if (bRead && file_version < 68)
+ case F_GB12_NOLONGERUSED:
+ case F_GB13_NOLONGERUSED:
+ case F_GB14_NOLONGERUSED:
+ // Implicit solvent parameters can still be read, but never used
+ if (bRead)
{
- gmx_fio_do_real(fio, rdum);
- gmx_fio_do_real(fio, rdum);
- gmx_fio_do_real(fio, rdum);
- gmx_fio_do_real(fio, rdum);
+ if (file_version < 68)
+ {
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ }
+ if (file_version < tpxv_RemoveImplicitSolvation)
+ {
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ gmx_fio_do_real(fio, rdum);
+ }
}
- gmx_fio_do_real(fio, iparams->gb.sar);
- gmx_fio_do_real(fio, iparams->gb.st);
- gmx_fio_do_real(fio, iparams->gb.pi);
- gmx_fio_do_real(fio, iparams->gb.gbr);
- gmx_fio_do_real(fio, iparams->gb.bmlt);
break;
case F_CMAP:
gmx_fio_do_int(fio, iparams->cmap.cmapA);
}
}
-static void do_ilist(t_fileio *fio, t_ilist *ilist, gmx_bool bRead, int file_version)
+static void do_ilist(t_fileio *fio, InteractionList *ilist, gmx_bool bRead)
{
- int i, idum;
-
- if (file_version < 44)
- {
- for (i = 0; i < MAXNODES; i++)
- {
- gmx_fio_do_int(fio, idum);
- }
- }
- gmx_fio_do_int(fio, ilist->nr);
+ int nr = ilist->size();
+ gmx_fio_do_int(fio, nr);
if (bRead)
{
- snew(ilist->iatoms, ilist->nr);
+ ilist->iatoms.resize(nr);
}
- gmx_fio_ndo_int(fio, ilist->iatoms, ilist->nr);
+ gmx_fio_ndo_int(fio, ilist->iatoms.data(), ilist->size());
}
static void do_ffparams(t_fileio *fio, gmx_ffparams_t *ffparams,
gmx_bool bRead, int file_version)
{
- int idum, i;
- unsigned int k;
-
gmx_fio_do_int(fio, ffparams->atnr);
- if (file_version < 57)
- {
- gmx_fio_do_int(fio, idum);
- }
- gmx_fio_do_int(fio, ffparams->ntypes);
+ int numTypes = ffparams->numTypes();
+ gmx_fio_do_int(fio, numTypes);
if (bRead)
{
- snew(ffparams->functype, ffparams->ntypes);
- snew(ffparams->iparams, ffparams->ntypes);
+ ffparams->functype.resize(numTypes);
+ ffparams->iparams.resize(numTypes);
}
/* Read/write all the function types */
- gmx_fio_ndo_int(fio, ffparams->functype, ffparams->ntypes);
+ gmx_fio_ndo_int(fio, ffparams->functype.data(), ffparams->functype.size());
if (file_version >= 66)
{
ffparams->reppow = 12.0;
}
- if (file_version >= 57)
- {
- gmx_fio_do_real(fio, ffparams->fudgeQQ);
- }
+ gmx_fio_do_real(fio, ffparams->fudgeQQ);
/* Check whether all these function types are supported by the code.
* In practice the code is backwards compatible, which means that the
* numbering may have to be altered from old numbering to new numbering
*/
- for (i = 0; (i < ffparams->ntypes); i++)
+ for (int i = 0; i < ffparams->numTypes(); i++)
{
if (bRead)
{
/* Loop over file versions */
- for (k = 0; (k < NFTUPD); k++)
+ for (int k = 0; k < NFTUPD; k++)
{
/* Compare the read file_version to the update table */
if ((file_version < ftupd[k].fvnr) &&
}
}
-static void add_settle_atoms(t_ilist *ilist)
+static void add_settle_atoms(InteractionList *ilist)
{
int i;
/* Settle used to only store the first atom: add the other two */
- srenew(ilist->iatoms, 2*ilist->nr);
- for (i = ilist->nr/2-1; i >= 0; i--)
+ ilist->iatoms.resize(2*ilist->size());
+ for (i = ilist->size()/4 - 1; i >= 0; i--)
{
ilist->iatoms[4*i+0] = ilist->iatoms[2*i+0];
ilist->iatoms[4*i+1] = ilist->iatoms[2*i+1];
ilist->iatoms[4*i+2] = ilist->iatoms[2*i+1] + 1;
ilist->iatoms[4*i+3] = ilist->iatoms[2*i+1] + 2;
}
- ilist->nr = 2*ilist->nr;
}
-static void do_ilists(t_fileio *fio, t_ilist *ilist, gmx_bool bRead,
+static void do_ilists(t_fileio *fio, InteractionLists *ilists, gmx_bool bRead,
int file_version)
{
- int j;
- gmx_bool bClear;
- unsigned int k;
+ GMX_RELEASE_ASSERT(ilists, "Need a valid ilists object");
+ GMX_RELEASE_ASSERT(ilists->size() == F_NRE, "The code needs to be in sync with InteractionLists");
- for (j = 0; (j < F_NRE); j++)
+ for (int j = 0; j < F_NRE; j++)
{
- bClear = FALSE;
+ InteractionList &ilist = (*ilists)[j];
+ gmx_bool bClear = FALSE;
if (bRead)
{
- for (k = 0; k < NFTUPD; k++)
+ for (int k = 0; k < NFTUPD; k++)
{
if ((file_version < ftupd[k].fvnr) && (j == ftupd[k].ftype))
{
}
if (bClear)
{
- ilist[j].nr = 0;
- ilist[j].iatoms = nullptr;
+ ilist.iatoms.clear();
}
else
{
- do_ilist(fio, &ilist[j], bRead, file_version);
- if (file_version < 78 && j == F_SETTLE && ilist[j].nr > 0)
+ do_ilist(fio, &ilist, bRead);
+ if (file_version < 78 && j == F_SETTLE && ilist.size() > 0)
{
- add_settle_atoms(&ilist[j]);
+ add_settle_atoms(&ilist);
}
}
}
}
-static void do_idef(t_fileio *fio, gmx_ffparams_t *ffparams, gmx_moltype_t *molt,
- gmx_bool bRead, int file_version)
+static void do_block(t_fileio *fio, t_block *block, gmx_bool bRead)
{
- do_ffparams(fio, ffparams, bRead, file_version);
-
- if (file_version >= 54)
- {
- gmx_fio_do_real(fio, ffparams->fudgeQQ);
- }
-
- do_ilists(fio, molt->ilist, bRead, file_version);
-}
-
-static void do_block(t_fileio *fio, t_block *block, gmx_bool bRead, int file_version)
-{
- int i, idum, dum_nra, *dum_a;
-
- if (file_version < 44)
- {
- for (i = 0; i < MAXNODES; i++)
- {
- gmx_fio_do_int(fio, idum);
- }
- }
gmx_fio_do_int(fio, block->nr);
- if (file_version < 51)
- {
- gmx_fio_do_int(fio, dum_nra);
- }
if (bRead)
{
if ((block->nalloc_index > 0) && (nullptr != block->index))
snew(block->index, block->nalloc_index);
}
gmx_fio_ndo_int(fio, block->index, block->nr+1);
-
- if (file_version < 51 && dum_nra > 0)
- {
- snew(dum_a, dum_nra);
- gmx_fio_ndo_int(fio, dum_a, dum_nra);
- sfree(dum_a);
- }
}
-static void do_blocka(t_fileio *fio, t_blocka *block, gmx_bool bRead,
- int file_version)
+static void do_blocka(t_fileio *fio, t_blocka *block, gmx_bool bRead)
{
- int i, idum;
-
- if (file_version < 44)
- {
- for (i = 0; i < MAXNODES; i++)
- {
- gmx_fio_do_int(fio, idum);
- }
- }
gmx_fio_do_int(fio, block->nr);
gmx_fio_do_int(fio, block->nra);
if (bRead)
}
-static void do_atom(t_fileio *fio, t_atom *atom, int ngrp, gmx_bool bRead,
- int file_version, gmx_groups_t *groups, int atnr)
+static void do_atom(t_fileio *fio, t_atom *atom, gmx_bool bRead)
{
- int i, myngrp;
-
gmx_fio_do_real(fio, atom->m);
gmx_fio_do_real(fio, atom->q);
gmx_fio_do_real(fio, atom->mB);
gmx_fio_do_ushort(fio, atom->typeB);
gmx_fio_do_int(fio, atom->ptype);
gmx_fio_do_int(fio, atom->resind);
- if (file_version >= 52)
- {
- gmx_fio_do_int(fio, atom->atomnumber);
- if (bRead)
- {
- /* Set element string from atomic number if present.
- * This routine returns an empty string if the name is not found.
- */
- std::strncpy(atom->elem, atomicnumber_to_element(atom->atomnumber), 4);
- /* avoid warnings about potentially unterminated string */
- atom->elem[3] = '\0';
- }
- }
- else if (bRead)
- {
- atom->atomnumber = 0;
- }
- if (file_version < 39)
- {
- myngrp = 9;
- }
- else
- {
- myngrp = ngrp;
- }
-
- if (file_version < 57)
+ gmx_fio_do_int(fio, atom->atomnumber);
+ if (bRead)
{
- unsigned char uchar[egcNR];
- gmx_fio_ndo_uchar(fio, uchar, myngrp);
- for (i = myngrp; (i < ngrp); i++)
- {
- uchar[i] = 0;
- }
- /* Copy the old data format to the groups struct */
- for (i = 0; i < ngrp; i++)
- {
- groups->grpnr[i][atnr] = uchar[i];
- }
+ /* Set element string from atomic number if present.
+ * This routine returns an empty string if the name is not found.
+ */
+ std::strncpy(atom->elem, atomicnumber_to_element(atom->atomnumber), 4);
+ /* avoid warnings about potentially unterminated string */
+ atom->elem[3] = '\0';
}
}
-static void do_grps(t_fileio *fio, int ngrp, t_grps grps[], gmx_bool bRead,
- int file_version)
+static void do_grps(t_fileio *fio, int ngrp, t_grps grps[], gmx_bool bRead)
{
- int j, myngrp;
-
- if (file_version < 39)
- {
- myngrp = 9;
- }
- else
- {
- myngrp = ngrp;
- }
-
- for (j = 0; (j < ngrp); j++)
+ for (int j = 0; j < ngrp; j++)
{
- if (j < myngrp)
- {
- gmx_fio_do_int(fio, grps[j].nr);
- if (bRead)
- {
- snew(grps[j].nm_ind, grps[j].nr);
- }
- gmx_fio_ndo_int(fio, grps[j].nm_ind, grps[j].nr);
- }
- else
+ gmx_fio_do_int(fio, grps[j].nr);
+ if (bRead)
{
- grps[j].nr = 1;
snew(grps[j].nm_ind, grps[j].nr);
}
+ gmx_fio_ndo_int(fio, grps[j].nm_ind, grps[j].nr);
}
}
}
static void do_atoms(t_fileio *fio, t_atoms *atoms, gmx_bool bRead, t_symtab *symtab,
- int file_version,
- gmx_groups_t *groups)
+ int file_version)
{
int i;
gmx_fio_do_int(fio, atoms->nr);
gmx_fio_do_int(fio, atoms->nres);
- if (file_version < 57)
- {
- gmx_fio_do_int(fio, groups->ngrpname);
- for (i = 0; i < egcNR; i++)
- {
- groups->ngrpnr[i] = atoms->nr;
- snew(groups->grpnr[i], groups->ngrpnr[i]);
- }
- }
if (bRead)
{
/* Since we have always written all t_atom properties in the tpr file
snew(atoms->atomtype, atoms->nr);
snew(atoms->atomtypeB, atoms->nr);
snew(atoms->resinfo, atoms->nres);
- if (file_version < 57)
- {
- snew(groups->grpname, groups->ngrpname);
- }
atoms->pdbinfo = nullptr;
}
else
}
for (i = 0; (i < atoms->nr); i++)
{
- do_atom(fio, &atoms->atom[i], egcNR, bRead, file_version, groups, i);
+ do_atom(fio, &atoms->atom[i], bRead);
}
do_strstr(fio, atoms->nr, atoms->atomname, bRead, symtab);
do_strstr(fio, atoms->nr, atoms->atomtype, bRead, symtab);
do_strstr(fio, atoms->nr, atoms->atomtypeB, bRead, symtab);
do_resinfo(fio, atoms->nres, atoms->resinfo, bRead, symtab, file_version);
-
- if (file_version < 57)
- {
- do_strstr(fio, groups->ngrpname, groups->grpname, bRead, symtab);
-
- do_grps(fio, egcNR, groups->grps, bRead, file_version);
- }
}
static void do_groups(t_fileio *fio, gmx_groups_t *groups,
- gmx_bool bRead, t_symtab *symtab,
- int file_version)
+ gmx_bool bRead, t_symtab *symtab)
{
int g;
- do_grps(fio, egcNR, groups->grps, bRead, file_version);
+ do_grps(fio, egcNR, groups->grps, bRead);
gmx_fio_do_int(fio, groups->ngrpname);
if (bRead)
{
j = atomtypes->nr;
if (bRead)
{
- snew(atomtypes->radius, j);
- snew(atomtypes->vol, j);
- snew(atomtypes->surftens, j);
snew(atomtypes->atomnumber, j);
- snew(atomtypes->gb_radius, j);
- snew(atomtypes->S_hct, j);
}
- gmx_fio_ndo_real(fio, atomtypes->radius, j);
- gmx_fio_ndo_real(fio, atomtypes->vol, j);
- gmx_fio_ndo_real(fio, atomtypes->surftens, j);
- if (file_version >= 40)
+ if (bRead && file_version < tpxv_RemoveImplicitSolvation)
{
- gmx_fio_ndo_int(fio, atomtypes->atomnumber, j);
+ std::vector<real> dummy(atomtypes->nr, 0);
+ gmx_fio_ndo_real(fio, dummy.data(), dummy.size());
+ gmx_fio_ndo_real(fio, dummy.data(), dummy.size());
+ gmx_fio_ndo_real(fio, dummy.data(), dummy.size());
}
- if (file_version >= 60)
+ gmx_fio_ndo_int(fio, atomtypes->atomnumber, j);
+
+ if (bRead && file_version >= 60 && file_version < tpxv_RemoveImplicitSolvation)
{
- gmx_fio_ndo_real(fio, atomtypes->gb_radius, j);
- gmx_fio_ndo_real(fio, atomtypes->S_hct, j);
+ std::vector<real> dummy(atomtypes->nr, 0);
+ gmx_fio_ndo_real(fio, dummy.data(), dummy.size());
+ gmx_fio_ndo_real(fio, dummy.data(), dummy.size());
}
}
static void do_cmap(t_fileio *fio, gmx_cmap_t *cmap_grid, gmx_bool bRead)
{
- int i, j, ngrid, gs, nelem;
- gmx_fio_do_int(fio, cmap_grid->ngrid);
+ int ngrid = cmap_grid->cmapdata.size();
+ gmx_fio_do_int(fio, ngrid);
gmx_fio_do_int(fio, cmap_grid->grid_spacing);
- ngrid = cmap_grid->ngrid;
- gs = cmap_grid->grid_spacing;
- nelem = gs * gs;
+ int gs = cmap_grid->grid_spacing;
+ int nelem = gs * gs;
if (bRead)
{
- snew(cmap_grid->cmapdata, ngrid);
+ cmap_grid->cmapdata.resize(ngrid);
- for (i = 0; i < cmap_grid->ngrid; i++)
+ for (int i = 0; i < ngrid; i++)
{
- snew(cmap_grid->cmapdata[i].cmap, 4*nelem);
+ cmap_grid->cmapdata[i].cmap.resize(4*nelem);
}
}
- for (i = 0; i < cmap_grid->ngrid; i++)
+ for (int i = 0; i < ngrid; i++)
{
- for (j = 0; j < nelem; j++)
+ for (int j = 0; j < nelem; j++)
{
gmx_fio_do_real(fio, cmap_grid->cmapdata[i].cmap[j*4]);
gmx_fio_do_real(fio, cmap_grid->cmapdata[i].cmap[j*4+1]);
static void do_moltype(t_fileio *fio, gmx_moltype_t *molt, gmx_bool bRead,
- t_symtab *symtab, int file_version,
- gmx_groups_t *groups)
+ t_symtab *symtab, int file_version)
{
- if (file_version >= 57)
- {
- do_symstr(fio, &(molt->name), bRead, symtab);
- }
+ do_symstr(fio, &(molt->name), bRead, symtab);
- do_atoms(fio, &molt->atoms, bRead, symtab, file_version, groups);
+ do_atoms(fio, &molt->atoms, bRead, symtab, file_version);
- if (file_version >= 57)
- {
- do_ilists(fio, molt->ilist, bRead, file_version);
+ do_ilists(fio, &molt->ilist, bRead, file_version);
- do_block(fio, &molt->cgs, bRead, file_version);
- }
+ do_block(fio, &molt->cgs, bRead);
/* This used to be in the atoms struct */
- do_blocka(fio, &molt->excls, bRead, file_version);
+ do_blocka(fio, &molt->excls, bRead);
}
-static void do_molblock(t_fileio *fio, gmx_molblock_t *molb, gmx_bool bRead)
+static void do_molblock(t_fileio *fio, gmx_molblock_t *molb,
+ int numAtomsPerMolecule,
+ gmx_bool bRead)
{
gmx_fio_do_int(fio, molb->type);
gmx_fio_do_int(fio, molb->nmol);
- gmx_fio_do_int(fio, molb->natoms_mol);
+ /* To maintain forward topology reading compatibility, we store #atoms.
+ * TODO: Change this to conditional reading of a dummy int when we
+ * increase tpx_generation.
+ */
+ gmx_fio_do_int(fio, numAtomsPerMolecule);
/* Position restraint coordinates */
- gmx_fio_do_int(fio, molb->nposres_xA);
- if (molb->nposres_xA > 0)
+ int numPosres_xA = molb->posres_xA.size();
+ gmx_fio_do_int(fio, numPosres_xA);
+ if (numPosres_xA > 0)
{
if (bRead)
{
- snew(molb->posres_xA, molb->nposres_xA);
+ molb->posres_xA.resize(numPosres_xA);
}
- gmx_fio_ndo_rvec(fio, molb->posres_xA, molb->nposres_xA);
+ gmx_fio_ndo_rvec(fio, as_rvec_array(molb->posres_xA.data()), numPosres_xA);
}
- gmx_fio_do_int(fio, molb->nposres_xB);
- if (molb->nposres_xB > 0)
+ int numPosres_xB = molb->posres_xB.size();
+ gmx_fio_do_int(fio, numPosres_xB);
+ if (numPosres_xB > 0)
{
if (bRead)
{
- snew(molb->posres_xB, molb->nposres_xB);
+ molb->posres_xB.resize(numPosres_xB);
}
- gmx_fio_ndo_rvec(fio, molb->posres_xB, molb->nposres_xB);
+ gmx_fio_ndo_rvec(fio, as_rvec_array(molb->posres_xB.data()), numPosres_xB);
}
}
-static t_block mtop_mols(gmx_mtop_t *mtop)
-{
- int mb, m, a, mol;
- t_block mols;
-
- mols.nr = 0;
- for (mb = 0; mb < mtop->nmolblock; mb++)
- {
- mols.nr += mtop->molblock[mb].nmol;
- }
- mols.nalloc_index = mols.nr + 1;
- snew(mols.index, mols.nalloc_index);
-
- a = 0;
- m = 0;
- mols.index[m] = a;
- for (mb = 0; mb < mtop->nmolblock; mb++)
- {
- for (mol = 0; mol < mtop->molblock[mb].nmol; mol++)
- {
- a += mtop->molblock[mb].natoms_mol;
- m++;
- mols.index[m] = a;
- }
- }
-
- return mols;
-}
-
-static void add_posres_molblock(gmx_mtop_t *mtop)
-{
- t_ilist *il, *ilfb;
- int am, i, mol, a;
- gmx_bool bFE;
- gmx_molblock_t *molb;
- t_iparams *ip;
-
- /* posres reference positions are stored in ip->posres (if present) and
- in ip->fbposres (if present). If normal and flat-bottomed posres are present,
- posres.pos0A are identical to fbposres.pos0. */
- il = &mtop->moltype[0].ilist[F_POSRES];
- ilfb = &mtop->moltype[0].ilist[F_FBPOSRES];
- if (il->nr == 0 && ilfb->nr == 0)
- {
- return;
- }
- am = 0;
- bFE = FALSE;
- for (i = 0; i < il->nr; i += 2)
- {
- ip = &mtop->ffparams.iparams[il->iatoms[i]];
- am = std::max(am, il->iatoms[i+1]);
- if (ip->posres.pos0B[XX] != ip->posres.pos0A[XX] ||
- ip->posres.pos0B[YY] != ip->posres.pos0A[YY] ||
- ip->posres.pos0B[ZZ] != ip->posres.pos0A[ZZ])
- {
- bFE = TRUE;
- }
- }
- /* This loop is required if we have only flat-bottomed posres:
- - set am
- - bFE == FALSE (no B-state for flat-bottomed posres) */
- if (il->nr == 0)
- {
- for (i = 0; i < ilfb->nr; i += 2)
- {
- am = std::max(am, ilfb->iatoms[i+1]);
- }
- }
- /* Make the posres coordinate block end at a molecule end */
- mol = 0;
- while (am >= mtop->mols.index[mol+1])
- {
- mol++;
- }
- molb = &mtop->molblock[0];
- molb->nposres_xA = mtop->mols.index[mol+1];
- snew(molb->posres_xA, molb->nposres_xA);
- if (bFE)
- {
- molb->nposres_xB = molb->nposres_xA;
- snew(molb->posres_xB, molb->nposres_xB);
- }
- else
- {
- molb->nposres_xB = 0;
- }
- for (i = 0; i < il->nr; i += 2)
- {
- ip = &mtop->ffparams.iparams[il->iatoms[i]];
- a = il->iatoms[i+1];
- molb->posres_xA[a][XX] = ip->posres.pos0A[XX];
- molb->posres_xA[a][YY] = ip->posres.pos0A[YY];
- molb->posres_xA[a][ZZ] = ip->posres.pos0A[ZZ];
- if (bFE)
- {
- molb->posres_xB[a][XX] = ip->posres.pos0B[XX];
- molb->posres_xB[a][YY] = ip->posres.pos0B[YY];
- molb->posres_xB[a][ZZ] = ip->posres.pos0B[ZZ];
- }
- }
- if (il->nr == 0)
- {
- /* If only flat-bottomed posres are present, take reference pos from them.
- Here: bFE == FALSE */
- for (i = 0; i < ilfb->nr; i += 2)
- {
- ip = &mtop->ffparams.iparams[ilfb->iatoms[i]];
- a = ilfb->iatoms[i+1];
- molb->posres_xA[a][XX] = ip->fbposres.pos0[XX];
- molb->posres_xA[a][YY] = ip->fbposres.pos0[YY];
- molb->posres_xA[a][ZZ] = ip->fbposres.pos0[ZZ];
- }
- }
-}
-
static void set_disres_npair(gmx_mtop_t *mtop)
{
- t_iparams *ip;
- gmx_mtop_ilistloop_t iloop;
- t_ilist *ilist, *il;
- int nmol, i, npair;
- t_iatom *a;
+ gmx_mtop_ilistloop_t iloop;
+ int nmol;
- ip = mtop->ffparams.iparams;
+ gmx::ArrayRef<t_iparams> ip = mtop->ffparams.iparams;
iloop = gmx_mtop_ilistloop_init(mtop);
- while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
+ while (const InteractionLists *ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
{
- il = &ilist[F_DISRES];
+ const InteractionList &il = (*ilist)[F_DISRES];
- if (il->nr > 0)
+ if (il.size() > 0)
{
- a = il->iatoms;
- npair = 0;
- for (i = 0; i < il->nr; i += 3)
+ gmx::ArrayRef<const int> a = il.iatoms;
+ int npair = 0;
+ for (int i = 0; i < il.size(); i += 3)
{
npair++;
- if (i+3 == il->nr || ip[a[i]].disres.label != ip[a[i+3]].disres.label)
+ if (i+3 == il.size() || ip[a[i]].disres.label != ip[a[i+3]].disres.label)
{
ip[a[i]].disres.npair = npair;
npair = 0;
static void do_mtop(t_fileio *fio, gmx_mtop_t *mtop, gmx_bool bRead,
int file_version)
{
- int mt, mb;
- t_blocka dumb;
-
- if (bRead)
- {
- init_mtop(mtop);
- }
do_symtab(fio, &(mtop->symtab), bRead);
do_symstr(fio, &(mtop->name), bRead, &(mtop->symtab));
- if (file_version >= 57)
- {
- do_ffparams(fio, &mtop->ffparams, bRead, file_version);
+ do_ffparams(fio, &mtop->ffparams, bRead, file_version);
- gmx_fio_do_int(fio, mtop->nmoltype);
- }
- else
- {
- mtop->nmoltype = 1;
- }
+ int nmoltype = mtop->moltype.size();
+ gmx_fio_do_int(fio, nmoltype);
if (bRead)
{
- snew(mtop->moltype, mtop->nmoltype);
- if (file_version < 57)
- {
- mtop->moltype[0].name = mtop->name;
- }
+ mtop->moltype.resize(nmoltype);
}
- for (mt = 0; mt < mtop->nmoltype; mt++)
+ for (gmx_moltype_t &moltype : mtop->moltype)
{
- do_moltype(fio, &mtop->moltype[mt], bRead, &mtop->symtab, file_version,
- &mtop->groups);
+ do_moltype(fio, &moltype, bRead, &mtop->symtab, file_version);
}
- if (file_version >= 57)
- {
- gmx_fio_do_int(fio, mtop->nmolblock);
- }
- else
- {
- mtop->nmolblock = 1;
- }
+ int nmolblock = mtop->molblock.size();
+ gmx_fio_do_int(fio, nmolblock);
if (bRead)
{
- snew(mtop->molblock, mtop->nmolblock);
- }
- if (file_version >= 57)
- {
- for (mb = 0; mb < mtop->nmolblock; mb++)
- {
- do_molblock(fio, &mtop->molblock[mb], bRead);
- }
- gmx_fio_do_int(fio, mtop->natoms);
+ mtop->molblock.resize(nmolblock);
}
- else
+ for (gmx_molblock_t &molblock : mtop->molblock)
{
- mtop->molblock[0].type = 0;
- mtop->molblock[0].nmol = 1;
- mtop->molblock[0].natoms_mol = mtop->moltype[0].atoms.nr;
- mtop->molblock[0].nposres_xA = 0;
- mtop->molblock[0].nposres_xB = 0;
+ int numAtomsPerMolecule = (bRead ? 0 : mtop->moltype[molblock.type].atoms.nr);
+ do_molblock(fio, &molblock, numAtomsPerMolecule, bRead);
}
+ gmx_fio_do_int(fio, mtop->natoms);
if (file_version >= tpxv_IntermolecularBondeds)
{
{
if (bRead)
{
- snew(mtop->intermolecular_ilist, F_NRE);
+ mtop->intermolecular_ilist = gmx::compat::make_unique<InteractionLists>();
}
- do_ilists(fio, mtop->intermolecular_ilist, bRead, file_version);
+ do_ilists(fio, mtop->intermolecular_ilist.get(), bRead, file_version);
}
}
else
do_atomtypes (fio, &(mtop->atomtypes), bRead, file_version);
- if (file_version < 57)
- {
- do_idef (fio, &mtop->ffparams, &mtop->moltype[0], bRead, file_version);
- mtop->natoms = mtop->moltype[0].atoms.nr;
- }
-
if (file_version >= 65)
{
do_cmap(fio, &mtop->ffparams.cmap_grid, bRead);
}
else
{
- mtop->ffparams.cmap_grid.ngrid = 0;
mtop->ffparams.cmap_grid.grid_spacing = 0;
- mtop->ffparams.cmap_grid.cmapdata = nullptr;
- }
-
- if (file_version >= 57)
- {
- do_groups(fio, &mtop->groups, bRead, &(mtop->symtab), file_version);
+ mtop->ffparams.cmap_grid.cmapdata.clear();
}
- if (file_version < 57)
- {
- do_block(fio, &mtop->moltype[0].cgs, bRead, file_version);
- do_block(fio, &mtop->mols, bRead, file_version);
- /* Add the posres coordinates to the molblock */
- add_posres_molblock(mtop);
- }
- if (bRead)
- {
- if (file_version >= 57)
- {
- done_block(&mtop->mols);
- mtop->mols = mtop_mols(mtop);
- }
- }
+ do_groups(fio, &mtop->groups, bRead, &(mtop->symtab));
- if (file_version < 51)
- {
- /* Here used to be the shake blocks */
- do_blocka(fio, &dumb, bRead, file_version);
- if (dumb.nr > 0)
- {
- sfree(dumb.index);
- }
- if (dumb.nra > 0)
- {
- sfree(dumb.a);
- }
- }
+ mtop->haveMoleculeIndices = true;
if (bRead)
{
if (bRead)
{
gmx_fio_do_string(fio, buf);
- if (std::strncmp(buf, "VERSION", 7))
+ if (std::strncmp(buf, "VERSION", 7) != 0)
{
gmx_fatal(FARGS, "Can not read file %s,\n"
" this file is from a GROMACS version which is older than 2.0\n"
if ((precision != sizeof(float)) && !bDouble)
{
gmx_fatal(FARGS, "Unknown precision in file %s: real is %d bytes "
- "instead of %d or %d",
+ "instead of %zu or %zu",
gmx_fio_getname(fio), precision, sizeof(float), sizeof(double));
}
gmx_fio_setprecision(fio, bDouble);
gmx_fio_do_int(fio, tpx->fep_state);
}
gmx_fio_do_real(fio, tpx->lambda);
- gmx_fio_do_int(fio, tpx->bIr);
- gmx_fio_do_int(fio, tpx->bTop);
- gmx_fio_do_int(fio, tpx->bX);
- gmx_fio_do_int(fio, tpx->bV);
- gmx_fio_do_int(fio, tpx->bF);
- gmx_fio_do_int(fio, tpx->bBox);
+ gmx_fio_do_gmx_bool(fio, tpx->bIr);
+ gmx_fio_do_gmx_bool(fio, tpx->bTop);
+ gmx_fio_do_gmx_bool(fio, tpx->bX);
+ gmx_fio_do_gmx_bool(fio, tpx->bV);
+ gmx_fio_do_gmx_bool(fio, tpx->bF);
+ gmx_fio_do_gmx_bool(fio, tpx->bBox);
if ((fileGeneration > tpx_generation))
{
gmx_mtop_t *mtop)
{
t_tpxheader tpx;
- gmx_mtop_t dum_top;
gmx_bool TopOnlyOK;
int ePBC;
gmx_bool bPeriodicMols;
tpx.ngtc = state->ngtc;
tpx.fep_state = state->fep_state;
tpx.lambda = state->lambda[efptFEP];
- tpx.bIr = (ir != nullptr);
- tpx.bTop = (mtop != nullptr);
- tpx.bX = (state->flags & (1 << estX));
- tpx.bV = (state->flags & (1 << estV));
+ tpx.bIr = ir != nullptr;
+ tpx.bTop = mtop != nullptr;
+ tpx.bX = (state->flags & (1 << estX)) != 0;
+ tpx.bV = (state->flags & (1 << estV)) != 0;
tpx.bF = FALSE;
tpx.bBox = TRUE;
}
if (x == nullptr)
{
- x = as_rvec_array(state->x.data());
- v = as_rvec_array(state->v.data());
+ x = state->x.rvec_array();
+ v = state->v.rvec_array();
}
-#define do_test(fio, b, p) if (bRead && (p != NULL) && !b) gmx_fatal(FARGS, "No %s in %s",#p, gmx_fio_getname(fio))
+#define do_test(fio, b, p) if (bRead && ((p) != NULL) && !(b)) gmx_fatal(FARGS, "No %s in %s",#p, gmx_fio_getname(fio))
do_test(fio, tpx.bBox, state->box);
if (tpx.bBox)
}
else
{
+ gmx_mtop_t dum_top;
do_mtop(fio, &dum_top, bRead, fileVersion);
- done_mtop(&dum_top);
}
}
do_test(fio, tpx.bX, x);
}
if (fileGeneration <= tpx_generation && ir)
{
- do_inputrec(fio, ir, bRead, fileVersion, mtop ? &mtop->ffparams.fudgeQQ : nullptr);
+ do_inputrec(fio, ir, bRead, fileVersion);
if (fileVersion < 51)
{
set_box_rel(ir, state);
{
if (fileVersion < 57)
{
- if (mtop->moltype[0].ilist[F_DISRES].nr > 0)
+ if (mtop->moltype[0].ilist[F_DISRES].size() > 0)
{
ir->eDisre = edrSimple;
}
fio = open_tpx(fn, "r");
ePBC = do_tpx(fio, TRUE, ir, &state, x, v, mtop);
close_tpx(fio);
- *natoms = mtop->natoms;
+ if (mtop != nullptr && natoms != nullptr)
+ {
+ *natoms = mtop->natoms;
+ }
if (box)
{
copy_mat(state.box, box);
struct t_tpxheader
{
- int bIr; /* Non zero if input_rec is present */
- int bBox; /* Non zero if a box is present */
- int bTop; /* Non zero if a topology is present */
- int bX; /* Non zero if coordinates are present */
- int bV; /* Non zero if velocities are present */
- int bF; /* Non zero if forces are present (no longer
+ bool bIr; /* Non zero if input_rec is present */
+ bool bBox; /* Non zero if a box is present */
+ bool bTop; /* Non zero if a topology is present */
+ bool bX; /* Non zero if coordinates are present */
+ bool bV; /* Non zero if velocities are present */
+ bool bF; /* Non zero if forces are present (no longer
supported, but retained so old .tpr can be read) */
int natoms; /* The total number of atoms */
void read_tpx_state(const char *fn,
t_inputrec *ir, t_state *state,
gmx_mtop_t *mtop);
+
+ /*! \brief
+ * Read a file and close it again.
+ *
+ * Reads a topology input file and populates the fields if the passed
+ * variables are valid. It is possible to pass \p ir, \p natoms,
+ * \p x, \p v or \p mtop as nullptr to the function. In those cases,
+ * the variables will not be populated from the input file. Passing both
+ * \p x and \p v as nullptr is not supported. If both \p natoms and
+ * \p mtop are passed as valid objects to the function, the total atom
+ * number from \p mtop will be set in \p natoms. Otherwise \p natoms
+ * will not be changed. If \p box is valid, the box will be set from
+ * the information read in from the file.
+ *
+ * \param[in] fn Input file name.
+ * \param[out] ir Input parameters to be set, or nullptr.
+ * \param[out] box Box matrix.
+ * \param[out] natoms Total atom numbers to be set, or nullptr.
+ * \param[out] x Positions to be filled from file, or nullptr.
+ * \param[out] v Velocities to be filled from file, or nullptr.
+ * \param[out] mtop Topology to be populated, or nullptr.
+ * \returns ir->ePBC if it was read from the file.
+ */
int read_tpx(const char *fn,
t_inputrec *ir, matrix box, int *natoms,
rvec *x, rvec *v, gmx_mtop_t *mtop);
- /* Read a file, and close it again.
- * When step, t or lambda are NULL they will not be stored.
- * Returns ir->ePBC, if it could be read from the file.
- */
int read_tpx_top(const char *fn,
t_inputrec *ir, matrix box, int *natoms,
"Last residue in helix" }
};
- typedef struct {
+ typedef struct { //NOLINT(clang-analyzer-optin.performance.Padding)
FILE *fp, *fp2;
gmx_bool bfp2;
const char *filenm;
/* Read reference frame from tpx file to compute helix length */
snew(xref, top->atoms.nr);
read_tpx(ftp2fn(efTPR, NFILE, fnm),
- nullptr, nullptr, &natoms, xref, nullptr, nullptr);
+ nullptr, nullptr, nullptr, xref, nullptr, nullptr);
calc_hxprops(nres, bb, xref);
do_start_end(nres, bb, &nbb, bbindex, &nca, caindex, bRange, rStart, rEnd);
sfree(xref);
static void periodic_dist(int ePBC,
- matrix box, rvec x[], int n, int index[],
+ matrix box, rvec x[], int n, const int index[],
real *rmin, real *rmax, int *min_ind)
{
#define NSHIFT_MAX 26
{
gmx_fatal(FARGS, "pbc = %s is not supported by g_mindist",
epbc_names[ePBC]);
- nsz = 0; /* Keep compilers quiet */
}
nshift = 0;
snew(leg, 1);
sprintf(buf, "Internal in %s", grpn[0]);
leg[0] = gmx_strdup(buf);
- xvgr_legend(dist, 0, (const char**)leg, oenv);
+ xvgr_legend(dist, 0, leg, oenv);
if (num)
{
- xvgr_legend(num, 0, (const char**)leg, oenv);
+ xvgr_legend(num, 0, leg, oenv);
}
}
else
leg[j] = gmx_strdup(buf);
}
}
- xvgr_legend(dist, j, (const char**)leg, oenv);
+ xvgr_legend(dist, j, leg, oenv);
if (num)
{
- xvgr_legend(num, j, (const char**)leg, oenv);
+ xvgr_legend(num, j, leg, oenv);
}
}
}
sprintf(buf, "%s-%s", grpn[0], grpn[i+1]);
leg[i] = gmx_strdup(buf);
}
- xvgr_legend(dist, ng-1, (const char**)leg, oenv);
+ xvgr_legend(dist, ng-1, leg, oenv);
if (num)
{
- xvgr_legend(num, ng-1, (const char**)leg, oenv);
+ xvgr_legend(num, ng-1, leg, oenv);
}
}
{
sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
respertime = xvgropen(rfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
- xvgr_legend(respertime, ng-1, (const char**)leg, oenv);
+ xvgr_legend(respertime, ng-1, leg, oenv);
if (bPrintResName && output_env_get_print_xvgr_codes(oenv) )
{
fprintf(respertime, "# ");
sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
res = xvgropen(rfile, buf, "Residue (#)", "Distance (nm)", oenv);
- xvgr_legend(res, ng-1, (const char**)leg, oenv);
+ xvgr_legend(res, ng-1, leg, oenv);
for (j = 0; j < nres; j++)
{
fprintf(res, "%4d", j+1);
sfree(x0);
}
-static int find_residues(const t_atoms *atoms, int n, int index[], int **resindex)
+static int find_residues(const t_atoms *atoms, int n, const int index[], int **resindex)
{
int i;
int nres = 0, resnr, presnr = 0;
dump_res(debug, nres, residues, index[0]);
}
}
+ else if (bEachResEachTime || bPrintResName)
+ {
+ gmx_fatal(FARGS, "Option -or needs to be set to print residues");
+ }
if (bPI)
{
#include "config.h"
-#include <assert.h>
-#include <string.h>
-
+#include <cassert>
#include <cmath>
+#include <cstring>
#include "gromacs/domdec/domdec.h"
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/listed-forces/listed-forces.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
+#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/forcerec-threading.h"
-#include "gromacs/mdlib/genborn.h"
#include "gromacs/mdlib/mdrun.h"
#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/qmmm.h"
+#include "gromacs/mdlib/rf_util.h"
+#include "gromacs/mdlib/wall.h"
#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/enerdata.h"
#include "gromacs/mdtypes/forceoutput.h"
+#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/pbcutil/ishift.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
-void ns(FILE *fp,
- t_forcerec *fr,
- matrix box,
- gmx_groups_t *groups,
- gmx_localtop_t *top,
- t_mdatoms *md,
- t_commrec *cr,
- t_nrnb *nrnb,
- gmx_bool bFillGrid)
+void ns(FILE *fp,
+ t_forcerec *fr,
+ matrix box,
+ const gmx_groups_t *groups,
+ gmx_localtop_t *top,
+ const t_mdatoms *md,
+ const t_commrec *cr,
+ t_nrnb *nrnb,
+ gmx_bool bFillGrid)
{
int nsearch;
}
}
-void do_force_lowlevel(t_forcerec *fr, t_inputrec *ir,
- t_idef *idef, t_commrec *cr,
- t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- t_mdatoms *md,
- rvec x[], history_t *hist,
- rvec *forceForUseWithShiftForces,
+void do_force_lowlevel(t_forcerec *fr,
+ const t_inputrec *ir,
+ const t_idef *idef,
+ const t_commrec *cr,
+ const gmx_multisim_t *ms,
+ t_nrnb *nrnb,
+ gmx_wallcycle_t wcycle,
+ const t_mdatoms *md,
+ rvec x[],
+ history_t *hist,
+ rvec *forceForUseWithShiftForces,
gmx::ForceWithVirial *forceWithVirial,
- gmx_enerdata_t *enerd,
- t_fcdata *fcd,
- gmx_localtop_t *top,
- gmx_genborn_t *born,
- gmx_bool bBornRadii,
- matrix box,
- t_lambda *fepvals,
- real *lambda,
- t_graph *graph,
- t_blocka *excl,
- rvec mu_tot[],
- int flags,
- float *cycles_pme)
+ gmx_enerdata_t *enerd,
+ t_fcdata *fcd,
+ matrix box,
+ t_lambda *fepvals,
+ real *lambda,
+ const t_graph *graph,
+ const t_blocka *excl,
+ rvec mu_tot[],
+ int flags,
+ float *cycles_pme)
{
int i, j;
int donb_flags;
if (ir->nwall)
{
/* foreign lambda component for walls */
- real dvdl_walls = do_walls(ir, fr, box, md, x, forceForUseWithShiftForces, lambda[efptVDW],
+ real dvdl_walls = do_walls(*ir, *fr, box, *md, x,
+ forceWithVirial, lambda[efptVDW],
enerd->grpp.ener[egLJSR], nrnb);
enerd->dvdl_lin[efptVDW] += dvdl_walls;
}
- /* If doing GB, reset dvda and calculate the Born radii */
- if (ir->implicit_solvent)
- {
- wallcycle_sub_start(wcycle, ewcsNONBONDED);
-
- for (i = 0; i < born->nr; i++)
- {
- fr->dvda[i] = 0;
- }
-
- if (bBornRadii)
- {
- calc_gb_rad(cr, fr, ir, top, x, fr->gblist, born, md, nrnb);
- }
-
- wallcycle_sub_stop(wcycle, ewcsNONBONDED);
- }
-
- where();
/* We only do non-bonded calculation with group scheme here, the verlet
* calls are done from do_force_cutsVERLET(). */
if (fr->cutoff_scheme == ecutsGROUP && (flags & GMX_FORCE_NONBONDED))
}
}
wallcycle_sub_stop(wcycle, ewcsNONBONDED);
- where();
- }
-
- /* If we are doing GB, calculate bonded forces and apply corrections
- * to the solvation forces */
- /* MRS: Eventually, many need to include free energy contribution here! */
- if (ir->implicit_solvent)
- {
- wallcycle_sub_start(wcycle, ewcsLISTED);
- calc_gb_forces(cr, md, born, top, x, forceForUseWithShiftForces, fr, idef,
- ir->gb_algorithm, ir->sa_algorithm, nrnb, &pbc, graph, enerd);
- wallcycle_sub_stop(wcycle, ewcsLISTED);
}
#if GMX_MPI
TRUE, box);
}
- do_force_listed(wcycle, box, ir->fepvals, cr,
- idef, (const rvec *) x, hist,
+ do_force_listed(wcycle, box, ir->fepvals, cr, ms,
+ idef, x, hist,
forceForUseWithShiftForces, forceWithVirial,
fr, &pbc, graph, enerd, nrnb, lambda, md, fcd,
- DOMAINDECOMP(cr) ? cr->dd->gatindex : nullptr,
+ DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr,
flags);
- where();
*cycles_pme = 0;
md->sqrt_c6A, md->sqrt_c6B,
md->sigmaA, md->sigmaB,
md->sigma3A, md->sigma3B,
- md->nChargePerturbed || md->nTypePerturbed,
+ (md->nChargePerturbed != 0) || (md->nTypePerturbed != 0),
ir->cutoff_scheme != ecutsVERLET,
excl, x, box, mu_tot,
ir->ewald_geometry,
{
real dvdl_rf_excl = 0;
enerd->term[F_RF_EXCL] =
- RF_excl_correction(fr, graph, md, excl, x, forceForUseWithShiftForces,
+ RF_excl_correction(fr, graph, md, excl, DOMAINDECOMP(cr),
+ x, forceForUseWithShiftForces,
fr->fshift, &pbc, lambda[efptCOUL], &dvdl_rf_excl);
enerd->dvdl_lin[efptCOUL] += dvdl_rf_excl;
}
}
- where();
if (debug)
{
}
}
-static real sum_v(int n, real v[])
+static real sum_v(int n, const real v[])
{
real t;
int i;
epot[F_LJ] = sum_v(grpp->nener, grpp->ener[egLJSR]);
epot[F_LJ14] = sum_v(grpp->nener, grpp->ener[egLJ14]);
epot[F_COUL14] = sum_v(grpp->nener, grpp->ener[egCOUL14]);
- /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
- epot[F_GBPOL] += sum_v(grpp->nener, grpp->ener[egGB]);
/* lattice part of LR doesnt belong to any group
* and has been added earlier
void sum_dhdl(gmx_enerdata_t *enerd, gmx::ArrayRef<const real> lambda, t_lambda *fepvals)
{
int index;
- double dlam;
enerd->dvdl_lin[efptVDW] += enerd->term[F_DVDL_VDW]; /* include dispersion correction */
enerd->term[F_DVDL] = 0.0;
}
}
- /* Notes on the foreign lambda free energy difference evaluation:
- * Adding the potential and ekin terms that depend linearly on lambda
- * as delta lam * dvdl to the energy differences is exact.
- * For the constraints this is not exact, but we have no other option
- * without literally changing the lengths and reevaluating the energies at each step.
- * (try to remedy this post 4.6 - MRS)
- */
if (fepvals->separate_dvdl[efptBONDED])
{
enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
{
enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
}
- enerd->term[F_DVDL_CONSTR] = 0;
for (int i = 0; i < fepvals->n_lambda; i++)
{
current lambda, because the contributions to the current
lambda are automatically zeroed */
- for (size_t j = 0; j < lambda.size(); j++)
+ double &enerpart_lambda = enerd->enerpart_lambda[i + 1];
+
+ for (gmx::index j = 0; j < lambda.size(); j++)
{
/* Note that this loop is over all dhdl components, not just the separated ones */
- dlam = (fepvals->all_lambda[j][i] - lambda[j]);
- enerd->enerpart_lambda[i+1] += dlam*enerd->dvdl_lin[j];
+ const double dlam = fepvals->all_lambda[j][i] - lambda[j];
+
+ enerpart_lambda += dlam*enerd->dvdl_lin[j];
+
+ /* Constraints can not be evaluated at foreign lambdas, so we add
+ * a linear extrapolation. This is an approximation, but usually
+ * quite accurate since constraints change little between lambdas.
+ */
+ if ((j == efptBONDED && fepvals->separate_dvdl[efptBONDED]) ||
+ (j == efptFEP && !fepvals->separate_dvdl[efptBONDED]))
+ {
+ enerpart_lambda += dlam*enerd->term[F_DVDL_CONSTR];
+ }
+
+ if (j == efptMASS)
+ {
+ enerpart_lambda += dlam*enerd->term[F_DKDL];
+ }
+
if (debug)
{
fprintf(debug, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
fepvals->all_lambda[j][i], efpt_names[j],
- (enerd->enerpart_lambda[i+1] - enerd->enerpart_lambda[0]),
+ enerpart_lambda - enerd->enerpart_lambda[0],
dlam, enerd->dvdl_lin[j]);
}
}
}
+
+ /* The constrain contribution is now included in other terms, so clear it */
+ enerd->term[F_DVDL_CONSTR] = 0;
}
--- /dev/null
- state->natoms -= n;
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2010,2011,2012,2013,2014,2015,2016,2017,2018, 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.
+ *
+ * GROMACS 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.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "membed.h"
+
+#include <csignal>
+#include <cstdlib>
+
+#include "gromacs/commandline/filenm.h"
+#include "gromacs/fileio/readinp.h"
+#include "gromacs/fileio/warninp.h"
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/state.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/topology/mtop_lookup.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/filestream.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
+
+/* information about scaling center */
+typedef struct {
+ rvec xmin; /* smallest coordinates of all embedded molecules */
+ rvec xmax; /* largest coordinates of all embedded molecules */
+ rvec *geom_cent; /* scaling center of each independent molecule to embed */
+ int pieces; /* number of molecules to embed independently */
+ int *nidx; /* n atoms for every independent embedded molecule (index in subindex) */
+ int **subindex; /* atomids for independent molecule *
+ * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */
+} pos_ins_t;
+
+/* variables needed in do_md */
+struct gmx_membed_t {
+ int it_xy; /* number of iterations (steps) used to grow something in the xy-plane */
+ int it_z; /* same, but for z */
+ real xy_step; /* stepsize used during resize in xy-plane */
+ real z_step; /* same, but in z */
+ rvec fac; /* initial scaling of the molecule to grow into the membrane */
+ rvec *r_ins; /* final coordinates of the molecule to grow */
+ pos_ins_t *pos_ins; /* scaling center for each piece to embed */
+};
+
+/* membrane related variables */
+typedef struct {
+ char *name; /* name of index group to embed molecule into (usually membrane) */
+ t_block mem_at; /* list all atoms in membrane */
+ int nmol; /* number of membrane molecules overlapping with the molecule to embed */
+ int *mol_id; /* list of molecules in membrane that overlap with the molecule to embed */
+ real lip_area; /* average area per lipid in membrane (only correct for homogeneous bilayers)*/
+ real zmin; /* minimum z coordinate of membrane */
+ real zmax; /* maximum z coordinate of membrane */
+ real zmed; /* median z coordinate of membrane */
+} mem_t;
+
+/* Lists all molecules in the membrane that overlap with the molecule to be embedded. *
+ * These will then be removed from the system */
+typedef struct {
+ int nr; /* number of molecules to remove */
+ int *mol; /* list of molecule ids to remove */
+ int *block; /* id of the molblock that the molecule to remove is part of */
+} rm_t;
+
+/* Get the global molecule id, and the corresponding molecule type and id of the *
+ * molblock from the global atom nr. */
+static int get_mol_id(int at, gmx_mtop_t *mtop, int *type, int *block)
+{
+ int mol_id = 0;
+ int i;
+ int atnr_mol;
+
+ *block = 0;
+ mtopGetMolblockIndex(mtop, at, block, &mol_id, &atnr_mol);
+ for (i = 0; i < *block; i++)
+ {
+ mol_id += mtop->molblock[i].nmol;
+ }
+ *type = mtop->molblock[*block].type;
+
+ return mol_id;
+}
+
+/* Get the id of the molblock from a global molecule id */
+static int get_molblock(int mol_id, const std::vector<gmx_molblock_t> &mblock)
+{
+ int nmol = 0;
+
+ for (size_t i = 0; i < mblock.size(); i++)
+ {
+ nmol += mblock[i].nmol;
+ if (mol_id < nmol)
+ {
+ return i;
+ }
+ }
+
+ gmx_fatal(FARGS, "mol_id %d larger than total number of molecules %d.\n", mol_id, nmol);
+}
+
+/* Get a list of all the molecule types that are present in a group of atoms. *
+ * Because all interaction within the group to embed are removed on the topology *
+ * level, if the same molecule type is found in another part of the system, these *
+ * would also be affected. Therefore we have to check if the embedded and rest group *
+ * share common molecule types. If so, membed will stop with an error. */
+static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
+{
+ int i, j, nr;
+ int type = 0, block = 0;
+ gmx_bool bNEW;
+
+ nr = 0;
+ snew(tlist->index, at->nr);
+ for (i = 0; i < at->nr; i++)
+ {
+ bNEW = TRUE;
+ get_mol_id(at->index[i], mtop, &type, &block);
+ for (j = 0; j < nr; j++)
+ {
+ if (tlist->index[j] == type)
+ {
+ bNEW = FALSE;
+ }
+ }
+
+ if (bNEW)
+ {
+ tlist->index[nr] = type;
+ nr++;
+ }
+ }
+ srenew(tlist->index, nr);
+ return nr;
+}
+
+/* Do the actual check of the molecule types between embedded and rest group */
+static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop)
+{
+ t_block *ins_mtype, *rest_mtype;
+ int i, j;
+
+ snew(ins_mtype, 1);
+ snew(rest_mtype, 1);
+ ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype );
+ rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
+
+ for (i = 0; i < ins_mtype->nr; i++)
+ {
+ for (j = 0; j < rest_mtype->nr; j++)
+ {
+ if (ins_mtype->index[i] == rest_mtype->index[j])
+ {
+ gmx_fatal(FARGS, "Moleculetype %s is found both in the group to insert and the rest of the system.\n"
+ "1. Your *.ndx and *.top do not match\n"
+ "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
+ "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
+ "we need to exclude all interactions between the atoms in the group to\n"
+ "insert, the same moleculetype can not be used in both groups. Change the\n"
+ "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
+ "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name),
+ *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name));
+ }
+ }
+ }
+
+ done_block(ins_mtype);
+ done_block(rest_mtype);
+ sfree(ins_mtype);
+ sfree(rest_mtype);
+}
+
+static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
+ int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
+ int *pieces, gmx_bool *bALLOW_ASYMMETRY)
+{
+ warninp_t wi;
+ std::vector <t_inpfile> inp;
+
+ wi = init_warning(TRUE, 0);
+
+ {
+ gmx::TextInputFile stream(membed_input);
+ inp = read_inpfile(&stream, membed_input, wi);
+ stream.close();
+ }
+ *it_xy = get_eint(&inp, "nxy", 1000, wi);
+ *it_z = get_eint(&inp, "nz", 0, wi);
+ *xy_fac = get_ereal(&inp, "xyinit", 0.5, wi);
+ *xy_max = get_ereal(&inp, "xyend", 1.0, wi);
+ *z_fac = get_ereal(&inp, "zinit", 1.0, wi);
+ *z_max = get_ereal(&inp, "zend", 1.0, wi);
+ *probe_rad = get_ereal(&inp, "rad", 0.22, wi);
+ *low_up_rm = get_eint(&inp, "ndiff", 0, wi);
+ *maxwarn = get_eint(&inp, "maxwarn", 0, wi);
+ *pieces = get_eint(&inp, "pieces", 1, wi);
+ const char *yesno_names[BOOL_NR+1] =
+ {
+ "no", "yes", nullptr
+ };
+ *bALLOW_ASYMMETRY = (get_eeenum(&inp, "asymmetry", yesno_names, wi) != 0);
+
+ check_warning_error(wi, FARGS);
+ {
+ gmx::TextOutputFile stream(membed_input);
+ write_inpfile(&stream, membed_input, &inp, FALSE, WriteMdpHeader::yes, wi);
+ stream.close();
+ }
+ done_warning(wi, FARGS);
+}
+
+/* Obtain the maximum and minimum coordinates of the group to be embedded */
+static int init_ins_at(t_block *ins_at, t_block *rest_at, t_state *state, pos_ins_t *pos_ins,
+ gmx_groups_t *groups, int ins_grp_id, real xy_max)
+{
+ int i, gid, c = 0;
+ real xmin, xmax, ymin, ymax, zmin, zmax;
+ const real min_memthick = 6.0; /* minimum thickness of the bilayer that will be used to *
+ * determine the overlap between molecule to embed and membrane */
+ const real fac_inp_size = 1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */
+ snew(rest_at->index, state->natoms);
+ auto x = makeArrayRef(state->x);
+
+ xmin = xmax = x[ins_at->index[0]][XX];
+ ymin = ymax = x[ins_at->index[0]][YY];
+ zmin = zmax = x[ins_at->index[0]][ZZ];
+
+ for (i = 0; i < state->natoms; i++)
+ {
+ gid = groups->grpnr[egcFREEZE][i];
+ if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id)
+ {
+ xmin = std::min(xmin, x[i][XX]);
+ xmax = std::max(xmax, x[i][XX]);
+ ymin = std::min(ymin, x[i][YY]);
+ ymax = std::max(ymax, x[i][YY]);
+ zmin = std::min(zmin, x[i][ZZ]);
+ zmax = std::max(zmax, x[i][ZZ]);
+ }
+ else
+ {
+ rest_at->index[c] = i;
+ c++;
+ }
+ }
+
+ rest_at->nr = c;
+ srenew(rest_at->index, c);
+
+ if (xy_max > fac_inp_size)
+ {
+ pos_ins->xmin[XX] = xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
+ pos_ins->xmin[YY] = ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
+
+ pos_ins->xmax[XX] = xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
+ pos_ins->xmax[YY] = ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
+ }
+ else
+ {
+ pos_ins->xmin[XX] = xmin;
+ pos_ins->xmin[YY] = ymin;
+
+ pos_ins->xmax[XX] = xmax;
+ pos_ins->xmax[YY] = ymax;
+ }
+
+ if ( (zmax-zmin) < min_memthick)
+ {
+ pos_ins->xmin[ZZ] = zmin+(zmax-zmin)/2.0-0.5*min_memthick;
+ pos_ins->xmax[ZZ] = zmin+(zmax-zmin)/2.0+0.5*min_memthick;
+ }
+ else
+ {
+ pos_ins->xmin[ZZ] = zmin;
+ pos_ins->xmax[ZZ] = zmax;
+ }
+
+ return c;
+}
+
+/* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the *
+ * xy-plane and counting the number of occupied grid points */
+static real est_prot_area(pos_ins_t *pos_ins, rvec *r, t_block *ins_at, mem_t *mem_p)
+{
+ real x, y, dx = 0.15, dy = 0.15, area = 0.0;
+ real add, memmin, memmax;
+ int c, at;
+
+ /* min and max membrane coordinate are altered to reduce the influence of the *
+ * boundary region */
+ memmin = mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin);
+ memmax = mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin);
+
+ //NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
+ for (x = pos_ins->xmin[XX]; x < pos_ins->xmax[XX]; x += dx)
+ {
+ //NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
+ for (y = pos_ins->xmin[YY]; y < pos_ins->xmax[YY]; y += dy)
+ {
+ c = 0;
+ add = 0.0;
+ do
+ {
+ at = ins_at->index[c];
+ if ( (r[at][XX] >= x) && (r[at][XX] < x+dx) &&
+ (r[at][YY] >= y) && (r[at][YY] < y+dy) &&
+ (r[at][ZZ] > memmin) && (r[at][ZZ] < memmax) )
+ {
+ add = 1.0;
+ }
+ c++;
+ }
+ while ( (c < ins_at->nr) && (add < 0.5) );
+ area += add;
+ }
+ }
+ area = area*dx*dy;
+
+ return area;
+}
+
+static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
+{
+ int i, j, at, mol, nmol, nmolbox, count;
+ t_block *mem_a;
+ real z, zmin, zmax, mem_area;
+ gmx_bool bNew;
+ int *mol_id;
+ int type = 0, block = 0;
+
+ nmol = count = 0;
+ mem_a = &(mem_p->mem_at);
+ snew(mol_id, mem_a->nr);
+ zmin = pos_ins->xmax[ZZ];
+ zmax = pos_ins->xmin[ZZ];
+ for (i = 0; i < mem_a->nr; i++)
+ {
+ at = mem_a->index[i];
+ if ( (r[at][XX] > pos_ins->xmin[XX]) && (r[at][XX] < pos_ins->xmax[XX]) &&
+ (r[at][YY] > pos_ins->xmin[YY]) && (r[at][YY] < pos_ins->xmax[YY]) &&
+ (r[at][ZZ] > pos_ins->xmin[ZZ]) && (r[at][ZZ] < pos_ins->xmax[ZZ]) )
+ {
+ mol = get_mol_id(at, mtop, &type, &block);
+ bNew = TRUE;
+ for (j = 0; j < nmol; j++)
+ {
+ if (mol == mol_id[j])
+ {
+ bNew = FALSE;
+ }
+ }
+
+ if (bNew)
+ {
+ mol_id[nmol] = mol;
+ nmol++;
+ }
+
+ z = r[at][ZZ];
+ if (z < zmin)
+ {
+ zmin = z;
+ }
+
+ if (z > zmax)
+ {
+ zmax = z;
+ }
+
+ count++;
+ }
+ }
+
+ mem_p->nmol = nmol;
+ srenew(mol_id, nmol);
+ mem_p->mol_id = mol_id;
+
+ if ((zmax-zmin) > (box[ZZ][ZZ]-0.5))
+ {
+ gmx_fatal(FARGS, "Something is wrong with your membrane. Max and min z values are %f and %f.\n"
+ "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
+ "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
+ "your water layer is not thick enough.\n", zmax, zmin);
+ }
+ mem_p->zmin = zmin;
+ mem_p->zmax = zmax;
+ mem_p->zmed = (zmax-zmin)/2+zmin;
+
+ /*number of membrane molecules in protein box*/
+ nmolbox = count/mtop->moltype[mtop->molblock[block].type].atoms.nr;
+ /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/
+ mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
+ /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/
+ mem_p->lip_area = 2.0*mem_area/static_cast<double>(nmolbox);
+
+ return mem_p->mem_at.nr;
+}
+
+static void init_resize(t_block *ins_at, rvec *r_ins, pos_ins_t *pos_ins, mem_t *mem_p, rvec *r,
+ gmx_bool bALLOW_ASYMMETRY)
+{
+ int i, j, at, c, outsidesum, gctr = 0;
+ int idxsum = 0;
+
+ /*sanity check*/
+ for (i = 0; i < pos_ins->pieces; i++)
+ {
+ idxsum += pos_ins->nidx[i];
+ }
+
+ if (idxsum != ins_at->nr)
+ {
+ gmx_fatal(FARGS, "Piecewise sum of inserted atoms not same as size of group selected to insert.");
+ }
+
+ snew(pos_ins->geom_cent, pos_ins->pieces);
+ for (i = 0; i < pos_ins->pieces; i++)
+ {
+ c = 0;
+ outsidesum = 0;
+ for (j = 0; j < DIM; j++)
+ {
+ pos_ins->geom_cent[i][j] = 0;
+ }
+
+ for (j = 0; j < pos_ins->nidx[i]; j++)
+ {
+ at = pos_ins->subindex[i][j];
+ copy_rvec(r[at], r_ins[gctr]);
+ if ( (r_ins[gctr][ZZ] < mem_p->zmax) && (r_ins[gctr][ZZ] > mem_p->zmin) )
+ {
+ rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]);
+ c++;
+ }
+ else
+ {
+ outsidesum++;
+ }
+ gctr++;
+ }
+
+ if (c > 0)
+ {
+ svmul(1/static_cast<double>(c), pos_ins->geom_cent[i], pos_ins->geom_cent[i]);
+ }
+
+ if (!bALLOW_ASYMMETRY)
+ {
+ pos_ins->geom_cent[i][ZZ] = mem_p->zmed;
+ }
+
+ fprintf(stderr, "Embedding piece %d with center of geometry: %f %f %f\n",
+ i, pos_ins->geom_cent[i][XX], pos_ins->geom_cent[i][YY], pos_ins->geom_cent[i][ZZ]);
+ }
+ fprintf(stderr, "\n");
+}
+
+/* resize performed in the md loop */
+static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, const rvec fac)
+{
+ int i, j, k, at, c = 0;
+ for (k = 0; k < pos_ins->pieces; k++)
+ {
+ for (i = 0; i < pos_ins->nidx[k]; i++)
+ {
+ at = pos_ins->subindex[k][i];
+ for (j = 0; j < DIM; j++)
+ {
+ r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
+ }
+ c++;
+ }
+ }
+}
+
+/* generate the list of membrane molecules that overlap with the molecule to be embedded. *
+ * The molecule to be embedded is already reduced in size. */
+static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc, gmx_mtop_t *mtop,
+ rvec *r, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad,
+ int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
+{
+ int i, j, k, l, at, at2, mol_id;
+ int type = 0, block = 0;
+ int nrm, nupper, nlower;
+ real r_min_rad, z_lip, min_norm;
+ gmx_bool bRM;
+ rvec dr, dr_tmp;
+ real *dist;
+ int *order;
+
+ r_min_rad = probe_rad*probe_rad;
+ gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop);
+ snew(rm_p->block, molecules.numBlocks());
+ nrm = nupper = 0;
+ nlower = low_up_rm;
+ for (i = 0; i < ins_at->nr; i++)
+ {
+ at = ins_at->index[i];
+ for (j = 0; j < rest_at->nr; j++)
+ {
+ at2 = rest_at->index[j];
+ pbc_dx(pbc, r[at], r[at2], dr);
+
+ if (norm2(dr) < r_min_rad)
+ {
+ mol_id = get_mol_id(at2, mtop, &type, &block);
+ bRM = TRUE;
+ for (l = 0; l < nrm; l++)
+ {
+ if (rm_p->mol[l] == mol_id)
+ {
+ bRM = FALSE;
+ }
+ }
+
+ if (bRM)
+ {
+ rm_p->mol[nrm] = mol_id;
+ rm_p->block[nrm] = block;
+ nrm++;
+ z_lip = 0.0;
+ for (l = 0; l < mem_p->nmol; l++)
+ {
+ if (mol_id == mem_p->mol_id[l])
+ {
+ for (int k : molecules.block(mol_id))
+ {
+ z_lip += r[k][ZZ];
+ }
+ z_lip /= molecules.block(mol_id).size();
+ if (z_lip < mem_p->zmed)
+ {
+ nlower++;
+ }
+ else
+ {
+ nupper++;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /*make sure equal number of lipids from upper and lower layer are removed */
+ if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) )
+ {
+ snew(dist, mem_p->nmol);
+ snew(order, mem_p->nmol);
+ for (i = 0; i < mem_p->nmol; i++)
+ {
+ at = molecules.block(mem_p->mol_id[i]).begin();
+ pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr);
+ if (pos_ins->pieces > 1)
+ {
+ /*minimum dr value*/
+ min_norm = norm2(dr);
+ for (k = 1; k < pos_ins->pieces; k++)
+ {
+ pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp);
+ if (norm2(dr_tmp) < min_norm)
+ {
+ min_norm = norm2(dr_tmp);
+ copy_rvec(dr_tmp, dr);
+ }
+ }
+ }
+ dist[i] = dr[XX]*dr[XX]+dr[YY]*dr[YY];
+ j = i-1;
+ while (j >= 0 && dist[i] < dist[order[j]])
+ {
+ order[j+1] = order[j];
+ j--;
+ }
+ order[j+1] = i;
+ }
+
+ i = 0;
+ while (nupper != nlower)
+ {
+ mol_id = mem_p->mol_id[order[i]];
+ block = get_molblock(mol_id, mtop->molblock);
+ bRM = TRUE;
+ for (l = 0; l < nrm; l++)
+ {
+ if (rm_p->mol[l] == mol_id)
+ {
+ bRM = FALSE;
+ }
+ }
+
+ if (bRM)
+ {
+ z_lip = 0;
+ for (int k : molecules.block(mol_id))
+ {
+ z_lip += r[k][ZZ];
+ }
+ z_lip /= molecules.block(mol_id).size();
+ if (nupper > nlower && z_lip < mem_p->zmed)
+ {
+ rm_p->mol[nrm] = mol_id;
+ rm_p->block[nrm] = block;
+ nrm++;
+ nlower++;
+ }
+ else if (nupper < nlower && z_lip > mem_p->zmed)
+ {
+ rm_p->mol[nrm] = mol_id;
+ rm_p->block[nrm] = block;
+ nrm++;
+ nupper++;
+ }
+ }
+ i++;
+
+ if (i > mem_p->nmol)
+ {
+ gmx_fatal(FARGS, "Trying to remove more lipid molecules than there are in the membrane");
+ }
+ }
+ sfree(dist);
+ sfree(order);
+ }
+
+ rm_p->nr = nrm;
+ srenew(rm_p->mol, nrm);
+ srenew(rm_p->block, nrm);
+
+ return nupper+nlower;
+}
+
+/*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/
+static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state,
+ t_block *ins_at, pos_ins_t *pos_ins)
+{
+ int j, k, n, rm, mol_id, at, block;
+ rvec *x_tmp, *v_tmp;
+ int *list;
+ unsigned char *new_egrp[egcNR];
+ gmx_bool bRM;
+ int RMmolblock;
+
+ /* Construct the molecule range information */
+ gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop);
+
+ snew(list, state->natoms);
+ n = 0;
+ for (int i = 0; i < rm_p->nr; i++)
+ {
+ mol_id = rm_p->mol[i];
+ at = molecules.block(mol_id).size();
+ block = rm_p->block[i];
+ mtop->molblock[block].nmol--;
+ for (j = 0; j < mtop->moltype[mtop->molblock[block].type].atoms.nr; j++)
+ {
+ list[n] = at+j;
+ n++;
+ }
+ }
+
+ mtop->natoms -= n;
++ state_change_natoms(state, state->natoms - n);
+ snew(x_tmp, state->natoms);
+ snew(v_tmp, state->natoms);
+
+ for (int i = 0; i < egcNR; i++)
+ {
+ if (groups->grpnr[i] != nullptr)
+ {
+ groups->ngrpnr[i] = state->natoms;
+ snew(new_egrp[i], state->natoms);
+ }
+ }
+
+ auto x = makeArrayRef(state->x);
+ auto v = makeArrayRef(state->v);
+ rm = 0;
+ for (int i = 0; i < state->natoms+n; i++)
+ {
+ bRM = FALSE;
+ for (j = 0; j < n; j++)
+ {
+ if (i == list[j])
+ {
+ bRM = TRUE;
+ rm++;
+ }
+ }
+
+ if (!bRM)
+ {
+ for (j = 0; j < egcNR; j++)
+ {
+ if (groups->grpnr[j] != nullptr)
+ {
+ new_egrp[j][i-rm] = groups->grpnr[j][i];
+ }
+ }
+ copy_rvec(x[i], x_tmp[i-rm]);
+ copy_rvec(v[i], v_tmp[i-rm]);
+ for (j = 0; j < ins_at->nr; j++)
+ {
+ if (i == ins_at->index[j])
+ {
+ ins_at->index[j] = i-rm;
+ }
+ }
+
+ for (j = 0; j < pos_ins->pieces; j++)
+ {
+ for (k = 0; k < pos_ins->nidx[j]; k++)
+ {
+ if (i == pos_ins->subindex[j][k])
+ {
+ pos_ins->subindex[j][k] = i-rm;
+ }
+ }
+ }
+ }
+ }
+ for (int i = 0; i < state->natoms; i++)
+ {
+ copy_rvec(x_tmp[i], x[i]);
+ }
+ sfree(x_tmp);
+ for (int i = 0; i < state->natoms; i++)
+ {
+ copy_rvec(v_tmp[i], v[i]);
+ }
+ sfree(v_tmp);
+
+ for (int i = 0; i < egcNR; i++)
+ {
+ if (groups->grpnr[i] != nullptr)
+ {
+ sfree(groups->grpnr[i]);
+ groups->grpnr[i] = new_egrp[i];
+ }
+ }
+
+ /* remove empty molblocks */
+ RMmolblock = 0;
+ for (size_t i = 0; i < mtop->molblock.size(); i++)
+ {
+ if (mtop->molblock[i].nmol == 0)
+ {
+ RMmolblock++;
+ }
+ else
+ {
+ mtop->molblock[i-RMmolblock] = mtop->molblock[i];
+ }
+ }
+ mtop->molblock.resize(mtop->molblock.size() - RMmolblock);
+}
+
+/* remove al bonded interactions from mtop for the molecule to be embedded */
+static int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
+{
+ int j, m;
+ int type, natom, nmol, at, atom1 = 0, rm_at = 0;
+ gmx_bool *bRM, bINS;
+ /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
+ /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make *
+ * sure that g_membed exits with a warning when there are molecules of the same type not in the *
+ * ins_at index group. MGWolf 050710 */
+
+
+ snew(bRM, mtop->moltype.size());
+ for (size_t i = 0; i < mtop->moltype.size(); i++)
+ {
+ bRM[i] = TRUE;
+ }
+
+ for (size_t i = 0; i < mtop->molblock.size(); i++)
+ {
+ /*loop over molecule blocks*/
+ type = mtop->molblock[i].type;
+ natom = mtop->moltype[type].atoms.nr;
+ nmol = mtop->molblock[i].nmol;
+
+ for (j = 0; j < natom*nmol && bRM[type]; j++)
+ {
+ /*loop over atoms in the block*/
+ at = j+atom1; /*atom index = block index + offset*/
+ bINS = FALSE;
+
+ for (m = 0; (m < ins_at->nr) && (!bINS); m++)
+ {
+ /*loop over atoms in insertion index group to determine if we're inserting one*/
+ if (at == ins_at->index[m])
+ {
+ bINS = TRUE;
+ }
+ }
+ bRM[type] = bINS;
+ }
+ atom1 += natom*nmol; /*update offset*/
+ if (bRM[type])
+ {
+ rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/
+ }
+ }
+
+ for (size_t i = 0; i < mtop->moltype.size(); i++)
+ {
+ if (bRM[i])
+ {
+ for (j = 0; j < F_LJ; j++)
+ {
+ mtop->moltype[i].ilist[j].iatoms.clear();
+ }
+
+ for (j = F_POSRES; j <= F_VSITEN; j++)
+ {
+ mtop->moltype[i].ilist[j].iatoms.clear();
+ }
+ }
+ }
+ sfree(bRM);
+
+ return rm_at;
+}
+
+/* Write a topology where the number of molecules is correct for the system after embedding */
+static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop)
+{
+ int bMolecules = 0;
+ FILE *fpin, *fpout;
+ char buf[STRLEN], buf2[STRLEN], *temp;
+ int i, *nmol_rm, nmol, line;
+ char temporary_filename[STRLEN];
+
+ fpin = gmx_ffopen(topfile, "r");
+ strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
+ gmx_tmpnam(temporary_filename);
+ fpout = gmx_ffopen(temporary_filename, "w");
+
+ snew(nmol_rm, mtop->moltype.size());
+ for (i = 0; i < rm_p->nr; i++)
+ {
+ nmol_rm[rm_p->block[i]]++;
+ }
+
+ line = 0;
+ while (fgets(buf, STRLEN, fpin))
+ {
+ line++;
+ if (buf[0] != ';')
+ {
+ strcpy(buf2, buf);
+ if ((temp = strchr(buf2, '\n')) != nullptr)
+ {
+ temp[0] = '\0';
+ }
+ ltrim(buf2);
+ if (buf2[0] == '[')
+ {
+ buf2[0] = ' ';
+ if ((temp = strchr(buf2, '\n')) != nullptr)
+ {
+ temp[0] = '\0';
+ }
+ rtrim(buf2);
+ if (buf2[strlen(buf2)-1] == ']')
+ {
+ buf2[strlen(buf2)-1] = '\0';
+ ltrim(buf2);
+ rtrim(buf2);
+ if (gmx_strcasecmp(buf2, "molecules") == 0)
+ {
+ bMolecules = 1;
+ }
+ }
+ fprintf(fpout, "%s", buf);
+ }
+ else if (bMolecules == 1)
+ {
+ for (size_t i = 0; i < mtop->molblock.size(); i++)
+ {
+ nmol = mtop->molblock[i].nmol;
+ sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol);
+ fprintf(fpout, "%s", buf);
+ }
+ bMolecules = 2;
+ }
+ else if (bMolecules == 2)
+ {
+ /* print nothing */
+ }
+ else
+ {
+ fprintf(fpout, "%s", buf);
+ }
+ }
+ else
+ {
+ fprintf(fpout, "%s", buf);
+ }
+ }
+
+ gmx_ffclose(fpout);
+ /* use gmx_ffopen to generate backup of topinout */
+ fpout = gmx_ffopen(topfile, "w");
+ gmx_ffclose(fpout);
+ rename(temporary_filename, topfile);
+}
+
+void rescale_membed(int step_rel, gmx_membed_t *membed, rvec *x)
+{
+ /* Set new positions for the group to embed */
+ if (step_rel <= membed->it_xy)
+ {
+ membed->fac[0] += membed->xy_step;
+ membed->fac[1] += membed->xy_step;
+ }
+ else if (step_rel <= (membed->it_xy+membed->it_z))
+ {
+ membed->fac[2] += membed->z_step;
+ }
+ resize(membed->r_ins, x, membed->pos_ins, membed->fac);
+}
+
+/* We would like gn to be const as well, but C doesn't allow this */
+/* TODO this is utility functionality (search for the index of a
+ string in a collection), so should be refactored and located more
+ centrally. Also, it nearly duplicates the same string in readir.c */
+static int search_string(const char *s, int ng, char *gn[])
+{
+ int i;
+
+ for (i = 0; (i < ng); i++)
+ {
+ if (gmx_strcasecmp(s, gn[i]) == 0)
+ {
+ return i;
+ }
+ }
+
+ gmx_fatal(FARGS,
+ "Group %s selected for embedding was not found in the index file.\n"
+ "Group names must match either [moleculetype] names or custom index group\n"
+ "names, in which case you must supply an index file to the '-n' option\n"
+ "of grompp.",
+ s);
+}
+
+gmx_membed_t *init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
+ t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt)
+{
+ char *ins, **gnames;
+ int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0;
+ int ng, j, max_lip_rm, ins_grp_id, ntype, lip_rm;
+ real prot_area;
+ rvec *r_ins = nullptr;
+ t_block *ins_at, *rest_at;
+ pos_ins_t *pos_ins;
+ mem_t *mem_p;
+ rm_t *rm_p;
+ gmx_groups_t *groups;
+ gmx_bool bExcl = FALSE;
+ t_atoms atoms;
+ t_pbc *pbc;
+ char **piecename = nullptr;
+ gmx_membed_t *membed = nullptr;
+
+ /* input variables */
+ real xy_fac = 0.5;
+ real xy_max = 1.0;
+ real z_fac = 1.0;
+ real z_max = 1.0;
+ int it_xy = 1000;
+ int it_z = 0;
+ real probe_rad = 0.22;
+ int low_up_rm = 0;
+ int maxwarn = 0;
+ int pieces = 1;
+ gmx_bool bALLOW_ASYMMETRY = FALSE;
+
+ /* sanity check constants */ /* Issue a warning when: */
+ const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule *
+ * and rest smaller than this value is probably too small */
+ const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */
+ const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */
+ const int min_it_z = 100; /* the number of steps to embed in z is smaller */
+ const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */
+ const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */
+
+ snew(membed, 1);
+ snew(ins_at, 1);
+ snew(pos_ins, 1);
+
+ if (MASTER(cr))
+ {
+ fprintf(fplog, "Note: it is expected that in future gmx mdrun -membed will not be the "
+ "way to access this feature, perhaps changing to e.g. gmx membed.");
+ /* get input data out membed file */
+ try
+ {
+ get_input(opt2fn("-membed", nfile, fnm),
+ &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm,
+ &maxwarn, &pieces, &bALLOW_ASYMMETRY);
+ }
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+
+ if (!EI_DYNAMICS(inputrec->eI) )
+ {
+ gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
+ }
+
+ if (PAR(cr))
+ {
+ gmx_input("Sorry, parallel membed is not yet fully functional.");
+ }
+
+ if (*cpt >= 0)
+ {
+ fprintf(stderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
+ *cpt = -1;
+ }
+ groups = &(mtop->groups);
+ snew(gnames, groups->ngrpname);
+ for (i = 0; i < groups->ngrpname; i++)
+ {
+ gnames[i] = *(groups->grpname[i]);
+ }
+
+ atoms = gmx_mtop_global_atoms(mtop);
+ snew(mem_p, 1);
+ fprintf(stderr, "\nSelect a group to embed in the membrane:\n");
+ get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins);
+ ins_grp_id = search_string(ins, groups->ngrpname, gnames);
+ fprintf(stderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins);
+ get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name));
+
+ pos_ins->pieces = pieces;
+ snew(pos_ins->nidx, pieces);
+ snew(pos_ins->subindex, pieces);
+ snew(piecename, pieces);
+ if (pieces > 1)
+ {
+ fprintf(stderr, "\nSelect pieces to embed:\n");
+ get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename);
+ }
+ else
+ {
+ /*use whole embedded group*/
+ snew(pos_ins->nidx, 1);
+ snew(pos_ins->subindex, 1);
+ pos_ins->nidx[0] = ins_at->nr;
+ pos_ins->subindex[0] = ins_at->index;
+ }
+
+ if (probe_rad < min_probe_rad)
+ {
+ warn++;
+ fprintf(stderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result "
+ "in overlap between waters and the group to embed, which will result "
+ "in Lincs errors etc.\n\n", warn);
+ }
+
+ if (xy_fac < min_xy_init)
+ {
+ warn++;
+ fprintf(stderr, "\nWarning %d:\nThe initial size of %s is probably too small.\n\n", warn, ins);
+ }
+
+ if (it_xy < min_it_xy)
+ {
+ warn++;
+ fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)"
+ " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy);
+ }
+
+ if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) )
+ {
+ warn++;
+ fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)"
+ " is probably too small.\nIncrease -nz or the maxwarn setting in the membed input file.\n\n", warn, ins, it_z);
+ }
+
+ if (it_xy+it_z > inputrec->nsteps)
+ {
+ warn++;
+ fprintf(stderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the "
+ "number of steps in the tpr.\n"
+ "(increase maxwarn in the membed input file to override)\n\n", warn);
+ }
+
+ fr_id = -1;
+ if (inputrec->opts.ngfrz == 1)
+ {
+ gmx_fatal(FARGS, "You did not specify \"%s\" as a freezegroup.", ins);
+ }
+
+ for (i = 0; i < inputrec->opts.ngfrz; i++)
+ {
+ tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
+ if (ins_grp_id == tmp_id)
+ {
+ fr_id = tmp_id;
+ fr_i = i;
+ }
+ }
+
+ if (fr_id == -1)
+ {
+ gmx_fatal(FARGS, "\"%s\" not as freezegroup defined in the mdp-file.", ins);
+ }
+
+ for (i = 0; i < DIM; i++)
+ {
+ if (inputrec->opts.nFreeze[fr_i][i] != 1)
+ {
+ gmx_fatal(FARGS, "freeze dimensions for %s are not Y Y Y\n", ins);
+ }
+ }
+
+ ng = groups->grps[egcENER].nr;
+ if (ng == 1)
+ {
+ gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
+ }
+
+ for (i = 0; i < ng; i++)
+ {
+ for (j = 0; j < ng; j++)
+ {
+ if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
+ {
+ bExcl = TRUE;
+ if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) ||
+ (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
+ {
+ gmx_fatal(FARGS, "Energy exclusions \"%s\" and \"%s\" do not match the group "
+ "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]],
+ *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins);
+ }
+ }
+ }
+ }
+
+ if (!bExcl)
+ {
+ gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "
+ "the freeze group");
+ }
+
+ /* Obtain the maximum and minimum coordinates of the group to be embedded */
+ snew(rest_at, 1);
+ init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max);
+ /* Check that moleculetypes in insertion group are not part of the rest of the system */
+ check_types(ins_at, rest_at, mtop);
+
+ init_mem_at(mem_p, mtop, state->x.rvec_array(), state->box, pos_ins);
+
+ prot_area = est_prot_area(pos_ins, state->x.rvec_array(), ins_at, mem_p);
+ if ( (prot_area > prot_vs_box) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX]) < box_vs_prot) )
+ {
+ warn++;
+ fprintf(stderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
+ "This might cause pressure problems during the growth phase. Just try with\n"
+ "current setup and increase 'maxwarn' in your membed settings file, but lower the\n"
+ "compressibility in the mdp-file or disable pressure coupling if problems occur.\n\n", warn);
+ }
+
+ if (warn > maxwarn)
+ {
+ gmx_fatal(FARGS, "Too many warnings (override by setting maxwarn in the membed input file)\n");
+ }
+
+ printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area);
+ printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n"
+ "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area);
+
+ /* Maximum number of lipids to be removed*/
+ max_lip_rm = static_cast<int>(2*prot_area/mem_p->lip_area);
+ printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm);
+
+ printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
+ "This resizing will be done with respect to the geometrical center of all protein atoms\n"
+ "that span the membrane region, i.e. z between %.3f and %.3f\n\n",
+ xy_fac, z_fac, mem_p->zmin, mem_p->zmax);
+
+ /* resize the protein by xy and by z if necessary*/
+ snew(r_ins, ins_at->nr);
+ init_resize(ins_at, r_ins, pos_ins, mem_p, state->x.rvec_array(), bALLOW_ASYMMETRY);
+ membed->fac[0] = membed->fac[1] = xy_fac;
+ membed->fac[2] = z_fac;
+
+ membed->xy_step = (xy_max-xy_fac)/static_cast<double>(it_xy);
+ membed->z_step = (z_max-z_fac)/static_cast<double>(it_z-1);
+
+ resize(r_ins, state->x.rvec_array(), pos_ins, membed->fac);
+
+ /* remove overlapping lipids and water from the membrane box*/
+ /*mark molecules to be removed*/
+ snew(pbc, 1);
+ set_pbc(pbc, inputrec->ePBC, state->box);
+
+ snew(rm_p, 1);
+ lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x.rvec_array(), mem_p, pos_ins,
+ probe_rad, low_up_rm, bALLOW_ASYMMETRY);
+ lip_rm -= low_up_rm;
+
+ if (fplog)
+ {
+ for (i = 0; i < rm_p->nr; i++)
+ {
+ fprintf(fplog, "rm mol %d\n", rm_p->mol[i]);
+ }
+ }
+
+ for (size_t i = 0; i < mtop->molblock.size(); i++)
+ {
+ ntype = 0;
+ for (j = 0; j < rm_p->nr; j++)
+ {
+ if (rm_p->block[j] == static_cast<int>(i))
+ {
+ ntype++;
+ }
+ }
+ printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name));
+ }
+
+ if (lip_rm > max_lip_rm)
+ {
+ warn++;
+ fprintf(stderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated "
+ "protein area\nTry making the -xyinit resize factor smaller or increase "
+ "maxwarn in the membed input file.\n\n", warn);
+ }
+
+ /*remove all lipids and waters overlapping and update all important structures*/
+ rm_group(groups, mtop, rm_p, state, ins_at, pos_ins);
+
+ rm_bonded_at = rm_bonded(ins_at, mtop);
+ if (rm_bonded_at != ins_at->nr)
+ {
+ fprintf(stderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, "
+ "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
+ "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr);
+ }
+
+ if (warn > maxwarn)
+ {
+ gmx_fatal(FARGS, "Too many warnings.\nIf you are sure these warnings are harmless,\n"
+ "you can increase the maxwarn setting in the membed input file.");
+ }
+
+ // Re-establish the invariants of the derived values within
+ // mtop.
+ gmx_mtop_finalize(mtop);
+
+ if (ftp2bSet(efTOP, nfile, fnm))
+ {
+ top_update(opt2fn("-mp", nfile, fnm), rm_p, mtop);
+ }
+
+ sfree(pbc);
+ sfree(rest_at);
+ if (pieces > 1)
+ {
+ sfree(piecename);
+ }
+
+ membed->it_xy = it_xy;
+ membed->it_z = it_z;
+ membed->pos_ins = pos_ins;
+ membed->r_ins = r_ins;
+ }
+
+ return membed;
+}
+
+void free_membed(gmx_membed_t *membed)
+{
+ sfree(membed);
+}
--- /dev/null
- double *enerpart_lambda; /* Partial energy for lambda and flambda[] */
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, 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.
+ *
+ * GROMACS 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.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#ifndef GMX_MDTYPES_TYPES_ENERDATA_H
+#define GMX_MDTYPES_TYPES_ENERDATA_H
+
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/utility/real.h"
+
+enum {
+ egCOULSR, egLJSR, egBHAMSR,
+ egCOUL14, egLJ14, egNR
+};
+
+struct gmx_grppairener_t
+{
+ int nener; /* The number of energy group pairs */
+ real *ener[egNR]; /* Energy terms for each pair of groups */
+};
+
+struct gmx_enerdata_t
+{
+ real term[F_NRE]; /* The energies for all different interaction types */
+ struct gmx_grppairener_t grpp;
+ double dvdl_lin[efptNR]; /* Contributions to dvdl with linear lam-dependence */
+ double dvdl_nonlin[efptNR]; /* Idem, but non-linear dependence */
++ /* The idea is that dvdl terms with linear lambda dependence will be added
++ * automatically to enerpart_lambda. Terms with non-linear lambda dependence
++ * should explicitly determine the energies at foreign lambda points
++ * when n_lambda > 0. */
++
+ int n_lambda;
+ int fep_state; /*current fep state -- just for printing */
++ double *enerpart_lambda; /* Partial Hamiltonian for lambda and flambda[], includes at least all perturbed terms */
+ real foreign_term[F_NRE]; /* alternate array for storing foreign lambda energies */
+ struct gmx_grppairener_t foreign_grpp; /* alternate array for storing foreign lambda energies */
+};
+
+#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff