#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2014,2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
+# Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
- # Copyright (c) 2019,2020, by the GROMACS development team, led by
++# Copyright (c) 2019,2020,2021, 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.
# of configure time, because 1) some of the version variables are only
# available during build time, and 2) we don't want to do all the Sphinx setup
# during configuration to save some time when not building the content.
-# All the generated values get put into conf-vars.py (generated from
-# conf-vars.py.cmakein), which in turn is included by the Sphinx configuration
-# file conf.py.
+# All of the generated values get put into conf.py (generated from
+# conf.cmakein.py).
set(SOURCE_MD5SUM "unknown" CACHE STRING
"MD5 sum of the source tarball, normally used only for the pre-release webpage build")
set(SPHINX_SOURCE_FILES
index.rst
download.rst
- conf.py
links.dat
dev-manual/build-system.rst
dev-manual/change-management.rst
dev-manual/commitstyle.rst
+ dev-manual/containers.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/gitlab.rst
dev-manual/gmxtree.rst
dev-manual/includestyle.rst
dev-manual/index.rst
+ dev-manual/infrastructure.rst
dev-manual/jenkins.rst
dev-manual/known-issues.rst
dev-manual/language-features.rst
how-to/visualize.rst
install-guide/index.rst
release-notes/index.rst
+ release-notes/2021/major/highlights.rst
+ release-notes/2021/major/features.rst
+ release-notes/2021/major/performance.rst
+ release-notes/2021/major/tools.rst
+ release-notes/2021/major/bugs-fixed.rst
+ release-notes/2021/major/removed-functionality.rst
+ release-notes/2021/major/deprecated-functionality.rst
+ release-notes/2021/major/portability.rst
+ release-notes/2021/major/miscellaneous.rst
release-notes/2020/2020.1.rst
release-notes/2020/2020.2.rst
release-notes/2020/2020.3.rst
release-notes/2020/2020.4.rst
release-notes/2020/2020.5.rst
+ release-notes/2020/2020.6.rst
release-notes/2020/major/highlights.rst
release-notes/2020/major/features.rst
release-notes/2020/major/performance.rst
set(IMAGE_CONVERT_STRING "impossible")
endif()
- set(SPHINX_CONFIG_VARS_FILE ${SPHINX_INPUT_DIR}/conf-vars.py)
+ set(SPHINX_CONFIG_FILE ${SPHINX_INPUT_DIR}/conf.py)
if (GMX_PYTHON_PACKAGE)
set(GMXAPI_PYTHON_STAGING_DIR ${CMAKE_BINARY_DIR}/python_packaging/src/gmxapi_staging)
# TODO: Resolve circular reference. We would like to get the CMake build-time directory for
# in this context?
endif ()
- gmx_configure_version_file(conf-vars.py.cmakein ${SPHINX_CONFIG_VARS_FILE}
+ gmx_configure_version_file(
+ conf.cmakein.py ${SPHINX_CONFIG_FILE}
EXTRA_VARS
- SPHINX_EXTENSION_PATH RELENG_PATH
- IMAGE_CONVERT_STRING
+ CMAKE_MINIMUM_REQUIRED_VERSION
EXPECTED_DOXYGEN_VERSION
EXPECTED_SPHINX_VERSION
- CMAKE_MINIMUM_REQUIRED_VERSION REQUIRED_CUDA_VERSION
- REQUIRED_OPENCL_MIN_VERSION
- REQUIRED_CUDA_COMPUTE_CAPABILITY REGRESSIONTEST_VERSION
- SOURCE_MD5SUM REGRESSIONTEST_MD5SUM_STRING
- GMX_TNG_MINIMUM_REQUIRED_VERSION
+ GMX_ADMIN_DIR
GMX_LMFIT_REQUIRED_VERSION
GMX_MANUAL_DOI_STRING
+ GMX_TNG_MINIMUM_REQUIRED_VERSION
GMX_SOURCE_DOI_STRING
GMXAPI_PYTHON_STAGING_DIR
+ IMAGE_CONVERT_STRING
+ REGRESSIONTEST_VERSION
+ REQUIRED_CUDA_COMPUTE_CAPABILITY
+ REQUIRED_CUDA_VERSION
+ REQUIRED_OPENCL_MIN_VERSION
+ REGRESSIONTEST_MD5SUM_STRING
+ RELENG_PATH
+ SOURCE_MD5SUM
+ SPHINX_EXTENSION_PATH
COMMENT "Configuring Sphinx configuration file")
- gmx_add_sphinx_input_file(${SPHINX_CONFIG_VARS_FILE})
+ gmx_add_sphinx_input_file(${SPHINX_CONFIG_FILE})
gmx_add_sphinx_source_files(FILES ${SPHINX_SOURCE_FILES})
if (EXISTS ${RELENG_PATH}/docs/FileList.cmake)
include(${RELENG_PATH}/docs/FileList.cmake)
)
endif ()
+ gmx_add_sphinx_source_files(
+ FILES
+ nblib/index.rst
+ nblib/guide-to-writing-MD-programs.rst
+ )
+
gmx_add_sphinx_source_files(
FILES
${REFERENCEMANUAL_SPHINX_FILES_GENERAL})
# Sphinx cache with pickled ReST documents
set(SPHINX_CACHE_DIR "${CMAKE_CURRENT_BINARY_DIR}/_doctrees")
+ set(SPHINX_CONFIG_OVERRIDES "")
+ if (GMX_DEVELOPER_BUILD)
+ set(SPHINX_CONFIG_OVERRIDES "-Dtodo_include_todos=1")
+ endif()
add_custom_target(webpage-sphinx
DEPENDS sphinx-programs
DEPENDS sphinx-input
-q -b html
-w sphinx-html.log
-d "${SPHINX_CACHE_DIR}"
+ ${SPHINX_CONFIG_OVERRIDES}
"${SPHINX_INPUT_DIR}"
"${HTML_OUTPUT_DIR}"
WORKING_DIRECTORY
set(HTML_BUILD_WARNINGS)
# Next, turn it off if any of the preconditions are unsatisified
-if (NOT PythonInterp_FOUND)
+if (NOT Python3_Interpreter_FOUND)
set(HTML_BUILD_IS_POSSIBLE OFF)
set(HTML_BUILD_NOT_POSSIBLE_REASON "Python is required")
elseif (NOT SPHINX_FOUND)
--- /dev/null
- # Copyright (c) 2020, by the GROMACS development team, led by
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2015,2016,2017,2018,2019 by the GROMACS development team.
- .. |gmx-source-package-ftp| replace:: As ftp ftp://ftp.gromacs.org/pub/gromacs/gromacs-{gmx_version_string}.tar.gz
- .. |gmx-source-package-http| replace:: As http http://ftp.gromacs.org/pub/gromacs/gromacs-{gmx_version_string}.tar.gz
- .. |gmx-regressiontests-package| replace:: http://ftp.gromacs.org/pub/regressiontests/regressiontests-{regressiontest_version}.tar.gz
- .. _up-to-date installation instructions: http://manual.gromacs.org/documentation/current/install-guide/index.html
++# Copyright (c) 2020,2021, 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.
+
+# -*- coding: utf-8 -*-
+#
+# GROMACS documentation build configuration file, created by
+# sphinx-quickstart on Tue Jan 13 14:28:44 2015.
+#
+# This file is execfile()d with the current directory set to its
+# containing dir.
+#
+# Note that not all possible configuration values are present in this
+# autogenerated file.
+#
+# All configuration values have a default; values that are commented out
+# serve to show the default.
+
+import datetime
+import os
+import sys
+
+# The following definitions are completed via CMake machinery.
+gmx_containers_path = '@GMX_ADMIN_DIR@/containers'
+gmx_sphinx_extension_path = '@SPHINX_EXTENSION_PATH@'
+gmxapi_staging_path = '@GMXAPI_PYTHON_STAGING_DIR@'
+releng_path = '@RELENG_PATH@'
+gmx_version_string = '@GMX_VERSION_STRING@'
+gmx_version_string_full = '@GMX_VERSION_STRING_FULL@'
+regressiontest_version = '@REGRESSIONTEST_VERSION@'
+gmx_min_sphinx = '@EXPECTED_SPHINX_VERSION@'
+gmx_image_convert = '@IMAGE_CONVERT_STRING@'
+variables = [
+ ('EXPECTED_DOXYGEN_VERSION', '@EXPECTED_DOXYGEN_VERSION@'),
+ ('EXPECTED_SPHINX_VERSION', '@EXPECTED_SPHINX_VERSION@'),
+ ('CMAKE_MINIMUM_REQUIRED_VERSION', '@CMAKE_MINIMUM_REQUIRED_VERSION@'),
+ ('REQUIRED_CUDA_VERSION', '@REQUIRED_CUDA_VERSION@'),
+ ('REQUIRED_CUDA_COMPUTE_CAPABILITY', '@REQUIRED_CUDA_COMPUTE_CAPABILITY@'),
+ ('REQUIRED_OPENCL_MIN_VERSION', '@REQUIRED_OPENCL_MIN_VERSION@'),
+ ('SOURCE_MD5SUM', '@SOURCE_MD5SUM@'),
+ ('REGRESSIONTEST_MD5SUM', '@REGRESSIONTEST_MD5SUM_STRING@'),
+ ('GMX_TNG_MINIMUM_REQUIRED_VERSION', '@GMX_TNG_MINIMUM_REQUIRED_VERSION@'),
+ ('GMX_LMFIT_REQUIRED_VERSION', '@GMX_LMFIT_REQUIRED_VERSION@'),
+ ('GMX_MANUAL_DOI_STRING', '@GMX_MANUAL_DOI_STRING@'),
+ ('GMX_SOURCE_DOI_STRING', '@GMX_SOURCE_DOI_STRING@')
+]
+# End of build-time substitutions.
+
+sys.path.append(gmx_containers_path)
+sys.path.append(gmx_sphinx_extension_path)
+if releng_path and os.path.isdir(releng_path):
+ sys.path.append(releng_path)
+if gmxapi_staging_path and os.path.isdir(gmxapi_staging_path):
+ sys.path.append(gmxapi_staging_path)
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#sys.path.insert(0, os.path.abspath('.'))
+
+# -- General configuration ------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+# gmx_min_sphinx is set from the expected minimum version of Sphinx
+# in CMakeLists.txt
+needs_sphinx = gmx_min_sphinx
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+ 'sphinx.ext.autodoc',
+ 'sphinx.ext.graphviz',
+ 'sphinx.ext.extlinks',
+ 'sphinx.ext.napoleon',
+ 'sphinx.ext.doctest',
+ 'sphinx.ext.ifconfig',
+ 'sphinx.ext.intersphinx',
+ 'sphinx.ext.mathjax',
+ 'sphinx.ext.napoleon',
+ 'sphinx.ext.todo',
+ # The viewcode extension requires extra configuration or consideration to
+ # avoid collecting sources for external packages (through intersphninx) or
+ # generating broken links for compiled extension modules that do not have
+ # inspectable 'code' member data.
+ # Ref: https://www.sphinx-doc.org/en/master/usage/extensions/viewcode.html
+ # Please run linkchecker and inspect the generated
+ # docs/html/_modules/index.html page before committing a change that enables
+ # 'sphinx.ext.viewcode',
+ 'gmxsphinx'
+]
+extlinks = {'issue': ('https://gitlab.com/gromacs/gromacs/-/issues/%s',
+ 'Issue ')}
+
+# Add any paths that contain templates here, relative to this directory.
+# templates_path = ['_templates']
+
+# The suffix of source filenames.
+source_suffix = '.rst'
+
+# The encoding of source files.
+#source_encoding = 'utf-8-sig'
+
+# The master toctree document.
+master_doc = 'index'
+
+# General information about the project.
+project = u'GROMACS'
+copyright = str(datetime.datetime.now().year) + u', GROMACS development team'
+thisyear_string = str(datetime.datetime.now().year)
+# The version info for the project you're documenting, acts as replacement for
+# |version| and |release|, also used in various other places throughout the
+# built documents.
+#
+# The short X.Y version.
+version = gmx_version_string
+# The full version, including alpha/beta/rc tags.
+release = gmx_version_string_full
+# default file extension for plots
+plotext = u'.eps'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#language = None
+
+# There are two options for replacing |today|: either, you set today to some
+# non-false value, then it is used:
+#today = ''
+# Else, today_fmt is used as the format for a strftime call.
+#today_fmt = '%B %d, %Y'
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+exclude_patterns = ['fragments']
+if not tags.has('do_man'):
+ exclude_patterns += ['man']
+
+# Set variable if documentation build can convert images or not
+# to selectively include files in reference manual
+def setup(app):
+ app.add_config_value('gmx_image_convert', 'impossible', 'env')
+
+# The reST default role (used for this markup: `text`) to use for all
+# documents.
+default_role = 'any'
+
+# If true, '()' will be appended to :func: etc. cross-reference text.
+#add_function_parentheses = True
+
+# If true, the current module name will be prepended to all description
+# unit titles (such as .. function::).
+#add_module_names = True
+
+# If true, sectionauthor and moduleauthor directives will be shown in the
+# output. They are ignored by default.
+#show_authors = False
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+# A list of ignored prefixes for module index sorting.
+#modindex_common_prefix = []
+
+# If true, keep warnings as "system message" paragraphs in the built documents.
+#keep_warnings = False
+
+# Configure the values for all the variables that might later configure any of the .rst files.
+substitutions = ['.. |{0}| replace:: {1}'.format(*x) for x in variables if x[1]]
+substitutions.extend(['.. |{0}| replace:: unknown'.format(x[0]) for x in variables if x[1] == ''])
+rst_epilog = "\n".join(substitutions)
+rst_epilog += """
+.. |Gromacs| replace:: GROMACS
+.. _gmx-manual: manual-{gmx_version_string}.pdf
+.. _gmx-manual-parent-dir: ../manual-{gmx_version_string}.pdf
- .. _ftp site: ftp://ftp.gromacs.org/pub/gromacs/
- .. _tutorials: http://www.gromacs.org/Documentation/Tutorials
++.. |gmx-source-package-ftp| replace:: As ftp ftp://ftp.gromacs.org/gromacs/gromacs-{gmx_version_string}.tar.gz
++.. |gmx-source-package-http| replace:: As https https://ftp.gromacs.org/gromacs/gromacs-{gmx_version_string}.tar.gz
++.. |gmx-regressiontests-package| replace:: https://ftp.gromacs.org/regressiontests/regressiontests-{regressiontest_version}.tar.gz
++.. _up-to-date installation instructions: https://manual.gromacs.org/documentation/current/install-guide/index.html
+.. _CUDA: http://www.nvidia.com/object/cuda_home_new.html
+.. _OpenCL: https://www.khronos.org/opencl/
+.. _OpenMPI: http://www.open-mpi.org
+.. _MPICH: http://www.mpich.org
+.. _LAM-MPI: http://www.lam-mpi.org
+.. _OpenMP: http://en.wikipedia.org/wiki/OpenMP
+.. _CMake installation page: http://www.cmake.org/install/
+.. _Ubuntu toolchain ppa page: https://launchpad.net/~ubuntu-toolchain-r/+archive/ubuntu/test
+.. _EPEL page: https://fedoraproject.org/wiki/EPEL
+.. _running CMake: http://www.cmake.org/runningcmake/
+.. _CMake environment variables: http://cmake.org/Wiki/CMake_Useful_Variables#Environment_Variables
+.. _FFTW: http://www.fftw.org
+.. _FFTW installation guide: http://www.fftw.org/doc/Installation-and-Customization.html#Installation-and-Customization
+.. _MKL: https://software.intel.com/en-us/intel-mkl
+.. _VMD: http://www.ks.uiuc.edu/Research/vmd/
+.. _PyMOL: http://www.pymol.org
+.. _webpage: http://www.gromacs.org
++.. _ftp site: ftp://ftp.gromacs.org/gromacs/
++.. _tutorials: http://www.mdtutorials.com/gmx/
+.. _issue tracker: https://gitlab.com/gromacs/gromacs/-/issues/
+.. _gitlab: https://gitlab.com/gromacs/gromacs/
+.. _download: ../download.html
+.. |thisyear| replace:: {thisyear_string}
+""".format(gmx_version_string=gmx_version_string, regressiontest_version=regressiontest_version, thisyear_string=thisyear_string)
+
+# -- Options for HTML output ----------------------------------------------
+
+# The theme to use for HTML and HTML Help pages. See the documentation for
+# a list of builtin themes.
+html_theme = 'classic'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further. For a list of options available for each theme, see the
+# documentation.
+#html_theme_options = {}
+
+# Add any paths that contain custom themes here, relative to this directory.
+html_theme_path = ['']
+
+# The name for this set of Sphinx documents. If None, it defaults to
+# "<project> v<release> documentation".
+html_title = u'GROMACS ' + release + ' documentation'
+
+# A shorter title for the navigation bar. Default is the same as html_title.
+html_short_title = u'GROMACS ' + version
+
+# The name of an image file (relative to this directory) to place at the top
+# of the sidebar.
+#html_logo = None
+
+# The name of an image file (within the static path) to use as favicon of the
+# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
+# pixels large.
+#html_favicon = None
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# Add any extra paths that contain custom files (such as robots.txt or
+# .htaccess) here, relative to this directory. These files are copied
+# directly to the root of the documentation.
+#html_extra_path = []
+
+# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
+# using the given strftime format.
+#html_last_updated_fmt = '%b %d, %Y'
+
+# If true, SmartyPants will be used to convert quotes and dashes to
+# typographically correct entities.
+#html_use_smartypants = True
+
+# Custom sidebar templates, maps document names to template names.
+#html_sidebars = {}
+
+# Additional templates that should be rendered to pages, maps page names to
+# template names.
+#html_additional_pages = {}
+
+# If false, no module index is generated.
+#html_domain_indices = True
+
+# If false, no index is generated.
+#html_use_index = True
+
+# If true, the index is split into individual pages for each letter.
+#html_split_index = False
+
+# If true, links to the reST sources are added to the pages.
+#html_show_sourcelink = True
+
+# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
+#html_show_sphinx = True
+
+# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
+#html_show_copyright = True
+
+# If true, an OpenSearch description file will be output, and all pages will
+# contain a <link> tag referring to it. The value of this option must be the
+# base URL from which the finished HTML is served.
+#html_use_opensearch = ''
+
+# This is the file name suffix for HTML files (e.g. ".xhtml").
+#html_file_suffix = None
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'Gromacsdoc'
+
+
+# -- Options for LaTeX output ---------------------------------------------
+
+latex_elements = {
+# The paper size ('letterpaper' or 'a4paper').
+'papersize': 'a4paper',
+
+
+# The font size ('10pt', '11pt' or '12pt').
+'pointsize': '11',
+
+# Additional stuff for the LaTeX preamble.
+# The tocdepth setting is needed to overwrite the default value set by Sphinx
+# to get a more detailed toctree in the pdf version of the manual.
+ 'preamble': r'''
+ \usepackage{here}
+ \usepackage{picins}
+ \usepackage{underscore}
+ \usepackage{tabularx}
+ \usepackage{multicol}
+ \usepackage{dcolumn}
+ \usepackage{makeidx}
+ \usepackage{times}
+ \usepackage{ifthen}
+ \usepackage{enumitem}
+ \usepackage{longtable}
+ \usepackage{pdflscape}
+ \pagenumbering{roman}
+ \usepackage{array}
+ \setcounter{tocdepth}{2}
+ ''',
+# Format output to avoid empty pages
+ 'classoptions': ',openany,oneside'
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title,
+# author, documentclass [howto, manual, or own class]).
+latex_documents = [
+ ('index', 'gromacs.tex', u'GROMACS Documentation',
+ u'GROMACS development team', 'manual'),
+]
+
+# The name of an image file (relative to this directory) to place at the top of
+# the title page.
+latex_logo = 'reference-manual/plots/peregrine.png'
+
+# For "manual" documents, if this is true, then toplevel headings are parts,
+# not chapters.
+latex_use_parts = True
+
+# If true, show page references after internal links.
+latex_show_pagerefs = True
+
+# If true, show URL addresses after external links.
+#latex_show_urls = False
+
+# Documents to append as an appendix to all manuals.
+#latex_appendices = []
+
+# If false, no module index is generated.
+#latex_domain_indices = True
+
+
+# -- Options for manual page output ---------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+if tags.has('do_man'):
+ exec(open('conf-man.py').read())
+
+# If true, show URL addresses after external links.
+#man_show_urls = False
+
+
+# -- Options for Texinfo output -------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+# dir menu entry, description, category)
+texinfo_documents = [
+ ('index', 'GROMACS', u'GROMACS Documentation',
+ u'GROMACS development team', 'GROMACS', 'One line description of project.',
+ 'Miscellaneous'),
+]
+
+# Documents to append as an appendix to all manuals.
+#texinfo_appendices = []
+
+# If false, no module index is generated.
+#texinfo_domain_indices = True
+
+# How to display URL addresses: 'footnote', 'no', or 'inline'.
+#texinfo_show_urls = 'footnote'
+
+# If true, do not generate a @detailmenu in the "Top" node's menu.
+#texinfo_no_detailmenu = False
+
+# Make it possible to use numbered labels for figures and tables
+numfig = True
+
+# -- Options for autodoc extension ----------------------------------------
+
+# http://www.sphinx-doc.org/en/master/usage/extensions/autodoc.html#confval-autodoc-mock-imports
+autodoc_mock_imports = ['hpccm', 'hpccm.config']
+
+# -- Options for intersphinx extension ------------------------------------
+
+intersphinx_mapping = {'python': ('https://docs.python.org/3', None)}
+intersphinx_cache_limit = -1
+intersphinx_timeout = 10
functionality supported, whereas patch releases contain only fixes for
issues identified in the corresponding major releases.
-Two versions of |Gromacs| are under active maintenance, the 2020
-series and the 2019 series. In the latter, only highly conservative
+Two versions of |Gromacs| are under active maintenance, the 2021
+series and the 2020 series. In the latter, only highly conservative
fixes will be made, and only to address issues that affect scientific
correctness. Naturally, some of those releases will be made after the
-year 2019 ends, but we keep 2019 in the name so users understand how
+year 2020 ends, but we keep 2019 in the name so users understand how
up to date their version is. Such fixes will also be incorporated into
-the 2020 release series, as appropriate. Around the time the 2021
-release is made, the 2019 series will no longer be maintained.
+the 2021 release series, as appropriate. Around the time the 2022
+release is made, the 2020 series will no longer be maintained.
Where issue numbers are reported in these release notes, more details
-can be found at https://gitlab.com/gromacs/gromacs/issues at that issue number.
+can be found on the `issue tracker`_ at that issue number.
+
+|Gromacs| 2021 series
+---------------------
+
+.. todolist::
+
+Major release
+^^^^^^^^^^^^^
+
+.. toctree::
+ :maxdepth: 1
+
+ 2021/major/highlights
+ 2021/major/features
+ 2021/major/performance
+ 2021/major/tools
+ 2021/major/bugs-fixed
+ 2021/major/deprecated-functionality
+ 2021/major/removed-functionality
+ 2021/major/portability
+ 2021/major/miscellaneous
|Gromacs| 2020 series
---------------------
.. toctree::
:maxdepth: 1
+ 2020/2020.6
2020/2020.5
2020/2020.4
2020/2020.3
--- /dev/null
- * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
- buffer[i] = pointState[i].inTargetRegion() ? std::exp(-pointState[i].logPmfSum()) : 0;
++ * Copyright (c) 2015,2016,2017,2018,2019, The GROMACS development team.
++ * Copyright (c) 2020,2021, 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.
+ */
+
+/*! \internal \file
+ * \brief
+ * Implements the BiasState class.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_awh
+ */
+
+#include "gmxpre.h"
+
+#include "biasstate.h"
+
+#include <cassert>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
+#include <optional>
+
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/mdrunutility/multisim.h"
+#include "gromacs/mdtypes/awh_history.h"
+#include "gromacs/mdtypes/awh_params.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/simd/simd.h"
+#include "gromacs/simd/simd_math.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringutil.h"
+
+#include "biasgrid.h"
+#include "pointstate.h"
+
+namespace gmx
+{
+
+void BiasState::getPmf(gmx::ArrayRef<float> pmf) const
+{
+ GMX_ASSERT(pmf.size() == points_.size(), "pmf should have the size of the bias grid");
+
+ /* The PMF is just the negative of the log of the sampled PMF histogram.
+ * Points with zero target weight are ignored, they will mostly contain noise.
+ */
+ for (size_t i = 0; i < points_.size(); i++)
+ {
+ pmf[i] = points_[i].inTargetRegion() ? -points_[i].logPmfSum() : GMX_FLOAT_MAX;
+ }
+}
+
+namespace
+{
+
+/*! \brief
+ * Sum an array over all simulations on the master rank of each simulation.
+ *
+ * \param[in,out] arrayRef The data to sum.
+ * \param[in] multiSimComm Struct for multi-simulation communication.
+ */
+void sumOverSimulations(gmx::ArrayRef<int> arrayRef, const gmx_multisim_t* multiSimComm)
+{
+ gmx_sumi_sim(arrayRef.size(), arrayRef.data(), multiSimComm);
+}
+
+/*! \brief
+ * Sum an array over all simulations on the master rank of each simulation.
+ *
+ * \param[in,out] arrayRef The data to sum.
+ * \param[in] multiSimComm Struct for multi-simulation communication.
+ */
+void sumOverSimulations(gmx::ArrayRef<double> arrayRef, const gmx_multisim_t* multiSimComm)
+{
+ gmx_sumd_sim(arrayRef.size(), arrayRef.data(), multiSimComm);
+}
+
+/*! \brief
+ * Sum an array over all simulations on all ranks of each simulation.
+ *
+ * This assumes the data is identical on all ranks within each simulation.
+ *
+ * \param[in,out] arrayRef The data to sum.
+ * \param[in] commRecord Struct for intra-simulation communication.
+ * \param[in] multiSimComm Struct for multi-simulation communication.
+ */
+template<typename T>
+void sumOverSimulations(gmx::ArrayRef<T> arrayRef, const t_commrec* commRecord, const gmx_multisim_t* multiSimComm)
+{
+ if (MASTER(commRecord))
+ {
+ sumOverSimulations(arrayRef, multiSimComm);
+ }
+ if (commRecord->nnodes > 1)
+ {
+ gmx_bcast(arrayRef.size() * sizeof(T), arrayRef.data(), commRecord->mpi_comm_mygroup);
+ }
+}
+
+/*! \brief
+ * Sum PMF over multiple simulations, when requested.
+ *
+ * \param[in,out] pointState The state of the points in the bias.
+ * \param[in] numSharedUpdate The number of biases sharing the histogram.
+ * \param[in] commRecord Struct for intra-simulation communication.
+ * \param[in] multiSimComm Struct for multi-simulation communication.
+ */
+void sumPmf(gmx::ArrayRef<PointState> pointState,
+ int numSharedUpdate,
+ const t_commrec* commRecord,
+ const gmx_multisim_t* multiSimComm)
+{
+ if (numSharedUpdate == 1)
+ {
+ return;
+ }
+ GMX_ASSERT(multiSimComm != nullptr && numSharedUpdate % multiSimComm->numSimulations_ == 0,
+ "numSharedUpdate should be a multiple of multiSimComm->numSimulations_");
+ GMX_ASSERT(numSharedUpdate == multiSimComm->numSimulations_,
+ "Sharing within a simulation is not implemented (yet)");
+
+ std::vector<double> buffer(pointState.size());
+
+ /* Need to temporarily exponentiate the log weights to sum over simulations */
+ for (size_t i = 0; i < buffer.size(); i++)
+ {
- pointState[i].setLogPmfSum(-std::log(buffer[i] * normFac));
++ buffer[i] = pointState[i].inTargetRegion() ? std::exp(pointState[i].logPmfSum()) : 0;
+ }
+
+ sumOverSimulations(gmx::ArrayRef<double>(buffer), commRecord, multiSimComm);
+
+ /* Take log again to get (non-normalized) PMF */
+ double normFac = 1.0 / numSharedUpdate;
+ for (gmx::index i = 0; i < pointState.ssize(); i++)
+ {
+ if (pointState[i].inTargetRegion())
+ {
++ pointState[i].setLogPmfSum(std::log(buffer[i] * normFac));
+ }
+ }
+}
+
+/*! \brief
+ * Find the minimum free energy value.
+ *
+ * \param[in] pointState The state of the points.
+ * \returns the minimum free energy value.
+ */
+double freeEnergyMinimumValue(gmx::ArrayRef<const PointState> pointState)
+{
+ double fMin = GMX_FLOAT_MAX;
+
+ for (auto const& ps : pointState)
+ {
+ if (ps.inTargetRegion() && ps.freeEnergy() < fMin)
+ {
+ fMin = ps.freeEnergy();
+ }
+ }
+
+ return fMin;
+}
+
+/*! \brief
+ * Find and return the log of the probability weight of a point given a coordinate value.
+ *
+ * The unnormalized weight is given by
+ * w(point|value) = exp(bias(point) - U(value,point)),
+ * where U is a harmonic umbrella potential.
+ *
+ * \param[in] dimParams The bias dimensions parameters
+ * \param[in] points The point state.
+ * \param[in] grid The grid.
+ * \param[in] pointIndex Point to evaluate probability weight for.
+ * \param[in] pointBias Bias for the point (as a log weight).
+ * \param[in] value Coordinate value.
+ * \param[in] neighborLambdaEnergies The energy of the system in neighboring lambdas states. Can be
+ * empty when there are no free energy lambda state dimensions.
+ * \param[in] gridpointIndex The index of the current grid point.
+ * \returns the log of the biased probability weight.
+ */
+double biasedLogWeightFromPoint(const std::vector<DimParams>& dimParams,
+ const std::vector<PointState>& points,
+ const BiasGrid& grid,
+ int pointIndex,
+ double pointBias,
+ const awh_dvec value,
+ gmx::ArrayRef<const double> neighborLambdaEnergies,
+ int gridpointIndex)
+{
+ double logWeight = detail::c_largeNegativeExponent;
+
+ /* Only points in the target region have non-zero weight */
+ if (points[pointIndex].inTargetRegion())
+ {
+ logWeight = pointBias;
+
+ /* Add potential for all parameter dimensions */
+ for (size_t d = 0; d < dimParams.size(); d++)
+ {
+ if (dimParams[d].isFepLambdaDimension())
+ {
+ /* If this is not a sampling step or if this function is called from
+ * calcConvolvedBias(), when writing energy subblocks, neighborLambdaEnergies will
+ * be empty. No convolution is required along the lambda dimension. */
+ if (!neighborLambdaEnergies.empty())
+ {
+ const int pointLambdaIndex = grid.point(pointIndex).coordValue[d];
+ const int gridpointLambdaIndex = grid.point(gridpointIndex).coordValue[d];
+ logWeight -= dimParams[d].fepDimParams().beta
+ * (neighborLambdaEnergies[pointLambdaIndex]
+ - neighborLambdaEnergies[gridpointLambdaIndex]);
+ }
+ }
+ else
+ {
+ double dev = getDeviationFromPointAlongGridAxis(grid, d, pointIndex, value[d]);
+ logWeight -= 0.5 * dimParams[d].pullDimParams().betak * dev * dev;
+ }
+ }
+ }
+ return logWeight;
+}
+
+/*! \brief
+ * Calculates the marginal distribution (marginal probability) for each value along
+ * a free energy lambda axis.
+ * The marginal distribution of one coordinate dimension value is the sum of the probability
+ * distribution of all values (herein all neighbor values) with the same value in the dimension
+ * of interest.
+ * \param[in] grid The bias grid.
+ * \param[in] neighbors The points to use for the calculation of the marginal distribution.
+ * \param[in] probWeightNeighbor Probability weights of the neighbors.
+ * \returns The calculated marginal distribution in a 1D array with
+ * as many elements as there are points along the axis of interest.
+ */
+std::vector<double> calculateFELambdaMarginalDistribution(const BiasGrid& grid,
+ gmx::ArrayRef<const int> neighbors,
+ gmx::ArrayRef<const double> probWeightNeighbor)
+{
+ const std::optional<int> lambdaAxisIndex = grid.lambdaAxisIndex();
+ GMX_RELEASE_ASSERT(lambdaAxisIndex,
+ "There must be a free energy lambda axis in order to calculate the free "
+ "energy lambda marginal distribution.");
+ const int numFepLambdaStates = grid.numFepLambdaStates();
+ std::vector<double> lambdaMarginalDistribution(numFepLambdaStates, 0);
+
+ for (size_t i = 0; i < neighbors.size(); i++)
+ {
+ const int neighbor = neighbors[i];
+ const int lambdaState = grid.point(neighbor).coordValue[lambdaAxisIndex.value()];
+ lambdaMarginalDistribution[lambdaState] += probWeightNeighbor[i];
+ }
+ return lambdaMarginalDistribution;
+}
+
+} // namespace
+
+void BiasState::calcConvolvedPmf(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ std::vector<float>* convolvedPmf) const
+{
+ size_t numPoints = grid.numPoints();
+
+ convolvedPmf->resize(numPoints);
+
+ /* Get the PMF to convolve. */
+ std::vector<float> pmf(numPoints);
+ getPmf(pmf);
+
+ for (size_t m = 0; m < numPoints; m++)
+ {
+ double freeEnergyWeights = 0;
+ const GridPoint& point = grid.point(m);
+ for (auto& neighbor : point.neighbor)
+ {
+ /* Do not convolve the bias along a lambda axis - only use the pmf from the current point */
+ if (!pointsHaveDifferentLambda(grid, m, neighbor))
+ {
+ /* The negative PMF is a positive bias. */
+ double biasNeighbor = -pmf[neighbor];
+
+ /* Add the convolved PMF weights for the neighbors of this point.
+ Note that this function only adds point within the target > 0 region.
+ Sum weights, take the logarithm last to get the free energy. */
+ double logWeight = biasedLogWeightFromPoint(dimParams, points_, grid, neighbor,
+ biasNeighbor, point.coordValue, {}, m);
+ freeEnergyWeights += std::exp(logWeight);
+ }
+ }
+
+ GMX_RELEASE_ASSERT(freeEnergyWeights > 0,
+ "Attempting to do log(<= 0) in AWH convolved PMF calculation.");
+ (*convolvedPmf)[m] = -std::log(static_cast<float>(freeEnergyWeights));
+ }
+}
+
+namespace
+{
+
+/*! \brief
+ * Updates the target distribution for all points.
+ *
+ * The target distribution is always updated for all points
+ * at the same time.
+ *
+ * \param[in,out] pointState The state of all points.
+ * \param[in] params The bias parameters.
+ */
+void updateTargetDistribution(gmx::ArrayRef<PointState> pointState, const BiasParams& params)
+{
+ double freeEnergyCutoff = 0;
+ if (params.eTarget == eawhtargetCUTOFF)
+ {
+ freeEnergyCutoff = freeEnergyMinimumValue(pointState) + params.freeEnergyCutoffInKT;
+ }
+
+ double sumTarget = 0;
+ for (PointState& ps : pointState)
+ {
+ sumTarget += ps.updateTargetWeight(params, freeEnergyCutoff);
+ }
+ GMX_RELEASE_ASSERT(sumTarget > 0, "We should have a non-zero distribution");
+
+ /* Normalize to 1 */
+ double invSum = 1.0 / sumTarget;
+ for (PointState& ps : pointState)
+ {
+ ps.scaleTarget(invSum);
+ }
+}
+
+/*! \brief
+ * Puts together a string describing a grid point.
+ *
+ * \param[in] grid The grid.
+ * \param[in] point BiasGrid point index.
+ * \returns a string for the point.
+ */
+std::string gridPointValueString(const BiasGrid& grid, int point)
+{
+ std::string pointString;
+
+ pointString += "(";
+
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ pointString += gmx::formatString("%g", grid.point(point).coordValue[d]);
+ if (d < grid.numDimensions() - 1)
+ {
+ pointString += ",";
+ }
+ else
+ {
+ pointString += ")";
+ }
+ }
+
+ return pointString;
+}
+
+} // namespace
+
+int BiasState::warnForHistogramAnomalies(const BiasGrid& grid, int biasIndex, double t, FILE* fplog, int maxNumWarnings) const
+{
+ GMX_ASSERT(fplog != nullptr, "Warnings can only be issued if there is log file.");
+ const double maxHistogramRatio = 0.5; /* Tolerance for printing a warning about the histogram ratios */
+
+ /* Sum up the histograms and get their normalization */
+ double sumVisits = 0;
+ double sumWeights = 0;
+ for (auto& pointState : points_)
+ {
+ if (pointState.inTargetRegion())
+ {
+ sumVisits += pointState.numVisitsTot();
+ sumWeights += pointState.weightSumTot();
+ }
+ }
+ GMX_RELEASE_ASSERT(sumVisits > 0, "We should have visits");
+ GMX_RELEASE_ASSERT(sumWeights > 0, "We should have weight");
+ double invNormVisits = 1.0 / sumVisits;
+ double invNormWeight = 1.0 / sumWeights;
+
+ /* Check all points for warnings */
+ int numWarnings = 0;
+ size_t numPoints = grid.numPoints();
+ for (size_t m = 0; m < numPoints; m++)
+ {
+ /* Skip points close to boundary or non-target region */
+ const GridPoint& gridPoint = grid.point(m);
+ bool skipPoint = false;
+ for (size_t n = 0; (n < gridPoint.neighbor.size()) && !skipPoint; n++)
+ {
+ int neighbor = gridPoint.neighbor[n];
+ skipPoint = !points_[neighbor].inTargetRegion();
+ for (int d = 0; (d < grid.numDimensions()) && !skipPoint; d++)
+ {
+ const GridPoint& neighborPoint = grid.point(neighbor);
+ skipPoint = neighborPoint.index[d] == 0
+ || neighborPoint.index[d] == grid.axis(d).numPoints() - 1;
+ }
+ }
+
+ /* Warn if the coordinate distribution is less than the target distribution with a certain fraction somewhere */
+ const double relativeWeight = points_[m].weightSumTot() * invNormWeight;
+ const double relativeVisits = points_[m].numVisitsTot() * invNormVisits;
+ if (!skipPoint && relativeVisits < relativeWeight * maxHistogramRatio)
+ {
+ std::string pointValueString = gridPointValueString(grid, m);
+ std::string warningMessage = gmx::formatString(
+ "\nawh%d warning: "
+ "at t = %g ps the obtained coordinate distribution at coordinate value %s "
+ "is less than a fraction %g of the reference distribution at that point. "
+ "If you are not certain about your settings you might want to increase your "
+ "pull force constant or "
+ "modify your sampling region.\n",
+ biasIndex + 1, t, pointValueString.c_str(), maxHistogramRatio);
+ gmx::TextLineWrapper wrapper;
+ wrapper.settings().setLineLength(c_linewidth);
+ fprintf(fplog, "%s", wrapper.wrapToString(warningMessage).c_str());
+
+ numWarnings++;
+ }
+ if (numWarnings >= maxNumWarnings)
+ {
+ break;
+ }
+ }
+
+ return numWarnings;
+}
+
+double BiasState::calcUmbrellaForceAndPotential(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ int point,
+ ArrayRef<const double> neighborLambdaDhdl,
+ gmx::ArrayRef<double> force) const
+{
+ double potential = 0;
+ for (size_t d = 0; d < dimParams.size(); d++)
+ {
+ if (dimParams[d].isFepLambdaDimension())
+ {
+ if (!neighborLambdaDhdl.empty())
+ {
+ const int coordpointLambdaIndex = grid.point(point).coordValue[d];
+ force[d] = neighborLambdaDhdl[coordpointLambdaIndex];
+ /* The potential should not be affected by the lambda dimension. */
+ }
+ }
+ else
+ {
+ double deviation =
+ getDeviationFromPointAlongGridAxis(grid, d, point, coordState_.coordValue()[d]);
+ double k = dimParams[d].pullDimParams().k;
+
+ /* Force from harmonic potential 0.5*k*dev^2 */
+ force[d] = -k * deviation;
+ potential += 0.5 * k * deviation * deviation;
+ }
+ }
+
+ return potential;
+}
+
+void BiasState::calcConvolvedForce(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ gmx::ArrayRef<const double> probWeightNeighbor,
+ ArrayRef<const double> neighborLambdaDhdl,
+ gmx::ArrayRef<double> forceWorkBuffer,
+ gmx::ArrayRef<double> force) const
+{
+ for (size_t d = 0; d < dimParams.size(); d++)
+ {
+ force[d] = 0;
+ }
+
+ /* Only neighboring points have non-negligible contribution. */
+ const std::vector<int>& neighbor = grid.point(coordState_.gridpointIndex()).neighbor;
+ gmx::ArrayRef<double> forceFromNeighbor = forceWorkBuffer;
+ for (size_t n = 0; n < neighbor.size(); n++)
+ {
+ double weightNeighbor = probWeightNeighbor[n];
+ int indexNeighbor = neighbor[n];
+
+ /* Get the umbrella force from this point. The returned potential is ignored here. */
+ calcUmbrellaForceAndPotential(dimParams, grid, indexNeighbor, neighborLambdaDhdl, forceFromNeighbor);
+
+ /* Add the weighted umbrella force to the convolved force. */
+ for (size_t d = 0; d < dimParams.size(); d++)
+ {
+ force[d] += forceFromNeighbor[d] * weightNeighbor;
+ }
+ }
+}
+
+double BiasState::moveUmbrella(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ gmx::ArrayRef<const double> probWeightNeighbor,
+ ArrayRef<const double> neighborLambdaDhdl,
+ gmx::ArrayRef<double> biasForce,
+ int64_t step,
+ int64_t seed,
+ int indexSeed,
+ bool onlySampleUmbrellaGridpoint)
+{
+ /* Generate and set a new coordinate reference value */
+ coordState_.sampleUmbrellaGridpoint(grid, coordState_.gridpointIndex(), probWeightNeighbor,
+ step, seed, indexSeed);
+
+ if (onlySampleUmbrellaGridpoint)
+ {
+ return 0;
+ }
+
+ std::vector<double> newForce(dimParams.size());
+ double newPotential = calcUmbrellaForceAndPotential(
+ dimParams, grid, coordState_.umbrellaGridpoint(), neighborLambdaDhdl, newForce);
+
+ /* A modification of the reference value at time t will lead to a different
+ force over t-dt/2 to t and over t to t+dt/2. For high switching rates
+ this means the force and velocity will change signs roughly as often.
+ To avoid any issues we take the average of the previous and new force
+ at steps when the reference value has been moved. E.g. if the ref. value
+ is set every step to (coord dvalue +/- delta) would give zero force.
+ */
+ for (gmx::index d = 0; d < biasForce.ssize(); d++)
+ {
+ /* Average of the current and new force */
+ biasForce[d] = 0.5 * (biasForce[d] + newForce[d]);
+ }
+
+ return newPotential;
+}
+
+namespace
+{
+
+/*! \brief
+ * Sets the histogram rescaling factors needed to control the histogram size.
+ *
+ * For sake of robustness, the reference weight histogram can grow at a rate
+ * different from the actual sampling rate. Typically this happens for a limited
+ * initial time, alternatively growth is scaled down by a constant factor for all
+ * times. Since the size of the reference histogram sets the size of the free
+ * energy update this should be reflected also in the PMF. Thus the PMF histogram
+ * needs to be rescaled too.
+ *
+ * This function should only be called by the bias update function or wrapped by a function that
+ * knows what scale factors should be applied when, e.g,
+ * getSkippedUpdateHistogramScaleFactors().
+ *
+ * \param[in] params The bias parameters.
+ * \param[in] newHistogramSize New reference weight histogram size.
+ * \param[in] oldHistogramSize Previous reference weight histogram size (before adding new samples).
+ * \param[out] weightHistScaling Scaling factor for the reference weight histogram.
+ * \param[out] logPmfSumScaling Log of the scaling factor for the PMF histogram.
+ */
+void setHistogramUpdateScaleFactors(const BiasParams& params,
+ double newHistogramSize,
+ double oldHistogramSize,
+ double* weightHistScaling,
+ double* logPmfSumScaling)
+{
+
+ /* The two scaling factors below are slightly different (ignoring the log factor) because the
+ reference and the PMF histogram apply weight scaling differently. The weight histogram
+ applies is locally, i.e. each sample is scaled down meaning all samples get equal weight.
+ It is done this way because that is what target type local Boltzmann (for which
+ target = weight histogram) needs. In contrast, the PMF histogram is rescaled globally
+ by repeatedly scaling down the whole histogram. The reasons for doing it this way are:
+ 1) empirically this is necessary for converging the PMF; 2) since the extraction of
+ the PMF is theoretically only valid for a constant bias, new samples should get more
+ weight than old ones for which the bias is fluctuating more. */
+ *weightHistScaling =
+ newHistogramSize / (oldHistogramSize + params.updateWeight * params.localWeightScaling);
+ *logPmfSumScaling = std::log(newHistogramSize / (oldHistogramSize + params.updateWeight));
+}
+
+} // namespace
+
+void BiasState::getSkippedUpdateHistogramScaleFactors(const BiasParams& params,
+ double* weightHistScaling,
+ double* logPmfSumScaling) const
+{
+ GMX_ASSERT(params.skipUpdates(),
+ "Calling function for skipped updates when skipping updates is not allowed");
+
+ if (inInitialStage())
+ {
+ /* In between global updates the reference histogram size is kept constant so we trivially
+ know what the histogram size was at the time of the skipped update. */
+ double histogramSize = histogramSize_.histogramSize();
+ setHistogramUpdateScaleFactors(params, histogramSize, histogramSize, weightHistScaling,
+ logPmfSumScaling);
+ }
+ else
+ {
+ /* In the final stage, the reference histogram grows at the sampling rate which gives trivial scale factors. */
+ *weightHistScaling = 1;
+ *logPmfSumScaling = 0;
+ }
+}
+
+void BiasState::doSkippedUpdatesForAllPoints(const BiasParams& params)
+{
+ double weightHistScaling;
+ double logPmfsumScaling;
+
+ getSkippedUpdateHistogramScaleFactors(params, &weightHistScaling, &logPmfsumScaling);
+
+ for (auto& pointState : points_)
+ {
+ bool didUpdate = pointState.performPreviouslySkippedUpdates(
+ params, histogramSize_.numUpdates(), weightHistScaling, logPmfsumScaling);
+
+ /* Update the bias for this point only if there were skipped updates in the past to avoid calculating the log unneccessarily */
+ if (didUpdate)
+ {
+ pointState.updateBias();
+ }
+ }
+}
+
+void BiasState::doSkippedUpdatesInNeighborhood(const BiasParams& params, const BiasGrid& grid)
+{
+ double weightHistScaling;
+ double logPmfsumScaling;
+
+ getSkippedUpdateHistogramScaleFactors(params, &weightHistScaling, &logPmfsumScaling);
+
+ /* For each neighbor point of the center point, refresh its state by adding the results of all past, skipped updates. */
+ const std::vector<int>& neighbors = grid.point(coordState_.gridpointIndex()).neighbor;
+ for (auto& neighbor : neighbors)
+ {
+ bool didUpdate = points_[neighbor].performPreviouslySkippedUpdates(
+ params, histogramSize_.numUpdates(), weightHistScaling, logPmfsumScaling);
+
+ if (didUpdate)
+ {
+ points_[neighbor].updateBias();
+ }
+ }
+}
+
+namespace
+{
+
+/*! \brief
+ * Merge update lists from multiple sharing simulations.
+ *
+ * \param[in,out] updateList Update list for this simulation (assumed >= npoints long).
+ * \param[in] numPoints Total number of points.
+ * \param[in] commRecord Struct for intra-simulation communication.
+ * \param[in] multiSimComm Struct for multi-simulation communication.
+ */
+void mergeSharedUpdateLists(std::vector<int>* updateList,
+ int numPoints,
+ const t_commrec* commRecord,
+ const gmx_multisim_t* multiSimComm)
+{
+ std::vector<int> numUpdatesOfPoint;
+
+ /* Flag the update points of this sim.
+ TODO: we can probably avoid allocating this array and just use the input array. */
+ numUpdatesOfPoint.resize(numPoints, 0);
+ for (auto& pointIndex : *updateList)
+ {
+ numUpdatesOfPoint[pointIndex] = 1;
+ }
+
+ /* Sum over the sims to get all the flagged points */
+ sumOverSimulations(arrayRefFromArray(numUpdatesOfPoint.data(), numPoints), commRecord, multiSimComm);
+
+ /* Collect the indices of the flagged points in place. The resulting array will be the merged update list.*/
+ updateList->clear();
+ for (int m = 0; m < numPoints; m++)
+ {
+ if (numUpdatesOfPoint[m] > 0)
+ {
+ updateList->push_back(m);
+ }
+ }
+}
+
+/*! \brief
+ * Generate an update list of points sampled since the last update.
+ *
+ * \param[in] grid The AWH bias.
+ * \param[in] points The point state.
+ * \param[in] originUpdatelist The origin of the rectangular region that has been sampled since
+ * last update. \param[in] endUpdatelist The end of the rectangular that has been sampled since
+ * last update. \param[in,out] updateList Local update list to set (assumed >= npoints long).
+ */
+void makeLocalUpdateList(const BiasGrid& grid,
+ const std::vector<PointState>& points,
+ const awh_ivec originUpdatelist,
+ const awh_ivec endUpdatelist,
+ std::vector<int>* updateList)
+{
+ awh_ivec origin;
+ awh_ivec numPoints;
+
+ /* Define the update search grid */
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ origin[d] = originUpdatelist[d];
+ numPoints[d] = endUpdatelist[d] - originUpdatelist[d] + 1;
+
+ /* Because the end_updatelist is unwrapped it can be > (npoints - 1) so that numPoints can be > npoints in grid.
+ This helps for calculating the distance/number of points but should be removed and fixed when the way of
+ updating origin/end updatelist is changed (see sampleProbabilityWeights). */
+ numPoints[d] = std::min(grid.axis(d).numPoints(), numPoints[d]);
+ }
+
+ /* Make the update list */
+ updateList->clear();
+ int pointIndex = -1;
+ bool pointExists = true;
+ while (pointExists)
+ {
+ pointExists = advancePointInSubgrid(grid, origin, numPoints, &pointIndex);
+
+ if (pointExists && points[pointIndex].inTargetRegion())
+ {
+ updateList->push_back(pointIndex);
+ }
+ }
+}
+
+} // namespace
+
+void BiasState::resetLocalUpdateRange(const BiasGrid& grid)
+{
+ const int gridpointIndex = coordState_.gridpointIndex();
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ /* This gives the minimum range consisting only of the current closest point. */
+ originUpdatelist_[d] = grid.point(gridpointIndex).index[d];
+ endUpdatelist_[d] = grid.point(gridpointIndex).index[d];
+ }
+}
+
+namespace
+{
+
+/*! \brief
+ * Add partial histograms (accumulating between updates) to accumulating histograms.
+ *
+ * \param[in,out] pointState The state of the points in the bias.
+ * \param[in,out] weightSumCovering The weights for checking covering.
+ * \param[in] numSharedUpdate The number of biases sharing the histrogram.
+ * \param[in] commRecord Struct for intra-simulation communication.
+ * \param[in] multiSimComm Struct for multi-simulation communication.
+ * \param[in] localUpdateList List of points with data.
+ */
+void sumHistograms(gmx::ArrayRef<PointState> pointState,
+ gmx::ArrayRef<double> weightSumCovering,
+ int numSharedUpdate,
+ const t_commrec* commRecord,
+ const gmx_multisim_t* multiSimComm,
+ const std::vector<int>& localUpdateList)
+{
+ /* The covering checking histograms are added before summing over simulations, so that the
+ weights from different simulations are kept distinguishable. */
+ for (int globalIndex : localUpdateList)
+ {
+ weightSumCovering[globalIndex] += pointState[globalIndex].weightSumIteration();
+ }
+
+ /* Sum histograms over multiple simulations if needed. */
+ if (numSharedUpdate > 1)
+ {
+ GMX_ASSERT(numSharedUpdate == multiSimComm->numSimulations_,
+ "Sharing within a simulation is not implemented (yet)");
+
+ /* Collect the weights and counts in linear arrays to be able to use gmx_sumd_sim. */
+ std::vector<double> weightSum;
+ std::vector<double> coordVisits;
+
+ weightSum.resize(localUpdateList.size());
+ coordVisits.resize(localUpdateList.size());
+
+ for (size_t localIndex = 0; localIndex < localUpdateList.size(); localIndex++)
+ {
+ const PointState& ps = pointState[localUpdateList[localIndex]];
+
+ weightSum[localIndex] = ps.weightSumIteration();
+ coordVisits[localIndex] = ps.numVisitsIteration();
+ }
+
+ sumOverSimulations(gmx::ArrayRef<double>(weightSum), commRecord, multiSimComm);
+ sumOverSimulations(gmx::ArrayRef<double>(coordVisits), commRecord, multiSimComm);
+
+ /* Transfer back the result */
+ for (size_t localIndex = 0; localIndex < localUpdateList.size(); localIndex++)
+ {
+ PointState& ps = pointState[localUpdateList[localIndex]];
+
+ ps.setPartialWeightAndCount(weightSum[localIndex], coordVisits[localIndex]);
+ }
+ }
+
+ /* Now add the partial counts and weights to the accumulating histograms.
+ Note: we still need to use the weights for the update so we wait
+ with resetting them until the end of the update. */
+ for (int globalIndex : localUpdateList)
+ {
+ pointState[globalIndex].addPartialWeightAndCount();
+ }
+}
+
+/*! \brief
+ * Label points along an axis as covered or not.
+ *
+ * A point is covered if it is surrounded by visited points up to a radius = coverRadius.
+ *
+ * \param[in] visited Visited? For each point.
+ * \param[in] checkCovering Check for covering? For each point.
+ * \param[in] numPoints The number of grid points along this dimension.
+ * \param[in] period Period in number of points.
+ * \param[in] coverRadius Cover radius, in points, needed for defining a point as covered.
+ * \param[in,out] covered In this array elements are 1 for covered points and 0 for
+ * non-covered points, this routine assumes that \p covered has at least size \p numPoints.
+ */
+void labelCoveredPoints(const std::vector<bool>& visited,
+ const std::vector<bool>& checkCovering,
+ int numPoints,
+ int period,
+ int coverRadius,
+ gmx::ArrayRef<int> covered)
+{
+ GMX_ASSERT(covered.ssize() >= numPoints, "covered should be at least as large as the grid");
+
+ bool haveFirstNotVisited = false;
+ int firstNotVisited = -1;
+ int notVisitedLow = -1;
+ int notVisitedHigh = -1;
+
+ for (int n = 0; n < numPoints; n++)
+ {
+ if (checkCovering[n] && !visited[n])
+ {
+ if (!haveFirstNotVisited)
+ {
+ notVisitedLow = n;
+ firstNotVisited = n;
+ haveFirstNotVisited = true;
+ }
+ else
+ {
+ notVisitedHigh = n;
+
+ /* Have now an interval I = [notVisitedLow,notVisitedHigh] of visited points bounded
+ by unvisited points. The unvisted end points affect the coveredness of the
+ visited with a reach equal to the cover radius. */
+ int notCoveredLow = notVisitedLow + coverRadius;
+ int notCoveredHigh = notVisitedHigh - coverRadius;
+ for (int i = notVisitedLow; i <= notVisitedHigh; i++)
+ {
+ covered[i] = static_cast<int>((i > notCoveredLow) && (i < notCoveredHigh));
+ }
+
+ /* Find a new interval to set covering for. Make the notVisitedHigh of this interval
+ the notVisitedLow of the next. */
+ notVisitedLow = notVisitedHigh;
+ }
+ }
+ }
+
+ /* Have labelled all the internal points. Now take care of the boundary regions. */
+ if (!haveFirstNotVisited)
+ {
+ /* No non-visited points <=> all points visited => all points covered. */
+
+ for (int n = 0; n < numPoints; n++)
+ {
+ covered[n] = 1;
+ }
+ }
+ else
+ {
+ int lastNotVisited = notVisitedLow;
+
+ /* For periodic boundaries, non-visited points can influence points
+ on the other side of the boundary so we need to wrap around. */
+
+ /* Lower end. For periodic boundaries the last upper end not visited point becomes the low-end not visited point.
+ For non-periodic boundaries there is no lower end point so a dummy value is used. */
+ int notVisitedHigh = firstNotVisited;
+ int notVisitedLow = period > 0 ? (lastNotVisited - period) : -(coverRadius + 1);
+
+ int notCoveredLow = notVisitedLow + coverRadius;
+ int notCoveredHigh = notVisitedHigh - coverRadius;
+
+ for (int i = 0; i <= notVisitedHigh; i++)
+ {
+ /* For non-periodic boundaries notCoveredLow = -1 will impose no restriction. */
+ covered[i] = static_cast<int>((i > notCoveredLow) && (i < notCoveredHigh));
+ }
+
+ /* Upper end. Same as for lower end but in the other direction. */
+ notVisitedHigh = period > 0 ? (firstNotVisited + period) : (numPoints + coverRadius);
+ notVisitedLow = lastNotVisited;
+
+ notCoveredLow = notVisitedLow + coverRadius;
+ notCoveredHigh = notVisitedHigh - coverRadius;
+
+ for (int i = notVisitedLow; i <= numPoints - 1; i++)
+ {
+ /* For non-periodic boundaries notCoveredHigh = numPoints will impose no restriction. */
+ covered[i] = static_cast<int>((i > notCoveredLow) && (i < notCoveredHigh));
+ }
+ }
+}
+
+} // namespace
+
+bool BiasState::isSamplingRegionCovered(const BiasParams& params,
+ const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ const t_commrec* commRecord,
+ const gmx_multisim_t* multiSimComm) const
+{
+ /* Allocate and initialize arrays: one for checking visits along each dimension,
+ one for keeping track of which points to check and one for the covered points.
+ Possibly these could be kept as AWH variables to avoid these allocations. */
+ struct CheckDim
+ {
+ std::vector<bool> visited;
+ std::vector<bool> checkCovering;
+ // We use int for the covering array since we might use gmx_sumi_sim.
+ std::vector<int> covered;
+ };
+
+ std::vector<CheckDim> checkDim;
+ checkDim.resize(grid.numDimensions());
+
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ const size_t numPoints = grid.axis(d).numPoints();
+ checkDim[d].visited.resize(numPoints, false);
+ checkDim[d].checkCovering.resize(numPoints, false);
+ checkDim[d].covered.resize(numPoints, 0);
+ }
+
+ /* Set visited points along each dimension and which points should be checked for covering.
+ Specifically, points above the free energy cutoff (if there is one) or points outside
+ of the target region are ignored. */
+
+ /* Set the free energy cutoff */
+ double maxFreeEnergy = GMX_FLOAT_MAX;
+
+ if (params.eTarget == eawhtargetCUTOFF)
+ {
+ maxFreeEnergy = freeEnergyMinimumValue(points_) + params.freeEnergyCutoffInKT;
+ }
+
+ /* Set the threshold weight for a point to be considered visited. */
+ double weightThreshold = 1;
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ if (grid.axis(d).isFepLambdaAxis())
+ {
+ /* Do not modify the weight threshold based on a FEP lambda axis. The spread
+ * of the sampling weights is not depending on a Gaussian distribution (like
+ * below). */
+ weightThreshold *= 1.0;
+ }
+ else
+ {
+ /* The spacing is proportional to 1/sqrt(betak). The weight threshold will be
+ * approximately (given that the spacing can be modified if the dimension is periodic)
+ * proportional to sqrt(1/(2*pi)). */
+ weightThreshold *= grid.axis(d).spacing()
+ * std::sqrt(dimParams[d].pullDimParams().betak * 0.5 * M_1_PI);
+ }
+ }
+
+ /* Project the sampling weights onto each dimension */
+ for (size_t m = 0; m < grid.numPoints(); m++)
+ {
+ const PointState& pointState = points_[m];
+
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ int n = grid.point(m).index[d];
+
+ /* Is visited if it was already visited or if there is enough weight at the current point */
+ checkDim[d].visited[n] = checkDim[d].visited[n] || (weightSumCovering_[m] > weightThreshold);
+
+ /* Check for covering if there is at least point in this slice that is in the target region and within the cutoff */
+ checkDim[d].checkCovering[n] =
+ checkDim[d].checkCovering[n]
+ || (pointState.inTargetRegion() && pointState.freeEnergy() < maxFreeEnergy);
+ }
+ }
+
+ /* Label each point along each dimension as covered or not. */
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ labelCoveredPoints(checkDim[d].visited, checkDim[d].checkCovering, grid.axis(d).numPoints(),
+ grid.axis(d).numPointsInPeriod(), params.coverRadius()[d], checkDim[d].covered);
+ }
+
+ /* Now check for global covering. Each dimension needs to be covered separately.
+ A dimension is covered if each point is covered. Multiple simulations collectively
+ cover the points, i.e. a point is covered if any of the simulations covered it.
+ However, visited points are not shared, i.e. if a point is covered or not is
+ determined by the visits of a single simulation. In general the covering criterion is
+ all points covered => all points are surrounded by visited points up to a radius = coverRadius.
+ For 1 simulation, all points covered <=> all points visited. For multiple simulations
+ however, all points visited collectively !=> all points covered, except for coverRadius = 0.
+ In the limit of large coverRadius, all points covered => all points visited by at least one
+ simulation (since no point will be covered until all points have been visited by a
+ single simulation). Basically coverRadius sets how much "connectedness" (or mixing) a point
+ needs with surrounding points before sharing covering information with other simulations. */
+
+ /* Communicate the covered points between sharing simulations if needed. */
+ if (params.numSharedUpdate > 1)
+ {
+ /* For multiple dimensions this may not be the best way to do it. */
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ sumOverSimulations(
+ gmx::arrayRefFromArray(checkDim[d].covered.data(), grid.axis(d).numPoints()),
+ commRecord, multiSimComm);
+ }
+ }
+
+ /* Now check if for each dimension all points are covered. Break if not true. */
+ bool allPointsCovered = true;
+ for (int d = 0; d < grid.numDimensions() && allPointsCovered; d++)
+ {
+ for (int n = 0; n < grid.axis(d).numPoints() && allPointsCovered; n++)
+ {
+ allPointsCovered = (checkDim[d].covered[n] != 0);
+ }
+ }
+
+ return allPointsCovered;
+}
+
+/*! \brief
+ * Normalizes the free energy and PMF sum.
+ *
+ * \param[in] pointState The state of the points.
+ */
+static void normalizeFreeEnergyAndPmfSum(std::vector<PointState>* pointState)
+{
+ double minF = freeEnergyMinimumValue(*pointState);
+
+ for (PointState& ps : *pointState)
+ {
+ ps.normalizeFreeEnergyAndPmfSum(minF);
+ }
+}
+
+void BiasState::updateFreeEnergyAndAddSamplesToHistogram(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ const BiasParams& params,
+ const t_commrec* commRecord,
+ const gmx_multisim_t* multiSimComm,
+ double t,
+ int64_t step,
+ FILE* fplog,
+ std::vector<int>* updateList)
+{
+ /* Note hat updateList is only used in this scope and is always
+ * re-initialized. We do not use a local vector, because that would
+ * cause reallocation every time this funtion is called and the vector
+ * can be the size of the whole grid.
+ */
+
+ /* Make a list of all local points, i.e. those that could have been touched since
+ the last update. These are the points needed for summing histograms below
+ (non-local points only add zeros). For local updates, this will also be the
+ final update list. */
+ makeLocalUpdateList(grid, points_, originUpdatelist_, endUpdatelist_, updateList);
+ if (params.numSharedUpdate > 1)
+ {
+ mergeSharedUpdateLists(updateList, points_.size(), commRecord, multiSimComm);
+ }
+
+ /* Reset the range for the next update */
+ resetLocalUpdateRange(grid);
+
+ /* Add samples to histograms for all local points and sync simulations if needed */
+ sumHistograms(points_, weightSumCovering_, params.numSharedUpdate, commRecord, multiSimComm, *updateList);
+
+ sumPmf(points_, params.numSharedUpdate, commRecord, multiSimComm);
+
+ /* Renormalize the free energy if values are too large. */
+ bool needToNormalizeFreeEnergy = false;
+ for (int& globalIndex : *updateList)
+ {
+ /* We want to keep the absolute value of the free energies to be less
+ c_largePositiveExponent to be able to safely pass these values to exp(). The check below
+ ensures this as long as the free energy values grow less than 0.5*c_largePositiveExponent
+ in a return time to this neighborhood. For reasonable update sizes it's unlikely that
+ this requirement would be broken. */
+ if (std::abs(points_[globalIndex].freeEnergy()) > 0.5 * detail::c_largePositiveExponent)
+ {
+ needToNormalizeFreeEnergy = true;
+ break;
+ }
+ }
+
+ /* Update target distribution? */
+ bool needToUpdateTargetDistribution =
+ (params.eTarget != eawhtargetCONSTANT && params.isUpdateTargetStep(step));
+
+ /* In the initial stage, the histogram grows dynamically as a function of the number of coverings. */
+ bool detectedCovering = false;
+ if (inInitialStage())
+ {
+ detectedCovering =
+ (params.isCheckCoveringStep(step)
+ && isSamplingRegionCovered(params, dimParams, grid, commRecord, multiSimComm));
+ }
+
+ /* The weighthistogram size after this update. */
+ double newHistogramSize = histogramSize_.newHistogramSize(params, t, detectedCovering, points_,
+ weightSumCovering_, fplog);
+
+ /* Make the update list. Usually we try to only update local points,
+ * but if the update has non-trivial or non-deterministic effects
+ * on non-local points a global update is needed. This is the case when:
+ * 1) a covering occurred in the initial stage, leading to non-trivial
+ * histogram rescaling factors; or
+ * 2) the target distribution will be updated, since we don't make any
+ * assumption on its form; or
+ * 3) the AWH parameters are such that we never attempt to skip non-local
+ * updates; or
+ * 4) the free energy values have grown so large that a renormalization
+ * is needed.
+ */
+ if (needToUpdateTargetDistribution || detectedCovering || !params.skipUpdates() || needToNormalizeFreeEnergy)
+ {
+ /* Global update, just add all points. */
+ updateList->clear();
+ for (size_t m = 0; m < points_.size(); m++)
+ {
+ if (points_[m].inTargetRegion())
+ {
+ updateList->push_back(m);
+ }
+ }
+ }
+
+ /* Set histogram scale factors. */
+ double weightHistScalingSkipped = 0;
+ double logPmfsumScalingSkipped = 0;
+ if (params.skipUpdates())
+ {
+ getSkippedUpdateHistogramScaleFactors(params, &weightHistScalingSkipped, &logPmfsumScalingSkipped);
+ }
+ double weightHistScalingNew;
+ double logPmfsumScalingNew;
+ setHistogramUpdateScaleFactors(params, newHistogramSize, histogramSize_.histogramSize(),
+ &weightHistScalingNew, &logPmfsumScalingNew);
+
+ /* Update free energy and reference weight histogram for points in the update list. */
+ for (int pointIndex : *updateList)
+ {
+ PointState* pointStateToUpdate = &points_[pointIndex];
+
+ /* Do updates from previous update steps that were skipped because this point was at that time non-local. */
+ if (params.skipUpdates())
+ {
+ pointStateToUpdate->performPreviouslySkippedUpdates(params, histogramSize_.numUpdates(),
+ weightHistScalingSkipped,
+ logPmfsumScalingSkipped);
+ }
+
+ /* Now do an update with new sampling data. */
+ pointStateToUpdate->updateWithNewSampling(params, histogramSize_.numUpdates(),
+ weightHistScalingNew, logPmfsumScalingNew);
+ }
+
+ /* Only update the histogram size after we are done with the local point updates */
+ histogramSize_.setHistogramSize(newHistogramSize, weightHistScalingNew);
+
+ if (needToNormalizeFreeEnergy)
+ {
+ normalizeFreeEnergyAndPmfSum(&points_);
+ }
+
+ if (needToUpdateTargetDistribution)
+ {
+ /* The target distribution is always updated for all points at once. */
+ updateTargetDistribution(points_, params);
+ }
+
+ /* Update the bias. The bias is updated separately and last since it simply a function of
+ the free energy and the target distribution and we want to avoid doing extra work. */
+ for (int pointIndex : *updateList)
+ {
+ points_[pointIndex].updateBias();
+ }
+
+ /* Increase the update counter. */
+ histogramSize_.incrementNumUpdates();
+}
+
+double BiasState::updateProbabilityWeightsAndConvolvedBias(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ gmx::ArrayRef<const double> neighborLambdaEnergies,
+ std::vector<double, AlignedAllocator<double>>* weight) const
+{
+ /* Only neighbors of the current coordinate value will have a non-negligible chance of getting sampled */
+ const std::vector<int>& neighbors = grid.point(coordState_.gridpointIndex()).neighbor;
+
+#if GMX_SIMD_HAVE_DOUBLE
+ typedef SimdDouble PackType;
+ constexpr int packSize = GMX_SIMD_DOUBLE_WIDTH;
+#else
+ typedef double PackType;
+ constexpr int packSize = 1;
+#endif
+ /* Round the size of the weight array up to packSize */
+ const int weightSize = ((neighbors.size() + packSize - 1) / packSize) * packSize;
+ weight->resize(weightSize);
+
+ double* gmx_restrict weightData = weight->data();
+ PackType weightSumPack(0.0);
+ for (size_t i = 0; i < neighbors.size(); i += packSize)
+ {
+ for (size_t n = i; n < i + packSize; n++)
+ {
+ if (n < neighbors.size())
+ {
+ const int neighbor = neighbors[n];
+ (*weight)[n] = biasedLogWeightFromPoint(
+ dimParams, points_, grid, neighbor, points_[neighbor].bias(),
+ coordState_.coordValue(), neighborLambdaEnergies, coordState_.gridpointIndex());
+ }
+ else
+ {
+ /* Pad with values that don't affect the result */
+ (*weight)[n] = detail::c_largeNegativeExponent;
+ }
+ }
+ PackType weightPack = load<PackType>(weightData + i);
+ weightPack = gmx::exp(weightPack);
+ weightSumPack = weightSumPack + weightPack;
+ store(weightData + i, weightPack);
+ }
+ /* Sum of probability weights */
+ double weightSum = reduce(weightSumPack);
+ GMX_RELEASE_ASSERT(weightSum > 0,
+ "zero probability weight when updating AWH probability weights.");
+
+ /* Normalize probabilities to sum to 1 */
+ double invWeightSum = 1 / weightSum;
+
+ /* When there is a free energy lambda state axis remove the convolved contributions along that
+ * axis from the total bias. This must be done after calculating invWeightSum (since weightSum
+ * will be modified), but before normalizing the weights (below). */
+ if (grid.hasLambdaAxis())
+ {
+ /* If there is only one axis the bias will not be convolved in any dimension. */
+ if (grid.axis().size() == 1)
+ {
+ weightSum = gmx::exp(points_[coordState_.gridpointIndex()].bias());
+ }
+ else
+ {
+ for (size_t i = 0; i < neighbors.size(); i++)
+ {
+ const int neighbor = neighbors[i];
+ if (pointsHaveDifferentLambda(grid, coordState_.gridpointIndex(), neighbor))
+ {
+ weightSum -= weightData[i];
+ }
+ }
+ }
+ }
+
+ for (double& w : *weight)
+ {
+ w *= invWeightSum;
+ }
+
+ /* Return the convolved bias */
+ return std::log(weightSum);
+}
+
+double BiasState::calcConvolvedBias(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ const awh_dvec& coordValue) const
+{
+ int point = grid.nearestIndex(coordValue);
+ const GridPoint& gridPoint = grid.point(point);
+
+ /* Sum the probability weights from the neighborhood of the given point */
+ double weightSum = 0;
+ for (int neighbor : gridPoint.neighbor)
+ {
+ /* No convolution is required along the lambda dimension. */
+ if (pointsHaveDifferentLambda(grid, point, neighbor))
+ {
+ continue;
+ }
+ double logWeight = biasedLogWeightFromPoint(dimParams, points_, grid, neighbor,
+ points_[neighbor].bias(), coordValue, {}, point);
+ weightSum += std::exp(logWeight);
+ }
+
+ /* Returns -GMX_FLOAT_MAX if no neighboring points were in the target region. */
+ return (weightSum > 0) ? std::log(weightSum) : -GMX_FLOAT_MAX;
+}
+
+void BiasState::sampleProbabilityWeights(const BiasGrid& grid, gmx::ArrayRef<const double> probWeightNeighbor)
+{
+ const std::vector<int>& neighbor = grid.point(coordState_.gridpointIndex()).neighbor;
+
+ /* Save weights for next update */
+ for (size_t n = 0; n < neighbor.size(); n++)
+ {
+ points_[neighbor[n]].increaseWeightSumIteration(probWeightNeighbor[n]);
+ }
+
+ /* Update the local update range. Two corner points define this rectangular
+ * domain. We need to choose two new corner points such that the new domain
+ * contains both the old update range and the current neighborhood.
+ * In the simplest case when an update is performed every sample,
+ * the update range would simply equal the current neighborhood.
+ */
+ int neighborStart = neighbor[0];
+ int neighborLast = neighbor[neighbor.size() - 1];
+ for (int d = 0; d < grid.numDimensions(); d++)
+ {
+ int origin = grid.point(neighborStart).index[d];
+ int last = grid.point(neighborLast).index[d];
+
+ if (origin > last)
+ {
+ /* Unwrap if wrapped around the boundary (only happens for periodic
+ * boundaries). This has been already for the stored index interval.
+ */
+ /* TODO: what we want to do is to find the smallest the update
+ * interval that contains all points that need to be updated.
+ * This amounts to combining two intervals, the current
+ * [origin, end] update interval and the new touched neighborhood
+ * into a new interval that contains all points from both the old
+ * intervals.
+ *
+ * For periodic boundaries it becomes slightly more complicated
+ * than for closed boundaries because then it needs not be
+ * true that origin < end (so one can't simply relate the origin/end
+ * in the min()/max() below). The strategy here is to choose the
+ * origin closest to a reference point (index 0) and then unwrap
+ * the end index if needed and choose the largest end index.
+ * This ensures that both intervals are in the new interval
+ * but it's not necessarily the smallest.
+ * Currently we solve this by going through each possibility
+ * and checking them.
+ */
+ last += grid.axis(d).numPointsInPeriod();
+ }
+
+ originUpdatelist_[d] = std::min(originUpdatelist_[d], origin);
+ endUpdatelist_[d] = std::max(endUpdatelist_[d], last);
+ }
+}
+
+void BiasState::sampleCoordAndPmf(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ gmx::ArrayRef<const double> probWeightNeighbor,
+ double convolvedBias)
+{
+ /* Sampling-based deconvolution extracting the PMF.
+ * Update the PMF histogram with the current coordinate value.
+ *
+ * Because of the finite width of the harmonic potential, the free energy
+ * defined for each coordinate point does not exactly equal that of the
+ * actual coordinate, the PMF. However, the PMF can be estimated by applying
+ * the relation exp(-PMF) = exp(-bias_convolved)*P_biased/Z, i.e. by keeping a
+ * reweighted histogram of the coordinate value. Strictly, this relies on
+ * the unknown normalization constant Z being either constant or known. Here,
+ * neither is true except in the long simulation time limit. Empirically however,
+ * it works (mainly because how the PMF histogram is rescaled).
+ */
+
+ const int gridPointIndex = coordState_.gridpointIndex();
+ const std::optional<int> lambdaAxisIndex = grid.lambdaAxisIndex();
+
+ /* Update the PMF of points along a lambda axis with their bias. */
+ if (lambdaAxisIndex)
+ {
+ const std::vector<int>& neighbors = grid.point(gridPointIndex).neighbor;
+
+ std::vector<double> lambdaMarginalDistribution =
+ calculateFELambdaMarginalDistribution(grid, neighbors, probWeightNeighbor);
+
+ awh_dvec coordValueAlongLambda = { coordState_.coordValue()[0], coordState_.coordValue()[1],
+ coordState_.coordValue()[2], coordState_.coordValue()[3] };
+ for (size_t i = 0; i < neighbors.size(); i++)
+ {
+ const int neighbor = neighbors[i];
+ double bias;
+ if (pointsAlongLambdaAxis(grid, gridPointIndex, neighbor))
+ {
+ const double neighborLambda = grid.point(neighbor).coordValue[lambdaAxisIndex.value()];
+ if (neighbor == gridPointIndex)
+ {
+ bias = convolvedBias;
+ }
+ else
+ {
+ coordValueAlongLambda[lambdaAxisIndex.value()] = neighborLambda;
+ bias = calcConvolvedBias(dimParams, grid, coordValueAlongLambda);
+ }
+
+ const double probWeight = lambdaMarginalDistribution[neighborLambda];
+ const double weightedBias = bias - std::log(std::max(probWeight, GMX_DOUBLE_MIN)); // avoid log(0)
+
+ if (neighbor == gridPointIndex && grid.covers(coordState_.coordValue()))
+ {
+ points_[neighbor].samplePmf(weightedBias);
+ }
+ else
+ {
+ points_[neighbor].updatePmfUnvisited(weightedBias);
+ }
+ }
+ }
+ }
+ else
+ {
+ /* Only save coordinate data that is in range (the given index is always
+ * in range even if the coordinate value is not).
+ */
+ if (grid.covers(coordState_.coordValue()))
+ {
+ /* Save PMF sum and keep a histogram of the sampled coordinate values */
+ points_[gridPointIndex].samplePmf(convolvedBias);
+ }
+ }
+
+ /* Save probability weights for the update */
+ sampleProbabilityWeights(grid, probWeightNeighbor);
+}
+
+void BiasState::initHistoryFromState(AwhBiasHistory* biasHistory) const
+{
+ biasHistory->pointState.resize(points_.size());
+}
+
+void BiasState::updateHistory(AwhBiasHistory* biasHistory, const BiasGrid& grid) const
+{
+ GMX_RELEASE_ASSERT(biasHistory->pointState.size() == points_.size(),
+ "The AWH history setup does not match the AWH state.");
+
+ AwhBiasStateHistory* stateHistory = &biasHistory->state;
+ stateHistory->umbrellaGridpoint = coordState_.umbrellaGridpoint();
+
+ for (size_t m = 0; m < biasHistory->pointState.size(); m++)
+ {
+ AwhPointStateHistory* psh = &biasHistory->pointState[m];
+
+ points_[m].storeState(psh);
+
+ psh->weightsum_covering = weightSumCovering_[m];
+ }
+
+ histogramSize_.storeState(stateHistory);
+
+ stateHistory->origin_index_updatelist = multiDimGridIndexToLinear(grid, originUpdatelist_);
+ stateHistory->end_index_updatelist = multiDimGridIndexToLinear(grid, endUpdatelist_);
+}
+
+void BiasState::restoreFromHistory(const AwhBiasHistory& biasHistory, const BiasGrid& grid)
+{
+ const AwhBiasStateHistory& stateHistory = biasHistory.state;
+
+ coordState_.restoreFromHistory(stateHistory);
+
+ if (biasHistory.pointState.size() != points_.size())
+ {
+ GMX_THROW(
+ InvalidInputError("Bias grid size in checkpoint and simulation do not match. "
+ "Likely you provided a checkpoint from a different simulation."));
+ }
+ for (size_t m = 0; m < points_.size(); m++)
+ {
+ points_[m].setFromHistory(biasHistory.pointState[m]);
+ }
+
+ for (size_t m = 0; m < weightSumCovering_.size(); m++)
+ {
+ weightSumCovering_[m] = biasHistory.pointState[m].weightsum_covering;
+ }
+
+ histogramSize_.restoreFromHistory(stateHistory);
+
+ linearGridindexToMultiDim(grid, stateHistory.origin_index_updatelist, originUpdatelist_);
+ linearGridindexToMultiDim(grid, stateHistory.end_index_updatelist, endUpdatelist_);
+}
+
+void BiasState::broadcast(const t_commrec* commRecord)
+{
+ gmx_bcast(sizeof(coordState_), &coordState_, commRecord->mpi_comm_mygroup);
+
+ gmx_bcast(points_.size() * sizeof(PointState), points_.data(), commRecord->mpi_comm_mygroup);
+
+ gmx_bcast(weightSumCovering_.size() * sizeof(double), weightSumCovering_.data(),
+ commRecord->mpi_comm_mygroup);
+
+ gmx_bcast(sizeof(histogramSize_), &histogramSize_, commRecord->mpi_comm_mygroup);
+}
+
+void BiasState::setFreeEnergyToConvolvedPmf(const std::vector<DimParams>& dimParams, const BiasGrid& grid)
+{
+ std::vector<float> convolvedPmf;
+
+ calcConvolvedPmf(dimParams, grid, &convolvedPmf);
+
+ for (size_t m = 0; m < points_.size(); m++)
+ {
+ points_[m].setFreeEnergy(convolvedPmf[m]);
+ }
+}
+
+/*! \brief
+ * Count trailing data rows containing only zeros.
+ *
+ * \param[in] data 2D data array.
+ * \param[in] numRows Number of rows in array.
+ * \param[in] numColumns Number of cols in array.
+ * \returns the number of trailing zero rows.
+ */
+static int countTrailingZeroRows(const double* const* data, int numRows, int numColumns)
+{
+ int numZeroRows = 0;
+ for (int m = numRows - 1; m >= 0; m--)
+ {
+ bool rowIsZero = true;
+ for (int d = 0; d < numColumns; d++)
+ {
+ if (data[d][m] != 0)
+ {
+ rowIsZero = false;
+ break;
+ }
+ }
+
+ if (!rowIsZero)
+ {
+ /* At a row with non-zero data */
+ break;
+ }
+ else
+ {
+ /* Still at a zero data row, keep checking rows higher up. */
+ numZeroRows++;
+ }
+ }
+
+ return numZeroRows;
+}
+
+/*! \brief
+ * Initializes the PMF and target with data read from an input table.
+ *
+ * \param[in] dimParams The dimension parameters.
+ * \param[in] grid The grid.
+ * \param[in] filename The filename to read PMF and target from.
+ * \param[in] numBias Number of biases.
+ * \param[in] biasIndex The index of the bias.
+ * \param[in,out] pointState The state of the points in this bias.
+ */
+static void readUserPmfAndTargetDistribution(const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ const std::string& filename,
+ int numBias,
+ int biasIndex,
+ std::vector<PointState>* pointState)
+{
+ /* Read the PMF and target distribution.
+ From the PMF, the convolved PMF, or the reference value free energy, can be calculated
+ base on the force constant. The free energy and target together determine the bias.
+ */
+ std::string filenameModified(filename);
+ if (numBias > 1)
+ {
+ size_t n = filenameModified.rfind('.');
+ GMX_RELEASE_ASSERT(n != std::string::npos,
+ "The filename should contain an extension starting with .");
+ filenameModified.insert(n, formatString("%d", biasIndex));
+ }
+
+ std::string correctFormatMessage = formatString(
+ "%s is expected in the following format. "
+ "The first ndim column(s) should contain the coordinate values for each point, "
+ "each column containing values of one dimension (in ascending order). "
+ "For a multidimensional coordinate, points should be listed "
+ "in the order obtained by traversing lower dimensions first. "
+ "E.g. for two-dimensional grid of size nxn: "
+ "(1, 1), (1, 2),..., (1, n), (2, 1), (2, 2), ..., , (n, n - 1), (n, n). "
+ "Column ndim + 1 should contain the PMF value for each coordinate value. "
+ "The target distribution values should be in column ndim + 2 or column ndim + 5. "
+ "Make sure the input file ends with a new line but has no trailing new lines.",
+ filename.c_str());
+ gmx::TextLineWrapper wrapper;
+ wrapper.settings().setLineLength(c_linewidth);
+ correctFormatMessage = wrapper.wrapToString(correctFormatMessage);
+
+ double** data;
+ int numColumns;
+ int numRows = read_xvg(filenameModified.c_str(), &data, &numColumns);
+
+ /* Check basic data properties here. BiasGrid takes care of more complicated things. */
+
+ if (numRows <= 0)
+ {
+ std::string mesg = gmx::formatString("%s is empty!.\n\n%s", filename.c_str(),
+ correctFormatMessage.c_str());
+ GMX_THROW(InvalidInputError(mesg));
+ }
+
+ /* Less than 2 points is not useful for PMF or target. */
+ if (numRows < 2)
+ {
+ std::string mesg = gmx::formatString(
+ "%s contains too few data points (%d)."
+ "The minimum number of points is 2.",
+ filename.c_str(), numRows);
+ GMX_THROW(InvalidInputError(mesg));
+ }
+
+ /* Make sure there are enough columns of data.
+
+ Two formats are allowed. Either with columns {coords, PMF, target} or
+ {coords, PMF, x, y, z, target, ...}. The latter format is allowed since that
+ is how AWH output is written (x, y, z being other AWH variables). For this format,
+ trailing columns are ignored.
+ */
+ int columnIndexTarget;
+ int numColumnsMin = dimParams.size() + 2;
+ int columnIndexPmf = dimParams.size();
+ if (numColumns == numColumnsMin)
+ {
+ columnIndexTarget = columnIndexPmf + 1;
+ }
+ else
+ {
+ columnIndexTarget = columnIndexPmf + 4;
+ }
+
+ if (numColumns < numColumnsMin)
+ {
+ std::string mesg = gmx::formatString(
+ "The number of columns in %s should be at least %d."
+ "\n\n%s",
+ filename.c_str(), numColumnsMin, correctFormatMessage.c_str());
+ GMX_THROW(InvalidInputError(mesg));
+ }
+
+ /* read_xvg can give trailing zero data rows for trailing new lines in the input. We allow 1 zero row,
+ since this could be real data. But multiple trailing zero rows cannot correspond to valid data. */
+ int numZeroRows = countTrailingZeroRows(data, numRows, numColumns);
+ if (numZeroRows > 1)
+ {
+ std::string mesg = gmx::formatString(
+ "Found %d trailing zero data rows in %s. Please remove trailing empty lines and "
+ "try again.",
+ numZeroRows, filename.c_str());
+ GMX_THROW(InvalidInputError(mesg));
+ }
+
+ /* Convert from user units to internal units before sending the data of to grid. */
+ for (size_t d = 0; d < dimParams.size(); d++)
+ {
+ double scalingFactor = dimParams[d].scaleUserInputToInternal(1);
+ if (scalingFactor == 1)
+ {
+ continue;
+ }
+ for (size_t m = 0; m < pointState->size(); m++)
+ {
+ data[d][m] *= scalingFactor;
+ }
+ }
+
+ /* Get a data point for each AWH grid point so that they all get data. */
+ std::vector<int> gridIndexToDataIndex(grid.numPoints());
+ mapGridToDataGrid(&gridIndexToDataIndex, data, numRows, filename, grid, correctFormatMessage);
+
+ /* Extract the data for each grid point.
+ * We check if the target distribution is zero for all points.
+ */
+ bool targetDistributionIsZero = true;
+ for (size_t m = 0; m < pointState->size(); m++)
+ {
+ (*pointState)[m].setLogPmfSum(-data[columnIndexPmf][gridIndexToDataIndex[m]]);
+ double target = data[columnIndexTarget][gridIndexToDataIndex[m]];
+
+ /* Check if the values are allowed. */
+ if (target < 0)
+ {
+ std::string mesg = gmx::formatString(
+ "Target distribution weight at point %zu (%g) in %s is negative.", m, target,
+ filename.c_str());
+ GMX_THROW(InvalidInputError(mesg));
+ }
+ if (target > 0)
+ {
+ targetDistributionIsZero = false;
+ }
+ (*pointState)[m].setTargetConstantWeight(target);
+ }
+
+ if (targetDistributionIsZero)
+ {
+ std::string mesg =
+ gmx::formatString("The target weights given in column %d in %s are all 0",
+ columnIndexTarget, filename.c_str());
+ GMX_THROW(InvalidInputError(mesg));
+ }
+
+ /* Free the arrays. */
+ for (int m = 0; m < numColumns; m++)
+ {
+ sfree(data[m]);
+ }
+ sfree(data);
+}
+
+void BiasState::normalizePmf(int numSharingSims)
+{
+ /* The normalization of the PMF estimate matters because it determines how big effect the next sample has.
+ Approximately (for large enough force constant) we should have:
+ sum_x(exp(-pmf(x)) = nsamples*sum_xref(exp(-f(xref)).
+ */
+
+ /* Calculate the normalization factor, i.e. divide by the pmf sum, multiply by the number of samples and the f sum */
+ double expSumPmf = 0;
+ double expSumF = 0;
+ for (const PointState& pointState : points_)
+ {
+ if (pointState.inTargetRegion())
+ {
+ expSumPmf += std::exp(pointState.logPmfSum());
+ expSumF += std::exp(-pointState.freeEnergy());
+ }
+ }
+ double numSamples = histogramSize_.histogramSize() / numSharingSims;
+
+ /* Renormalize */
+ double logRenorm = std::log(numSamples * expSumF / expSumPmf);
+ for (PointState& pointState : points_)
+ {
+ if (pointState.inTargetRegion())
+ {
+ pointState.setLogPmfSum(pointState.logPmfSum() + logRenorm);
+ }
+ }
+}
+
+void BiasState::initGridPointState(const AwhBiasParams& awhBiasParams,
+ const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid,
+ const BiasParams& params,
+ const std::string& filename,
+ int numBias)
+{
+ /* Modify PMF, free energy and the constant target distribution factor
+ * to user input values if there is data given.
+ */
+ if (awhBiasParams.bUserData)
+ {
+ readUserPmfAndTargetDistribution(dimParams, grid, filename, numBias, params.biasIndex, &points_);
+ setFreeEnergyToConvolvedPmf(dimParams, grid);
+ }
+
+ /* The local Boltzmann distribution is special because the target distribution is updated as a function of the reference weighthistogram. */
+ GMX_RELEASE_ASSERT(params.eTarget != eawhtargetLOCALBOLTZMANN || points_[0].weightSumRef() != 0,
+ "AWH reference weight histogram not initialized properly with local "
+ "Boltzmann target distribution.");
+
+ updateTargetDistribution(points_, params);
+
+ for (PointState& pointState : points_)
+ {
+ if (pointState.inTargetRegion())
+ {
+ pointState.updateBias();
+ }
+ else
+ {
+ /* Note that for zero target this is a value that represents -infinity but should not be used for biasing. */
+ pointState.setTargetToZero();
+ }
+ }
+
+ /* Set the initial reference weighthistogram. */
+ const double histogramSize = histogramSize_.histogramSize();
+ for (auto& pointState : points_)
+ {
+ pointState.setInitialReferenceWeightHistogram(histogramSize);
+ }
+
+ /* Make sure the pmf is normalized consistently with the histogram size.
+ Note: the target distribution and free energy need to be set here. */
+ normalizePmf(params.numSharedUpdate);
+}
+
+BiasState::BiasState(const AwhBiasParams& awhBiasParams,
+ double histogramSizeInitial,
+ const std::vector<DimParams>& dimParams,
+ const BiasGrid& grid) :
+ coordState_(awhBiasParams, dimParams, grid),
+ points_(grid.numPoints()),
+ weightSumCovering_(grid.numPoints()),
+ histogramSize_(awhBiasParams, histogramSizeInitial)
+{
+ /* The minimum and maximum multidimensional point indices that are affected by the next update */
+ for (size_t d = 0; d < dimParams.size(); d++)
+ {
+ int index = grid.point(coordState_.gridpointIndex()).index[d];
+ originUpdatelist_[d] = index;
+ endUpdatelist_[d] = index;
+ }
+}
+
+} // namespace gmx
/*
* This file is part of the GROMACS molecular simulation package.
*
-- * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
++ * Copyright (c) 2015,2016,2017,2018,2019, The GROMACS development team.
++ * Copyright (c) 2020,2021, 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.
void calculateForces(const ForceProviderInput& forceProviderInput,
ForceProviderOutput* forceProviderOutput) override;
+ void subscribeToSimulationSetupNotifications(MdModulesNotifier* /* notifier */) override {}
+ void subscribeToPreProcessingNotifications(MdModulesNotifier* /* notifier */) override {}
+
private:
//! Return whether or not to apply a field
bool isActive() const;
if (fieldStrength != 0)
{
// TODO: Check parallellism
- for (index i = 0; i != ssize(f); ++i)
+ for (int i = 0; i < mdatoms.homenr; ++i)
{
// NOTE: Not correct with perturbed charges
f[i][m] += mdatoms.chargeA[i] * fieldStrength;
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013-2020, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017, The GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
++ * Copyright (c) 2018,2019,2020,2021, 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.
#include <cstdlib>
#include <algorithm>
+#include <memory>
#include <string>
-#include "gromacs/awh/read_params.h"
+#include "gromacs/applied_forces/awh/read_params.h"
#include "gromacs/fileio/readinp.h"
#include "gromacs/fileio/warninp.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/mdrun/mdmodules.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/multipletimestepping.h"
#include "gromacs/mdtypes/pull_params.h"
#include "gromacs/options/options.h"
#include "gromacs/options/treesupport.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/symtab.h"
#include "gromacs/topology/topology.h"
+ #include "gromacs/utility/arrayref.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
* message.
*/
-typedef struct t_inputrec_strings
+struct gmx_inputrec_strings
{
char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], acc[STRLEN], accgrps[STRLEN], freeze[STRLEN],
frdim[STRLEN], energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN],
x_compressed_groups[STRLEN], couple_moltype[STRLEN], orirefitgrp[STRLEN],
egptable[STRLEN], egpexcl[STRLEN], wall_atomtype[STRLEN], wall_density[STRLEN],
deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
- char fep_lambda[efptNR][STRLEN];
- char lambda_weights[STRLEN];
- char** pull_grp;
- char** rot_grp;
- char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
- char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN], bSH[STRLEN],
- CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN], SAoff[STRLEN], SAsteps[STRLEN];
-
-} gmx_inputrec_strings;
+ char fep_lambda[efptNR][STRLEN];
+ char lambda_weights[STRLEN];
+ std::vector<std::string> pullGroupNames;
+ std::vector<std::string> rotateGroupNames;
+ char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
+};
-static gmx_inputrec_strings* is = nullptr;
+static gmx_inputrec_strings* inputrecStrings = nullptr;
void init_inputrec_strings()
{
- if (is)
+ if (inputrecStrings)
{
gmx_incons(
"Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
"Only one inputrec (i.e. .mdp file) can be parsed at a time.");
}
- snew(is, 1);
+ inputrecStrings = new gmx_inputrec_strings();
}
void done_inputrec_strings()
{
- sfree(is);
- is = nullptr;
+ delete inputrecStrings;
+ inputrecStrings = nullptr;
}
}
}
+static void checkMtsRequirement(const t_inputrec& ir, const char* param, const int nstValue, warninp_t wi)
+{
+ GMX_RELEASE_ASSERT(ir.mtsLevels.size() >= 2, "Need at least two levels for MTS");
+ const int mtsFactor = ir.mtsLevels.back().stepFactor;
+ if (nstValue % mtsFactor != 0)
+ {
+ auto message = gmx::formatString(
+ "With MTS, %s = %d should be a multiple of mts-factor = %d", param, nstValue, mtsFactor);
+ warning_error(wi, message.c_str());
+ }
+}
+
+static void setupMtsLevels(gmx::ArrayRef<gmx::MtsLevel> mtsLevels,
+ const t_inputrec& ir,
+ const t_gromppopts& opts,
+ warninp_t wi)
+{
+ /* MD-VV has no MTS support yet.
+ * SD1 needs different scaling coefficients for the different MTS forces
+ * and the different forces are currently not available in ForceBuffers.
+ */
+ if (ir.eI != eiMD)
+ {
+ auto message = gmx::formatString(
+ "Multiple time stepping is only supported with integrator %s", ei_names[eiMD]);
+ warning_error(wi, message.c_str());
+ }
+ if (opts.numMtsLevels != 2)
+ {
+ warning_error(wi, "Only mts-levels = 2 is supported");
+ }
+ else
+ {
+ const std::vector<std::string> inputForceGroups = gmx::splitString(opts.mtsLevel2Forces);
+ auto& forceGroups = mtsLevels[1].forceGroups;
+ for (const auto& inputForceGroup : inputForceGroups)
+ {
+ bool found = false;
+ int nameIndex = 0;
+ for (const auto& forceGroupName : gmx::mtsForceGroupNames)
+ {
+ if (gmx::equalCaseInsensitive(inputForceGroup, forceGroupName))
+ {
+ forceGroups.set(nameIndex);
+ found = true;
+ }
+ nameIndex++;
+ }
+ if (!found)
+ {
+ auto message =
+ gmx::formatString("Unknown MTS force group '%s'", inputForceGroup.c_str());
+ warning_error(wi, message.c_str());
+ }
+ }
+
+ if (mtsLevels[1].stepFactor <= 1)
+ {
+ gmx_fatal(FARGS, "mts-factor should be larger than 1");
+ }
+
+ // Make the level 0 use the complement of the force groups of group 1
+ mtsLevels[0].forceGroups = ~mtsLevels[1].forceGroups;
+ mtsLevels[0].stepFactor = 1;
+
+ if ((EEL_FULL(ir.coulombtype) || EVDW_PME(ir.vdwtype))
+ && !mtsLevels[1].forceGroups[static_cast<int>(gmx::MtsForceGroups::LongrangeNonbonded)])
+ {
+ warning_error(wi,
+ "With long-range electrostatics and/or LJ treatment, the long-range part "
+ "has to be part of the mts-level2-forces");
+ }
+
+ if (ir.nstcalcenergy > 0)
+ {
+ checkMtsRequirement(ir, "nstcalcenergy", ir.nstcalcenergy, wi);
+ }
+ checkMtsRequirement(ir, "nstenergy", ir.nstenergy, wi);
+ checkMtsRequirement(ir, "nstlog", ir.nstlog, wi);
+ if (ir.efep != efepNO)
+ {
+ checkMtsRequirement(ir, "nstdhdl", ir.fepvals->nstdhdl, wi);
+ }
+
+ if (ir.bPull)
+ {
+ const int pullMtsLevel = gmx::forceGroupMtsLevel(ir.mtsLevels, gmx::MtsForceGroups::Pull);
+ if (ir.pull->nstxout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
+ {
+ warning_error(wi, "pull-nstxout should be a multiple of mts-factor");
+ }
+ if (ir.pull->nstfout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
+ {
+ warning_error(wi, "pull-nstfout should be a multiple of mts-factor");
+ }
+ }
+ }
+}
+
void check_ir(const char* mdparin,
const gmx::MdModulesNotifier& mdModulesNotifier,
t_inputrec* ir,
{
ir->nstpcouple = ir_optimal_nstpcouple(ir);
}
+ if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
+ {
+ warning_error(wi,
+ "With multiple time stepping, nstpcouple should be a mutiple of "
+ "mts-factor");
+ }
}
if (ir->nstcalcenergy > 0)
// Inquire all MdModules, if their parameters match with the energy
// calculation frequency
gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
- mdModulesNotifier.notifier_.notify(&energyCalculationFrequencyErrors);
+ mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
// Emit all errors from the energy calculation frequency checks
for (const std::string& energyFrequencyErrorMessage :
/* TPI STUFF */
if (EI_TPI(ir->eI))
{
- sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
- CHECK(ir->ePBC != epbcXYZ);
+ sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
+ CHECK(ir->pbcType != PbcType::Xyz);
sprintf(err_buf, "with TPI nstlist should be larger than zero");
CHECK(ir->nstlist <= 0);
sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
- sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
+ sprintf(err_buf,
+ "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
+ "supported.)",
static_cast<int>(fep->sc_r_power));
- CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
+ CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
sprintf(err_buf,
"Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
}
/* PBC/WALLS */
- sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
- CHECK(ir->nwall && ir->ePBC != epbcXY);
+ sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
+ CHECK(ir->nwall && ir->pbcType != PbcType::XY);
/* VACUUM STUFF */
- if (ir->ePBC != epbcXYZ && ir->nwall != 2)
+ if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
{
- if (ir->ePBC == epbcNONE)
+ if (ir->pbcType == PbcType::No)
{
if (ir->epc != epcNO)
{
}
else
{
- sprintf(err_buf, "Can not have pressure coupling with pbc=%s", epbc_names[ir->ePBC]);
+ sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
+ c_pbcTypeNames[ir->pbcType].c_str());
CHECK(ir->epc != epcNO);
}
- sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
+ sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
CHECK(EEL_FULL(ir->coulombtype));
- sprintf(err_buf, "Can not have dispersion correction with pbc=%s", epbc_names[ir->ePBC]);
+ sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
+ c_pbcTypeNames[ir->pbcType].c_str());
CHECK(ir->eDispCorr != edispcNO);
}
"with coulombtype = %s or coulombtype = %s\n"
"without periodic boundary conditions (pbc = %s) and\n"
"rcoulomb and rvdw set to zero",
- eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
+ eel_names[eelCUT], eel_names[eelUSER], c_pbcTypeNames[PbcType::No].c_str());
CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
- || (ir->ePBC != epbcNONE) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
+ || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
if (ir->nstlist > 0)
{
"Can not remove the rotation around the center of mass with periodic "
"molecules");
CHECK(ir->bPeriodicMols);
- if (ir->ePBC != epbcNONE)
+ if (ir->pbcType != PbcType::No)
{
warning(wi,
"Removing the rotation around the center of mass in a periodic system, "
}
}
- if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
+ if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
{
sprintf(warn_buf,
"Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
{
if (ir->ewald_geometry == eewg3D)
{
- sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s", epbc_names[ir->ePBC],
- eewg_names[eewg3DC]);
+ sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
+ c_pbcTypeNames[ir->pbcType].c_str(), eewg_names[eewg3DC]);
warning(wi, warn_buf);
}
/* This check avoids extra pbc coding for exclusion corrections */
sprintf(err_buf, "wall-ewald-zfac should be >= 2");
CHECK(ir->wall_ewald_zfac < 2);
}
- if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) && EEL_FULL(ir->coulombtype))
+ if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
{
sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
- eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
+ eel_names[ir->coulombtype], eewg_names[eewg3DC], c_pbcTypeNames[PbcType::XY].c_str());
warning(wi, warn_buf);
}
if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
if (ir->bQMMM)
{
- warning_error(wi, "QMMM is currently not supported");
- if (!EI_DYNAMICS(ir->eI))
- {
- char buf[STRLEN];
- sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
- warning_error(wi, buf);
- }
+ warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
}
if (ir->bAdress)
}
- /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
- if (fep->sc_r_power == 48)
- {
- if (fep->sc_alpha > 0.1)
- {
- gmx_fatal(FARGS,
- "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004",
- fep->sc_alpha);
- }
- }
-
/* now read in the weights */
parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
if (nweights == 0)
GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
}
-static void convertYesNos(warninp_t /*wi*/,
- gmx::ArrayRef<const std::string> inputs,
- const char* /*name*/,
- gmx_bool* outputs)
-{
- int i = 0;
- for (const auto& input : inputs)
- {
- outputs[i] = gmx::equalCaseInsensitive(input, "Y", 1);
- ++i;
- }
-}
-
template<typename T>
void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
{
gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu", ir->nwall,
wallAtomTypes.size());
}
+ GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
for (int i = 0; i < ir->nwall; i++)
{
opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
printStringNoNewline(
&inp, "Part index is updated automatically on checkpointing (keeps files separate)");
ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
+ printStringNoNewline(&inp, "Multiple time-stepping");
+ ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
+ if (ir->useMts)
+ {
+ opts->numMtsLevels = get_eint(&inp, "mts-levels", 2, wi);
+ ir->mtsLevels.resize(2);
+ gmx::MtsLevel& mtsLevel = ir->mtsLevels[1];
+ opts->mtsLevel2Forces = setStringEntry(&inp, "mts-level2-forces",
+ "longrange-nonbonded nonbonded pair dihedral");
+ mtsLevel.stepFactor = get_eint(&inp, "mts-level2-factor", 2, wi);
+
+ // We clear after reading without dynamics to not force the user to remove MTS mdp options
+ if (!EI_DYNAMICS(ir->eI))
+ {
+ ir->useMts = false;
+ ir->mtsLevels.clear();
+ }
+ }
printStringNoNewline(&inp, "mode for center of mass motion removal");
ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
printStringNoNewline(&inp, "number of steps for center of mass motion removal");
ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
printStringNoNewline(&inp, "group(s) for center of mass motion removal");
- setStringEntry(&inp, "comm-grps", is->vcm, nullptr);
+ setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
printStringNoNewline(&inp, "default, all atoms will be written.");
- setStringEntry(&inp, "compressed-x-grps", is->x_compressed_groups, nullptr);
+ setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
printStringNoNewline(&inp, "Selection of energy groups");
- setStringEntry(&inp, "energygrps", is->energy, nullptr);
+ setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
/* Neighbor searching */
printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
printStringNoNewline(&inp, "nblist update frequency");
ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
- ir->ePBC = get_eeenum(&inp, "pbc", epbc_names, wi);
+ // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
+ std::vector<const char*> pbcTypesNamesChar;
+ for (const auto& pbcTypeName : c_pbcTypeNames)
+ {
+ pbcTypesNamesChar.push_back(pbcTypeName.c_str());
+ }
+ ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
printStringNoNewline(&inp,
"Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
printStringNoNewline(&inp, "Separate tables between energy group pairs");
- setStringEntry(&inp, "energygrp-table", is->egptable, nullptr);
+ setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
printStringNoNewline(&inp, "Groups to couple separately");
- setStringEntry(&inp, "tc-grps", is->tcgrps, nullptr);
+ setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
- setStringEntry(&inp, "tau-t", is->tau_t, nullptr);
- setStringEntry(&inp, "ref-t", is->ref_t, nullptr);
+ setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
+ setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
printStringNoNewline(&inp, "pressure coupling");
ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
/* QMMM */
printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
- printStringNoNewline(&inp, "Groups treated Quantum Mechanically");
- setStringEntry(&inp, "QMMM-grps", is->QMMM, nullptr);
- printStringNoNewline(&inp, "QM method");
- setStringEntry(&inp, "QMmethod", is->QMmethod, nullptr);
- printStringNoNewline(&inp, "QMMM scheme");
- ir->QMMMscheme = get_eeenum(&inp, "QMMMscheme", eQMMMscheme_names, wi);
- printStringNoNewline(&inp, "QM basisset");
- setStringEntry(&inp, "QMbasis", is->QMbasis, nullptr);
- printStringNoNewline(&inp, "QM charge");
- setStringEntry(&inp, "QMcharge", is->QMcharge, nullptr);
- printStringNoNewline(&inp, "QM multiplicity");
- setStringEntry(&inp, "QMmult", is->QMmult, nullptr);
- printStringNoNewline(&inp, "Surface Hopping");
- setStringEntry(&inp, "SH", is->bSH, nullptr);
- printStringNoNewline(&inp, "CAS space options");
- setStringEntry(&inp, "CASorbitals", is->CASorbitals, nullptr);
- setStringEntry(&inp, "CASelectrons", is->CASelectrons, nullptr);
- setStringEntry(&inp, "SAon", is->SAon, nullptr);
- setStringEntry(&inp, "SAoff", is->SAoff, nullptr);
- setStringEntry(&inp, "SAsteps", is->SAsteps, nullptr);
- printStringNoNewline(&inp, "Scale factor for MM charges");
- ir->scalefactor = get_ereal(&inp, "MMChargeScaleFactor", 1.0, wi);
+ printStringNoNewline(&inp, "Groups treated with MiMiC");
+ setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
/* Simulated annealing */
printStringNewline(&inp, "SIMULATED ANNEALING");
printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
- setStringEntry(&inp, "annealing", is->anneal, nullptr);
+ setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
printStringNoNewline(&inp,
"Number of time points to use for specifying annealing in each group");
- setStringEntry(&inp, "annealing-npoints", is->anneal_npoints, nullptr);
+ setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
printStringNoNewline(&inp, "List of times at the annealing points for each group");
- setStringEntry(&inp, "annealing-time", is->anneal_time, nullptr);
+ setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
- setStringEntry(&inp, "annealing-temp", is->anneal_temp, nullptr);
+ setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
/* Startup run */
printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
printStringNoNewline(
&inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
- setStringEntry(&inp, "energygrp-excl", is->egpexcl, nullptr);
+ setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
/* Walls */
printStringNewline(&inp, "WALLS");
ir->nwall = get_eint(&inp, "nwall", 0, wi);
ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
- setStringEntry(&inp, "wall-atomtype", is->wall_atomtype, nullptr);
- setStringEntry(&inp, "wall-density", is->wall_density, nullptr);
+ setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
+ setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
/* COM pulling */
ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
if (ir->bPull)
{
- snew(ir->pull, 1);
- is->pull_grp = read_pullparams(&inp, ir->pull, wi);
+ ir->pull = std::make_unique<pull_params_t>();
+ inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
+
+ if (ir->useMts)
+ {
+ for (int c = 0; c < ir->pull->ncoord; c++)
+ {
+ if (ir->pull->coord[c].eType == epullCONSTRAINT)
+ {
+ warning_error(wi,
+ "Constraint COM pulling is not supported in combination with "
+ "multiple time stepping");
+ break;
+ }
+ }
+ }
}
/* AWH biasing
- NOTE: needs COM pulling input */
+ NOTE: needs COM pulling or free energy input */
printStringNewline(&inp, "AWH biasing");
ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
if (ir->bDoAwh)
{
- if (ir->bPull)
- {
- ir->awhParams = gmx::readAndCheckAwhParams(&inp, ir, wi);
- }
- else
- {
- gmx_fatal(FARGS, "AWH biasing is only compatible with COM pulling turned on");
- }
+ ir->awhParams = gmx::readAwhParams(&inp, wi);
}
/* Enforced rotation */
if (ir->bRot)
{
snew(ir->rot, 1);
- is->rot_grp = read_rotparams(&inp, ir->rot, wi);
+ inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
}
/* Interactive MD */
ir->bIMD = FALSE;
printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
- setStringEntry(&inp, "IMD-group", is->imd_grp, nullptr);
- if (is->imd_grp[0] != '\0')
+ setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
+ if (inputrecStrings->imd_grp[0] != '\0')
{
snew(ir->imd, 1);
ir->bIMD = TRUE;
printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
- setStringEntry(&inp, "orire-fitgrp", is->orirefitgrp, nullptr);
+ setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
/* free energy variables */
printStringNewline(&inp, "Free energy variables");
ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
- setStringEntry(&inp, "couple-moltype", is->couple_moltype, nullptr);
+ setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
- setStringEntry(&inp, "fep-lambdas", is->fep_lambda[efptFEP], nullptr);
- setStringEntry(&inp, "mass-lambdas", is->fep_lambda[efptMASS], nullptr);
- setStringEntry(&inp, "coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
- setStringEntry(&inp, "vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
- setStringEntry(&inp, "bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
- setStringEntry(&inp, "restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
- setStringEntry(&inp, "temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
+ setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
+ setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
+ setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
+ setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
+ setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
+ setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
+ setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
- setStringEntry(&inp, "init-lambda-weights", is->lambda_weights, nullptr);
+ setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
/* Non-equilibrium MD stuff */
printStringNewline(&inp, "Non-equilibrium MD stuff");
- setStringEntry(&inp, "acc-grps", is->accgrps, nullptr);
- setStringEntry(&inp, "accelerate", is->acc, nullptr);
- setStringEntry(&inp, "freezegrps", is->freeze, nullptr);
- setStringEntry(&inp, "freezedim", is->frdim, nullptr);
+ setStringEntry(&inp, "acc-grps", inputrecStrings->accgrps, nullptr);
+ setStringEntry(&inp, "accelerate", inputrecStrings->acc, nullptr);
+ setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
+ setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
- setStringEntry(&inp, "deform", is->deform, nullptr);
+ setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
/* simulated tempering variables */
printStringNewline(&inp, "simulated tempering variables");
/* User defined thingies */
printStringNewline(&inp, "User defined thingies");
- setStringEntry(&inp, "user1-grps", is->user1, nullptr);
- setStringEntry(&inp, "user2-grps", is->user2, nullptr);
+ setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
+ setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
ir->userint1 = get_eint(&inp, "userint1", 0, wi);
ir->userint2 = get_eint(&inp, "userint2", 0, wi);
ir->userint3 = get_eint(&inp, "userint3", 0, wi);
}
opts->couple_moltype = nullptr;
- if (strlen(is->couple_moltype) > 0)
+ if (strlen(inputrecStrings->couple_moltype) > 0)
{
if (ir->efep != efepNO)
{
- opts->couple_moltype = gmx_strdup(is->couple_moltype);
+ opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
if (opts->couple_lam0 == opts->couple_lam1)
{
warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
/* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
if (ir->efep != efepNO)
{
- if (fep->delta_lambda > 0)
+ if (fep->delta_lambda != 0)
{
ir->efep = efepSLOWGROWTH;
}
{
ir->bExpanded = TRUE;
}
- do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
+ do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
if (ir->bSimTemp) /* done after fep params */
{
do_simtemp_params(ir);
/* WALL PARAMETERS */
- do_wall_params(ir, is->wall_atomtype, is->wall_density, opts, wi);
+ do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
/* ORIENTATION RESTRAINT PARAMETERS */
- if (opts->bOrire && gmx::splitString(is->orirefitgrp).size() != 1)
+ if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
{
warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
}
}
double gmx_unused canary;
- int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]), &(dumdub[0][1]),
- &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
+ int ndeform = sscanf(inputrecStrings->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]),
+ &(dumdub[0][1]), &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]),
+ &(dumdub[0][5]), &canary);
- if (strlen(is->deform) > 0 && ndeform != 6)
+ if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
{
- warning_error(
- wi, gmx::formatString(
- "Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform)
- .c_str());
+ warning_error(wi,
+ gmx::formatString(
+ "Cannot parse exactly 6 box deformation velocities from string '%s'",
+ inputrecStrings->deform)
+ .c_str());
}
for (i = 0; i < 3; i++)
{
}
}
- sfree(dumstr[0]);
- sfree(dumstr[1]);
-}
-
-static int search_QMstring(const char* s, int ng, const char* gn[])
-{
- /* same as normal search_string, but this one searches QM strings */
- int i;
+ /* Set up MTS levels, this needs to happen before checking AWH parameters */
+ if (ir->useMts)
+ {
+ setupMtsLevels(ir->mtsLevels, *ir, *opts, wi);
+ }
- for (i = 0; (i < ng); i++)
+ if (ir->bDoAwh)
{
- if (gmx_strcasecmp(s, gn[i]) == 0)
- {
- return i;
- }
+ gmx::checkAwhParams(ir->awhParams, ir, wi);
}
- gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
-} /* search_QMstring */
+ sfree(dumstr[0]);
+ sfree(dumstr[1]);
+}
/* 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
s);
}
+ static void atomGroupRangeValidation(int natoms, int groupIndex, const t_blocka& block)
+ {
+ /* Now go over the atoms in the group */
+ for (int j = block.index[groupIndex]; (j < block.index[groupIndex + 1]); j++)
+ {
+ int aj = block.a[j];
+
+ /* Range checking */
+ if ((aj < 0) || (aj >= natoms))
+ {
+ gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
+ }
+ }
+ }
+
static void do_numbering(int natoms,
SimulationGroups* groups,
gmx::ArrayRef<std::string> groupsFromMdpFile,
{
unsigned short* cbuf;
AtomGroupIndices* grps = &(groups->groups[gtype]);
- int j, gid, aj, ognr, ntot = 0;
+ int ntot = 0;
const char* title;
char warn_buf[STRLEN];
for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
{
/* Lookup the group name in the block structure */
- gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
+ const int gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
if ((grptp != egrptpONE) || (i == 0))
{
grps->emplace_back(gid);
}
-
+ GMX_ASSERT(block, "Can't have a nullptr block");
+ atomGroupRangeValidation(natoms, gid, *block);
/* Now go over the atoms in the group */
- for (j = block->index[gid]; (j < block->index[gid + 1]); j++)
+ for (int j = block->index[gid]; (j < block->index[gid + 1]); j++)
{
-
- aj = block->a[j];
-
- /* Range checking */
- if ((aj < 0) || (aj >= natoms))
- {
- gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
- }
+ const int aj = block->a[j];
/* Lookup up the old group number */
- ognr = cbuf[aj];
+ const int ognr = cbuf[aj];
if (ognr != NOGID)
{
gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
warning_note(wi, warn_buf);
}
/* Assign all atoms currently unassigned to a rest group */
- for (j = 0; (j < natoms); j++)
+ for (int j = 0; (j < natoms); j++)
{
if (cbuf[j] == NOGID)
{
grps->emplace_back(restnm);
/* Assign the rest name to all atoms not currently assigned to a group */
- for (j = 0; (j < natoms); j++)
+ for (int j = 0; (j < natoms); j++)
{
if (cbuf[j] == NOGID)
{
* belong to different TC or VCM groups it is anyhow difficult
* to determine the optimal nrdf assignment.
*/
- pull = ir->pull;
+ pull = ir->pull.get();
for (int i = 0; i < pull->ncoord; i++)
{
pgrp = &pull->group[pull->coord[i].group[j]];
- if (pgrp->nat > 0)
+ if (!pgrp->ind.empty())
{
/* Subtract 1/2 dof from each group */
int ai = pgrp->ind[0];
set_warning_line(wi, mdparin, -1);
- auto temperatureCouplingTauValues = gmx::splitString(is->tau_t);
- auto temperatureCouplingReferenceValues = gmx::splitString(is->ref_t);
- auto temperatureCouplingGroupNames = gmx::splitString(is->tcgrps);
+ auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
+ auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
+ auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
|| temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
{
}
/* Simulated annealing for each group. There are nr groups */
- auto simulatedAnnealingGroupNames = gmx::splitString(is->anneal);
+ auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
if (simulatedAnnealingGroupNames.size() == 1
&& gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
{
if (bAnneal)
{
/* Read the other fields too */
- auto simulatedAnnealingPoints = gmx::splitString(is->anneal_npoints);
+ auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
{
gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
}
- auto simulatedAnnealingTimes = gmx::splitString(is->anneal_time);
+ auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
{
gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
}
- auto simulatedAnnealingTemperatures = gmx::splitString(is->anneal_temp);
+ auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
{
gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
if (ir->bPull)
{
- const int gid = search_string(is->pull_grp[i], defaultIndexGroups->nr, gnames);
+ for (int i = 1; i < ir->pull->ngroup; i++)
+ {
- make_pull_groups(ir->pull, is->pull_grp, defaultIndexGroups, gnames);
++ const int gid = search_string(inputrecStrings->pullGroupNames[i].c_str(),
++ defaultIndexGroups->nr, gnames);
+ GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
+ atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
+ }
+
+ process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
- make_pull_coords(ir->pull);
+ checkPullCoords(ir->pull->group, ir->pull->coord);
}
if (ir->bRot)
{
- make_rotation_groups(ir->rot, is->rot_grp, defaultIndexGroups, gnames);
+ make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
}
if (ir->eSwapCoords != eswapNO)
/* Make indices for IMD session */
if (ir->bIMD)
{
- make_IMD_group(ir->imd, is->imd_grp, defaultIndexGroups, gnames);
+ make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
}
gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
*defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
- notifier.notifier_.notify(defaultIndexGroupsAndNames);
+ notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
- auto accelerations = gmx::splitString(is->acc);
- auto accelerationGroupNames = gmx::splitString(is->accgrps);
+ auto accelerations = gmx::splitString(inputrecStrings->acc);
+ auto accelerationGroupNames = gmx::splitString(inputrecStrings->accgrps);
if (accelerationGroupNames.size() * DIM != accelerations.size())
{
gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
- auto freezeDims = gmx::splitString(is->frdim);
- auto freezeGroupNames = gmx::splitString(is->freeze);
+ auto freezeDims = gmx::splitString(inputrecStrings->frdim);
+ auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
if (freezeDims.size() != DIM * freezeGroupNames.size())
{
gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
}
}
- auto energyGroupNames = gmx::splitString(is->energy);
+ auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::EnergyOutput, restnm, egrptpALL_GENREST, bVerbose, wi);
add_wall_energrps(groups, ir->nwall, symtab);
ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
- auto vcmGroupNames = gmx::splitString(is->vcm);
+ auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
/* Now we have filled the freeze struct, so we can calculate NRDF */
calc_nrdf(mtop, ir, gnames);
- auto user1GroupNames = gmx::splitString(is->user1);
+ auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::User1, restnm, egrptpALL_GENREST, bVerbose, wi);
- auto user2GroupNames = gmx::splitString(is->user2);
+ auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::User2, restnm, egrptpALL_GENREST, bVerbose, wi);
- auto compressedXGroupNames = gmx::splitString(is->x_compressed_groups);
+ auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::CompressedPositionOutput, restnm, egrptpONE, bVerbose, wi);
- auto orirefFitGroupNames = gmx::splitString(is->orirefitgrp);
+ auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
SimulationAtomGroupType::OrientationRestraintsFit, restnm, egrptpALL_GENREST,
bVerbose, wi);
- /* QMMM input processing */
- auto qmGroupNames = gmx::splitString(is->QMMM);
- auto qmMethods = gmx::splitString(is->QMmethod);
- auto qmBasisSets = gmx::splitString(is->QMbasis);
- if (ir->eI != eiMimic)
- {
- if (qmMethods.size() != qmGroupNames.size() || qmBasisSets.size() != qmGroupNames.size())
- {
- gmx_fatal(FARGS,
- "Invalid QMMM input: %zu groups %zu basissets"
- " and %zu methods\n",
- qmGroupNames.size(), qmBasisSets.size(), qmMethods.size());
- }
- /* group rest, if any, is always MM! */
- do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
- nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
- ir->opts.ngQM = qmGroupNames.size();
- snew(ir->opts.QMmethod, nr);
- snew(ir->opts.QMbasis, nr);
- for (i = 0; i < nr; i++)
- {
- /* input consists of strings: RHF CASSCF PM3 .. These need to be
- * converted to the corresponding enum in names.c
- */
- ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(), eQMmethodNR, eQMmethod_names);
- ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(), eQMbasisNR, eQMbasis_names);
- }
- auto qmMultiplicities = gmx::splitString(is->QMmult);
- auto qmCharges = gmx::splitString(is->QMcharge);
- auto qmbSH = gmx::splitString(is->bSH);
- snew(ir->opts.QMmult, nr);
- snew(ir->opts.QMcharge, nr);
- snew(ir->opts.bSH, nr);
- convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult);
- convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge);
- convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH);
-
- auto CASelectrons = gmx::splitString(is->CASelectrons);
- auto CASorbitals = gmx::splitString(is->CASorbitals);
- snew(ir->opts.CASelectrons, nr);
- snew(ir->opts.CASorbitals, nr);
- convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons);
- convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals);
-
- auto SAon = gmx::splitString(is->SAon);
- auto SAoff = gmx::splitString(is->SAoff);
- auto SAsteps = gmx::splitString(is->SAsteps);
- snew(ir->opts.SAon, nr);
- snew(ir->opts.SAoff, nr);
- snew(ir->opts.SAsteps, nr);
- convertInts(wi, SAon, "SAon", ir->opts.SAon);
- convertInts(wi, SAoff, "SAoff", ir->opts.SAoff);
- convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps);
- }
- else
+ /* MiMiC QMMM input processing */
+ auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
+ if (qmGroupNames.size() > 1)
{
- /* MiMiC */
- if (qmGroupNames.size() > 1)
- {
- gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
- }
- /* group rest, if any, is always MM! */
- do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
- SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
-
- ir->opts.ngQM = qmGroupNames.size();
+ gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
}
+ /* group rest, if any, is always MM! */
+ do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
+ SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
+ ir->opts.ngQM = qmGroupNames.size();
- /* end of QMMM input */
+ /* end of MiMiC QMMM input */
if (bVerbose)
{
nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
snew(ir->opts.egp_flags, nr * nr);
- bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
+ bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
if (bExcl && ir->cutoff_scheme == ecutsVERLET)
{
warning_error(wi, "Energy group exclusions are currently not supported");
warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
}
- bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
+ bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
&& !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
{
void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
{
+ // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
+ gmx::assertMtsRequirements(*ir);
+
char err_buf[STRLEN];
int i, m, c, nmol;
bool bCharge, bAcc;
char warn_buf[STRLEN];
const char* ptr;
- ptr = check_box(ir->ePBC, box);
+ ptr = check_box(ir->pbcType, box);
if (ptr)
{
warning_error(wi, ptr);
ir->LincsWarnAngle = 90.0;
}
- if (ir->ePBC != epbcNONE)
+ if (ir->pbcType != PbcType::No)
{
if (ir->nstlist == 0)
{
"With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
"atoms might cause the simulation to crash.");
}
- if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
+ if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
{
sprintf(warn_buf,
"ERROR: The cut-off length is longer than half the shortest box vector or "
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
++ * Copyright (c) 2018,2019,2020,2021, 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.
*/
#include "gmxpre.h"
+#include "coupling.h"
+
#include <cassert>
#include <cmath>
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
+#include "gromacs/mdlib/boxdeformation.h"
#include "gromacs/mdlib/expanded.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/stat.h"
static const double* sy_const[] = { nullptr, sy_const_1, nullptr, sy_const_3, nullptr, sy_const_5 };
+
+void update_tcouple(int64_t step,
+ const t_inputrec* inputrec,
+ t_state* state,
+ gmx_ekindata_t* ekind,
+ const t_extmass* MassQ,
+ const t_mdatoms* md)
+
+{
+ // This condition was explicitly checked in previous version, but should have never been satisfied
+ GMX_ASSERT(!(EI_VV(inputrec->eI)
+ && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec)
+ || inputrecNphTrotter(inputrec))),
+ "Temperature coupling was requested with velocity verlet and trotter");
+
+ bool doTemperatureCoupling = false;
+
+ // For VV temperature coupling parameters are updated on the current
+ // step, for the others - one step before.
+ if (inputrec->etc == etcNO)
+ {
+ doTemperatureCoupling = false;
+ }
+ else if (EI_VV(inputrec->eI))
+ {
+ doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple);
+ }
+ else
+ {
+ doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple);
+ }
+
+ if (doTemperatureCoupling)
+ {
+ real dttc = inputrec->nsttcouple * inputrec->delta_t;
+
+ // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
+ // temperature coupling parameters, which should be reflected in the name of these
+ // subroutines
+ switch (inputrec->etc)
+ {
+ case etcNO: break;
+ case etcBERENDSEN:
+ berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral);
+ break;
+ case etcNOSEHOOVER:
+ nosehoover_tcoupl(&(inputrec->opts), ekind, dttc, state->nosehoover_xi.data(),
+ state->nosehoover_vxi.data(), MassQ);
+ break;
+ case etcVRESCALE:
+ vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data());
+ break;
+ }
+ /* rescale in place here */
+ if (EI_VV(inputrec->eI))
+ {
+ rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array());
+ }
+ }
+ else
+ {
+ // Set the T scaling lambda to 1 to have no scaling
+ // TODO: Do we have to do it on every non-t-couple step?
+ for (int i = 0; (i < inputrec->opts.ngtc); i++)
+ {
+ ekind->tcstat[i].lambda = 1.0;
+ }
+ }
+}
+
+void update_pcouple_before_coordinates(FILE* fplog,
+ int64_t step,
+ const t_inputrec* inputrec,
+ t_state* state,
+ matrix parrinellorahmanMu,
+ matrix M,
+ gmx_bool bInitStep)
+{
+ /* Berendsen P-coupling is completely handled after the coordinate update.
+ * Trotter P-coupling is handled by separate calls to trotter_update().
+ */
+ if (inputrec->epc == epcPARRINELLORAHMAN
+ && do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
+ {
+ real dtpc = inputrec->nstpcouple * inputrec->delta_t;
+
+ parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box,
+ state->box_rel, state->boxv, M, parrinellorahmanMu, bInitStep);
+ }
+}
+
+void update_pcouple_after_coordinates(FILE* fplog,
+ int64_t step,
+ const t_inputrec* inputrec,
+ const t_mdatoms* md,
+ const matrix pressure,
+ const matrix forceVirial,
+ const matrix constraintVirial,
+ matrix pressureCouplingMu,
+ t_state* state,
+ t_nrnb* nrnb,
+ gmx::BoxDeformation* boxDeformation,
+ const bool scaleCoordinates)
+{
+ int start = 0;
+ int homenr = md->homenr;
+
+ /* Cast to real for faster code, no loss in precision (see comment above) */
+ real dt = inputrec->delta_t;
+
+
+ /* now update boxes */
+ switch (inputrec->epc)
+ {
+ case (epcNO): break;
+ case (epcBERENDSEN):
+ if (do_per_step(step, inputrec->nstpcouple))
+ {
+ real dtpc = inputrec->nstpcouple * dt;
+ berendsen_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial,
+ constraintVirial, pressureCouplingMu, &state->baros_integral);
+ berendsen_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start,
+ homenr, state->x.rvec_array(), md->cFREEZE, nrnb, scaleCoordinates);
+ }
+ break;
+ case (epcCRESCALE):
+ if (do_per_step(step, inputrec->nstpcouple))
+ {
+ real dtpc = inputrec->nstpcouple * dt;
+ crescale_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial,
+ constraintVirial, pressureCouplingMu, &state->baros_integral);
+ crescale_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start,
+ homenr, state->x.rvec_array(), state->v.rvec_array(), md->cFREEZE,
+ nrnb, scaleCoordinates);
+ }
+ break;
+ case (epcPARRINELLORAHMAN):
+ if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
+ {
+ /* The box velocities were updated in do_pr_pcoupl,
+ * but we dont change the box vectors until we get here
+ * since we need to be able to shift/unshift above.
+ */
+ real dtpc = inputrec->nstpcouple * dt;
+ for (int i = 0; i < DIM; i++)
+ {
+ for (int m = 0; m <= i; m++)
+ {
+ state->box[i][m] += dtpc * state->boxv[i][m];
+ }
+ }
+ preserve_box_shape(inputrec, state->box_rel, state->box);
+
+ /* Scale the coordinates */
+ if (scaleCoordinates)
+ {
+ auto x = state->x.rvec_array();
+ for (int n = start; n < start + homenr; n++)
+ {
+ tmvmul_ur0(pressureCouplingMu, x[n], x[n]);
+ }
+ }
+ }
+ break;
+ case (epcMTTK):
+ switch (inputrec->epct)
+ {
+ case (epctISOTROPIC):
+ /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
+ ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
+ Side length scales as exp(veta*dt) */
+
+ msmul(state->box, std::exp(state->veta * dt), state->box);
+
+ /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
+ o If we assume isotropic scaling, and box length scaling
+ factor L, then V = L^DIM (det(M)). So dV/dt = DIM
+ L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
+ determinant of B is L^DIM det(M), and the determinant
+ of dB/dt is (dL/dT)^DIM det (M). veta will be
+ (det(dB/dT)/det(B))^(1/3). Then since M =
+ B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
+
+ msmul(state->box, state->veta, state->boxv);
+ break;
+ default: break;
+ }
+ break;
+ default: break;
+ }
+
+ if (boxDeformation)
+ {
+ auto localX = makeArrayRef(state->x).subArray(start, homenr);
+ boxDeformation->apply(localX, state->box, step);
+ }
+}
+
+extern gmx_bool update_randomize_velocities(const t_inputrec* ir,
+ int64_t step,
+ const t_commrec* cr,
+ const t_mdatoms* md,
+ gmx::ArrayRef<gmx::RVec> v,
+ const gmx::Update* upd,
+ const gmx::Constraints* constr)
+{
+
+ real rate = (ir->delta_t) / ir->opts.tau_t[0];
+
+ if (ir->etc == etcANDERSEN && constr != nullptr)
+ {
+ /* Currently, Andersen thermostat does not support constrained
+ systems. Functionality exists in the andersen_tcoupl
+ function in GROMACS 4.5.7 to allow this combination. That
+ code could be ported to the current random-number
+ generation approach, but has not yet been done because of
+ lack of time and resources. */
+ gmx_fatal(FARGS,
+ "Normal Andersen is currently not supported with constraints, use massive "
+ "Andersen instead");
+ }
+
+ /* proceed with andersen if 1) it's fixed probability per
+ particle andersen or 2) it's massive andersen and it's tau_t/dt */
+ if ((ir->etc == etcANDERSEN) || do_per_step(step, gmx::roundToInt(1.0 / rate)))
+ {
+ andersen_tcoupl(ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(),
+ upd->getBoltzmanFactor());
+ return TRUE;
+ }
+ return FALSE;
+}
+
/*
static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
{},
/* for now, we use Elr = 0, because if you want to get it right, you
really should be using PME. Maybe print a warning? */
- pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres) + pcorr;
+ pscal = calc_pres(ir->pbcType, nwall, box, ekinmod, vir, localpres) + pcorr;
vol = det(box);
GW = (vol * (MassQ->Winv / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
*
*/
-real calc_pres(int ePBC, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
+real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
{
int n, m;
real fac;
- if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
+ if (pbcType == PbcType::No || (pbcType == PbcType::XY && nwall != 2))
{
clear_mat(pres);
}
matrix mu,
double* baros_integral)
{
- int d, n;
real scalar_pressure, xy_pressure, p_corr_z;
char buf[STRLEN];
*/
scalar_pressure = 0;
xy_pressure = 0;
- for (d = 0; d < DIM; d++)
+ for (int d = 0; d < DIM; d++)
{
scalar_pressure += pres[d][d] / DIM;
if (d != ZZ)
switch (ir->epct)
{
case epctISOTROPIC:
- for (d = 0; d < DIM; d++)
+ for (int d = 0; d < DIM; d++)
{
mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
}
break;
case epctSEMIISOTROPIC:
- for (d = 0; d < ZZ; d++)
+ for (int d = 0; d < ZZ; d++)
{
mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - xy_pressure) / DIM;
}
mu[ZZ][ZZ] = 1.0 - factor(ZZ, ZZ) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
break;
case epctANISOTROPIC:
- for (d = 0; d < DIM; d++)
+ for (int d = 0; d < DIM; d++)
{
- for (n = 0; n < DIM; n++)
+ for (int n = 0; n < DIM; n++)
{
mu[d][n] = (d == n ? 1.0 : 0.0) - factor(d, n) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
}
p_corr_z = 0;
}
mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
- for (d = 0; d < DIM - 1; d++)
+ for (int d = 0; d < DIM - 1; d++)
{
mu[d][d] = 1.0
+ factor(d, d)
}
}
+void crescale_pcoupl(FILE* fplog,
+ int64_t step,
+ const t_inputrec* ir,
+ real dt,
+ const tensor pres,
+ const matrix box,
+ const matrix force_vir,
+ const matrix constraint_vir,
+ matrix mu,
+ double* baros_integral)
+{
+ /*
+ * Calculate the scaling matrix mu
+ */
+ real scalar_pressure = 0;
+ real xy_pressure = 0;
+ for (int d = 0; d < DIM; d++)
+ {
+ scalar_pressure += pres[d][d] / DIM;
+ if (d != ZZ)
+ {
+ xy_pressure += pres[d][d] / (DIM - 1);
+ }
+ }
+ clear_mat(mu);
+
+ gmx::ThreeFry2x64<64> rng(ir->ld_seed, gmx::RandomDomain::Barostat);
+ gmx::NormalDistribution<real> normalDist;
+ rng.restart(step, 0);
+ real vol = 1.0;
+ for (int d = 0; d < DIM; d++)
+ {
+ vol *= box[d][d];
+ }
+ real gauss;
+ real gauss2;
+ real kt = ir->opts.ref_t[0] * BOLTZ;
+ if (kt < 0.0)
+ {
+ kt = 0.0;
+ }
+
+ switch (ir->epct)
+ {
+ case epctISOTROPIC:
+ gauss = normalDist(rng);
+ for (int d = 0; d < DIM; d++)
+ {
+ const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
+ mu[d][d] = std::exp(-compressibilityFactor * (ir->ref_p[d][d] - scalar_pressure) / DIM
+ + std::sqrt(2.0 * kt * compressibilityFactor * PRESFAC / vol)
+ * gauss / DIM);
+ }
+ break;
+ case epctSEMIISOTROPIC:
+ gauss = normalDist(rng);
+ gauss2 = normalDist(rng);
+ for (int d = 0; d < ZZ; d++)
+ {
+ const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
+ mu[d][d] = std::exp(
+ -compressibilityFactor * (ir->ref_p[d][d] - xy_pressure) / DIM
+ + std::sqrt((DIM - 1) * 2.0 * kt * compressibilityFactor * PRESFAC / vol / DIM)
+ / (DIM - 1) * gauss);
+ }
+ {
+ const real compressibilityFactor = ir->compress[ZZ][ZZ] * dt / ir->tau_p;
+ mu[ZZ][ZZ] = std::exp(
+ -compressibilityFactor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
+ + std::sqrt(2.0 * kt * compressibilityFactor * PRESFAC / vol / DIM) * gauss2);
+ }
+ break;
+ case epctSURFACETENSION:
+ gauss = normalDist(rng);
+ gauss2 = normalDist(rng);
+ for (int d = 0; d < ZZ; d++)
+ {
+ const real compressibilityFactor = ir->compress[d][d] * dt / ir->tau_p;
+ /* Notice: we here use ref_p[ZZ][ZZ] as isotropic pressure and ir->ref_p[d][d] as surface tension */
+ mu[d][d] = std::exp(
+ -compressibilityFactor
+ * (ir->ref_p[ZZ][ZZ] - ir->ref_p[d][d] / box[ZZ][ZZ] - xy_pressure) / DIM
+ + std::sqrt(4.0 / 3.0 * kt * compressibilityFactor * PRESFAC / vol)
+ / (DIM - 1) * gauss);
+ }
+ {
+ const real compressibilityFactor = ir->compress[ZZ][ZZ] * dt / ir->tau_p;
+ mu[ZZ][ZZ] = std::exp(
+ -compressibilityFactor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
+ + std::sqrt(2.0 / 3.0 * kt * compressibilityFactor * PRESFAC / vol) * gauss2);
+ }
+ break;
+ default:
+ gmx_fatal(FARGS, "C-rescale pressure coupling type %s not supported yet\n",
+ EPCOUPLTYPETYPE(ir->epct));
+ }
+ /* To fullfill the orientation restrictions on triclinic boxes
+ * we will set mu_yx, mu_zx and mu_zy to 0 and correct
+ * the other elements of mu to first order.
+ */
+ mu[YY][XX] += mu[XX][YY];
+ mu[ZZ][XX] += mu[XX][ZZ];
+ mu[ZZ][YY] += mu[YY][ZZ];
+ mu[XX][YY] = 0;
+ mu[XX][ZZ] = 0;
+ mu[YY][ZZ] = 0;
+
+ /* Keep track of the work the barostat applies on the system.
+ * Without constraints force_vir tells us how Epot changes when scaling.
+ * With constraints constraint_vir gives us the constraint contribution
+ * to both Epot and Ekin. Although we are not scaling velocities, scaling
+ * the coordinates leads to scaling of distances involved in constraints.
+ * This in turn changes the angular momentum (even if the constrained
+ * distances are corrected at the next step). The kinetic component
+ * of the constraint virial captures the angular momentum change.
+ */
+ for (int d = 0; d < DIM; d++)
+ {
+ for (int n = 0; n <= d; n++)
+ {
+ *baros_integral -=
+ 2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
+ }
+ }
+
+ if (debug)
+ {
+ pr_rvecs(debug, 0, "PC: pres ", pres, 3);
+ pr_rvecs(debug, 0, "PC: mu ", mu, 3);
+ }
+
+ if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
+ || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
+ {
+ char buf[STRLEN];
+ char buf2[22];
+ sprintf(buf,
+ "\nStep %s Warning: pressure scaling more than 1%%, "
+ "mu: %g %g %g\n",
+ gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
+ if (fplog)
+ {
+ fprintf(fplog, "%s", buf);
+ }
+ fprintf(stderr, "%s", buf);
+ }
+}
+
+void crescale_pscale(const t_inputrec* ir,
+ const matrix mu,
+ matrix box,
+ matrix box_rel,
+ int start,
+ int nr_atoms,
+ rvec x[],
+ rvec v[],
+ const unsigned short cFREEZE[],
+ t_nrnb* nrnb,
+ const bool scaleCoordinates)
+{
+ ivec* nFreeze = ir->opts.nFreeze;
+ int nthreads gmx_unused;
+ matrix inv_mu;
+
+#ifndef __clang_analyzer__
+ nthreads = gmx_omp_nthreads_get(emntUpdate);
+#endif
+
+ gmx::invertBoxMatrix(mu, inv_mu);
+
+ /* Scale the positions and the velocities */
+ if (scaleCoordinates)
+ {
+#pragma omp parallel for num_threads(nthreads) schedule(static)
+ for (int n = start; n < start + nr_atoms; n++)
+ {
+ // Trivial OpenMP region that does not throw
+ int g;
+
+ if (cFREEZE == nullptr)
+ {
+ g = 0;
+ }
+ else
+ {
+ g = cFREEZE[n];
+ }
+
+ if (!nFreeze[g][XX])
+ {
+ x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
+ v[n][XX] = inv_mu[XX][XX] * v[n][XX] + inv_mu[YY][XX] * v[n][YY]
+ + inv_mu[ZZ][XX] * v[n][ZZ];
+ }
+ if (!nFreeze[g][YY])
+ {
+ x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
+ v[n][YY] = inv_mu[YY][YY] * v[n][YY] + inv_mu[ZZ][YY] * v[n][ZZ];
+ }
+ if (!nFreeze[g][ZZ])
+ {
+ x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
+ v[n][ZZ] = inv_mu[ZZ][ZZ] * v[n][ZZ];
+ }
+ }
+ }
+ /* compute final boxlengths */
+ for (int d = 0; d < DIM; d++)
+ {
+ box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
+ box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
+ box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
+ }
+
+ preserve_box_shape(ir, box_rel, box);
+
+ /* (un)shifting should NOT be done after this,
+ * since the box vectors might have changed
+ */
+ inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
+}
+
void berendsen_pscale(const t_inputrec* ir,
const matrix mu,
matrix box,
const double* ivxi = &state->nosehoover_vxi[i * nh];
const double* iQinv = &(MassQ->Qinv[i * nh]);
- int nd = static_cast<int>(ir->opts.nrdf[i]);
+ real nd = ir->opts.nrdf[i];
real reft = std::max<real>(ir->opts.ref_t[i], 0);
real kT = BOLTZ * reft;
{
energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
/* contribution from the thermal variable of the NH chain */
- int ndj;
+ real ndj = 0;
if (j == 0)
{
ndj = nd;
energyNPT += energyPressureMTTK(ir, state, MassQ);
}
break;
- case epcBERENDSEN: energyNPT += state->baros_integral; break;
+ case epcBERENDSEN:
+ case epcCRESCALE: energyNPT += state->baros_integral; break;
default:
GMX_RELEASE_ASSERT(
false,
}
}
- update_temperature_constants(upd->sd(), ir);
+ upd->update_temperature_constants(*ir);
}
void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
{
please_cite(fplog, "Bussi2007a");
}
+ if (ir.epc == epcCRESCALE)
+ {
+ please_cite(fplog, "Bernetti2020");
+ }
// TODO this is actually an integrator, not a coupling algorithm
if (ir.eI == eiSD1)
{
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2012-2018, The GROMACS development team.
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020,2021, 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.
#include "gromacs/utility/gmxmpi.h"
#include "gromacs/utility/smalloc.h"
+ #include "expanded_internal.h"
+
static void init_df_history_weights(df_history_t* dfhist, const t_expanded* expand, int nlim)
{
int i;
int64_t step)
{
gmx_bool bSufficientSamples;
+ real acceptanceWeight;
int i;
- int n0, np1, nm1, nval, min_nvalm, min_nvalp, maxc;
- real omega_m1_0, omega_p1_0, clam_osum;
- real de, de_function;
- real cnval, zero_sum_weights;
+ int min_nvalm, min_nvalp, maxc;
+ real omega_m1_0, omega_p1_0;
+ real zero_sum_weights;
real *omegam_array, *weightsm_array, *omegap_array, *weightsp_array, *varm_array, *varp_array,
*dwp_array, *dwm_array;
- real clam_varm, clam_varp, clam_weightsm, clam_weightsp, clam_minvar;
+ real clam_varm, clam_varp, clam_osum, clam_weightsm, clam_weightsp, clam_minvar;
real * lam_variance, *lam_dg;
double* p_k;
double pks = 0;
- real chi_m1_0, chi_p1_0, chi_m2_0, chi_p2_0, chi_p1_m1, chi_p2_m1, chi_m1_p1, chi_m2_p1;
- /* if we have equilibrated the weights, exit now */
+ /* Future potential todos for this function (see #3848):
+ * - Update the names in the dhist structure to be clearer. Not done for now since this
+ * a bugfix update and we are mininizing other code changes.
+ * - Modularize the code some more.
+ * - potentially merge with accelerated weight histogram functionality, since it's very similar.
+ */
+ /* if we have equilibrated the expanded ensemble weights, we are not updating them, so exit now */
if (dfhist->bEquil)
{
return FALSE;
if (EWL(expand->elamstats))
{
- if (expand->elamstats == elamstatsWL) /* Standard Wang-Landau */
+ if (expand->elamstats == elamstatsWL) /* Using standard Wang-Landau for weight updates */
{
dfhist->sum_weights[fep_state] -= dfhist->wl_delta;
dfhist->wl_histo[fep_state] += 1.0;
}
- else if (expand->elamstats == elamstatsWWL) /* Weighted Wang-Landau */
+ else if (expand->elamstats == elamstatsWWL)
+ /* Using weighted Wang-Landau for weight updates.
+ * Very closly equivalent to accelerated weight histogram approach
+ * applied to expanded ensemble. */
{
snew(p_k, nlim);
if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMETROPOLIS
|| expand->elamstats == elamstatsMINVAR)
{
-
- de_function = 0; /* to get rid of warnings, but this value will not be used because of the logic */
maxc = 2 * expand->c_range + 1;
snew(lam_dg, nlim);
snew(varm_array, maxc);
snew(dwm_array, maxc);
- /* unpack the current lambdas -- we will only update 2 of these */
+ /* unpack the values of the free energy differences and the
+ * variance in their estimates between nearby lambdas. We will
+ * only actually update 2 of these, the state we are currently
+ * at and the one we end up moving to
+ */
for (i = 0; i < nlim - 1; i++)
{ /* only through the second to last */
gmx::square(dfhist->sum_variance[i + 1]) - gmx::square(dfhist->sum_variance[i]);
}
- /* accumulate running averages */
- for (nval = 0; nval < maxc; nval++)
+ /* accumulate running averages of thermodynamic averages for Bennett Acceptance Ratio-based
+ * estimates of the free energy .
+ * Rather than peforming self-consistent estimation of the free energies at each step,
+ * we keep track of an array of possible different free energies (cnvals),
+ * and we self-consistently choose the best one. The one that leads to a free energy estimate
+ * that is closest to itself is the best estimate of the free energy. It is essentially a
+ * parallellized version of self-consistent iteration. maxc is the number of these constants. */
+
+ for (int nval = 0; nval < maxc; nval++)
{
- /* constants for later use */
- cnval = static_cast<real>(nval - expand->c_range);
- /* actually, should be able to rewrite it w/o exponential, for better numerical stability */
+ const real cnval = static_cast<real>(nval - expand->c_range);
+
+ /* Compute acceptance criterion weight to the state below this one for use in averages.
+ * Note we do not have to have just moved from that state to use this free energy
+ * estimate; these are essentially "virtual" moves. */
+
if (fep_state > 0)
{
- de = std::exp(cnval - (scaled_lamee[fep_state] - scaled_lamee[fep_state - 1]));
- if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR)
- {
- de_function = 1.0 / (1.0 + de);
- }
- else if (expand->elamstats == elamstatsMETROPOLIS)
- {
- if (de < 1.0)
- {
- de_function = 1.0;
- }
- else
- {
- de_function = 1.0 / de;
- }
- }
- dfhist->accum_m[fep_state][nval] += de_function;
- dfhist->accum_m2[fep_state][nval] += de_function * de_function;
+ const auto lambdaEnergyDifference =
+ cnval - (scaled_lamee[fep_state] - scaled_lamee[fep_state - 1]);
+ acceptanceWeight =
+ gmx::calculateAcceptanceWeight(expand->elamstats, lambdaEnergyDifference);
+ dfhist->accum_m[fep_state][nval] += acceptanceWeight;
+ dfhist->accum_m2[fep_state][nval] += acceptanceWeight * acceptanceWeight;
}
+ // Compute acceptance criterion weight to transition to the next state
if (fep_state < nlim - 1)
{
- de = std::exp(-cnval + (scaled_lamee[fep_state + 1] - scaled_lamee[fep_state]));
- if (expand->elamstats == elamstatsBARKER || expand->elamstats == elamstatsMINVAR)
- {
- de_function = 1.0 / (1.0 + de);
- }
- else if (expand->elamstats == elamstatsMETROPOLIS)
- {
- if (de < 1.0)
- {
- de_function = 1.0;
- }
- else
- {
- de_function = 1.0 / de;
- }
- }
- dfhist->accum_p[fep_state][nval] += de_function;
- dfhist->accum_p2[fep_state][nval] += de_function * de_function;
+ const auto lambdaEnergyDifference =
+ -cnval + (scaled_lamee[fep_state + 1] - scaled_lamee[fep_state]);
+ acceptanceWeight =
+ gmx::calculateAcceptanceWeight(expand->elamstats, lambdaEnergyDifference);
+ dfhist->accum_p[fep_state][nval] += acceptanceWeight;
+ dfhist->accum_p2[fep_state][nval] += acceptanceWeight * acceptanceWeight;
}
- /* Metropolis transition and Barker transition (unoptimized Bennett) acceptance weight determination */
+ /* Determination of Metropolis transition and Barker transition weights */
- n0 = dfhist->n_at_lam[fep_state];
+ int numObservationsCurrentState = dfhist->n_at_lam[fep_state];
+ /* determine the number of observations above and below the current state */
+ int numObservationsLowerState = 0;
if (fep_state > 0)
{
- nm1 = dfhist->n_at_lam[fep_state - 1];
- }
- else
- {
- nm1 = 0;
+ numObservationsLowerState = dfhist->n_at_lam[fep_state - 1];
}
+ int numObservationsHigherState = 0;
if (fep_state < nlim - 1)
{
- np1 = dfhist->n_at_lam[fep_state + 1];
- }
- else
- {
- np1 = 0;
+ numObservationsHigherState = dfhist->n_at_lam[fep_state + 1];
}
- /* logic SHOULD keep these all set correctly whatever the logic, but apparently it can't figure it out. */
- chi_m1_0 = chi_p1_0 = chi_m2_0 = chi_p2_0 = chi_p1_m1 = chi_p2_m1 = chi_m1_p1 = chi_m2_p1 = 0;
+ /* Calculate the biases for each expanded ensemble state that minimize the total
+ * variance, as implemented in Martinez-Veracoechea and Escobedo,
+ * J. Phys. Chem. B 2008, 112, 8120-8128
+ *
+ * The variance associated with the free energy estimate between two states i and j
+ * is calculated as
+ * Var(i,j) = {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} / numObservations(i->j)
+ * + {avg[xi(j->i)^2] / avg[xi(j->i)]^2 - 1} / numObservations(j->i)
+ * where xi(i->j) is the acceptance factor / weight associated with moving from state i to j
+ * As we are calculating the acceptance factor to the neighbors every time we're visiting
+ * a state, numObservations(i->j) == numObservations(i) and numObservations(j->i) == numObservations(j)
+ */
- if (n0 > 0)
+ /* Accumulation of acceptance weight averages between the current state and the
+ * states +1 (p1) and -1 (m1), averaged at current state (0)
+ */
+ real avgAcceptanceCurrentToLower = 0;
+ real avgAcceptanceCurrentToHigher = 0;
+ /* Accumulation of acceptance weight averages quantities between states 0
+ * and states +1 and -1, squared
+ */
+ real avgAcceptanceCurrentToLowerSquared = 0;
+ real avgAcceptanceCurrentToHigherSquared = 0;
+ /* Accumulation of free energy quantities from lower state (m1) to current state (0) and squared */
+ real avgAcceptanceLowerToCurrent = 0;
+ real avgAcceptanceLowerToCurrentSquared = 0;
+ /* Accumulation of free energy quantities from upper state (p1) to current state (0) and squared */
+ real avgAcceptanceHigherToCurrent = 0;
+ real avgAcceptanceHigherToCurrentSquared = 0;
+
+ if (numObservationsCurrentState > 0)
{
- chi_m1_0 = dfhist->accum_m[fep_state][nval] / n0;
- chi_p1_0 = dfhist->accum_p[fep_state][nval] / n0;
- chi_m2_0 = dfhist->accum_m2[fep_state][nval] / n0;
- chi_p2_0 = dfhist->accum_p2[fep_state][nval] / n0;
+ avgAcceptanceCurrentToLower = dfhist->accum_m[fep_state][nval] / numObservationsCurrentState;
+ avgAcceptanceCurrentToHigher =
+ dfhist->accum_p[fep_state][nval] / numObservationsCurrentState;
+ avgAcceptanceCurrentToLowerSquared =
+ dfhist->accum_m2[fep_state][nval] / numObservationsCurrentState;
+ avgAcceptanceCurrentToHigherSquared =
+ dfhist->accum_p2[fep_state][nval] / numObservationsCurrentState;
}
- if ((fep_state > 0) && (nm1 > 0))
+ if ((fep_state > 0) && (numObservationsLowerState > 0))
{
- chi_p1_m1 = dfhist->accum_p[fep_state - 1][nval] / nm1;
- chi_p2_m1 = dfhist->accum_p2[fep_state - 1][nval] / nm1;
+ avgAcceptanceLowerToCurrent =
+ dfhist->accum_p[fep_state - 1][nval] / numObservationsLowerState;
+ avgAcceptanceLowerToCurrentSquared =
+ dfhist->accum_p2[fep_state - 1][nval] / numObservationsLowerState;
}
- if ((fep_state < nlim - 1) && (np1 > 0))
+ if ((fep_state < nlim - 1) && (numObservationsHigherState > 0))
{
- chi_m1_p1 = dfhist->accum_m[fep_state + 1][nval] / np1;
- chi_m2_p1 = dfhist->accum_m2[fep_state + 1][nval] / np1;
+ avgAcceptanceHigherToCurrent =
+ dfhist->accum_m[fep_state + 1][nval] / numObservationsHigherState;
+ avgAcceptanceHigherToCurrentSquared =
+ dfhist->accum_m2[fep_state + 1][nval] / numObservationsHigherState;
}
-
- omega_m1_0 = 0;
- omega_p1_0 = 0;
- clam_weightsm = 0;
- clam_weightsp = 0;
- clam_varm = 0;
- clam_varp = 0;
+ /* These are accumulation of positive values (see definition of acceptance functions
+ * above), or of squares of positive values.
+ * We're taking this for granted in the following calculation, so make sure
+ * here that nothing weird happened. Although technically all values should be positive,
+ * because of floating point precisions, they might be numerically zero. */
+ GMX_RELEASE_ASSERT(
+ avgAcceptanceCurrentToLower >= 0 && avgAcceptanceCurrentToLowerSquared >= 0
+ && avgAcceptanceCurrentToHigher >= 0
+ && avgAcceptanceCurrentToHigherSquared >= 0 && avgAcceptanceLowerToCurrent >= 0
+ && avgAcceptanceLowerToCurrentSquared >= 0 && avgAcceptanceHigherToCurrent >= 0
+ && avgAcceptanceHigherToCurrentSquared >= 0,
+ "By definition, the acceptance factors should all be nonnegative.");
+
+ real varianceCurrentToLower = 0;
+ real varianceCurrentToHigher = 0;
+ real weightDifferenceToLower = 0;
+ real weightDifferenceToHigher = 0;
+ real varianceToLower = 0;
+ real varianceToHigher = 0;
if (fep_state > 0)
{
- if (n0 > 0)
+ if (numObservationsCurrentState > 0)
{
- omega_m1_0 = chi_m2_0 / (chi_m1_0 * chi_m1_0) - 1.0;
- if (nm1 > 0)
+ /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
+ *
+ * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
+ * acceptances are all positive!), and hence
+ * {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0 for avg[xi(i->j)] -> 0
+ * We're catching that case explicitly to avoid numerical
+ * problems dividing by zero when the overlap between states is small (#3304)
+ */
+ if (avgAcceptanceCurrentToLower > 0)
{
- real omega_p1_m1 = chi_p2_m1 / (chi_p1_m1 * chi_p1_m1) - 1.0;
- clam_weightsm = (std::log(chi_m1_0) - std::log(chi_p1_m1)) + cnval;
- clam_varm = (1.0 / n0) * (omega_m1_0) + (1.0 / nm1) * (omega_p1_m1);
+ varianceCurrentToLower =
+ avgAcceptanceCurrentToLowerSquared
+ / (avgAcceptanceCurrentToLower * avgAcceptanceCurrentToLower)
+ - 1.0;
+ }
+ if (numObservationsLowerState > 0)
+ {
+ /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
+ *
+ * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
+ * acceptances are all positive!), and hence
+ * {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0 for avg[xi(i->j)] -> 0
+ * We're catching that case explicitly to avoid numerical
+ * problems dividing by zero when the overlap between states is small (#3304)
+ */
+ real varianceLowerToCurrent = 0;
+ if (avgAcceptanceLowerToCurrent > 0)
+ {
+ varianceLowerToCurrent =
+ avgAcceptanceLowerToCurrentSquared
+ / (avgAcceptanceLowerToCurrent * avgAcceptanceLowerToCurrent)
+ - 1.0;
+ }
+ /* Free energy difference to the state one state lower */
+ /* if these either of these quantities are zero, the energies are */
+ /* way too large for the dynamic range. We need an alternate guesstimate */
+ if ((avgAcceptanceCurrentToLower == 0) || (avgAcceptanceLowerToCurrent == 0))
+ {
+ weightDifferenceToLower =
+ (scaled_lamee[fep_state] - scaled_lamee[fep_state - 1]);
+ }
+ else
+ {
+ weightDifferenceToLower = (std::log(avgAcceptanceCurrentToLower)
+ - std::log(avgAcceptanceLowerToCurrent))
+ + cnval;
+ }
+ /* Variance of the free energy difference to the one state lower */
+ varianceToLower =
+ (1.0 / numObservationsCurrentState) * (varianceCurrentToLower)
+ + (1.0 / numObservationsLowerState) * (varianceLowerToCurrent);
}
}
}
if (fep_state < nlim - 1)
{
- if (n0 > 0)
+ if (numObservationsCurrentState > 0)
{
- omega_p1_0 = chi_p2_0 / (chi_p1_0 * chi_p1_0) - 1.0;
- if (np1 > 0)
+ /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
+ *
+ * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
+ * acceptances are all positive!), and hence
+ * {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0 for avg[xi(i->j)] -> 0
+ * We're catching that case explicitly to avoid numerical
+ * problems dividing by zero when the overlap between states is small (#3304)
+ */
+
+ if (avgAcceptanceCurrentToHigher < 0)
+ {
+ varianceCurrentToHigher =
+ avgAcceptanceCurrentToHigherSquared
+ / (avgAcceptanceCurrentToHigher * avgAcceptanceCurrentToHigher)
+ - 1.0;
+ }
+ if (numObservationsHigherState > 0)
{
- real omega_m1_p1 = chi_m2_p1 / (chi_m1_p1 * chi_m1_p1) - 1.0;
- clam_weightsp = (std::log(chi_m1_p1) - std::log(chi_p1_0)) + cnval;
- clam_varp = (1.0 / np1) * (omega_m1_p1) + (1.0 / n0) * (omega_p1_0);
+ /* Calculate {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1}
+ *
+ * Note that if avg[xi(i->j)] == 0, also avg[xi(i->j)^2] == 0 (since the
+ * acceptances are all positive!), and hence
+ * {avg[xi(i->j)^2] / avg[xi(i->j)]^2 - 1} -> 0 for avg[xi(i->j)] -> 0
+ * We're catching that case explicitly to avoid numerical
+ * problems dividing by zero when the overlap between states is small (#3304)
+ */
+ real varianceHigherToCurrent = 0;
+ if (avgAcceptanceHigherToCurrent > 0)
+ {
+ varianceHigherToCurrent =
+ avgAcceptanceHigherToCurrentSquared
+ / (avgAcceptanceHigherToCurrent * avgAcceptanceHigherToCurrent)
+ - 1.0;
+ }
+ /* Free energy difference to the state one state higher */
+ /* if these either of these quantities are zero, the energies are */
+ /* way too large for the dynamic range. We need an alternate guesstimate */
+ if ((avgAcceptanceHigherToCurrent == 0) || (avgAcceptanceCurrentToHigher == 0))
+ {
+ weightDifferenceToHigher =
+ (scaled_lamee[fep_state + 1] - scaled_lamee[fep_state]);
+ }
+ else
+ {
+ weightDifferenceToHigher = (std::log(avgAcceptanceHigherToCurrent)
+ - std::log(avgAcceptanceCurrentToHigher))
+ + cnval;
+ }
+ /* Variance of the free energy difference to the one state higher */
+ varianceToHigher =
+ (1.0 / numObservationsHigherState) * (varianceHigherToCurrent)
+ + (1.0 / numObservationsCurrentState) * (varianceCurrentToHigher);
}
}
}
- if (n0 > 0)
+ if (numObservationsCurrentState > 0)
{
- omegam_array[nval] = omega_m1_0;
+ omegam_array[nval] = varianceCurrentToLower;
}
else
{
omegam_array[nval] = 0;
}
- weightsm_array[nval] = clam_weightsm;
- varm_array[nval] = clam_varm;
- if (nm1 > 0)
+ weightsm_array[nval] = weightDifferenceToLower;
+ varm_array[nval] = varianceToLower;
+ if (numObservationsLowerState > 0)
{
- dwm_array[nval] = fabs((cnval + std::log((1.0 * n0) / nm1)) - lam_dg[fep_state - 1]);
+ dwm_array[nval] =
+ fabs((cnval + std::log((1.0 * numObservationsCurrentState) / numObservationsLowerState))
+ - lam_dg[fep_state - 1]);
}
else
{
dwm_array[nval] = std::fabs(cnval - lam_dg[fep_state - 1]);
}
- if (n0 > 0)
+ if (numObservationsCurrentState > 0)
{
- omegap_array[nval] = omega_p1_0;
+ omegap_array[nval] = varianceCurrentToHigher;
}
else
{
omegap_array[nval] = 0;
}
- weightsp_array[nval] = clam_weightsp;
- varp_array[nval] = clam_varp;
- if ((np1 > 0) && (n0 > 0))
+ weightsp_array[nval] = weightDifferenceToHigher;
+ varp_array[nval] = varianceToHigher;
+ if ((numObservationsHigherState > 0) && (numObservationsCurrentState > 0))
{
- dwp_array[nval] = fabs((cnval + std::log((1.0 * np1) / n0)) - lam_dg[fep_state]);
+ dwp_array[nval] =
+ fabs((cnval + std::log((1.0 * numObservationsHigherState) / numObservationsCurrentState))
+ - lam_dg[fep_state]);
}
else
{
}
}
- /* find the C's closest to the old weights value */
+ /* find the free energy estimate closest to the guessed weight's value */
min_nvalm = FindMinimum(dwm_array, maxc);
omega_m1_0 = omegam_array[min_nvalm];
if (expand->elamstats == elamstatsMINVAR)
{
bSufficientSamples = TRUE;
- /* make sure they are all past a threshold */
+ /* make sure the number of samples in each state are all
+ * past a user-specified threshold
+ */
for (i = 0; i < nlim; i++)
{
if (dfhist->n_at_lam[i] < expand->minvarmin)
de = weighted_lamee[lamtrial] - weighted_lamee[fep_state];
if (expand->elmcmove == elmcmoveMETROPOLIS)
{
- tprob = 1.0;
- trialprob = std::exp(de);
- if (trialprob < tprob)
+ tprob = 1.0;
+ if (de < 0)
{
- tprob = trialprob;
+ tprob = std::exp(de);
}
propose[fep_state] = 0;
propose[lamtrial] = 1.0; /* note that this overwrites the above line if fep_state = ntrial, which only occurs at the ends */
}
else if (expand->elmcmove == elmcmoveBARKER)
{
- tprob = 1.0 / (1.0 + std::exp(-de));
-
+ if (de > 0) /* Numerically stable version */
+ {
+ tprob = 1.0 / (1.0 + std::exp(-de));
+ }
+ else if (de < 0)
+ {
+ tprob = std::exp(de) / (std::exp(de) + 1.0);
+ }
propose[fep_state] = (1 - tprob);
propose[lamtrial] +=
tprob; /* we add, to account for the fact that at the end, they might be the same point */
if (ir->bSimTemp)
{
/* Note -- this assumes no mass changes, since kinetic energy is not added . . . */
- scaled_lamee[i] = (enerd->enerpart_lambda[i + 1] - enerd->enerpart_lambda[0])
- / (simtemp->temperatures[i] * BOLTZ)
+ scaled_lamee[i] = enerd->foreignLambdaTerms.deltaH(i) / (simtemp->temperatures[i] * BOLTZ)
+ enerd->term[F_EPOT]
* (1.0 / (simtemp->temperatures[i])
- 1.0 / (simtemp->temperatures[fep_state]))
}
else
{
- scaled_lamee[i] = (enerd->enerpart_lambda[i + 1] - enerd->enerpart_lambda[0])
- / (expand->mc_temp * BOLTZ);
+ scaled_lamee[i] = enerd->foreignLambdaTerms.deltaH(i) / (expand->mc_temp * BOLTZ);
/* mc_temp is currently set to the system reft unless otherwise defined */
}
#
# This file is part of the GROMACS molecular simulation package.
#
- # Copyright (c) 2014,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+ # Copyright (c) 2014,2016,2017,2018,2019,2020,2021, 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.
# To help us fund GROMACS development, we humbly ask that you cite
# the research papers on the package. Check out http://www.gromacs.org.
-gmx_add_unit_test(MdlibUnitTest mdlib-test
- calc_verletbuf.cpp
- constr.cpp
- constrtestdata.cpp
- constrtestrunners.cpp
- ebin.cpp
- energyoutput.cpp
- expanded.cpp
- leapfrog.cpp
- leapfrogtestdata.cpp
- leapfrogtestrunners.cpp
- settle.cpp
- settletestdata.cpp
- settletestrunners.cpp
- shake.cpp
- simulationsignal.cpp
- updategroups.cpp
- updategroupscog.cpp)
-
-# TODO: Make CUDA source to compile inside the testing framework
-if(GMX_USE_CUDA)
- gmx_add_libgromacs_sources(constrtestrunners.cu
- leapfrogtestrunners.cu
- settletestrunners.cu)
-endif()
+gmx_add_unit_test(MdlibUnitTest mdlib-test HARDWARE_DETECTION
+ CPP_SOURCE_FILES
+ calc_verletbuf.cpp
+ constr.cpp
+ constrtestdata.cpp
+ constrtestrunners.cpp
+ ebin.cpp
+ energydrifttracker.cpp
+ energyoutput.cpp
++ expanded.cpp
+ freeenergyparameters.cpp
+ leapfrog.cpp
+ leapfrogtestdata.cpp
+ leapfrogtestrunners.cpp
+ settle.cpp
+ settletestdata.cpp
+ settletestrunners.cpp
+ shake.cpp
+ simulationsignal.cpp
+ updategroups.cpp
+ updategroupscog.cpp
+ GPU_CPP_SOURCE_FILES
+ leapfrogtestrunners_gpu.cpp
+ CUDA_CU_SOURCE_FILES
+ constrtestrunners.cu
+ settletestrunners.cu
+ )
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
-- * Copyright (c) 2011-2019,2020, by the GROMACS development team, led by
++ * Copyright (c) 2011-2019,2020,2021, 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.
#include <algorithm>
#include <memory>
+#include <numeric>
-#include "gromacs/awh/awh.h"
+#include "gromacs/applied_forces/awh/awh.h"
#include "gromacs/commandline/filenm.h"
#include "gromacs/domdec/collect.h"
#include "gromacs/domdec/dlbtiming.h"
#include "gromacs/domdec/domdec.h"
#include "gromacs/domdec/domdec_network.h"
#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/domdec/gpuhaloexchange.h"
#include "gromacs/domdec/mdsetup.h"
#include "gromacs/domdec/partition.h"
#include "gromacs/essentialdynamics/edsam.h"
-#include "gromacs/ewald/pme.h"
#include "gromacs/ewald/pme_load_balancing.h"
+#include "gromacs/ewald/pme_pp.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/gpu_utils/device_stream_manager.h"
#include "gromacs/gpu_utils/gpu_utils.h"
#include "gromacs/imd/imd.h"
-#include "gromacs/listed_forces/manage_threading.h"
+#include "gromacs/listed_forces/listed_forces.h"
#include "gromacs/math/functions.h"
-#include "gromacs/math/utilities.h"
+#include "gromacs/math/invertmatrix.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/mdlib/checkpointhandler.h"
#include "gromacs/mdlib/compute_io.h"
#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/coupling.h"
#include "gromacs/mdlib/ebin.h"
#include "gromacs/mdlib/enerdata_utils.h"
#include "gromacs/mdlib/energyoutput.h"
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/freeenergyparameters.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
#include "gromacs/mdlib/tgroup.h"
#include "gromacs/mdlib/trajectory_writing.h"
#include "gromacs/mdlib/update.h"
-#include "gromacs/mdlib/update_constrain_cuda.h"
+#include "gromacs/mdlib/update_constrain_gpu.h"
#include "gromacs/mdlib/vcm.h"
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdrunutility/handlerestart.h"
#include "gromacs/mdtypes/df_history.h"
#include "gromacs/mdtypes/energyhistory.h"
#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forcebuffers.h"
#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/group.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/mdrunoptions.h"
+#include "gromacs/mdtypes/multipletimestepping.h"
#include "gromacs/mdtypes/observableshistory.h"
#include "gromacs/mdtypes/pullhistory.h"
#include "gromacs/mdtypes/simulation_workload.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/mdtypes/state_propagator_data_gpu.h"
-#include "gromacs/modularsimulator/energyelement.h"
+#include "gromacs/modularsimulator/energydata.h"
#include "gromacs/nbnxm/gpu_data_mgmt.h"
#include "gromacs/nbnxm/nbnxm.h"
-#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/output.h"
#include "gromacs/pulling/pull.h"
// will go away eventually.
t_inputrec* ir = inputrec;
int64_t step, step_rel;
- double t, t0 = ir->init_t, lam0[efptNR];
+ double t, t0 = ir->init_t;
gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
gmx_bool bNS = FALSE, bNStList, bStopCM, bFirstStep, bInitStep, bLastStep = FALSE;
gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
gmx_bool do_ene, do_log, do_verbose;
gmx_bool bMasterState;
unsigned int force_flags;
- tensor force_vir = { { 0 } }, shake_vir = { { 0 } }, total_vir = { { 0 } }, tmp_vir = { { 0 } },
- pres = { { 0 } };
- int i, m;
- rvec mu_tot;
- matrix pressureCouplingMu, M;
- gmx_repl_ex_t repl_ex = nullptr;
- gmx_localtop_t top;
- PaddedHostVector<gmx::RVec> f{};
- gmx_global_stat_t gstat;
- t_graph* graph = nullptr;
- gmx_shellfc_t* shellfc;
- gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
- gmx_bool bTemp, bPres, bTrotter;
- real dvdl_constr;
- std::vector<RVec> cbuf;
- matrix lastbox;
- int lamnew = 0;
+ tensor force_vir = { { 0 } }, shake_vir = { { 0 } }, total_vir = { { 0 } }, pres = { { 0 } };
+ int i, m;
+ rvec mu_tot;
+ matrix pressureCouplingMu, M;
+ gmx_repl_ex_t repl_ex = nullptr;
+ gmx_global_stat_t gstat;
+ gmx_shellfc_t* shellfc;
+ gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
+ gmx_bool bTemp, bPres, bTrotter;
+ real dvdl_constr;
+ std::vector<RVec> cbuf;
+ matrix lastbox;
+ int lamnew = 0;
/* for FEP */
int nstfep = 0;
double cycles;
int nstglobalcomm = computeGlobalCommunicationPeriod(mdlog, ir, cr);
bGStatEveryStep = (nstglobalcomm == 1);
- SimulationGroups* groups = &top_global->groups;
+ const SimulationGroups* groups = &top_global->groups;
std::unique_ptr<EssentialDynamics> ed = nullptr;
if (opt2bSet("-ei", nfile, fnm))
"Either specify the -ei option to mdrun, or do not use this checkpoint file.");
}
- initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda, lam0);
- Update upd(ir, deform);
+ int* fep_state = MASTER(cr) ? &state_global->fep_state : nullptr;
+ gmx::ArrayRef<real> lambda = MASTER(cr) ? state_global->lambda : gmx::ArrayRef<real>();
+ initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda);
+ Update upd(*ir, deform);
const bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
+ const t_fcdata& fcdata = *fr->fcdata;
+
bool simulationsShareState = false;
int nstSignalComm = nstglobalcomm;
{
// TODO This implementation of ensemble orientation restraints is nasty because
// a user can't just do multi-sim with single-sim orientation restraints.
bool usingEnsembleRestraints =
- (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
+ (fcdata.disres->nsystems > 1) || ((ms != nullptr) && (fcdata.orires->nr != 0));
bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
// Replica exchange, ensemble restraints and AWH need all
init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, ir,
top_global, oenv, wcycle, startingBehavior, simulationsShareState, ms);
gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work,
- mdoutf_get_fp_dhdl(outf), false, startingBehavior, mdModulesNotifier);
+ mdoutf_get_fp_dhdl(outf), false, startingBehavior,
+ simulationsShareState, mdModulesNotifier);
gstat = global_stat_init(ir);
+ const auto& simulationWork = runScheduleWork->simulationWork;
+ const bool useGpuForPme = simulationWork.useGpuPme;
+ const bool useGpuForNonbonded = simulationWork.useGpuNonbonded;
+ const bool useGpuForBufferOps = simulationWork.useGpuBufferOps;
+ const bool useGpuForUpdate = simulationWork.useGpuUpdate;
+
/* Check for polarizable models and flexible constraints */
shellfc = init_shell_flexcon(fplog, top_global, constr ? constr->numFlexibleConstraints() : 0,
- ir->nstcalcenergy, DOMAINDECOMP(cr));
+ ir->nstcalcenergy, DOMAINDECOMP(cr), useGpuForPme);
{
double io = compute_io(ir, top_global->natoms, *groups, energyOutput.numEnergyTerms(), 1);
std::unique_ptr<t_state> stateInstance;
t_state* state;
+ gmx_localtop_t top(top_global->ffparams);
auto mdatoms = mdAtoms->mdatoms();
- std::unique_ptr<UpdateConstrainCuda> integrator;
-
+ ForceBuffers f(fr->useMts, ((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
+ ? PinningPolicy::PinnedIfSupported
+ : PinningPolicy::CannotBePinned);
if (DOMAINDECOMP(cr))
{
- dd_init_local_top(*top_global, &top);
-
stateInstance = std::make_unique<t_state>();
state = stateInstance.get();
dd_init_local_state(cr->dd, state_global, state);
else
{
state_change_natoms(state_global, state_global->natoms);
- f.resizeWithPadding(state_global->natoms);
/* Copy the pointer to the global state */
state = state_global;
/* Generate and initialize new topology */
- mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &graph, mdAtoms, constr, vsite, shellfc);
+ mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &f, mdAtoms, constr, vsite, shellfc);
upd.setNumAtoms(state->natoms);
}
- const auto& simulationWork = runScheduleWork->simulationWork;
- const bool useGpuForPme = simulationWork.useGpuPme;
- const bool useGpuForNonbonded = simulationWork.useGpuNonbonded;
- const bool useGpuForBufferOps = simulationWork.useGpuBufferOps;
- const bool useGpuForUpdate = simulationWork.useGpuUpdate;
+ std::unique_ptr<UpdateConstrainGpu> integrator;
StatePropagatorDataGpu* stateGpu = fr->stateGpu;
+ // TODO: the assertions below should be handled by UpdateConstraintsBuilder.
if (useGpuForUpdate)
{
GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr) || ddUsesUpdateGroups(*cr->dd) || constr == nullptr
GMX_RELEASE_ASSERT(
ir->etc != etcNOSEHOOVER,
"Nose-Hoover temperature coupling is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(ir->epc == epcNO || ir->epc == epcPARRINELLORAHMAN || ir->epc == epcBERENDSEN,
- "Only Parrinello-Rahman and Berendsen pressure coupling are supported "
- "with the GPU update.\n");
+ GMX_RELEASE_ASSERT(
+ ir->epc == epcNO || ir->epc == epcPARRINELLORAHMAN || ir->epc == epcBERENDSEN
+ || ir->epc == epcCRESCALE,
+ "Only Parrinello-Rahman, Berendsen, and C-rescale pressure coupling are supported "
+ "with the GPU update.\n");
GMX_RELEASE_ASSERT(!mdatoms->haveVsites,
"Virtual sites are not supported with the GPU update.\n");
GMX_RELEASE_ASSERT(ed == nullptr,
"Essential dynamics is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(!ir->bPull || !pull_have_constraint(ir->pull),
+ GMX_RELEASE_ASSERT(!ir->bPull || !pull_have_constraint(*ir->pull),
"Constraints pulling is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(fcd->orires.nr == 0,
+ GMX_RELEASE_ASSERT(fcdata.orires->nr == 0,
"Orientation restraints are not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(ir->efep == efepNO,
- "Free energy perturbations are not supported with the GPU update.");
- GMX_RELEASE_ASSERT(graph == nullptr, "The graph is not supported with GPU update.");
+ GMX_RELEASE_ASSERT(
+ ir->efep == efepNO
+ || (!haveFepPerturbedMasses(*top_global) && !havePerturbedConstraints(*top_global)),
+ "Free energy perturbation of masses and constraints are not supported with the GPU "
+ "update.");
if (constr != nullptr && constr->numConstraintsTotal() > 0)
{
{
GMX_LOG(mdlog.info).asParagraph().appendText("Updating coordinates on the GPU.");
}
- integrator = std::make_unique<UpdateConstrainCuda>(
- *ir, *top_global, stateGpu->getUpdateStream(), stateGpu->xUpdatedOnDevice());
-
- t_pbc pbc;
- set_pbc(&pbc, epbcXYZ, state->box);
- integrator->setPbc(&pbc);
+ GMX_RELEASE_ASSERT(fr->deviceStreamManager != nullptr,
+ "Device stream manager should be initialized in order to use GPU "
+ "update-constraints.");
+ GMX_RELEASE_ASSERT(
+ fr->deviceStreamManager->streamIsValid(gmx::DeviceStreamType::UpdateAndConstraints),
+ "Update stream should be initialized in order to use GPU "
+ "update-constraints.");
+ integrator = std::make_unique<UpdateConstrainGpu>(
+ *ir, *top_global, fr->deviceStreamManager->context(),
+ fr->deviceStreamManager->stream(gmx::DeviceStreamType::UpdateAndConstraints),
+ stateGpu->xUpdatedOnDevice(), wcycle);
+
+ integrator->setPbc(PbcType::Xyz, state->box);
}
if (useGpuForPme || (useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
{
changePinningPolicy(&state->x, PinningPolicy::PinnedIfSupported);
}
- if ((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
- {
- changePinningPolicy(&f, PinningPolicy::PinnedIfSupported);
- }
if (useGpuForUpdate)
{
changePinningPolicy(&state->v, PinningPolicy::PinnedIfSupported);
if (MASTER(cr))
{
- EnergyElement::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput);
+ EnergyData::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput);
}
- preparePrevStepPullCom(ir, pull_work, mdatoms, state, state_global, cr,
+ preparePrevStepPullCom(ir, pull_work, mdatoms->massT, state, state_global, cr,
startingBehavior != StartingBehavior::NewSimulation);
// TODO: Remove this by converting AWH into a ForceProvider
if (constr)
{
/* Constrain the initial coordinates and velocities */
- do_constrain_first(fplog, constr, ir, mdatoms, state->natoms, state->x.arrayRefWithPadding(),
- state->v.arrayRefWithPadding(), state->box, state->lambda[efptBONDED]);
+ do_constrain_first(fplog, constr, ir, mdatoms->nr, mdatoms->homenr,
+ state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(),
+ state->box, state->lambda[efptBONDED]);
}
if (vsite)
{
/* Construct the virtual sites for the initial configuration */
- construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr, top.idef.iparams,
- top.idef.il, fr->ePBC, fr->bMolPBC, cr, state->box);
+ vsite->construct(state->x, ir->delta_t, {}, state->box);
}
}
nstfep = ir->fepvals->nstdhdl;
if (ir->bExpanded)
{
- nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
+ nstfep = std::gcd(ir->expandedvals->nstexpanded, nstfep);
}
if (useReplicaExchange)
{
- nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
+ nstfep = std::gcd(replExParams.exchangeInterval, nstfep);
+ }
+ if (ir->bDoAwh)
+ {
+ nstfep = std::gcd(ir->awhParams->nstSampleCoord, nstfep);
}
}
bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
if (PAR(cr))
{
- gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr);
+ gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr->mpi_comm_mygroup);
}
if (hasReadEkinState)
{
cglo_flags_iteration |= CGLO_STOPCM;
cglo_flags_iteration &= ~CGLO_TEMPERATURE;
}
- compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
- state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, nullptr, enerd,
- force_vir, shake_vir, total_vir, pres, mu_tot, constr, &nullSignaller,
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, nullptr,
+ enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
cglo_flags_iteration
| (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
/* At initialization, do not pass x with acceleration-correction mode
* to avoid (incorrect) correction of the initial coordinates.
*/
- rvec* xPtr = nullptr;
- if (vcm.mode != ecmLINEAR_ACCELERATION_CORRECTION)
- {
- xPtr = state->x.rvec_array();
- }
- process_and_stopcm_grp(fplog, &vcm, *mdatoms, xPtr, state->v.rvec_array());
+ auto x = (vcm.mode == ecmLINEAR_ACCELERATION_CORRECTION) ? ArrayRef<RVec>()
+ : makeArrayRef(state->x);
+ process_and_stopcm_grp(fplog, &vcm, *mdatoms, x, makeArrayRef(state->v));
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
}
}
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global, &top,
- state->x.rvec_array(), state->box,
+ makeConstArrayRef(state->x), state->box,
&shouldCheckNumberOfBondedInteractions);
if (ir->eI == eiVVAK)
{
kinetic energy calculation. This minimized excess variables, but
perhaps loses some logic?*/
- compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
- state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, nullptr, enerd,
- force_vir, shake_vir, total_vir, pres, mu_tot, constr, &nullSignaller,
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, nullptr,
+ enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
state->box, nullptr, &bSumEkinhOld, cglo_flags & ~CGLO_PRESSURE);
}
const DDBalanceRegionHandler ddBalanceRegionHandler(cr);
- // TODO extract this to new multi-simulation module
if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
{
- if (!multisim_int_all_are_equal(ms, ir->nsteps))
- {
- GMX_LOG(mdlog.warning)
- .appendText(
- "Note: The number of steps is not consistent across multi "
- "simulations,\n"
- "but we are proceeding anyway!");
- }
- if (!multisim_int_all_are_equal(ms, ir->init_step))
- {
- if (simulationsShareState)
- {
- if (MASTER(cr))
- {
- gmx_fatal(FARGS,
- "The initial step is not consistent across multi simulations which "
- "share the state");
- }
- gmx_barrier(cr);
- }
- else
- {
- GMX_LOG(mdlog.warning)
- .appendText(
- "Note: The initial step is not consistent across multi "
- "simulations,\n"
- "but we are proceeding anyway!");
- }
- }
+ logInitialMultisimStatus(ms, cr, mdlog, simulationsShareState, ir->nsteps, ir->init_step);
}
/* and stop now if we should */
if (ir->efep != efepNO || ir->bSimTemp)
{
/* find and set the current lambdas */
- setCurrentLambdasLocal(step, ir->fepvals, lam0, state->lambda, state->fep_state);
+ state->lambda = currentLambdas(step, *(ir->fepvals), state->fep_state);
bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
/* Correct the new box if it is too skewed */
if (inputrecDynamicBox(ir))
{
- if (correct_box(fplog, step, state->box, graph))
+ if (correct_box(fplog, step, state->box))
{
bMasterState = TRUE;
// If update is offloaded, it should be informed about the box size change
if (useGpuForUpdate)
{
- t_pbc pbc;
- set_pbc(&pbc, epbcXYZ, state->box);
- integrator->setPbc(&pbc);
+ integrator->setPbc(PbcType::Xyz, state->box);
}
}
}
}
}
+ // Allocate or re-size GPU halo exchange object, if necessary
+ if (bNS && havePPDomainDecomposition(cr) && simulationWork.useGpuHaloExchange)
+ {
+ GMX_RELEASE_ASSERT(fr->deviceStreamManager != nullptr,
+ "GPU device manager has to be initialized to use GPU "
+ "version of halo exchange.");
+ constructGpuHaloExchange(mdlog, *cr, *fr->deviceStreamManager, wcycle);
+ }
+
if (MASTER(cr) && do_log)
{
- energyOutput.printHeader(fplog, step, t); /* can we improve the information printed here? */
+ gmx::EnergyOutput::printHeader(fplog, step,
+ t); /* can we improve the information printed here? */
}
if (ir->efep != efepNO)
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
- state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, wcycle, enerd,
- nullptr, nullptr, nullptr, nullptr, mu_tot, constr, &nullSignaller,
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
+ enerd, nullptr, nullptr, nullptr, nullptr, constr, &nullSignaller,
state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global,
- &top, state->x.rvec_array(), state->box,
+ &top, makeConstArrayRef(state->x), state->box,
&shouldCheckNumberOfBondedInteractions);
}
clear_mat(force_vir);
force_flags = (GMX_FORCE_STATECHANGED | ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0)
| GMX_FORCE_ALLFORCES | (bCalcVir ? GMX_FORCE_VIRIAL : 0)
| (bCalcEner ? GMX_FORCE_ENERGY : 0) | (bDoFEP ? GMX_FORCE_DHDL : 0));
+ if (fr->useMts && !do_per_step(step, ir->nstfout))
+ {
+ force_flags |= GMX_FORCE_DO_NOT_NEED_NORMAL_FORCE;
+ }
if (shellfc)
{
/* Now is the time to relax the shells */
relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir,
- imdSession, pull_work, bNS, force_flags, &top, constr, enerd, fcd,
+ imdSession, pull_work, bNS, force_flags, &top, constr, enerd,
state->natoms, state->x.arrayRefWithPadding(),
- state->v.arrayRefWithPadding(), state->box, state->lambda, &state->hist,
- f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, graph,
- shellfc, fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler);
+ state->v.arrayRefWithPadding(), state->box, state->lambda,
+ &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc,
+ fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler);
}
else
{
*/
do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession, pull_work, step,
nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist,
- f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, graph,
- fr, runScheduleWork, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
+ &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork,
+ vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
(bNS ? GMX_FORCE_NS : 0) | force_flags, ddBalanceRegionHandler);
}
trotter_seq, ettTSEQ1);
}
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
- etrtVELOCITY1, cr, constr);
+ upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
+ M, etrtVELOCITY1, cr, constr != nullptr);
wallcycle_stop(wcycle, ewcUPDATE);
- constrain_velocities(step, nullptr, state, shake_vir, constr, bCalcVir, do_log, do_ene);
+ constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir);
wallcycle_start(wcycle, ewcUPDATE);
/* if VV, compute the pressure and constraints */
/* For VV2, we strictly only need this if using pressure
if (bGStat || do_per_step(step - 1, nstglobalcomm))
{
wallcycle_stop(wcycle, ewcUPDATE);
- compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
- state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, wcycle, enerd,
- force_vir, shake_vir, total_vir, pres, mu_tot, constr, &nullSignaller,
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
+ enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
(bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
| (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
EkinAveVel because it's needed for the pressure */
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
- top_global, &top, state->x.rvec_array(), state->box,
- &shouldCheckNumberOfBondedInteractions);
+ top_global, &top, makeConstArrayRef(state->x),
+ state->box, &shouldCheckNumberOfBondedInteractions);
if (bStopCM)
{
- process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(),
- state->v.rvec_array());
+ process_and_stopcm_grp(fplog, &vcm, *mdatoms, makeArrayRef(state->x),
+ makeArrayRef(state->v));
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
}
wallcycle_start(wcycle, ewcUPDATE);
copy_mat(shake_vir, state->svir_prev);
copy_mat(force_vir, state->fvir_prev);
}
- if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
+ if ((inputrecNptTrotter(ir) || inputrecNvtTrotter(ir)) && ir->eI == eiVV)
{
/* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
enerd->term[F_TEMP] =
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(),
- state->v.rvec_array(), state->box, state->lambda[efptVDW],
- mdatoms, nrnb, &vcm, wcycle, enerd, nullptr, nullptr, nullptr,
- nullptr, mu_tot, constr, &nullSignaller, state->box, nullptr,
- &bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
+ enerd, nullptr, nullptr, nullptr, nullptr, constr, &nullSignaller,
+ state->box, nullptr, &bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
wallcycle_start(wcycle, ewcUPDATE);
}
}
{
saved_conserved_quantity -= enerd->term[F_DISPCORR];
}
- /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
+ /* sum up the foreign kinetic energy and dK/dl terms for vv. currently done every step so that dhdl is correct in the .edr */
if (ir->efep != efepNO)
{
- sum_dhdl(enerd, state->lambda, *ir->fepvals);
+ accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
}
}
if (runScheduleWork->stepWork.useGpuFBufferOps && (simulationWork.useGpuUpdate && !vsite)
&& do_per_step(step, ir->nstfout))
{
- stateGpu->copyForcesFromGpu(ArrayRef<RVec>(f), AtomLocality::Local);
+ stateGpu->copyForcesFromGpu(f.view().force(), AtomLocality::Local);
stateGpu->waitForcesReadyOnHost(AtomLocality::Local);
}
/* Now we have the energies and forces corresponding to the
* the update.
*/
do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, state_global,
- observablesHistory, top_global, fr, outf, energyOutput, ekind, f,
- checkpointHandler->isCheckpointingStep(), bRerunMD, bLastStep,
- mdrunOptions.writeConfout, bSumEkinhOld);
+ observablesHistory, top_global, fr, outf, energyOutput, ekind,
+ f.view().force(), checkpointHandler->isCheckpointingStep(),
+ bRerunMD, bLastStep, mdrunOptions.writeConfout, bSumEkinhOld);
/* Check if IMD step and do IMD communication, if bIMD is TRUE. */
bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x.rvec_array(), t);
/* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
if (constr && bIfRandomize)
{
- constrain_velocities(step, nullptr, state, tmp_vir, constr, bCalcVir, do_log, do_ene);
+ constrain_velocities(constr, do_log, do_ene, step, state, nullptr, false, nullptr);
}
}
/* Box is changed in update() when we do pressure coupling,
dvdl_constr = 0;
- wallcycle_start(wcycle, ewcUPDATE);
+ if (!useGpuForUpdate)
+ {
+ wallcycle_start(wcycle, ewcUPDATE);
+ }
/* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
if (bTrotter)
{
if (EI_VV(ir->eI))
{
/* velocity half-step update */
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
- etrtVELOCITY2, cr, constr);
+ upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
+ M, etrtVELOCITY2, cr, constr != nullptr);
}
/* Above, initialize just copies ekinh into ekin,
stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
}
- // If the buffer ops were not offloaded this step, the forces are on the host and have to be copied
- if (!runScheduleWork->stepWork.useGpuFBufferOps)
+ if (simulationWork.useGpuPme && !runScheduleWork->simulationWork.useGpuPmePpCommunication
+ && !thisRankHasDuty(cr, DUTY_PME))
+ {
+ // The PME forces were recieved to the host, so have to be copied
+ stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::All);
+ }
+ else if (!runScheduleWork->stepWork.useGpuFBufferOps)
{
- stateGpu->copyForcesToGpu(ArrayRef<RVec>(f), AtomLocality::Local);
+ // The buffer ops were not offloaded this step, so the forces are on the
+ // host and have to be copied
+ stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::Local);
}
const bool doTemperatureScaling =
}
else
{
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
- etrtPOSITION, cr, constr);
+ /* With multiple time stepping we need to do an additional normal
+ * update step to obtain the virial, as the actual MTS integration
+ * using an acceleration where the slow forces are multiplied by mtsFactor.
+ * Using that acceleration would result in a virial with the slow
+ * force contribution would be a factor mtsFactor too large.
+ */
+ if (fr->useMts && bCalcVir && constr != nullptr)
+ {
+ upd.update_for_constraint_virial(*ir, *mdatoms, *state, f.view().forceWithPadding(), *ekind);
+
+ constrain_coordinates(constr, do_log, do_ene, step, state,
+ upd.xp()->arrayRefWithPadding(), &dvdl_constr, bCalcVir, shake_vir);
+ }
+
+ ArrayRefWithPadding<const RVec> forceCombined =
+ (fr->useMts && step % ir->mtsLevels[1].stepFactor == 0)
+ ? f.view().forceMtsCombinedWithPadding()
+ : f.view().forceWithPadding();
+ upd.update_coords(*ir, step, mdatoms, state, forceCombined, fcdata, ekind, M,
+ etrtPOSITION, cr, constr != nullptr);
wallcycle_stop(wcycle, ewcUPDATE);
- constrain_coordinates(step, &dvdl_constr, state, shake_vir, &upd, constr, bCalcVir,
- do_log, do_ene);
+ constrain_coordinates(constr, do_log, do_ene, step, state, upd.xp()->arrayRefWithPadding(),
+ &dvdl_constr, bCalcVir && !fr->useMts, shake_vir);
- update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state, cr, nrnb, wcycle, &upd,
- constr, do_log, do_ene);
- finish_update(ir, mdatoms, state, graph, nrnb, wcycle, &upd, constr);
+ upd.update_sd_second_half(*ir, step, &dvdl_constr, mdatoms, state, cr, nrnb, wcycle,
+ constr, do_log, do_ene);
+ upd.finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
}
if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
{
/* erase F_EKIN and F_TEMP here? */
/* just compute the kinetic energy at the half step to perform a trotter step */
- compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
- state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, wcycle, enerd,
- force_vir, shake_vir, total_vir, pres, mu_tot, constr, &nullSignaller, lastbox,
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle, enerd,
+ force_vir, shake_vir, total_vir, pres, constr, &nullSignaller, lastbox,
nullptr, &bSumEkinhOld, (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE);
wallcycle_start(wcycle, ewcUPDATE);
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
/* now we know the scaling, we can compute the positions again */
std::copy(cbuf.begin(), cbuf.end(), state->x.begin());
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
- etrtPOSITION, cr, constr);
+ upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
+ M, etrtPOSITION, cr, constr != nullptr);
wallcycle_stop(wcycle, ewcUPDATE);
/* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
* to numerical errors, or are they important
* physically? I'm thinking they are just errors, but not completely sure.
* For now, will call without actually constraining, constr=NULL*/
- finish_update(ir, mdatoms, state, graph, nrnb, wcycle, &upd, nullptr);
+ upd.finish_update(*ir, mdatoms, state, wcycle, false);
}
if (EI_VV(ir->eI))
{
/* this factor or 2 correction is necessary
because half of the constraint force is removed
in the vv step, so we have to double it. See
- the Redmine issue #1255. It is not yet clear
+ the Issue #1255. It is not yet clear
if the factor of 2 is exact, or just a very
good approximation, and this will be
investigated. The next step is to see if this
if (vsite != nullptr)
{
wallcycle_start(wcycle, ewcVSITECONSTR);
- if (graph != nullptr)
- {
- shift_self(graph, state->box, state->x.rvec_array());
- }
- construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
- top.idef.iparams, top.idef.il, fr->ePBC, fr->bMolPBC, cr, state->box);
-
- if (graph != nullptr)
- {
- unshift_self(graph, state->box, state->x.rvec_array());
- }
+ vsite->construct(state->x, ir->delta_t, state->v, state->box);
wallcycle_stop(wcycle, ewcVSITECONSTR);
}
bool doIntraSimSignal = true;
SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
- compute_globals(
- gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
- state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, wcycle, enerd,
- force_vir, shake_vir, total_vir, pres, mu_tot, constr, &signaller, lastbox,
- &totalNumberOfBondedInteractions, &bSumEkinhOld,
- (bGStat ? CGLO_GSTAT : 0) | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
- | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
- | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
- | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0) | CGLO_CONSTRAINT
- | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
- : 0));
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm,
+ wcycle, enerd, force_vir, shake_vir, total_vir, pres, constr,
+ &signaller, lastbox, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ (bGStat ? CGLO_GSTAT : 0) | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
+ | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
+ | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
+ | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0) | CGLO_CONSTRAINT
+ | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
+ : 0));
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
- top_global, &top, state->x.rvec_array(), state->box,
- &shouldCheckNumberOfBondedInteractions);
+ top_global, &top, makeConstArrayRef(state->x),
+ state->box, &shouldCheckNumberOfBondedInteractions);
if (!EI_VV(ir->eI) && bStopCM)
{
- process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(),
- state->v.rvec_array());
+ process_and_stopcm_grp(fplog, &vcm, *mdatoms, makeArrayRef(state->x),
+ makeArrayRef(state->v));
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
// TODO: The special case of removing CM motion should be dealt more gracefully
if (ir->efep != efepNO && !EI_VV(ir->eI))
{
- /* Sum up the foreign energy and dhdl terms for md and sd.
- Currently done every step so that dhdl is correct in the .edr */
- sum_dhdl(enerd, state->lambda, *ir->fepvals);
+ /* Sum up the foreign energy and dK/dl terms for md and sd.
+ Currently done every step so that dH/dl is correct in the .edr */
+ accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
}
update_pcouple_after_coordinates(fplog, step, ir, mdatoms, pres, force_vir, shake_vir,
- pressureCouplingMu, state, nrnb, &upd, !useGpuForUpdate);
+ pressureCouplingMu, state, nrnb, upd.deform(), !useGpuForUpdate);
const bool doBerendsenPressureCoupling =
(inputrec->epc == epcBERENDSEN && do_per_step(step, inputrec->nstpcouple));
- if (useGpuForUpdate && (doBerendsenPressureCoupling || doParrinelloRahman))
+ const bool doCRescalePressureCoupling =
+ (inputrec->epc == epcCRESCALE && do_per_step(step, inputrec->nstpcouple));
+ if (useGpuForUpdate
+ && (doBerendsenPressureCoupling || doCRescalePressureCoupling || doParrinelloRahman))
{
integrator->scaleCoordinates(pressureCouplingMu);
- t_pbc pbc;
- set_pbc(&pbc, epbcXYZ, state->box);
- integrator->setPbc(&pbc);
+ if (doCRescalePressureCoupling)
+ {
+ matrix pressureCouplingInvMu;
+ gmx::invertBoxMatrix(pressureCouplingMu, pressureCouplingInvMu);
+ integrator->scaleVelocities(pressureCouplingInvMu);
+ }
+ integrator->setPbc(PbcType::Xyz, state->box);
}
/* ################# END UPDATE STEP 2 ################# */
}
if (bCalcEner)
{
- energyOutput.addDataAtEnergyStep(bDoDHDL, bCalcEnerStep, t, mdatoms->tmass, enerd, state,
- ir->fepvals, ir->expandedvals, lastbox, shake_vir,
- force_vir, total_vir, pres, ekind, mu_tot, constr);
+ energyOutput.addDataAtEnergyStep(
+ bDoDHDL, bCalcEnerStep, t, mdatoms->tmass, enerd, ir->fepvals,
+ ir->expandedvals, lastbox,
+ PTCouplingArrays{ state->boxv, state->nosehoover_xi, state->nosehoover_vxi,
+ state->nhpres_xi, state->nhpres_vxi },
+ state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr);
}
else
{
if (doSimulatedAnnealing)
{
- energyOutput.printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
+ gmx::EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups,
+ &(ir->opts));
}
if (do_log || do_ene || do_dr || do_or)
{
energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
- do_log ? fplog : nullptr, step, t, fcd, awh.get());
+ do_log ? fplog : nullptr, step, t,
+ fr->fcdata.get(), awh.get());
+ }
+ if (do_log && ir->bDoAwh && awh->hasFepLambdaDimension())
+ {
+ const bool isInitialOutput = false;
+ printLambdaStateToLog(fplog, state->lambda, isInitialOutput);
}
if (ir->bPull)
/* Gets written into the state at the beginning of next loop*/
state->fep_state = lamnew;
}
+ else if (ir->bDoAwh && awh->needForeignEnergyDifferences(step))
+ {
+ state->fep_state = awh->fepLambdaState();
+ }
/* Print the remaining wall clock time for the run */
if (isMasterSimMasterRank(ms, MASTER(cr)) && (do_verbose || gmx_got_usr_signal()) && !bPMETunePrinting)
{
{
if (ir->nstcalcenergy > 0)
{
- energyOutput.printAnnealingTemperatures(fplog, groups, &(ir->opts));
+ energyOutput.printEnergyConservation(fplog, ir->simulation_part, EI_MD(ir->eI));
+
+ gmx::EnergyOutput::printAnnealingTemperatures(fplog, groups, &(ir->opts));
energyOutput.printAverages(fplog, groups);
}
}
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
* Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
-- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
++ * Copyright (c) 2018,2019,2020,2021, 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.
#include "gromacs/mdtypes/forceoutput.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/mtop_lookup.h"
static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
{
- if (params.nat <= 1)
+ if (params.ind.size() <= 1)
{
/* no PBC required */
return epgrppbcNONE;
static void apply_forces_grp_part(const pull_group_work_t* pgrp,
int ind_start,
int ind_end,
- const t_mdatoms* md,
+ const real* masses,
const dvec f_pull,
int sign,
rvec* f)
for (int i = ind_start; i < ind_end; i++)
{
int ii = localAtomIndices[i];
- double wmass = md->massT[ii];
+ double wmass = masses[ii];
if (!pgrp->localWeights.empty())
{
wmass *= pgrp->localWeights[i];
/* Apply forces in a mass weighted fashion */
static void apply_forces_grp(const pull_group_work_t* pgrp,
- const t_mdatoms* md,
+ const real* masses,
const dvec f_pull,
int sign,
rvec* f,
{
auto localAtomIndices = pgrp->atomSet.localIndex();
- if (pgrp->params.nat == 1 && pgrp->atomSet.numAtomsLocal() == 1)
+ if (pgrp->params.ind.size() == 1 && pgrp->atomSet.numAtomsLocal() == 1)
{
/* Only one atom and our rank has this atom: we can skip
* the mass weighting, which means that this code also works
{
if (localAtomIndices.size() <= c_pullMaxNumLocalAtomsSingleThreaded)
{
- apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), md, f_pull, sign, f);
+ apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), masses, f_pull, sign, f);
}
else
{
{
int ind_start = (localAtomIndices.size() * (th + 0)) / nthreads;
int ind_end = (localAtomIndices.size() * (th + 1)) / nthreads;
- apply_forces_grp_part(pgrp, ind_start, ind_end, md, f_pull, sign, f);
+ apply_forces_grp_part(pgrp, ind_start, ind_end, masses, f_pull, sign, f);
}
}
}
/* Apply forces in a mass weighted fashion to a cylinder group */
static void apply_forces_cyl_grp(const pull_group_work_t* pgrp,
const double dv_corr,
- const t_mdatoms* md,
+ const real* masses,
const dvec f_pull,
double f_scal,
int sign,
continue;
}
int ii = localAtomIndices[i];
- double mass = md->massT[ii];
+ double mass = masses[ii];
/* The stored axial distance from the cylinder center (dv) needs
* to be corrected for an offset (dv_corr), which was unknown when
* we calculated dv.
*/
static void apply_forces_vec_torque(const struct pull_t* pull,
const pull_coord_work_t* pcrd,
- const t_mdatoms* md,
+ const real* masses,
rvec* f)
{
const PullCoordSpatialData& spatialData = pcrd->spatialData;
}
/* Apply the force to the groups defining the vector using opposite signs */
- apply_forces_grp(&pull->group[pcrd->params.group[2]], md, f_perp, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd->params.group[3]], md, f_perp, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd->params.group[2]], masses, f_perp, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd->params.group[3]], masses, f_perp, 1, f, pull->nthreads);
}
/* Apply forces in a mass weighted fashion */
static void apply_forces_coord(struct pull_t* pull,
int coord,
const PullCoordVectorForces& forces,
- const t_mdatoms* md,
+ const real* masses,
rvec* f)
{
/* Here it would be more efficient to use one large thread-parallel
if (pcrd.params.eGeom == epullgCYL)
{
- apply_forces_cyl_grp(&pull->dyna[coord], pcrd.spatialData.cyl_dev, md, forces.force01,
+ apply_forces_cyl_grp(&pull->dyna[coord], pcrd.spatialData.cyl_dev, masses, forces.force01,
pcrd.scalarForce, -1, f, pull->nthreads);
/* Sum the force along the vector and the radial force */
{
f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
}
- apply_forces_grp(&pull->group[pcrd.params.group[1]], md, f_tot, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, f_tot, 1, f, pull->nthreads);
}
else
{
/* We need to apply the torque forces to the pull groups
* that define the pull vector.
*/
- apply_forces_vec_torque(pull, &pcrd, md, f);
+ apply_forces_vec_torque(pull, &pcrd, masses, f);
}
- if (pull->group[pcrd.params.group[0]].params.nat > 0)
+ if (!pull->group[pcrd.params.group[0]].params.ind.empty())
{
- apply_forces_grp(&pull->group[pcrd.params.group[0]], md, forces.force01, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[0]], masses, forces.force01, -1, f,
+ pull->nthreads);
}
- apply_forces_grp(&pull->group[pcrd.params.group[1]], md, forces.force01, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, forces.force01, 1, f, pull->nthreads);
if (pcrd.params.ngroup >= 4)
{
- apply_forces_grp(&pull->group[pcrd.params.group[2]], md, forces.force23, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd.params.group[3]], md, forces.force23, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[2]], masses, forces.force23, -1, f,
+ pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[3]], masses, forces.force23, 1, f,
+ pull->nthreads);
}
if (pcrd.params.ngroup >= 6)
{
- apply_forces_grp(&pull->group[pcrd.params.group[4]], md, forces.force45, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd.params.group[5]], md, forces.force45, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[4]], masses, forces.force45, -1, f,
+ pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[5]], masses, forces.force45, 1, f,
+ pull->nthreads);
}
}
}
/* This function returns the distance based on coordinates xg and xref.
* Note that the pull coordinate struct pcrd is not modified.
+ *
+ * \param[in] pull The pull struct
+ * \param[in] pcrd The pull coordinate to compute a distance for
+ * \param[in] pbc The periodic boundary conditions
+ * \param[in] xg The coordinate of group 1
+ * \param[in] xref The coordinate of group 0
+ * \param[in] groupIndex0 The index of group 0 in the pcrd->params.group
+ * \param[in] groupIndex1 The index of group 1 in the pcrd->params.group
+ * \param[in] max_dist2 The maximum distance squared
+ * \param[out] dr The distance vector
*/
static void low_get_pull_coord_dr(const struct pull_t* pull,
const pull_coord_work_t* pcrd,
const t_pbc* pbc,
- dvec xg,
+ const dvec xg,
dvec xref,
- double max_dist2,
+ const int groupIndex0,
+ const int groupIndex1,
+ const double max_dist2,
dvec dr)
{
const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
/* Only the first group can be an absolute reference, in that case nat=0 */
- if (pgrp0->params.nat == 0)
+ if (pgrp0->params.ind.empty())
{
for (int m = 0; m < DIM; m++)
{
gmx_fatal(FARGS,
"Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
"box size (%f).\n%s",
- pcrd->params.group[0], pcrd->params.group[1], sqrt(dr2),
+ pcrd->params.group[groupIndex0], pcrd->params.group[groupIndex1], sqrt(dr2),
sqrt(0.98 * 0.98 * max_dist2), pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
}
pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
low_get_pull_coord_dr(pull, pcrd, pbc, pgrp1->x,
- pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, md2,
- spatialData.dr01);
+ pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, 0,
+ 1, md2, spatialData.dr01);
if (pcrd->params.ngroup >= 4)
{
pgrp2 = &pull->group[pcrd->params.group[2]];
pgrp3 = &pull->group[pcrd->params.group[3]];
- low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, md2, spatialData.dr23);
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, 2, 3, md2, spatialData.dr23);
}
if (pcrd->params.ngroup >= 6)
{
pgrp4 = &pull->group[pcrd->params.group[4]];
pgrp5 = &pull->group[pcrd->params.group[5]];
- low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, md2, spatialData.dr45);
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, 4, 5, md2, spatialData.dr45);
}
}
/* Get the current difference vector */
low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->params.group[1]],
- rnew[pcrd->params.group[0]], -1, unc_ij);
+ rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
if (debug)
{
g0 = pcrd->params.group[0];
g1 = pcrd->params.group[1];
- low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
- low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
+ low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
+ low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", rnew[g0][0],
rnew[g0][1], rnew[g0][2], rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
fprintf(debug, "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", "", "", "", "", "",
}
low_get_pull_coord_dr(pull, &coord, pbc, rnew[coord.params.group[1]],
- rnew[coord.params.group[0]], -1, unc_ij);
+ rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
switch (coord.params.eGeom)
{
provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
}
- GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr,
+ GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
"The external potential provider string for a pull coordinate is NULL");
- if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider) != 0)
+ if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
{
gmx_fatal(FARGS,
"Module '%s' attempted to register an external potential for pull coordinate %d "
"which expects the external potential to be provided by a module named '%s'",
- provider, coord_index + 1, pcrd->params.externalPotentialProvider);
+ provider, coord_index + 1, pcrd->params.externalPotentialProvider.c_str());
}
/* Lock to avoid (extremely unlikely) simultaneous reading and writing of
"pull coordinates. The first coordinate without provider is number %zu, which "
"expects a module named '%s' to provide the external potential.",
pull->numUnregisteredExternalPotentials, c + 1,
- pull->coord[c].params.externalPotentialProvider);
+ pull->coord[c].params.externalPotentialProvider.c_str());
}
}
void apply_external_pull_coord_force(struct pull_t* pull,
int coord_index,
double coord_force,
- const t_mdatoms* mdatoms,
+ const real* masses,
gmx::ForceWithVirial* forceWithVirial)
{
pull_coord_work_t* pcrd;
forceWithVirial->addVirialContribution(virial);
}
- apply_forces_coord(pull, coord_index, pullCoordForces, mdatoms,
+ apply_forces_coord(pull, coord_index, pullCoordForces, masses,
as_rvec_array(forceWithVirial->force_.data()));
}
}
real pull_potential(struct pull_t* pull,
- const t_mdatoms* md,
+ const real* masses,
t_pbc* pbc,
const t_commrec* cr,
double t,
{
real dVdl = 0;
- pull_calc_coms(cr, pull, md, pbc, t, x, nullptr);
+ pull_calc_coms(cr, pull, masses, pbc, t, x, nullptr);
rvec* f = as_rvec_array(force->force_.data());
matrix virial = { { 0 } };
pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
/* Distribute the force over the atoms in the pulled groups */
- apply_forces_coord(pull, c, pullCoordForces, md, f);
+ apply_forces_coord(pull, c, pullCoordForces, masses, f);
}
if (MASTER(cr))
}
void pull_constraint(struct pull_t* pull,
- const t_mdatoms* md,
+ const real* masses,
t_pbc* pbc,
const t_commrec* cr,
double dt,
if (pull->comm.bParticipate)
{
- pull_calc_coms(cr, pull, md, pbc, t, x, xp);
+ pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
}
"date prev. COM "
"to bcast here as well as to e.g. checkpointing");
- gmx_bcast(sizeof(dvec), group.x_prev_step, cr);
+ gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
}
}
}
* But we still want to have the correct mass-weighted COMs.
* So we store the real masses in the weights.
*/
- const bool setWeights = (pg->params.nweight > 0 || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
+ const bool setWeights =
+ (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
/* In parallel, store we need to extract localWeights from weights at DD time */
std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
double wmass = 0;
double wwmass = 0;
int molb = 0;
- for (int i = 0; i < pg->params.nat; i++)
+ for (int i = 0; i < int(pg->params.ind.size()); i++)
{
int ii = pg->params.ind[i];
if (bConstraint && ir->opts.nFreeze)
m = (1 - lambda) * atom.m + lambda * atom.mB;
}
real w;
- if (pg->params.nweight > 0)
+ if (!pg->params.weight.empty())
{
w = pg->params.weight[i];
}
/* We can have single atom groups with zero mass with potential pulling
* without cosine weighting.
*/
- if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
+ if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
{
/* With one atom the mass doesn't matter */
wwmass = 1;
else
{
gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
- pg->params.weight ? " weighted" : "", g);
+ !pg->params.weight.empty() ? " weighted" : "", g);
}
}
if (fplog)
{
- fprintf(fplog, "Pull group %d: %5d atoms, mass %9.3f", g, pg->params.nat, tmass);
- if (pg->params.weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
+ fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
+ if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
{
fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
}
int ndim = 0;
for (int d = 0; d < DIM; d++)
{
- ndim += pulldim_con[d] * pg->params.nat;
+ ndim += pulldim_con[d] * pg->params.ind.size();
}
if (fplog && nfrozen > 0 && nfrozen < ndim)
{
/* Copy the pull parameters */
pull->params = *pull_params;
- /* Avoid pointer copies */
- pull->params.group = nullptr;
- pull->params.coord = nullptr;
for (int i = 0; i < pull_params->ngroup; ++i)
{
- pull->group.emplace_back(pull_params->group[i],
- atomSets->add({ pull_params->group[i].ind,
- pull_params->group[i].ind + pull_params->group[i].nat }),
+ pull->group.emplace_back(pull_params->group[i], atomSets->add(pull_params->group[i].ind),
pull_params->bSetPbcRefToPrevStepCOM);
}
pull->bXOutAverage = pull_params->bXOutAverage;
pull->bFOutAverage = pull_params->bFOutAverage;
- GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0,
+ GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
"pull group 0 is an absolute reference group and should not contain atoms");
pull->numCoordinatesWithExternalPotential = 0;
epull_names[epullUMBRELLA]);
}
+ GMX_RELEASE_ASSERT(
+ !ir->useMts,
+ "Constraint pulling can not be combined with multiple time stepping");
+
pull->bConstraint = TRUE;
}
else
pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
- pull->ePBC = ir->ePBC;
- switch (pull->ePBC)
+ pull->pbcType = ir->pbcType;
+ switch (pull->pbcType)
{
- case epbcNONE: pull->npbcdim = 0; break;
- case epbcXY: pull->npbcdim = 2; break;
+ case PbcType::No: pull->npbcdim = 0; break;
+ case PbcType::XY: pull->npbcdim = 2; break;
default: pull->npbcdim = 3; break;
}
bAbs = FALSE;
for (const pull_coord_work_t& coord : pull->coord)
{
- if (pull->group[coord.params.group[0]].params.nat == 0
- || pull->group[coord.params.group[1]].params.nat == 0)
+ if (pull->group[coord.params.group[0]].params.ind.empty()
+ || pull->group[coord.params.group[1]].params.ind.empty())
{
bAbs = TRUE;
}
// Don't include the reference group 0 in loop
for (size_t g = 1; g < pull->group.size(); g++)
{
- if (pull->group[g].params.nat > 1 && pull->group[g].params.pbcatom < 0)
+ if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
{
/* We are using cosine weighting */
fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
pull_group_work_t* pgrp;
pgrp = &pull->group[g];
- if (pgrp->params.nat > 0)
+ if (!pgrp->params.ind.empty())
{
/* There is an issue when a group is used in multiple coordinates
* and constraints are applied in different dimensions with atoms
{
case epgrppbcREFAT: pull->bRefAt = TRUE; break;
case epgrppbcCOS:
- if (pgrp->params.weight != nullptr)
+ if (!pgrp->params.weight.empty())
{
gmx_fatal(FARGS,
"Pull groups can not have relative weights and cosine weighting "
{
if (coord.params.eGeom == epullgCYL)
{
- if (pull->group[coord.params.group[0]].params.nat == 0)
+ if (pull->group[coord.params.group[0]].params.ind.empty())
{
gmx_fatal(FARGS,
"A cylinder pull group is not supported when using absolute "
void preparePrevStepPullCom(const t_inputrec* ir,
pull_t* pull_work,
- const t_mdatoms* md,
+ const real* masses,
t_state* state,
const t_state* state_global,
const t_commrec* cr,
if (PAR(cr))
{
/* Only the master rank has the checkpointed COM from the previous step */
- gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(), &state->pull_com_prev_step[0], cr);
+ gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(),
+ &state->pull_com_prev_step[0], cr->mpi_comm_mygroup);
}
setPrevStepPullComFromState(pull_work, state);
}
else
{
t_pbc pbc;
- set_pbc(&pbc, ir->ePBC, state->box);
- initPullComFromPrevStep(cr, pull_work, md, &pbc, state->x.rvec_array());
+ set_pbc(&pbc, ir->pbcType, state->box);
+ initPullComFromPrevStep(cr, pull_work, masses, &pbc, state->x.rvec_array());
updatePrevStepPullCom(pull_work, state);
}
}
destroy_pull(pull);
}
-gmx_bool pull_have_potential(const struct pull_t* pull)
+bool pull_have_potential(const pull_t& pull)
{
- return pull->bPotential;
+ return pull.bPotential;
}
-gmx_bool pull_have_constraint(const struct pull_t* pull)
+bool pull_have_constraint(const pull_t& pull)
{
- return pull->bConstraint;
+ return pull.bConstraint;
}
-bool pull_have_constraint(const pull_params_t* pullParameters)
+bool pull_have_constraint(const pull_params_t& pullParameters)
{
- if (pullParameters == nullptr)
- {
- return false;
- }
- for (int c = 0; c < pullParameters->ncoord; c++)
+ for (int c = 0; c < pullParameters.ncoord; c++)
{
- if (pullParameters->coord[c].eType == epullCONSTRAINT)
+ if (pullParameters.coord[c].eType == epullCONSTRAINT)
{
return true;
}
# This file is part of the GROMACS molecular simulation package.
#
# Copyright (c) 2012,2013,2014,2015,2016, The GROMACS development team.
--# Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
++# Copyright (c) 2017,2018,2019,2020,2021, 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.
set(REGRESSIONTEST_URL https://gitlab.com/gromacs/gromacs-regressiontests/-/archive/${REGRESSIONTEST_BRANCH}/gromacs-regressiontests-${REGRESSIONTEST_BRANCH}.tar.gz)
# REGRESSIONTEST_PATH for dev trees is set later based on the dirname found in the tar
else()
- set(REGRESSIONTEST_URL http://ftp.gromacs.org/regressiontests/regressiontests-${REGRESSIONTEST_VERSION}.tar.gz)
+ set(REGRESSIONTEST_URL https://ftp.gromacs.org/regressiontests/regressiontests-${REGRESSIONTEST_VERSION}.tar.gz)
set(REGRESSIONTEST_PATH
"${CMAKE_CURRENT_BINARY_DIR}/regressiontests-${REGRESSIONTEST_VERSION}"
CACHE PATH "Path to auto-downloaded regressiontests" FORCE)
# End copied from regression tests.
#
- if (NOT PYTHONINTERP_FOUND)
+ if (NOT Python3_Interpreter_FOUND)
message(FATAL_ERROR
"Python not found. Physical validation requires python. \
- Install python, set PYTHON_EXECUTABLE to a valid python location, \
+ Install python, set Python3_ROOT_DIR or PYTHON_EXECUTABLE to a valid location, \
or set GMX_PHYSICAL_VALIDATION=OFF to disable the physical validation tests.")
endif()
#
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff