# Test AMD OpenCL
# Test external clFFT (for build speed)
+# Test AVX2_256 SIMD
# Test newest gcc at time of release
# Test hwloc-1 support
gcc-9 openmp simd=avx2_256 gpuhw=amd opencl-1.2 clFFT-2.14 hwloc libhwloc-1.11.2 cmake-3.14.5
# Test with newest CUDA at time of release
gcc-5 gpuhw=nvidia cuda-10.0 release
+# Test with OpenCL support
+gcc-8 simd=avx2_256 gpuhw=amd opencl-1.2 release
+
# TODO items
# Avoid specifying cmake versions just to move jobs away from bs_nix-amd
-# Add an OpenCL GPU build
# build the regressiontests tarball with all the right naming. The
# naming affects the md5sum that has to go here, and if it isn't right
# release workflow will report a failure.
-set(REGRESSIONTEST_MD5SUM "162edf7d9d520672f5fa5888aae60989" CACHE INTERNAL "MD5 sum of the regressiontests tarball for this GROMACS version")
+set(REGRESSIONTEST_MD5SUM "639ae20f00e0311dcec8e5af7876abdb" CACHE INTERNAL "MD5 sum of the regressiontests tarball for this GROMACS version")
math(EXPR GMX_VERSION_NUMERIC
"${GMX_VERSION_MAJOR}*10000 + ${GMX_VERSION_PATCH}")
source $HOME/gromacs-gmxapi/bin/GMXRC
+.. _gmxapi venv:
+
Set up a Python virtual environment
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
command line tool or the :py:func:`help` interactive Python function, or refer to
the :doc:`pythonreference`.
-Python does not wrap any command-line tool, so once installation is complete,
-there shouldn't be any additional configuration necessary, and any errors that
-occur should be caught at the Python level. Any Python *exception* raised by gmxapi
+Any Python *exception* raised by gmxapi
should be descended from (and catchable as) :class:`gmxapi.exceptions.Error`.
+Additional status messages can be acquired through the :ref:`gmxapi logging`
+facility.
+Unfortunately, some errors occurring in the GROMACS library are not yet
+recoverable at the Python level, and much of the standard GROMACS terminal
+output is not yet accessible through Python.
+If you find a particularly problematic scenario, please file a GROMACS bug report.
-As an exception to the preceding paragraph, we have a tool specifically for
-wrapping arbitrary (unintegrated) command line tools: See :class:`gmxapi.commandline_operation`.
+During installation, the *gmxapi* Python package becomes tied to a specific
+GROMACS installation.
+If you would like to access multiple GROMACS installations
+from Python, build and install *gmxapi* in separate
+:ref:`virtual environments <gmxapi venv>`.
+
+In some cases *gmxapi* still needs help finding infrastructure from the
+GROMACS installation.
+For instance, :py:func:`gmxapi.commandline_operation` is not a pure API utility,
+but a wrapper for command line tools.
+Make sure that the command line tools you intend to use are discoverable in
+your :envvar:`PATH`, such as by "source"ing your :file:`GMXRC` before launching
+a *gmxapi* script.
+
+.. todo:: Get relevant GROMACS paths in Python environment.
+
+ :py:class:`gmxapi.commandline_operation` relies on the environment :envvar:`PATH`
+ to locate executables, including the :command:`gmx` wrapper binary.
+ Relates to `#2961 <https://redmine.gromacs.org/issues/2961>`__.
+
+.. _parallelism:
+
+Notes on parallelism and MPI
+============================
+
+When launching a *gmxapi* script in an MPI environment,
+such as with :command:`mpiexec` or :command:`srun`,
+you must help *gmxapi* detect the MPI environment by ensuring that :py:mod:`mpi4py`
+is loaded.
+Refer to :ref:`mpi_requirements` for more on installing :py:mod:`mpi4py`.
+
+Assuming you use :command:`mpiexec` to launch MPI jobs in your environment,
+run a *gmxapi* script on two ranks with something like the following.
+Note that it can be helpful to provide :command:`mpiexec` with the full path to
+the intended Python interpreter since new process environments are being created.
+
+::
+
+ mpiexec -n 2 `which python` -m mpi4py myscript.py
+
+*gmxapi* 0.1 has limited parallelism, but future versions will include seamless
+acceleration as integration improves with the GROMACS library and computing
+environment runtime resources.
+Currently, *gmxapi* and the GROMACS library do not have an effective way to
+share an MPI environment.
+Therefore, if you intend to run more than one simulation at a time, in parallel,
+in a *gmxapi* script, you should build GROMACS with *thread-MPI* instead of a
+standard MPI library.
+I.e. configure GROMACS with the CMake flag ``-DGMX_THREAD_MPI=ON``.
+Then, launch your *gmxapi* script with one MPI rank per node, and *gmxapi* will
+assign each (non-MPI) simulation to its own node, while keeping the full MPI
+environment available for use via :py:mod:`mpi4py`.
+
+Running simple simulations
+==========================
+
+Once the ``gmxapi`` package is installed, running simulations is easy with
+:py:func:`gmxapi.read_tpr`.
+
+::
+
+ import gmxapi as gmx
+ simulation_input = gmx.read_tpr(tpr_filename)
+ md = gmx.mdrun(simulation_input)
+
+Note that this sets up the work you want to perform, but does not immediately
+trigger execution. You can explicitly trigger execution with::
+
+ md.run()
+
+or you can let gmxapi automatically launch work in response to the data you
+request.
+
+The :py:func:`gmxapi.mdrun` operation produces a simulation trajectory output.
+You can use ``md.output.trajectory`` as input to other operations,
+or you can get the output directly by calling ``md.output.trajectory.result()``.
+If the simulation has not been run yet when ``result()`` is called,
+the simulation will be run before the function returns.
+
+Running ensemble simulations
+============================
+
+To run a batch of simulations, just pass an array of inputs.::
+
+ md = gmx.read_tpr([tpr_filename1, tpr_filename2, ...])
+ md.run()
+
+Make sure to launch the script in an MPI environment with a sufficient number
+of ranks to allow one rank per simulation.
+
+For *gmxapi* 0.1, we recommend configuring the GROMACS build with
+``GMX_THREAD_MPI=ON`` and allowing one rank per node in order to allow each
+simulation ensemble member to run on a separate node.
+
+.. seealso:: :ref:`parallelism`
+
+.. _commandline:
+
+Accessing command line tools
+============================
+
+In *gmxapi* 0.1, most GROMACS tools are not yet exposed as *gmxapi* Python operations.
+:class:`gmxapi.commandline_operation` provides a way to convert a :command:`gmx`
+(or other) command line tool into an operation that can be used in a *gmxapi*
+script.
+
+In order to establish data dependencies, input and output files need to be
+indicated with the ``input_files`` and ``output_files`` parameters.
+``input_files`` and ``output_files`` key word arguments are dictionaries
+consisting of files keyed by command line flags.
+
+For example, you might create a :command:`gmx solvate` operation as::
+
+ solvate = gmx.commandline_operation('gmx',
+ arguments=['solvate', '-box', '5', '5', '5'],
+ input_files={'-cs': structurefile},
+ output_files={'-p': topfile,
+ '-o': structurefile,
+ }
+
+To check the status or error output of a command line operation, refer to the
+``returncode`` and ``erroroutput`` outputs.
+To access the results from the output file arguments, use the command line flags
+as keys in the ``file`` dictionary output.
+
+Example::
+
+ structurefile = solvate.output.file['-o'].result()
+ if solvate.output.returncode.result() != 0:
+ print(solvate.output.erroroutput.result())
+
+Preparing simulations
+=====================
+
+Continuing the previous example, the output of ``solvate`` may be used as the
+input for ``grompp``::
+
+ grompp = gmx.commandline_operation('gmx', 'grompp',
+ input_files={
+ '-f': mdpfile,
+ '-p': solvate.output.file['-p'],
+ '-c': solvate.output.file['-o'],
+ '-po': mdout_mdp,
+ },
+ output_files={'-o': tprfile})
+
+Then, ``grompp.output.file['-o']`` can be used as the input for :py:func:`gmxapi.read_tpr`.
+
+Simulation input can be modified with the :py:func:`gmxapi.modify_input` operation
+before being passed to :py:func:`gmxapi.mdrun`.
+For *gmxapi* 0.1, a subset of MDP parameters may be overridden using the
+dictionary passed with the ``parameters`` key word argument.
+
+Example::
+
+ simulation_input = gmx.read_tpr(grompp.output.file['-o'])
+ modified_input = gmx.modify_input(input=simulation_input, parameters={'nsteps': 1000})
+ md = gmx.mdrun(input=modified_input)
+ md.run()
+
+Using arbitrary Python functions
+================================
+
+Generally, a function in the *gmxapi* package returns an object that references
+a node in a work graph,
+representing an operation that will be run when the graph executes.
+The object has an ``output`` attribute providing access to data Futures that
+can be provided as inputs to other operations before computation has actually
+been performed.
+
+You can also provide native Python data as input to operations,
+or you can operate on native results retrieved from a Future's ``result()``
+method.
+However, it is trivial to convert most Python functions into *gmxapi* compatible
+operations with :py:func:`gmxapi.function_wrapper`.
+All function inputs and outputs must have a name and type.
+Additionally, functions should be stateless and importable
+(e.g. via Python ``from some.module import myfunction``)
+for future compatibility.
+
+Simple functions can just use :py:func:`return` to publish their output,
+as long as they are defined with a return value type annotation.
+Functions with multiple outputs can accept an ``output`` key word argument and
+assign values to named attributes on the received argument.
+
+Examples::
+
+ from gmxapi import function_wrapper
+
+ @function_wrapper(output={'data': float})
+ def add_float(a: float, b: float) -> float:
+ return a + b
+
+
+ @function_wrapper(output={'data': bool})
+ def less_than(lhs: float, rhs: float, output=None):
+ output.data = lhs < rhs
+
+.. seealso::
+
+ For more on Python type hinting with function annotations,
+ check out :pep:`3107`.
+
+Subgraphs
+=========
+
+Basic *gmxapi* work consists of a flow of data from operation outputs to
+operation inputs, forming a directed acyclic graph (DAG).
+In many cases, it can be useful to repeat execution of a subgraph with updated
+inputs.
+You may want a data reference that is not tied to the immutable result
+of a single node in the work graph, but which instead refers to the most recent
+result of a repeated operation.
+
+One or more operations can be staged in a :py:class:`gmxapi.operation.Subgraph`,
+a sort of meta-operation factory that can store input binding behavior so that
+instances can be created without providing input arguments.
+
+The subgraph *variables* serve as input, output, and mutable internal data
+references which can be updated by operations in the subgraph.
+Variables also allow state to be propagated between iterations when a subgraph
+is used in a *while* loop.
+
+Use :py:func:`gmxapi.subgraph` to create a new empty subgraph.
+The ``variables`` argument declares data handles that define the state of the
+subgraph when it is run.
+To initialize input to the subgraph, give each variable a name and a value.
+
+To populate a subgraph, enter a SubgraphContext by using a :py:func:`with` statement.
+Operations created in the *with* block will be captued by the SubgraphContext.
+Define the subgraph outputs by assigning operation outputs to subgraph variables
+within the *with* block.
+
+After exiting the *with* block, the subgraph may be used to create operation
+instances or may be executed repeatedly in a *while* loop.
+
+.. note::
+
+ The object returned by :py:func:`gmxapi.subgraph` is atypical of *gmxapi*
+ operations, and has some special behaviors. When used as a Python
+ `context manager <https://docs.python.org/3/reference/datamodel.html#context-managers>`__,
+ it enters a "builder" state that changes the behavior of its attribute
+ variables and of operaton instantiation. After exiting the :py:func:`with`
+ block, the subgraph variables are no longer assignable, and operation
+ references obtained within the block are no longer valid.
+
+Looping
+=======
+
+An operation can be executed an arbitrary number of times with a
+:py:func:`gmxapi.while_loop` by providing a factory function as the
+*operation* argument.
+When the loop operation is run, the *operation* is instantiated and run repeatedly
+until *condition* evaluates ``True``.
+
+:py:func:`gmxapi.while_loop` does not provide a direct way to provide *operation*
+arguments. Use a *subgraph* to define the data flow for iterative operations.
+
+When a *condition* is a subgraph variable, the variable is evaluated in the
+running subgraph instance at the beginning of an iteration.
+
+Example::
+
+ subgraph = gmx.subgraph(variables={'float_with_default': 1.0, 'bool_data': True})
+ with subgraph:
+ # Define the update for float_with_default to come from an add_float operation.
+ subgraph.float_with_default = add_float(subgraph.float_with_default, 1.).output.data
+ subgraph.bool_data = less_than(lhs=subgraph.float_with_default, rhs=6.).output.data
+ operation_instance = subgraph()
+ operation_instance.run()
+ assert operation_instance.values['float_with_default'] == 2.
+
+ loop = gmx.while_loop(operation=subgraph, condition=subgraph.bool_data)
+ handle = loop()
+ assert handle.output.float_with_default.result() == 6
+
+.. _gmxapi logging:
+
+Logging
+=======
+
+*gmxapi* uses the Python :py:mod:`logging` module to provide hierarchical
+logging, organized by submodule.
+You can access the logger at ``gmxapi.logger`` or, after importing *gmxapi*,
+through the Python logging framework::
+
+ import gmxapi as gmx
+ import logging
+
+ # Get the root gmxapi logger.
+ gmx_logger = logging.getLogger('gmxapi')
+ # Set a low default logging level
+ gmx_logger.setLevel(logging.WARNING)
+ # Make some tools very verbose
+ # by descending the hierarchy
+ gmx_logger.getChild('commandline').setLevel(logging.DEBUG)
+ # or by direct reference
+ logging.getLogger('gmxapi.mdrun').setLevel(logging.DEBUG)
+
+You may prefer to adjust the log format or manipulate the log handlers.
+For example, tag the log output with MPI rank::
+
+ try:
+ from mpi4py import MPI
+ rank_number = MPI.COMM_WORLD.Get_rank()
+ except ImportError:
+ rank_number = 0
+ rank_tag = ''
+ MPI = None
+ else:
+ rank_tag = 'rank{}:'.format(rank_number)
+
+ formatter = logging.Formatter(rank_tag + '%(name)s:%(levelname)s: %(message)s')
+
+ # For additional console logging, create and attach a stream handler.
+ ch = logging.StreamHandler()
+ ch.setFormatter(formatter)
+ logging.getLogger().addHandler(ch)
+
+For more information, refer to the Python `logging documentation <https://docs.python.org/3/library/logging.html>`__.
+
+More
+====
+
+Refer to the :doc:`pythonreference` for complete and granular documentation.
+
+For more information on writing or using pluggable simulation extension code,
+refer to https://redmine.gromacs.org/issues/3133.
+(For gmxapi 0.0.7 and GROMACS 2019, see https://github.com/kassonlab/sample_restraint)
+
+.. todo:: :issue:`3133`: Replace these links as resources for pluggable extension code become available.
.. TODO needs link to ref list
-When citing this document in any scientific publication
-please refer to it as:
-
-.. parsed-literal::
-
- M.J. Abraham, D. van der Spoel, E. Lindahl, B. Hess, and the GROMACS development team,
- |Gromacs| User Manual version |version|,
- `http://www.gromacs.org <http://www.gromacs.org>`__
+|GMX_MANUAL_DOI_STRING|
However, we prefer that you cite (some of) the |Gromacs|
papers:
import os
import shutil
import tempfile
+import warnings
+from contextlib import contextmanager
import pytest
-@pytest.fixture()
-def cleandir():
+def pytest_addoption(parser):
+ """Add a command-line user option for the pytest invocation."""
+ parser.addoption(
+ '--rm',
+ action='store',
+ default='always',
+ choices=['always', 'never', 'success'],
+ help='Remove temporary directories "always", "never", or on "success".'
+ )
+
+
+@pytest.fixture(scope='session')
+def remove_tempdir(request):
+ """pytest fixture to get access to the --rm CLI option."""
+ return request.config.getoption('--rm')
+
+
+@contextmanager
+def scoped_chdir(dir):
+ oldpath = os.getcwd()
+ os.chdir(dir)
+ try:
+ yield dir
+ # If the `with` block using scoped_chdir produces an exception, it will
+ # be raised at this point in this function. We want the exception to
+ # propagate out of the `with` block, but first we want to restore the
+ # original working directory, so we skip `except` but provide a `finally`.
+ finally:
+ os.chdir(oldpath)
+
+
+@contextmanager
+def _cleandir(remove_tempdir):
+ """Context manager for a clean temporary working directory.
+
+ Arguments:
+ remove_tempdir (str): whether to remove temporary directory "always",
+ "never", or on "success"
+
+ The context manager will issue a warning for each temporary directory that
+ is not removed.
+ """
+
+ newpath = tempfile.mkdtemp()
+
+ def remove():
+ shutil.rmtree(newpath)
+
+ def warn():
+ warnings.warn('Temporary directory not removed: {}'.format(newpath))
+
+ if remove_tempdir == 'always':
+ callback = remove
+ else:
+ callback = warn
+ try:
+ with scoped_chdir(newpath):
+ yield newpath
+ # If we get to this line, the `with` block using _cleandir did not throw.
+ # Clean up the temporary directory unless the user specified `--rm never`.
+ # I.e. If the user specified `--rm success`, then we need to toggle from `warn` to `remove`.
+ if remove_tempdir != 'never':
+ callback = remove
+ finally:
+ callback()
+
+
+@pytest.fixture
+def cleandir(remove_tempdir):
"""Provide a clean temporary working directory for a test.
Example usage:
Ref: https://docs.pytest.org/en/latest/fixture.html#using-fixtures-from-classes-modules-or-projects
"""
- oldpath = os.getcwd()
- newpath = tempfile.mkdtemp()
- os.chdir(newpath)
- yield newpath
- os.chdir(oldpath)
- # TODO: Allow override to prevent removal of the temporary directory
- shutil.rmtree(newpath)
+ with _cleandir(remove_tempdir) as newdir:
+ yield newdir
@pytest.fixture(scope='session')
# TODO: (#2896) Find a more canonical way to identify the GROMACS commandline wrapper binary.
# We should be able to get the GMXRC contents and related hints from a gmxapi
# package resource or from module attributes of a ``gromacs`` stub package.
- command = shutil.which('gmx')
- if command is None:
- gmxbindir = os.getenv('GMXBIN')
- if gmxbindir is None:
- gromacsdir = os.getenv('GROMACS_DIR')
- if gromacsdir is not None and gromacsdir != '':
- gmxbindir = os.path.join(gromacsdir, 'bin')
- if gmxbindir is None:
- gmxapidir = os.getenv('gmxapi_DIR')
- if gmxapidir is not None and gmxapidir != '':
- gmxbindir = os.path.join(gmxapidir, 'bin')
- if gmxbindir is not None:
- gmxbindir = os.path.abspath(gmxbindir)
- command = shutil.which('gmx', path=gmxbindir)
- else:
- message = "Tests need 'gmx' command line tool, but could not find it on the path."
- raise RuntimeError(message)
+ allowed_command_names = ['gmx', 'gmx_mpi']
+ command = None
+ for command_name in allowed_command_names:
+ if command is not None:
+ break
+ command = shutil.which(command_name)
+ if command is None:
+ gmxbindir = os.getenv('GMXBIN')
+ if gmxbindir is None:
+ gromacsdir = os.getenv('GROMACS_DIR')
+ if gromacsdir is not None and gromacsdir != '':
+ gmxbindir = os.path.join(gromacsdir, 'bin')
+ if gmxbindir is None:
+ gmxapidir = os.getenv('gmxapi_DIR')
+ if gmxapidir is not None and gmxapidir != '':
+ gmxbindir = os.path.join(gmxapidir, 'bin')
+ if gmxbindir is not None:
+ gmxbindir = os.path.abspath(gmxbindir)
+ command = shutil.which(command_name, path=gmxbindir)
if command is None:
message = "Tests need 'gmx' command line tool, but could not find it on the path."
raise RuntimeError(message)
@pytest.fixture(scope='class')
-def spc_water_box(gmxcli):
+def spc_water_box(gmxcli, remove_tempdir):
"""Provide a TPR input file for a simple simulation.
Prepare the MD input in a freshly created working directory.
# # Ref https://setuptools.readthedocs.io/en/latest/setuptools.html#including-data-files
# from gmx.data import tprfilename
- tempdir = tempfile.mkdtemp()
-
- testdir = os.path.dirname(__file__)
- with open(os.path.join(testdir, 'testdata.json'), 'r') as fh:
- testdata = json.load(fh)
-
- # TODO: (#2756) Don't rely on so many automagical behaviors (as described in comments below)
-
- structurefile = os.path.join(tempdir, 'structure.gro')
- # We let `gmx solvate` use the default solvent. Otherwise, we would do
- # gro_input = testdata['solvent_structure']
- # with open(structurefile, 'w') as fh:
- # fh.write('\n'.join(gro_input))
- # fh.write('\n')
-
- topfile = os.path.join(tempdir, 'topology.top')
- top_input = testdata['solvent_topology']
- # `gmx solvate` will append a line to the provided file with the molecule count,
- # so we strip the last line from the input topology.
- with open(topfile, 'w') as fh:
- fh.write('\n'.join(top_input[:-1]))
- fh.write('\n')
-
- solvate = gmx.commandline_operation(gmxcli,
- arguments=['solvate', '-box', '5', '5', '5'],
- # We use the default solvent instead of specifying one.
- # input_files={'-cs': structurefile},
- output_files={'-p': topfile,
- '-o': structurefile,
- }
- )
- if solvate.output.returncode.result() != 0:
- logging.debug(solvate.output.erroroutput.result())
- raise RuntimeError('solvate failed in spc_water_box testing fixture.')
-
- mdp_input = [('integrator', 'md'),
- ('cutoff-scheme', 'Verlet'),
- ('nsteps', 2),
- ('nstxout', 1),
- ('nstvout', 1),
- ('nstfout', 1),
- ('tcoupl', 'v-rescale'),
- ('tc-grps', 'System'),
- ('tau-t', 1),
- ('ref-t', 298)]
- mdp_input = '\n'.join([' = '.join([str(item) for item in kvpair]) for kvpair in mdp_input])
- mdpfile = os.path.join(tempdir, 'md.mdp')
- with open(mdpfile, 'w') as fh:
- fh.write(mdp_input)
- fh.write('\n')
- tprfile = os.path.join(tempdir, 'topol.tpr')
- # We don't use mdout_mdp, but if we don't specify it to grompp,
- # it will be created in the current working directory.
- mdout_mdp = os.path.join(tempdir, 'mdout.mdp')
-
- grompp = gmx.commandline_operation(gmxcli, 'grompp',
- input_files={
- '-f': mdpfile,
- '-p': solvate.output.file['-p'],
- '-c': solvate.output.file['-o'],
- '-po': mdout_mdp,
- },
- output_files={'-o': tprfile})
- tprfilename = grompp.output.file['-o'].result()
- if grompp.output.returncode.result() != 0:
- logging.debug(grompp.output.erroroutput.result())
- raise RuntimeError('grompp failed in spc_water_box testing fixture.')
-
- # TODO: more inspection of grompp errors...
- assert os.path.exists(tprfilename)
- yield tprfilename
-
- # Clean up.
- # Note: these lines are not executed (and tempdir is not removed) if there
- # are exceptions before `yield`
- # TODO: Allow override to prevent removal of the temporary directory
- shutil.rmtree(tempdir)
+ with _cleandir(remove_tempdir) as tempdir:
+
+ testdir = os.path.dirname(__file__)
+ with open(os.path.join(testdir, 'testdata.json'), 'r') as fh:
+ testdata = json.load(fh)
+
+ # TODO: (#2756) Don't rely on so many automagical behaviors (as described in comments below)
+
+ structurefile = os.path.join(tempdir, 'structure.gro')
+ # We let `gmx solvate` use the default solvent. Otherwise, we would do
+ # gro_input = testdata['solvent_structure']
+ # with open(structurefile, 'w') as fh:
+ # fh.write('\n'.join(gro_input))
+ # fh.write('\n')
+
+ topfile = os.path.join(tempdir, 'topology.top')
+ top_input = testdata['solvent_topology']
+ # `gmx solvate` will append a line to the provided file with the molecule count,
+ # so we strip the last line from the input topology.
+ with open(topfile, 'w') as fh:
+ fh.write('\n'.join(top_input[:-1]))
+ fh.write('\n')
+
+ assert os.path.exists(topfile)
+ solvate = gmx.commandline_operation(gmxcli,
+ arguments=['solvate', '-box', '5', '5', '5'],
+ # We use the default solvent instead of specifying one.
+ # input_files={'-cs': structurefile},
+ output_files={'-p': topfile,
+ '-o': structurefile,
+ }
+ )
+ assert os.path.exists(topfile)
+
+ if solvate.output.returncode.result() != 0:
+ logging.debug(solvate.output.erroroutput.result())
+ raise RuntimeError('solvate failed in spc_water_box testing fixture.')
+
+ # Choose an exactly representable dt of 2^-9 ps (approximately 0.002)
+ dt = 2.**-9.
+ mdp_input = [('integrator', 'md'),
+ ('dt', dt),
+ ('cutoff-scheme', 'Verlet'),
+ ('nsteps', 2),
+ ('nstxout', 1),
+ ('nstvout', 1),
+ ('nstfout', 1),
+ ('tcoupl', 'v-rescale'),
+ ('tc-grps', 'System'),
+ ('tau-t', 1),
+ ('ref-t', 298)]
+ mdp_input = '\n'.join([' = '.join([str(item) for item in kvpair]) for kvpair in mdp_input])
+ mdpfile = os.path.join(tempdir, 'md.mdp')
+ with open(mdpfile, 'w') as fh:
+ fh.write(mdp_input)
+ fh.write('\n')
+ tprfile = os.path.join(tempdir, 'topol.tpr')
+ # We don't use mdout_mdp, but if we don't specify it to grompp,
+ # it will be created in the current working directory.
+ mdout_mdp = os.path.join(tempdir, 'mdout.mdp')
+
+ grompp = gmx.commandline_operation(gmxcli, 'grompp',
+ input_files={
+ '-f': mdpfile,
+ '-p': solvate.output.file['-p'],
+ '-c': solvate.output.file['-o'],
+ '-po': mdout_mdp,
+ },
+ output_files={'-o': tprfile})
+ tprfilename = grompp.output.file['-o'].result()
+ if grompp.output.returncode.result() != 0:
+ logging.debug(grompp.output.erroroutput.result())
+ raise RuntimeError('grompp failed in spc_water_box testing fixture.')
+
+ # TODO: more inspection of grompp errors...
+ assert os.path.exists(tprfilename)
+ yield tprfilename
import os
import shutil
import tempfile
+import warnings
+from contextlib import contextmanager
import pytest
-@pytest.fixture()
-def cleandir():
+def pytest_addoption(parser):
+ """Add a command-line user option for the pytest invocation."""
+ parser.addoption(
+ '--rm',
+ action='store',
+ default='always',
+ choices=['always', 'never', 'success'],
+ help='Remove temporary directories "always", "never", or on "success".'
+ )
+
+
+@pytest.fixture(scope='session')
+def remove_tempdir(request):
+ """pytest fixture to get access to the --rm CLI option."""
+ return request.config.getoption('--rm')
+
+
+@contextmanager
+def scoped_chdir(dir):
+ oldpath = os.getcwd()
+ os.chdir(dir)
+ try:
+ yield dir
+ # If the `with` block using scoped_chdir produces an exception, it will
+ # be raised at this point in this function. We want the exception to
+ # propagate out of the `with` block, but first we want to restore the
+ # original working directory, so we skip `except` but provide a `finally`.
+ finally:
+ os.chdir(oldpath)
+
+
+@contextmanager
+def _cleandir(remove_tempdir):
+ """Context manager for a clean temporary working directory.
+
+ Arguments:
+ remove_tempdir (str): whether to remove temporary directory "always",
+ "never", or on "success"
+
+ The context manager will issue a warning for each temporary directory that
+ is not removed.
+ """
+
+ newpath = tempfile.mkdtemp()
+
+ def remove():
+ shutil.rmtree(newpath)
+
+ def warn():
+ warnings.warn('Temporary directory not removed: {}'.format(newpath))
+
+ if remove_tempdir == 'always':
+ callback = remove
+ else:
+ callback = warn
+ try:
+ with scoped_chdir(newpath):
+ yield newpath
+ # If we get to this line, the `with` block using _cleandir did not throw.
+ # Clean up the temporary directory unless the user specified `--rm never`.
+ # I.e. If the user specified `--rm success`, then we need to toggle from `warn` to `remove`.
+ if remove_tempdir != 'never':
+ callback = remove
+ finally:
+ callback()
+
+
+@pytest.fixture
+def cleandir(remove_tempdir):
"""Provide a clean temporary working directory for a test.
Example usage:
Ref: https://docs.pytest.org/en/latest/fixture.html#using-fixtures-from-classes-modules-or-projects
"""
- oldpath = os.getcwd()
- newpath = tempfile.mkdtemp()
- os.chdir(newpath)
- yield newpath
- os.chdir(oldpath)
- # TODO: Allow override to prevent removal of the temporary directory
- shutil.rmtree(newpath)
+ with _cleandir(remove_tempdir) as newdir:
+ yield newdir
@pytest.fixture(scope='session')
@pytest.fixture(scope='class')
-def spc_water_box(gmxcli):
+def spc_water_box(gmxcli, remove_tempdir):
"""Provide a TPR input file for a simple simulation.
Prepare the MD input in a freshly created working directory.
# # Ref https://setuptools.readthedocs.io/en/latest/setuptools.html#including-data-files
# from gmx.data import tprfilename
- tempdir = tempfile.mkdtemp()
-
- testdir = os.path.dirname(__file__)
- with open(os.path.join(testdir, 'testdata.json'), 'r') as fh:
- testdata = json.load(fh)
-
- # TODO: (#2756) Don't rely on so many automagical behaviors (as described in comments below)
-
- structurefile = os.path.join(tempdir, 'structure.gro')
- # We let `gmx solvate` use the default solvent. Otherwise, we would do
- # gro_input = testdata['solvent_structure']
- # with open(structurefile, 'w') as fh:
- # fh.write('\n'.join(gro_input))
- # fh.write('\n')
-
- topfile = os.path.join(tempdir, 'topology.top')
- top_input = testdata['solvent_topology']
- # `gmx solvate` will append a line to the provided file with the molecule count,
- # so we strip the last line from the input topology.
- with open(topfile, 'w') as fh:
- fh.write('\n'.join(top_input[:-1]))
- fh.write('\n')
-
- solvate = gmx.commandline_operation(gmxcli,
- arguments=['solvate', '-box', '5', '5', '5'],
- # We use the default solvent instead of specifying one.
- # input_files={'-cs': structurefile},
- output_files={'-p': topfile,
- '-o': structurefile,
- }
- )
- if solvate.output.returncode.result() != 0:
- logging.debug(solvate.output.erroroutput.result())
- raise RuntimeError('solvate failed in spc_water_box testing fixture.')
-
- # Choose an exactly representable dt of 2^-9 ps (approximately 0.002)
- dt = 2.**-9.
- mdp_input = [('integrator', 'md'),
- ('dt', dt),
- ('cutoff-scheme', 'Verlet'),
- ('nsteps', 2),
- ('nstxout', 1),
- ('nstvout', 1),
- ('nstfout', 1),
- ('tcoupl', 'v-rescale'),
- ('tc-grps', 'System'),
- ('tau-t', 1),
- ('ref-t', 298)]
- mdp_input = '\n'.join([' = '.join([str(item) for item in kvpair]) for kvpair in mdp_input])
- mdpfile = os.path.join(tempdir, 'md.mdp')
- with open(mdpfile, 'w') as fh:
- fh.write(mdp_input)
- fh.write('\n')
- tprfile = os.path.join(tempdir, 'topol.tpr')
- # We don't use mdout_mdp, but if we don't specify it to grompp,
- # it will be created in the current working directory.
- mdout_mdp = os.path.join(tempdir, 'mdout.mdp')
-
- grompp = gmx.commandline_operation(gmxcli, 'grompp',
- input_files={
- '-f': mdpfile,
- '-p': solvate.output.file['-p'],
- '-c': solvate.output.file['-o'],
- '-po': mdout_mdp,
- },
- output_files={'-o': tprfile})
- tprfilename = grompp.output.file['-o'].result()
- if grompp.output.returncode.result() != 0:
- logging.debug(grompp.output.erroroutput.result())
- raise RuntimeError('grompp failed in spc_water_box testing fixture.')
-
- # TODO: more inspection of grompp errors...
- assert os.path.exists(tprfilename)
- yield tprfilename
-
- # Clean up.
- # Note: these lines are not executed (and tempdir is not removed) if there
- # are exceptions before `yield`
- # TODO: Allow override to prevent removal of the temporary directory
- shutil.rmtree(tempdir)
+ with _cleandir(remove_tempdir) as tempdir:
+
+ testdir = os.path.dirname(__file__)
+ with open(os.path.join(testdir, 'testdata.json'), 'r') as fh:
+ testdata = json.load(fh)
+
+ # TODO: (#2756) Don't rely on so many automagical behaviors (as described in comments below)
+
+ structurefile = os.path.join(tempdir, 'structure.gro')
+ # We let `gmx solvate` use the default solvent. Otherwise, we would do
+ # gro_input = testdata['solvent_structure']
+ # with open(structurefile, 'w') as fh:
+ # fh.write('\n'.join(gro_input))
+ # fh.write('\n')
+
+ topfile = os.path.join(tempdir, 'topology.top')
+ top_input = testdata['solvent_topology']
+ # `gmx solvate` will append a line to the provided file with the molecule count,
+ # so we strip the last line from the input topology.
+ with open(topfile, 'w') as fh:
+ fh.write('\n'.join(top_input[:-1]))
+ fh.write('\n')
+
+ assert os.path.exists(topfile)
+ solvate = gmx.commandline_operation(gmxcli,
+ arguments=['solvate', '-box', '5', '5', '5'],
+ # We use the default solvent instead of specifying one.
+ # input_files={'-cs': structurefile},
+ output_files={'-p': topfile,
+ '-o': structurefile,
+ }
+ )
+ assert os.path.exists(topfile)
+
+ if solvate.output.returncode.result() != 0:
+ logging.debug(solvate.output.erroroutput.result())
+ raise RuntimeError('solvate failed in spc_water_box testing fixture.')
+
+ # Choose an exactly representable dt of 2^-9 ps (approximately 0.002)
+ dt = 2.**-9.
+ mdp_input = [('integrator', 'md'),
+ ('dt', dt),
+ ('cutoff-scheme', 'Verlet'),
+ ('nsteps', 2),
+ ('nstxout', 1),
+ ('nstvout', 1),
+ ('nstfout', 1),
+ ('tcoupl', 'v-rescale'),
+ ('tc-grps', 'System'),
+ ('tau-t', 1),
+ ('ref-t', 298)]
+ mdp_input = '\n'.join([' = '.join([str(item) for item in kvpair]) for kvpair in mdp_input])
+ mdpfile = os.path.join(tempdir, 'md.mdp')
+ with open(mdpfile, 'w') as fh:
+ fh.write(mdp_input)
+ fh.write('\n')
+ tprfile = os.path.join(tempdir, 'topol.tpr')
+ # We don't use mdout_mdp, but if we don't specify it to grompp,
+ # it will be created in the current working directory.
+ mdout_mdp = os.path.join(tempdir, 'mdout.mdp')
+
+ grompp = gmx.commandline_operation(gmxcli, 'grompp',
+ input_files={
+ '-f': mdpfile,
+ '-p': solvate.output.file['-p'],
+ '-c': solvate.output.file['-o'],
+ '-po': mdout_mdp,
+ },
+ output_files={'-o': tprfile})
+ tprfilename = grompp.output.file['-o'].result()
+ if grompp.output.returncode.result() != 0:
+ logging.debug(grompp.output.erroroutput.result())
+ raise RuntimeError('grompp failed in spc_water_box testing fixture.')
+
+ # TODO: more inspection of grompp errors...
+ assert os.path.exists(tprfilename)
+ yield tprfilename
@gmx.function_wrapper(output={'data': bool})
-def less_than(lhs: float, rhs: float) -> bool:
- return lhs < rhs
+def less_than(lhs: float, rhs: float, output=None):
+ output.data = lhs < rhs
def test_subgraph_function():
import os
import shutil
import tempfile
+import warnings
+from contextlib import contextmanager
import pytest
+def pytest_addoption(parser):
+ """Add a command-line user option for the pytest invocation."""
+ parser.addoption(
+ '--rm',
+ action='store',
+ default='always',
+ choices=['always', 'never', 'success'],
+ help='Remove temporary directories "always", "never", or on "success".'
+ )
-@pytest.fixture()
-def cleandir():
+
+@pytest.fixture(scope='session')
+def remove_tempdir(request):
+ """pytest fixture to get access to the --rm CLI option."""
+ return request.config.getoption('--rm')
+
+
+@contextmanager
+def scoped_chdir(dir):
+ oldpath = os.getcwd()
+ os.chdir(dir)
+ try:
+ yield dir
+ # If the `with` block using scoped_chdir produces an exception, it will
+ # be raised at this point in this function. We want the exception to
+ # propagate out of the `with` block, but first we want to restore the
+ # original working directory, so we skip `except` but provide a `finally`.
+ finally:
+ os.chdir(oldpath)
+
+
+@contextmanager
+def _cleandir(remove_tempdir):
+ """Context manager for a clean temporary working directory.
+
+ Arguments:
+ remove_tempdir (str): whether to remove temporary directory "always",
+ "never", or on "success"
+
+ The context manager will issue a warning for each temporary directory that
+ is not removed.
+ """
+
+ newpath = tempfile.mkdtemp()
+
+ def remove():
+ shutil.rmtree(newpath)
+
+ def warn():
+ warnings.warn('Temporary directory not removed: {}'.format(newpath))
+
+ if remove_tempdir == 'always':
+ callback = remove
+ else:
+ callback = warn
+ try:
+ with scoped_chdir(newpath):
+ yield newpath
+ # If we get to this line, the `with` block using _cleandir did not throw.
+ # Clean up the temporary directory unless the user specified `--rm never`.
+ # I.e. If the user specified `--rm success`, then we need to toggle from `warn` to `remove`.
+ if remove_tempdir != 'never':
+ callback = remove
+ finally:
+ callback()
+
+
+@pytest.fixture
+def cleandir(remove_tempdir):
"""Provide a clean temporary working directory for a test.
Example usage:
Ref: https://docs.pytest.org/en/latest/fixture.html#using-fixtures-from-classes-modules-or-projects
"""
- oldpath = os.getcwd()
- newpath = tempfile.mkdtemp()
- os.chdir(newpath)
- yield newpath
- os.chdir(oldpath)
- # TODO: Allow override to prevent removal of the temporary directory
- shutil.rmtree(newpath)
+ with _cleandir(remove_tempdir) as newdir:
+ yield newdir
@pytest.fixture(scope='session')
# TODO: (#2896) Find a more canonical way to identify the GROMACS commandline wrapper binary.
# We should be able to get the GMXRC contents and related hints from a gmxapi
# package resource or from module attributes of a ``gromacs`` stub package.
- command = shutil.which('gmx')
- if command is None:
- gmxbindir = os.getenv('GMXBIN')
- if gmxbindir is None:
- gromacsdir = os.getenv('GROMACS_DIR')
- if gromacsdir is not None and gromacsdir != '':
- gmxbindir = os.path.join(gromacsdir, 'bin')
- if gmxbindir is None:
- gmxapidir = os.getenv('gmxapi_DIR')
- if gmxapidir is not None and gmxapidir != '':
- gmxbindir = os.path.join(gmxapidir, 'bin')
- if gmxbindir is not None:
- gmxbindir = os.path.abspath(gmxbindir)
- command = shutil.which('gmx', path=gmxbindir)
- else:
- message = "Tests need 'gmx' command line tool, but could not find it on the path."
- raise RuntimeError(message)
+ allowed_command_names = ['gmx', 'gmx_mpi']
+ command = None
+ for command_name in allowed_command_names:
+ if command is not None:
+ break
+ command = shutil.which(command_name)
+ if command is None:
+ gmxbindir = os.getenv('GMXBIN')
+ if gmxbindir is None:
+ gromacsdir = os.getenv('GROMACS_DIR')
+ if gromacsdir is not None and gromacsdir != '':
+ gmxbindir = os.path.join(gromacsdir, 'bin')
+ if gmxbindir is None:
+ gmxapidir = os.getenv('gmxapi_DIR')
+ if gmxapidir is not None and gmxapidir != '':
+ gmxbindir = os.path.join(gmxapidir, 'bin')
+ if gmxbindir is not None:
+ gmxbindir = os.path.abspath(gmxbindir)
+ command = shutil.which(command_name, path=gmxbindir)
if command is None:
message = "Tests need 'gmx' command line tool, but could not find it on the path."
raise RuntimeError(message)
@pytest.fixture(scope='class')
-def spc_water_box(gmxcli):
+def spc_water_box(gmxcli, remove_tempdir):
"""Provide a TPR input file for a simple simulation.
Prepare the MD input in a freshly created working directory.
# # Ref https://setuptools.readthedocs.io/en/latest/setuptools.html#including-data-files
# from gmx.data import tprfilename
- tempdir = tempfile.mkdtemp()
-
- testdir = os.path.dirname(__file__)
- with open(os.path.join(testdir, 'testdata.json'), 'r') as fh:
- testdata = json.load(fh)
-
- # TODO: (#2756) Don't rely on so many automagical behaviors (as described in comments below)
-
- structurefile = os.path.join(tempdir, 'structure.gro')
- # We let `gmx solvate` use the default solvent. Otherwise, we would do
- # gro_input = testdata['solvent_structure']
- # with open(structurefile, 'w') as fh:
- # fh.write('\n'.join(gro_input))
- # fh.write('\n')
-
- topfile = os.path.join(tempdir, 'topology.top')
- top_input = testdata['solvent_topology']
- # `gmx solvate` will append a line to the provided file with the molecule count,
- # so we strip the last line from the input topology.
- with open(topfile, 'w') as fh:
- fh.write('\n'.join(top_input[:-1]))
- fh.write('\n')
-
- solvate = gmx.commandline_operation(gmxcli,
- arguments=['solvate', '-box', '5', '5', '5'],
- # We use the default solvent instead of specifying one.
- # input_files={'-cs': structurefile},
- output_files={'-p': topfile,
- '-o': structurefile,
- }
- )
- if solvate.output.returncode.result() != 0:
- logging.debug(solvate.output.erroroutput.result())
- raise RuntimeError('solvate failed in spc_water_box testing fixture.')
-
- mdp_input = [('integrator', 'md'),
- ('cutoff-scheme', 'Verlet'),
- ('nsteps', 2),
- ('nstxout', 1),
- ('nstvout', 1),
- ('nstfout', 1),
- ('tcoupl', 'v-rescale'),
- ('tc-grps', 'System'),
- ('tau-t', 1),
- ('ref-t', 298)]
- mdp_input = '\n'.join([' = '.join([str(item) for item in kvpair]) for kvpair in mdp_input])
- mdpfile = os.path.join(tempdir, 'md.mdp')
- with open(mdpfile, 'w') as fh:
- fh.write(mdp_input)
- fh.write('\n')
- tprfile = os.path.join(tempdir, 'topol.tpr')
- # We don't use mdout_mdp, but if we don't specify it to grompp,
- # it will be created in the current working directory.
- mdout_mdp = os.path.join(tempdir, 'mdout.mdp')
-
- grompp = gmx.commandline_operation(gmxcli, 'grompp',
- input_files={
- '-f': mdpfile,
- '-p': solvate.output.file['-p'],
- '-c': solvate.output.file['-o'],
- '-po': mdout_mdp,
- },
- output_files={'-o': tprfile})
- tprfilename = grompp.output.file['-o'].result()
- if grompp.output.returncode.result() != 0:
- logging.debug(grompp.output.erroroutput.result())
- raise RuntimeError('grompp failed in spc_water_box testing fixture.')
-
- # TODO: more inspection of grompp errors...
- assert os.path.exists(tprfilename)
- yield tprfilename
-
- # Clean up.
- # Note: these lines are not executed (and tempdir is not removed) if there
- # are exceptions before `yield`
- # TODO: Allow override to prevent removal of the temporary directory
- shutil.rmtree(tempdir)
+ with _cleandir(remove_tempdir) as tempdir:
+
+ testdir = os.path.dirname(__file__)
+ with open(os.path.join(testdir, 'testdata.json'), 'r') as fh:
+ testdata = json.load(fh)
+
+ # TODO: (#2756) Don't rely on so many automagical behaviors (as described in comments below)
+
+ structurefile = os.path.join(tempdir, 'structure.gro')
+ # We let `gmx solvate` use the default solvent. Otherwise, we would do
+ # gro_input = testdata['solvent_structure']
+ # with open(structurefile, 'w') as fh:
+ # fh.write('\n'.join(gro_input))
+ # fh.write('\n')
+
+ topfile = os.path.join(tempdir, 'topology.top')
+ top_input = testdata['solvent_topology']
+ # `gmx solvate` will append a line to the provided file with the molecule count,
+ # so we strip the last line from the input topology.
+ with open(topfile, 'w') as fh:
+ fh.write('\n'.join(top_input[:-1]))
+ fh.write('\n')
+
+ assert os.path.exists(topfile)
+ solvate = gmx.commandline_operation(gmxcli,
+ arguments=['solvate', '-box', '5', '5', '5'],
+ # We use the default solvent instead of specifying one.
+ # input_files={'-cs': structurefile},
+ output_files={'-p': topfile,
+ '-o': structurefile,
+ }
+ )
+ assert os.path.exists(topfile)
+
+ if solvate.output.returncode.result() != 0:
+ logging.debug(solvate.output.erroroutput.result())
+ raise RuntimeError('solvate failed in spc_water_box testing fixture.')
+
+ # Choose an exactly representable dt of 2^-9 ps (approximately 0.002)
+ dt = 2.**-9.
+ mdp_input = [('integrator', 'md'),
+ ('dt', dt),
+ ('cutoff-scheme', 'Verlet'),
+ ('nsteps', 2),
+ ('nstxout', 1),
+ ('nstvout', 1),
+ ('nstfout', 1),
+ ('tcoupl', 'v-rescale'),
+ ('tc-grps', 'System'),
+ ('tau-t', 1),
+ ('ref-t', 298)]
+ mdp_input = '\n'.join([' = '.join([str(item) for item in kvpair]) for kvpair in mdp_input])
+ mdpfile = os.path.join(tempdir, 'md.mdp')
+ with open(mdpfile, 'w') as fh:
+ fh.write(mdp_input)
+ fh.write('\n')
+ tprfile = os.path.join(tempdir, 'topol.tpr')
+ # We don't use mdout_mdp, but if we don't specify it to grompp,
+ # it will be created in the current working directory.
+ mdout_mdp = os.path.join(tempdir, 'mdout.mdp')
+
+ grompp = gmx.commandline_operation(gmxcli, 'grompp',
+ input_files={
+ '-f': mdpfile,
+ '-p': solvate.output.file['-p'],
+ '-c': solvate.output.file['-o'],
+ '-po': mdout_mdp,
+ },
+ output_files={'-o': tprfile})
+ tprfilename = grompp.output.file['-o'].result()
+ if grompp.output.returncode.result() != 0:
+ logging.debug(grompp.output.erroroutput.result())
+ raise RuntimeError('grompp failed in spc_water_box testing fixture.')
+
+ # TODO: more inspection of grompp errors...
+ assert os.path.exists(tprfilename)
+ yield tprfilename
they will not generate warnings from such headers included in GROMACS
test code.
+As gtest is not building libraries inside googletest/googletest/src,
+the link_directories command for that is removed.
+
The googlemock/CMakeLists.txt no longer builds
googletest/src/gtest-all.cc as part of the gmock target, because it is
already built as part of the gtest target, and this leads to duplicate
${gtest_SOURCE_DIR}/include
${gtest_SOURCE_DIR})
-# Where Google Test's libraries can be found.
-link_directories(${gtest_BINARY_DIR}/src)
-
# Summary of tuple support for Microsoft Visual Studio:
# Compiler version(MS) version(cmake) Support
# ---------- ----------- -------------- -----------------------------
#include "pme_internal.h"
#include "pme_solve.h"
-
-/*! \brief
- * CUDA only.
- * Controls if we should use order (i.e. 4) threads per atom for the GPU
- * or order*order (i.e. 16) threads per atom.
- */
-constexpr bool c_useOrderThreadsPerAtom = false;
/*! \brief
- * CUDA only.
- * Controls if we should recalculate the splines in the gather or
- * save the values in the spread and reload in the gather.
+ * CUDA only
+ * Atom limit above which it is advantageous to turn on the
+ * recalcuating of the splines in the gather and using less threads per atom in the spline and spread
*/
-constexpr bool c_recalculateSplines = false;
+constexpr int c_pmeGpuPerformanceAtomLimit = 23000;
/*! \internal \brief
* Wrapper for getting a pointer to the plain C++ part of the GPU kernel parameters structure.
int pme_gpu_get_atoms_per_warp(const PmeGpu* pmeGpu)
{
- if (c_useOrderThreadsPerAtom)
+ if (pmeGpu->settings.useOrderThreadsPerAtom)
{
return pmeGpu->programHandle_->impl_->warpSize / c_pmeSpreadGatherThreadsPerAtom4ThPerAtom;
}
pmeGpu->common->boxScaler = pme->boxScaler;
}
+/*! \libinternal \brief
+ * uses heuristics to select the best performing PME gather and scatter kernels
+ *
+ * \param[in,out] pmeGpu The PME GPU structure.
+ */
+static void pme_gpu_select_best_performing_pme_spreadgather_kernels(PmeGpu* pmeGpu)
+{
+ if (pmeGpu->kernelParams->atoms.nAtoms > c_pmeGpuPerformanceAtomLimit && (GMX_GPU == GMX_GPU_CUDA))
+ {
+ pmeGpu->settings.useOrderThreadsPerAtom = true;
+ pmeGpu->settings.recalculateSplines = true;
+ }
+ else
+ {
+ pmeGpu->settings.useOrderThreadsPerAtom = false;
+ pmeGpu->settings.recalculateSplines = false;
+ }
+}
+
+
/*! \libinternal \brief
* Initializes the PME GPU data at the beginning of the run.
* TODO: this should become PmeGpu::PmeGpu()
* update for mixed mode on grid switch. TODO: use shared recipbox field.
*/
std::memset(pme->gpu->common->previousBox, 0, sizeof(pme->gpu->common->previousBox));
+ pme_gpu_select_best_performing_pme_spreadgather_kernels(pme->gpu);
}
void pme_gpu_destroy(PmeGpu* pmeGpu)
const int order = pmeGpu->common->pme_order;
GMX_ASSERT(order == c_pmeGpuOrder, "Only PME order 4 is implemented");
- const bool writeGlobal = pmeGpu->settings.copyAllOutputs;
+ const bool writeGlobal = pmeGpu->settings.copyAllOutputs;
+ const bool useOrderThreadsPerAtom = pmeGpu->settings.useOrderThreadsPerAtom;
+ const bool recalculateSplines = pmeGpu->settings.recalculateSplines;
#if GMX_GPU == GMX_GPU_OPENCL
- GMX_ASSERT(!c_useOrderThreadsPerAtom, "Only 16 threads per atom supported in OpenCL");
- GMX_ASSERT(!c_recalculateSplines, "Recalculating splines not supported in OpenCL");
+ GMX_ASSERT(!useOrderThreadsPerAtom, "Only 16 threads per atom supported in OpenCL");
+ GMX_ASSERT(!recalculateSplines, "Recalculating splines not supported in OpenCL");
#endif
- const int atomsPerBlock = c_useOrderThreadsPerAtom ? blockSize / c_pmeSpreadGatherThreadsPerAtom4ThPerAtom
- : blockSize / c_pmeSpreadGatherThreadsPerAtom;
+ const int atomsPerBlock = useOrderThreadsPerAtom ? blockSize / c_pmeSpreadGatherThreadsPerAtom4ThPerAtom
+ : blockSize / c_pmeSpreadGatherThreadsPerAtom;
// TODO: pick smaller block size in runtime if needed
// (e.g. on 660 Ti where 50% occupancy is ~25% faster than 100% occupancy with RNAse (~17.8k atoms))
KernelLaunchConfig config;
config.blockSize[0] = order;
- config.blockSize[1] = c_useOrderThreadsPerAtom ? 1 : order;
+ config.blockSize[1] = useOrderThreadsPerAtom ? 1 : order;
config.blockSize[2] = atomsPerBlock;
config.gridSize[0] = dimGrid.first;
config.gridSize[1] = dimGrid.second;
if (spreadCharges)
{
timingId = gtPME_SPLINEANDSPREAD;
- kernelPtr = selectSplineAndSpreadKernelPtr(pmeGpu, c_useOrderThreadsPerAtom,
- writeGlobal || (!c_recalculateSplines));
+ kernelPtr = selectSplineAndSpreadKernelPtr(pmeGpu, useOrderThreadsPerAtom,
+ writeGlobal || (!recalculateSplines));
}
else
{
timingId = gtPME_SPLINE;
- kernelPtr = selectSplineKernelPtr(pmeGpu, c_useOrderThreadsPerAtom,
- writeGlobal || (!c_recalculateSplines));
+ kernelPtr = selectSplineKernelPtr(pmeGpu, useOrderThreadsPerAtom,
+ writeGlobal || (!recalculateSplines));
}
}
else
{
timingId = gtPME_SPREAD;
- kernelPtr = selectSpreadKernelPtr(pmeGpu, c_useOrderThreadsPerAtom,
- writeGlobal || (!c_recalculateSplines));
+ kernelPtr = selectSpreadKernelPtr(pmeGpu, useOrderThreadsPerAtom,
+ writeGlobal || (!recalculateSplines));
}
}
/* Set if we have unit tests */
- const bool readGlobal = pmeGpu->settings.copyAllOutputs;
- const size_t blockSize = pmeGpu->programHandle_->impl_->gatherWorkGroupSize;
+ const bool readGlobal = pmeGpu->settings.copyAllOutputs;
+ const size_t blockSize = pmeGpu->programHandle_->impl_->gatherWorkGroupSize;
+ const bool useOrderThreadsPerAtom = pmeGpu->settings.useOrderThreadsPerAtom;
+ const bool recalculateSplines = pmeGpu->settings.recalculateSplines;
#if GMX_GPU == GMX_GPU_OPENCL
- GMX_ASSERT(!c_useOrderThreadsPerAtom, "Only 16 threads per atom supported in OpenCL");
- GMX_ASSERT(!c_recalculateSplines, "Recalculating splines not supported in OpenCL");
+ GMX_ASSERT(!useOrderThreadsPerAtom, "Only 16 threads per atom supported in OpenCL");
+ GMX_ASSERT(!recalculateSplines, "Recalculating splines not supported in OpenCL");
#endif
- const int atomsPerBlock = c_useOrderThreadsPerAtom ? blockSize / c_pmeSpreadGatherThreadsPerAtom4ThPerAtom
- : blockSize / c_pmeSpreadGatherThreadsPerAtom;
+ const int atomsPerBlock = useOrderThreadsPerAtom ? blockSize / c_pmeSpreadGatherThreadsPerAtom4ThPerAtom
+ : blockSize / c_pmeSpreadGatherThreadsPerAtom;
GMX_ASSERT(!c_usePadding || !(c_pmeAtomDataAlignment % atomsPerBlock),
"inconsistent atom data padding vs. gathering block size");
KernelLaunchConfig config;
config.blockSize[0] = order;
- config.blockSize[1] = c_useOrderThreadsPerAtom ? 1 : order;
+ config.blockSize[1] = useOrderThreadsPerAtom ? 1 : order;
config.blockSize[2] = atomsPerBlock;
config.gridSize[0] = dimGrid.first;
config.gridSize[1] = dimGrid.second;
int timingId = gtPME_GATHER;
PmeGpuProgramImpl::PmeKernelHandle kernelPtr = selectGatherKernelPtr(
- pmeGpu, c_useOrderThreadsPerAtom, readGlobal || (!c_recalculateSplines), forceTreatment);
+ pmeGpu, useOrderThreadsPerAtom, readGlobal || (!recalculateSplines), forceTreatment);
// TODO design kernel selection getters and make PmeGpu a friend of PmeGpuProgramImpl
pme_gpu_start_timing(pmeGpu, timingId);
GpuApiCallBehavior transferKind;
/*! \brief Various flags for the current PME computation, corresponding to the GMX_PME_ flags in pme.h. */
int currentFlags;
+ /*! \brief
+ * Currently only supported by CUDA.
+ * Controls if we should use order (i.e. 4) threads per atom for the GPU
+ * or order*order (i.e. 16) threads per atom.
+ */
+ bool useOrderThreadsPerAtom;
+ /*! \brief
+ * Currently only supported by CUDA.
+ * Controls if we should recalculate the splines in the gather or
+ * save the values in the spread and reload in the gather.
+ */
+ bool recalculateSplines;
};
// TODO There's little value in computing the Coulomb and LJ virial
#include "pme_internal.h"
//! Calculate the slab indices and store in \p atc, store counts in \p count
-static void pme_calc_pidx(int start, int end, const matrix recipbox, const rvec x[], PmeAtomComm* atc, int* count)
+static void pme_calc_pidx(int start,
+ int end,
+ const matrix recipbox,
+ gmx::ArrayRef<const gmx::RVec> x,
+ PmeAtomComm* atc,
+ int* count)
{
int nslab, i;
int si;
try
{
const int natoms = x.ssize();
- pme_calc_pidx(natoms * thread / nthread, natoms * (thread + 1) / nthread, recipbox,
- as_rvec_array(x.data()), atc, atc->count_thread[thread].data());
+ pme_calc_pidx(natoms * thread / nthread, natoms * (thread + 1) / nthread, recipbox, x,
+ atc, atc->count_thread[thread].data());
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
/* Neighbor searching */
printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
- printStringNoNewline(
- &inp, "cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
+ printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
printStringNoNewline(&inp, "nblist update frequency");
ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
energygrps =
; NEIGHBORSEARCHING PARAMETERS
-; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+; cut-off scheme (Verlet: particle based cut-offs)
cutoff-scheme = Verlet
; nblist update frequency
nstlist = 10
return V;
}
+// Avoid gcc 386 -O3 code generation bug in this function (see Redmine
+// #3205 for more information)
+#if defined(__GNUC__) && defined(__i386__) && defined(__OPTIMIZE__)
+# pragma GCC push_options
+# pragma GCC optimize("O2")
+# define avoid_gcc_i386_o3_code_generation_bug
+#endif
+
template<BondedKernelFlavor flavor>
std::enable_if_t<flavor != BondedKernelFlavor::ForcesSimdWhenAvailable || !GMX_SIMD_HAVE_REAL, real>
angles(int nbonds,
return vtot;
}
+#ifdef avoid_gcc_i386_o3_code_generation_bug
+# pragma GCC pop_options
+# undef avoid_gcc_i386_o3_code_generation_bug
+#endif
+
#if GMX_SIMD_HAVE_REAL
/* As angles, but using SIMD to calculate many angles at once.
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/idef.h"
#include "gromacs/utility/strconvert.h"
+#include "gromacs/utility/stringstream.h"
+#include "gromacs/utility/textwriter.h"
#include "testutils/refdata.h"
#include "testutils/testasserts.h"
//! Interaction parameters
t_iparams iparams = { { 0 } };
+ friend std::ostream& operator<<(std::ostream& out, const iListInput& input);
+
//! Constructor
iListInput() {}
}
};
+//! Prints the interaction and parameters to a stream
+std::ostream& operator<<(std::ostream& out, const iListInput& input)
+{
+ using std::endl;
+ out << "Function type " << input.ftype << " called " << interaction_function[input.ftype].name
+ << " ie. labelled '" << interaction_function[input.ftype].longname << "' in an energy file"
+ << endl;
+
+ // Organize to print the legacy C union t_iparams, whose
+ // relevant contents vary with ftype.
+ StringOutputStream stream;
+ {
+ TextWriter writer(&stream);
+ printInteractionParameters(&writer, input.ftype, input.iparams);
+ }
+ out << "Function parameters " << stream.toString();
+ out << "Parameters trigger FEP? " << (input.fep ? "true" : "false") << endl;
+ return out;
+}
+
/*! \brief Utility to fill iatoms struct
*
* \param[in] ftype Function type
BondedKernelFlavor::ForcesAndVirialAndEnergy);
// Internal consistency test of both test input
// and bonded functions.
- EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0)));
+ EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0))) << "dvdlambda was " << output.dvdlambda;
checkOutput(checker, output);
}
void testIfunc()
MultiDimArray<std::vector<float>, dynamicExtents3D> comparedDensity(100, 100, 100);
std::iota(begin(comparedDensity), end(comparedDensity), -10000);
- const real expectedSimilarity = 1;
- EXPECT_REAL_EQ(expectedSimilarity, measure.similarity(comparedDensity.asConstView()));
+ const real expectedSimilarity = 1;
+ FloatingPointTolerance tolerance(relativeToleranceAsUlp(1.0, 100'000));
+ EXPECT_REAL_EQ_TOL(expectedSimilarity, measure.similarity(comparedDensity.asConstView()), tolerance);
}
TEST(DensitySimilarityTest, CrossCorrelationIsMinusOneWhenAntiCorrelated)
MultiDimArray<std::vector<float>, dynamicExtents3D> comparedDensity(100, 100, 100);
std::iota(begin(comparedDensity), end(comparedDensity), -10000);
- const real expectedSimilarity = -1;
- EXPECT_REAL_EQ(expectedSimilarity, measure.similarity(comparedDensity.asConstView()));
+ const real expectedSimilarity = -1;
+ FloatingPointTolerance tolerance(relativeToleranceAsUlp(-1.0, 100'000));
+ EXPECT_REAL_EQ_TOL(expectedSimilarity, measure.similarity(comparedDensity.asConstView()), tolerance);
}
TEST(DensitySimilarityTest, CrossCorrelationGradientIsZeroWhenCorrelated)
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pbcutil/pbc_aiuc_cuda.cuh"
#include "gromacs/topology/ifunc.h"
+#include "gromacs/topology/topology.h"
namespace gmx
{
}
}
+/*! \brief Constructs and returns an atom constraint adjacency list
+ *
+ * Each constraint will be represented as a tuple, containing index of the second
+ * constrained atom, index of the constraint and a sign that indicates the order of atoms in
+ * which they are listed. Sign is needed to compute the mass factors.
+ */
+static std::vector<std::vector<std::tuple<int, int, int>>>
+constructAtomsAdjacencyList(const int numAtoms, ArrayRef<const int> iatoms)
+{
+ const int stride = 1 + NRAL(F_CONSTR);
+ const int numConstraints = iatoms.ssize() / stride;
+ std::vector<std::vector<std::tuple<int, int, int>>> atomsAdjacencyList(numAtoms);
+ for (int c = 0; c < numConstraints; c++)
+ {
+ int a1 = iatoms[stride * c + 1];
+ int a2 = iatoms[stride * c + 2];
+
+ // Each constraint will be represented as a tuple, containing index of the second
+ // constrained atom, index of the constraint and a sign that indicates the order of atoms in
+ // which they are listed. Sign is needed to compute the mass factors.
+ atomsAdjacencyList[a1].push_back(std::make_tuple(a2, c, +1));
+ atomsAdjacencyList[a2].push_back(std::make_tuple(a1, c, -1));
+ }
+
+ return atomsAdjacencyList;
+}
+
/*! \brief Helper function to go through constraints recursively.
*
- * For each constraint, counts the number of coupled constraints and stores the value in numCoupledConstraints array.
+ * For each constraint, counts the number of coupled constraints and stores the value in \p numCoupledConstraints array.
* This information is used to split the array of constraints between thread blocks on a GPU so there is no
- * coupling between constraints from different thread blocks. After the 'numCoupledConstraints' array is filled, the
- * value numCoupledConstraints[c] should be equal to the number of constraints that are coupled to 'c' and located
+ * coupling between constraints from different thread blocks. After the \p numCoupledConstraints array is filled, the
+ * value \p numCoupledConstraints[c] should be equal to the number of constraints that are coupled to \p c and located
* after it in the constraints array.
*
* \param[in] a Atom index.
* \param[in] c Sequential index for constraint to consider adding.
* \param[in,out] currentMapIndex The rolling index for the constraints mapping.
*/
-inline void addWithCoupled(const t_iatom* iatoms,
+inline void addWithCoupled(ArrayRef<const int> iatoms,
const int stride,
ArrayRef<const std::vector<std::tuple<int, int, int>>> atomsAdjacencyList,
ArrayRef<int> splitMap,
}
}
+/*! \brief Computes and returns how many constraints are coupled to each constraint
+ *
+ * Needed to introduce splits in data so that all coupled constraints will be computed in a
+ * single GPU block. The position \p c of the vector \p numCoupledConstraints should have the number
+ * of constraints that are coupled to a constraint \p c and are after \p c in the vector. Only
+ * first index of the connected group of the constraints is needed later in the code, hence the
+ * numCoupledConstraints vector is also used to keep track if the constrain was already counted.
+ */
+static std::vector<int> countNumCoupledConstraints(ArrayRef<const int> iatoms,
+ ArrayRef<const std::vector<std::tuple<int, int, int>>> atomsAdjacencyList)
+{
+ const int stride = 1 + NRAL(F_CONSTR);
+ const int numConstraints = iatoms.ssize() / stride;
+ std::vector<int> numCoupledConstraints(numConstraints, -1);
+ for (int c = 0; c < numConstraints; c++)
+ {
+ const int a1 = iatoms[stride * c + 1];
+ const int a2 = iatoms[stride * c + 2];
+ if (numCoupledConstraints[c] == -1)
+ {
+ numCoupledConstraints[c] = countCoupled(a1, numCoupledConstraints, atomsAdjacencyList)
+ + countCoupled(a2, numCoupledConstraints, atomsAdjacencyList);
+ }
+ }
+
+ return numCoupledConstraints;
+}
+
+bool LincsCuda::isNumCoupledConstraintsSupported(const gmx_mtop_t& mtop)
+{
+ for (const gmx_moltype_t& molType : mtop.moltype)
+ {
+ ArrayRef<const int> iatoms = molType.ilist[F_CONSTR].iatoms;
+ const auto atomsAdjacencyList = constructAtomsAdjacencyList(molType.atoms.nr, iatoms);
+ // Compute, how many constraints are coupled to each constraint
+ const auto numCoupledConstraints = countNumCoupledConstraints(iatoms, atomsAdjacencyList);
+ for (const int numCoupled : numCoupledConstraints)
+ {
+ if (numCoupled > c_threadsPerBlock)
+ {
+ return false;
+ }
+ }
+ }
+
+ return true;
+}
+
void LincsCuda::set(const t_idef& idef, const t_mdatoms& md)
{
int numAtoms = md.nr;
std::vector<float> massFactorsHost;
// List of constrained atoms in local topology
- t_iatom* iatoms = idef.il[F_CONSTR].iatoms;
+ gmx::ArrayRef<const int> iatoms =
+ constArrayRefFromArray(idef.il[F_CONSTR].iatoms, idef.il[F_CONSTR].nr);
const int stride = NRAL(F_CONSTR) + 1;
const int numConstraints = idef.il[F_CONSTR].nr / stride;
return;
}
- // Constructing adjacency list --- usefull intermediate structure
- std::vector<std::vector<std::tuple<int, int, int>>> atomsAdjacencyList(numAtoms);
- for (int c = 0; c < numConstraints; c++)
- {
- int a1 = iatoms[stride * c + 1];
- int a2 = iatoms[stride * c + 2];
-
- // Each constraint will be represented as a tuple, containing index of the second
- // constrained atom, index of the constraint and a sign that indicates the order of atoms in
- // which they are listed. Sign is needed to compute the mass factors.
- atomsAdjacencyList.at(a1).push_back(std::make_tuple(a2, c, +1));
- atomsAdjacencyList.at(a2).push_back(std::make_tuple(a1, c, -1));
- }
+ // Construct the adjacency list, a useful intermediate structure
+ const auto atomsAdjacencyList = constructAtomsAdjacencyList(numAtoms, iatoms);
- // Compute, how many constraints are coupled to each constraint.
- // Needed to introduce splits in data so that all coupled constraints will be computed in a
- // single GPU block. The position 'c' of the vector numCoupledConstraints should have the number
- // of constraints that are coupled to a constraint 'c' and are after 'c' in the vector. Only
- // first index of the connected group of the constraints is needed later in the code, hence the
- // numCoupledConstraints vector is also used to keep track if the constrain was already counted.
- std::vector<int> numCoupledConstraints(numConstraints, -1);
- for (int c = 0; c < numConstraints; c++)
- {
- int a1 = iatoms[stride * c + 1];
- int a2 = iatoms[stride * c + 2];
- if (numCoupledConstraints.at(c) == -1)
- {
- numCoupledConstraints.at(c) = countCoupled(a1, numCoupledConstraints, atomsAdjacencyList)
- + countCoupled(a2, numCoupledConstraints, atomsAdjacencyList);
- }
- }
+ // Compute, how many constraints are coupled to each constraint
+ const auto numCoupledConstraints = countNumCoupledConstraints(iatoms, atomsAdjacencyList);
// Map of splits in the constraints data. For each 'old' constraint index gives 'new' which
// takes into account the empty spaces which might be needed in the end of each thread block.
* Applies LINCS to coordinates and velocities, stored on GPU.
* The results are not automatically copied back to the CPU memory.
* Method uses this class data structures which should be updated
- * when needed using set() and setPbc() method.
+ * when needed using set() method.
*
* \param[in] d_x Coordinates before timestep (in GPU memory)
* \param[in,out] d_xp Coordinates after timestep (in GPU memory). The
*/
void set(const t_idef& idef, const t_mdatoms& md);
+ /*! \brief
+ * Returns whether the maximum number of coupled constraints is supported
+ * by the CUDA LINCS code.
+ *
+ * \param[in] mtop The molecular topology
+ */
+ static bool isNumCoupledConstraintsSupported(const gmx_mtop_t& mtop);
+
private:
//! CUDA stream
CommandStream commandStream_;
*/
GpuEventSynchronizer* getCoordinatesReadySync();
+ /*! \brief
+ * Returns whether the maximum number of coupled constraints is supported
+ * by the CUDA LINCS code.
+ *
+ * \param[in] mtop The molecular topology
+ */
+ static bool isNumCoupledConstraintsSupported(const gmx_mtop_t& mtop);
+
private:
class Impl;
gmx::PrivateImplPointer<Impl> impl_;
return nullptr;
}
+bool UpdateConstrainCuda::isNumCoupledConstraintsSupported(const gmx_mtop_t& /* mtop */)
+{
+ return false;
+}
+
} // namespace gmx
#endif /* GMX_GPU != GMX_GPU_CUDA */
return impl_->getCoordinatesReadySync();
}
+bool UpdateConstrainCuda::isNumCoupledConstraintsSupported(const gmx_mtop_t& mtop)
+{
+ return LincsCuda::isNumCoupledConstraintsSupported(mtop);
+}
+
} // namespace gmx
*/
GpuEventSynchronizer* getCoordinatesReadySync();
+ /*! \brief
+ * Returns whether the maximum number of coupled constraints is supported
+ * by the CUDA LINCS code.
+ *
+ * \param[in] mtop The molecular topology
+ */
+ static bool isNumCoupledConstraintsSupported(const gmx_mtop_t& mtop);
+
private:
//! CUDA stream
CommandStream commandStream_ = nullptr;
//! True if the Buffer ops development feature is enabled
// TODO: when the trigger of the buffer ops offload is fully automated this should go away
bool enableGpuBufferOps = false;
- //! If true, forces 'mdrun -update auto' default to 'gpu'
- bool forceGpuUpdateDefaultOn = false;
//! True if the GPU halo exchange development feature is enabled
bool enableGpuHaloExchange = false;
//! True if the PME PP direct communication GPU development feature is enabled
#pragma GCC diagnostic ignored "-Wunused-result"
devFlags.enableGpuBufferOps = (getenv("GMX_USE_GPU_BUFFER_OPS") != nullptr)
&& (GMX_GPU == GMX_GPU_CUDA) && useGpuForNonbonded;
- devFlags.forceGpuUpdateDefaultOn = (getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") != nullptr);
devFlags.enableGpuHaloExchange =
(getenv("GMX_GPU_DD_COMMS") != nullptr && GMX_THREAD_MPI && (GMX_GPU == GMX_GPU_CUDA));
devFlags.enableGpuPmePPComm =
}
}
- if (devFlags.forceGpuUpdateDefaultOn)
- {
- GMX_LOG(mdlog.warning)
- .asParagraph()
- .appendTextFormatted(
- "NOTE: This run will default to '-update gpu' as requested by the "
- "GMX_FORCE_UPDATE_DEFAULT_GPU environment variable.");
- }
-
if (devFlags.enableGpuPmePPComm)
{
if (pmeRunMode == PmeRunMode::GPU)
try
{
useGpuForUpdate = decideWhetherToUseGpuForUpdate(
- devFlags.forceGpuUpdateDefaultOn, useDomainDecomposition, useGpuForPme,
- useGpuForNonbonded, updateTarget, gpusWereDetected, *inputrec,
- gmx_mtop_interaction_count(mtop, IF_VSITE) > 0, doEssentialDynamics,
+ useDomainDecomposition, useGpuForPme, useGpuForNonbonded, updateTarget,
+ gpusWereDetected, *inputrec, mtop, doEssentialDynamics,
gmx_mtop_ftype_count(mtop, F_ORIRES) > 0, replExParams.exchangeInterval > 0);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
#include "gromacs/hardware/hardwaretopology.h"
#include "gromacs/hardware/hw_info.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
+#include "gromacs/mdlib/update_constrain_cuda.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdrunoptions.h"
#include "gromacs/taskassignment/taskassignment.h"
+#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/baseversion.h"
#include "gromacs/utility/exceptions.h"
return gpusWereDetected && usingOurCpuForPmeOrEwald;
}
-bool decideWhetherToUseGpuForUpdate(const bool forceGpuUpdateDefaultOn,
- const bool isDomainDecomposition,
+bool decideWhetherToUseGpuForUpdate(const bool isDomainDecomposition,
const bool useGpuForPme,
const bool useGpuForNonbonded,
const TaskTarget updateTarget,
const bool gpusWereDetected,
const t_inputrec& inputrec,
- const bool haveVSites,
+ const gmx_mtop_t& mtop,
const bool useEssentialDynamics,
const bool doOrientationRestraints,
const bool useReplicaExchange)
// The graph is needed, but not supported
errorMessage += "Ewald surface correction is not supported.\n";
}
- if (haveVSites)
+ if (gmx_mtop_interaction_count(mtop, IF_VSITE) > 0)
{
errorMessage += "Virtual sites are not supported.\n";
}
errorMessage += "Swapping the coordinates is not supported.\n";
}
- // \todo Check for coupled constraint block size restriction needs to be added
- // when update auto chooses GPU in some cases. Currently exceeding the restriction
- // triggers a fatal error during LINCS setup.
+ // TODO: F_CONSTRNC is only unsupported, because isNumCoupledConstraintsSupported()
+ // does not support it, the actual CUDA LINCS code does support it
+ if (gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0)
+ {
+ errorMessage += "Non-connecting constraints are not supported";
+ }
+ if (!UpdateConstrainCuda::isNumCoupledConstraintsSupported(mtop))
+ {
+ errorMessage +=
+ "The number of coupled constraints is higher than supported in the CUDA LINCS "
+ "code.\n";
+ }
if (!errorMessage.empty())
{
return false;
}
- return ((forceGpuUpdateDefaultOn && updateTarget == TaskTarget::Auto)
- || (updateTarget == TaskTarget::Gpu));
+ return true;
}
} // namespace gmx
/*! \brief Decide whether to use GPU for update.
*
- * \param[in] forceGpuUpdateDefaultOn If the update should be offloaded by default.
* \param[in] isDomainDecomposition Whether there more than one domain.
* \param[in] useGpuForPme Whether GPUs will be used for PME interactions.
* \param[in] useGpuForNonbonded Whether GPUs will be used for nonbonded interactions.
* \param[in] updateTarget User choice for running simulation on GPU.
* \param[in] gpusWereDetected Whether compatible GPUs were detected on any node.
* \param[in] inputrec The user input.
- * \param[in] haveVSites If there are virtual sites in the system.
+ * \param[in] mtop The global topology.
* \param[in] useEssentialDynamics If essential dynamics is active.
* \param[in] doOrientationRestraints If orientation restraints are enabled.
* \param[in] useReplicaExchange If this is a REMD simulation.
* \throws std::bad_alloc If out of memory
* InconsistentInputError If the user requirements are inconsistent.
*/
-bool decideWhetherToUseGpuForUpdate(bool forceGpuUpdateDefaultOn,
- bool isDomainDecomposition,
+bool decideWhetherToUseGpuForUpdate(bool isDomainDecomposition,
bool useGpuForPme,
bool useGpuForNonbonded,
TaskTarget updateTarget,
bool gpusWereDetected,
const t_inputrec& inputrec,
- bool haveVSites,
+ const gmx_mtop_t& mtop,
bool useEssentialDynamics,
bool doOrientationRestraints,
bool useReplicaExchange);
// Avoid all ranks spamming the error stream
//
// TODO improve this so that unique errors on different ranks
- // are all reported.
+ // are all reported on one rank.
if (countOfExceptionsOnThisRank > 0 && physicalNodeComm.rank_ == 0)
{
try
std::rethrow_exception(exceptionPtr);
}
}
- catch (const std::exception& ex)
- {
- printFatalErrorMessage(stderr, ex);
- }
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
+ // TODO This implements a global barrier so that MPI runtimes can
+ // organize an orderly shutdown if one of the ranks has had to
+ // issue a fatal error above. When we have MPI-aware error
+ // handling and reporting, this should be improved (perhaps
+ // centralized there).
+ simulationBarrier(cr);
+ multiSimBarrier(ms);
+ simulationBarrier(cr);
if (countOfExceptionsOverAllRanks > 0)
{
gmx_fatal(FARGS,
- "Exiting because task assignment failed. If there is no descriptive "
- "error message above this, please report this failure as a bug.");
+ "Exiting because task assignment failed. If there is no descriptive error "
+ "message in the terminal output, please report this failure as a bug.");
}
// TODO There is no check that mdrun -nb gpu or -pme gpu or
pr_indent(fp, indent + INDENT);
fprintf(fp, "functype[%d]=%s, ", bShowNumbers ? i : -1,
interaction_function[ffparams->functype[i]].name);
- pr_iparams(fp, ffparams->functype[i], &ffparams->iparams[i]);
+ pr_iparams(fp, ffparams->functype[i], ffparams->iparams[i]);
}
pr_double(fp, indent, "reppow", ffparams->reppow);
pr_real(fp, indent, "fudgeQQ", ffparams->fudgeQQ);
#include "gromacs/topology/ifunc.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringstream.h"
+#include "gromacs/utility/textwriter.h"
#include "gromacs/utility/txtdump.h"
-static void pr_harm(FILE* fp, const t_iparams* iparams, const char* r, const char* kr)
+static void printHarmonicInteraction(gmx::TextWriter* writer,
+ const t_iparams& iparams,
+ const char* r,
+ const char* kr)
{
- fprintf(fp, "%sA=%12.5e, %sA=%12.5e, %sB=%12.5e, %sB=%12.5e\n", r, iparams->harmonic.rA, kr,
- iparams->harmonic.krA, r, iparams->harmonic.rB, kr, iparams->harmonic.krB);
+ writer->writeLineFormatted("%sA=%12.5e, %sA=%12.5e, %sB=%12.5e, %sB=%12.5e", r,
+ iparams.harmonic.rA, kr, iparams.harmonic.krA, r,
+ iparams.harmonic.rB, kr, iparams.harmonic.krB);
}
-void pr_iparams(FILE* fp, t_functype ftype, const t_iparams* iparams)
+void pr_iparams(FILE* fp, t_functype ftype, const t_iparams& iparams)
+{
+ gmx::StringOutputStream stream;
+ {
+ gmx::TextWriter writer(&stream);
+ printInteractionParameters(&writer, ftype, iparams);
+ }
+ fputs(stream.toString().c_str(), fp);
+}
+
+void printInteractionParameters(gmx::TextWriter* writer, t_functype ftype, const t_iparams& iparams)
{
switch (ftype)
{
case F_ANGLES:
- case F_G96ANGLES: pr_harm(fp, iparams, "th", "ct"); break;
+ case F_G96ANGLES: printHarmonicInteraction(writer, iparams, "th", "ct"); break;
case F_CROSS_BOND_BONDS:
- fprintf(fp, "r1e=%15.8e, r2e=%15.8e, krr=%15.8e\n", iparams->cross_bb.r1e,
- iparams->cross_bb.r2e, iparams->cross_bb.krr);
+ writer->writeLineFormatted("r1e=%15.8e, r2e=%15.8e, krr=%15.8e", iparams.cross_bb.r1e,
+ iparams.cross_bb.r2e, iparams.cross_bb.krr);
break;
case F_CROSS_BOND_ANGLES:
- fprintf(fp, "r1e=%15.8e, r1e=%15.8e, r3e=%15.8e, krt=%15.8e\n", iparams->cross_ba.r1e,
- iparams->cross_ba.r2e, iparams->cross_ba.r3e, iparams->cross_ba.krt);
+ writer->writeLineFormatted("r1e=%15.8e, r1e=%15.8e, r3e=%15.8e, krt=%15.8e",
+ iparams.cross_ba.r1e, iparams.cross_ba.r2e,
+ iparams.cross_ba.r3e, iparams.cross_ba.krt);
break;
case F_LINEAR_ANGLES:
- fprintf(fp, "klinA=%15.8e, aA=%15.8e, klinB=%15.8e, aB=%15.8e\n", iparams->linangle.klinA,
- iparams->linangle.aA, iparams->linangle.klinB, iparams->linangle.aB);
+ writer->writeLineFormatted("klinA=%15.8e, aA=%15.8e, klinB=%15.8e, aB=%15.8e",
+ iparams.linangle.klinA, iparams.linangle.aA,
+ iparams.linangle.klinB, iparams.linangle.aB);
break;
case F_UREY_BRADLEY:
- fprintf(fp,
+ writer->writeLineFormatted(
"thetaA=%15.8e, kthetaA=%15.8e, r13A=%15.8e, kUBA=%15.8e, thetaB=%15.8e, "
- "kthetaB=%15.8e, r13B=%15.8e, kUBB=%15.8e\n",
- iparams->u_b.thetaA, iparams->u_b.kthetaA, iparams->u_b.r13A, iparams->u_b.kUBA,
- iparams->u_b.thetaB, iparams->u_b.kthetaB, iparams->u_b.r13B, iparams->u_b.kUBB);
+ "kthetaB=%15.8e, r13B=%15.8e, kUBB=%15.8e",
+ iparams.u_b.thetaA, iparams.u_b.kthetaA, iparams.u_b.r13A, iparams.u_b.kUBA,
+ iparams.u_b.thetaB, iparams.u_b.kthetaB, iparams.u_b.r13B, iparams.u_b.kUBB);
break;
case F_QUARTIC_ANGLES:
- fprintf(fp, "theta=%15.8e", iparams->qangle.theta);
+ writer->writeStringFormatted("theta=%15.8e", iparams.qangle.theta);
for (int i = 0; i < 5; i++)
{
- fprintf(fp, ", c%c=%15.8e", '0' + i, iparams->qangle.c[i]);
+ writer->writeStringFormatted(", c%c=%15.8e", '0' + i, iparams.qangle.c[i]);
}
- fprintf(fp, "\n");
+ writer->ensureLineBreak();
break;
case F_BHAM:
- fprintf(fp, "a=%15.8e, b=%15.8e, c=%15.8e\n", iparams->bham.a, iparams->bham.b,
- iparams->bham.c);
+ writer->writeLineFormatted("a=%15.8e, b=%15.8e, c=%15.8e", iparams.bham.a,
+ iparams.bham.b, iparams.bham.c);
break;
case F_BONDS:
case F_G96BONDS:
- case F_HARMONIC: pr_harm(fp, iparams, "b0", "cb"); break;
- case F_IDIHS: pr_harm(fp, iparams, "xi", "cx"); break;
+ case F_HARMONIC: printHarmonicInteraction(writer, iparams, "b0", "cb"); break;
+ case F_IDIHS: printHarmonicInteraction(writer, iparams, "xi", "cx"); break;
case F_MORSE:
- fprintf(fp,
- "b0A=%15.8e, cbA=%15.8e, betaA=%15.8e, b0B=%15.8e, cbB=%15.8e, betaB=%15.8e\n",
- iparams->morse.b0A, iparams->morse.cbA, iparams->morse.betaA,
- iparams->morse.b0B, iparams->morse.cbB, iparams->morse.betaB);
+ writer->writeLineFormatted(
+ "b0A=%15.8e, cbA=%15.8e, betaA=%15.8e, b0B=%15.8e, cbB=%15.8e, betaB=%15.8e",
+ iparams.morse.b0A, iparams.morse.cbA, iparams.morse.betaA, iparams.morse.b0B,
+ iparams.morse.cbB, iparams.morse.betaB);
break;
case F_CUBICBONDS:
- fprintf(fp, "b0=%15.8e, kb=%15.8e, kcub=%15.8e\n", iparams->cubic.b0, iparams->cubic.kb,
- iparams->cubic.kcub);
+ writer->writeLineFormatted("b0=%15.8e, kb=%15.8e, kcub=%15.8e", iparams.cubic.b0,
+ iparams.cubic.kb, iparams.cubic.kcub);
break;
- case F_CONNBONDS: fprintf(fp, "\n"); break;
+ case F_CONNBONDS: writer->ensureEmptyLine(); break;
case F_FENEBONDS:
- fprintf(fp, "bm=%15.8e, kb=%15.8e\n", iparams->fene.bm, iparams->fene.kb);
+ writer->writeLineFormatted("bm=%15.8e, kb=%15.8e", iparams.fene.bm, iparams.fene.kb);
break;
case F_RESTRBONDS:
- fprintf(fp,
+ writer->writeLineFormatted(
"lowA=%15.8e, up1A=%15.8e, up2A=%15.8e, kA=%15.8e, lowB=%15.8e, up1B=%15.8e, "
- "up2B=%15.8e, kB=%15.8e,\n",
- iparams->restraint.lowA, iparams->restraint.up1A, iparams->restraint.up2A,
- iparams->restraint.kA, iparams->restraint.lowB, iparams->restraint.up1B,
- iparams->restraint.up2B, iparams->restraint.kB);
+ "up2B=%15.8e, kB=%15.8e,",
+ iparams.restraint.lowA, iparams.restraint.up1A, iparams.restraint.up2A,
+ iparams.restraint.kA, iparams.restraint.lowB, iparams.restraint.up1B,
+ iparams.restraint.up2B, iparams.restraint.kB);
break;
case F_TABBONDS:
case F_TABBONDSNC:
case F_TABANGLES:
case F_TABDIHS:
- fprintf(fp, "tab=%d, kA=%15.8e, kB=%15.8e\n", iparams->tab.table, iparams->tab.kA,
- iparams->tab.kB);
+ writer->writeLineFormatted("tab=%d, kA=%15.8e, kB=%15.8e", iparams.tab.table,
+ iparams.tab.kA, iparams.tab.kB);
+ break;
+ case F_POLARIZATION:
+ writer->writeLineFormatted("alpha=%15.8e", iparams.polarize.alpha);
break;
- case F_POLARIZATION: fprintf(fp, "alpha=%15.8e\n", iparams->polarize.alpha); break;
case F_ANHARM_POL:
- fprintf(fp, "alpha=%15.8e drcut=%15.8e khyp=%15.8e\n", iparams->anharm_polarize.alpha,
- iparams->anharm_polarize.drcut, iparams->anharm_polarize.khyp);
+ writer->writeLineFormatted("alpha=%15.8e drcut=%15.8e khyp=%15.8e",
+ iparams.anharm_polarize.alpha, iparams.anharm_polarize.drcut,
+ iparams.anharm_polarize.khyp);
break;
case F_THOLE_POL:
- fprintf(fp, "a=%15.8e, alpha1=%15.8e, alpha2=%15.8e, rfac=%15.8e\n", iparams->thole.a,
- iparams->thole.alpha1, iparams->thole.alpha2, iparams->thole.rfac);
+ writer->writeLineFormatted("a=%15.8e, alpha1=%15.8e, alpha2=%15.8e, rfac=%15.8e",
+ iparams.thole.a, iparams.thole.alpha1, iparams.thole.alpha2,
+ iparams.thole.rfac);
break;
case F_WATER_POL:
- fprintf(fp, "al_x=%15.8e, al_y=%15.8e, al_z=%15.8e, rOH=%9.6f, rHH=%9.6f, rOD=%9.6f\n",
- iparams->wpol.al_x, iparams->wpol.al_y, iparams->wpol.al_z, iparams->wpol.rOH,
- iparams->wpol.rHH, iparams->wpol.rOD);
+ writer->writeLineFormatted(
+ "al_x=%15.8e, al_y=%15.8e, al_z=%15.8e, rOH=%9.6f, rHH=%9.6f, rOD=%9.6f",
+ iparams.wpol.al_x, iparams.wpol.al_y, iparams.wpol.al_z, iparams.wpol.rOH,
+ iparams.wpol.rHH, iparams.wpol.rOD);
+ break;
+ case F_LJ:
+ writer->writeLineFormatted("c6=%15.8e, c12=%15.8e", iparams.lj.c6, iparams.lj.c12);
break;
- case F_LJ: fprintf(fp, "c6=%15.8e, c12=%15.8e\n", iparams->lj.c6, iparams->lj.c12); break;
case F_LJ14:
- fprintf(fp, "c6A=%15.8e, c12A=%15.8e, c6B=%15.8e, c12B=%15.8e\n", iparams->lj14.c6A,
- iparams->lj14.c12A, iparams->lj14.c6B, iparams->lj14.c12B);
+ writer->writeLineFormatted("c6A=%15.8e, c12A=%15.8e, c6B=%15.8e, c12B=%15.8e",
+ iparams.lj14.c6A, iparams.lj14.c12A, iparams.lj14.c6B,
+ iparams.lj14.c12B);
break;
case F_LJC14_Q:
- fprintf(fp, "fqq=%15.8e, qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e\n", iparams->ljc14.fqq,
- iparams->ljc14.qi, iparams->ljc14.qj, iparams->ljc14.c6, iparams->ljc14.c12);
+ writer->writeLineFormatted("fqq=%15.8e, qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e",
+ iparams.ljc14.fqq, iparams.ljc14.qi, iparams.ljc14.qj,
+ iparams.ljc14.c6, iparams.ljc14.c12);
break;
case F_LJC_PAIRS_NB:
- fprintf(fp, "qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e\n", iparams->ljcnb.qi,
- iparams->ljcnb.qj, iparams->ljcnb.c6, iparams->ljcnb.c12);
+ writer->writeLineFormatted("qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e", iparams.ljcnb.qi,
+ iparams.ljcnb.qj, iparams.ljcnb.c6, iparams.ljcnb.c12);
break;
case F_PDIHS:
case F_PIDIHS:
case F_ANGRES:
case F_ANGRESZ:
- fprintf(fp, "phiA=%15.8e, cpA=%15.8e, phiB=%15.8e, cpB=%15.8e, mult=%d\n",
- iparams->pdihs.phiA, iparams->pdihs.cpA, iparams->pdihs.phiB,
- iparams->pdihs.cpB, iparams->pdihs.mult);
+ writer->writeLineFormatted("phiA=%15.8e, cpA=%15.8e, phiB=%15.8e, cpB=%15.8e, mult=%d",
+ iparams.pdihs.phiA, iparams.pdihs.cpA, iparams.pdihs.phiB,
+ iparams.pdihs.cpB, iparams.pdihs.mult);
break;
case F_DISRES:
- fprintf(fp, "label=%4d, type=%1d, low=%15.8e, up1=%15.8e, up2=%15.8e, fac=%15.8e)\n",
- iparams->disres.label, iparams->disres.type, iparams->disres.low,
- iparams->disres.up1, iparams->disres.up2, iparams->disres.kfac);
+ writer->writeLineFormatted(
+ "label=%4d, type=%1d, low=%15.8e, up1=%15.8e, up2=%15.8e, fac=%15.8e)",
+ iparams.disres.label, iparams.disres.type, iparams.disres.low,
+ iparams.disres.up1, iparams.disres.up2, iparams.disres.kfac);
break;
case F_ORIRES:
- fprintf(fp, "ex=%4d, label=%d, power=%4d, c=%15.8e, obs=%15.8e, kfac=%15.8e)\n",
- iparams->orires.ex, iparams->orires.label, iparams->orires.power,
- iparams->orires.c, iparams->orires.obs, iparams->orires.kfac);
+ writer->writeLineFormatted(
+ "ex=%4d, label=%d, power=%4d, c=%15.8e, obs=%15.8e, kfac=%15.8e)",
+ iparams.orires.ex, iparams.orires.label, iparams.orires.power, iparams.orires.c,
+ iparams.orires.obs, iparams.orires.kfac);
break;
case F_DIHRES:
- fprintf(fp,
+ writer->writeLineFormatted(
"phiA=%15.8e, dphiA=%15.8e, kfacA=%15.8e, phiB=%15.8e, dphiB=%15.8e, "
- "kfacB=%15.8e\n",
- iparams->dihres.phiA, iparams->dihres.dphiA, iparams->dihres.kfacA,
- iparams->dihres.phiB, iparams->dihres.dphiB, iparams->dihres.kfacB);
+ "kfacB=%15.8e",
+ iparams.dihres.phiA, iparams.dihres.dphiA, iparams.dihres.kfacA,
+ iparams.dihres.phiB, iparams.dihres.dphiB, iparams.dihres.kfacB);
break;
case F_POSRES:
- fprintf(fp,
+ writer->writeLineFormatted(
"pos0A=(%15.8e,%15.8e,%15.8e), fcA=(%15.8e,%15.8e,%15.8e), "
- "pos0B=(%15.8e,%15.8e,%15.8e), fcB=(%15.8e,%15.8e,%15.8e)\n",
- iparams->posres.pos0A[XX], iparams->posres.pos0A[YY], iparams->posres.pos0A[ZZ],
- iparams->posres.fcA[XX], iparams->posres.fcA[YY], iparams->posres.fcA[ZZ],
- iparams->posres.pos0B[XX], iparams->posres.pos0B[YY], iparams->posres.pos0B[ZZ],
- iparams->posres.fcB[XX], iparams->posres.fcB[YY], iparams->posres.fcB[ZZ]);
+ "pos0B=(%15.8e,%15.8e,%15.8e), fcB=(%15.8e,%15.8e,%15.8e)",
+ iparams.posres.pos0A[XX], iparams.posres.pos0A[YY], iparams.posres.pos0A[ZZ],
+ iparams.posres.fcA[XX], iparams.posres.fcA[YY], iparams.posres.fcA[ZZ],
+ iparams.posres.pos0B[XX], iparams.posres.pos0B[YY], iparams.posres.pos0B[ZZ],
+ iparams.posres.fcB[XX], iparams.posres.fcB[YY], iparams.posres.fcB[ZZ]);
break;
case F_FBPOSRES:
- fprintf(fp, "pos0=(%15.8e,%15.8e,%15.8e), geometry=%d, r=%15.8e, k=%15.8e\n",
- iparams->fbposres.pos0[XX], iparams->fbposres.pos0[YY], iparams->fbposres.pos0[ZZ],
- iparams->fbposres.geom, iparams->fbposres.r, iparams->fbposres.k);
+ writer->writeLineFormatted(
+ "pos0=(%15.8e,%15.8e,%15.8e), geometry=%d, r=%15.8e, k=%15.8e",
+ iparams.fbposres.pos0[XX], iparams.fbposres.pos0[YY], iparams.fbposres.pos0[ZZ],
+ iparams.fbposres.geom, iparams.fbposres.r, iparams.fbposres.k);
break;
case F_RBDIHS:
for (int i = 0; i < NR_RBDIHS; i++)
{
- fprintf(fp, "%srbcA[%d]=%15.8e", i == 0 ? "" : ", ", i, iparams->rbdihs.rbcA[i]);
+ writer->writeStringFormatted("%srbcA[%d]=%15.8e", i == 0 ? "" : ", ", i,
+ iparams.rbdihs.rbcA[i]);
}
- fprintf(fp, "\n");
+ writer->ensureLineBreak();
for (int i = 0; i < NR_RBDIHS; i++)
{
- fprintf(fp, "%srbcB[%d]=%15.8e", i == 0 ? "" : ", ", i, iparams->rbdihs.rbcB[i]);
+ writer->writeStringFormatted("%srbcB[%d]=%15.8e", i == 0 ? "" : ", ", i,
+ iparams.rbdihs.rbcB[i]);
}
- fprintf(fp, "\n");
+ writer->ensureLineBreak();
break;
case F_FOURDIHS:
{
/* Use the OPLS -> Ryckaert-Bellemans formula backwards to get
* the OPLS potential constants back.
*/
- const real* rbcA = iparams->rbdihs.rbcA;
- const real* rbcB = iparams->rbdihs.rbcB;
+ const real* rbcA = iparams.rbdihs.rbcA;
+ const real* rbcB = iparams.rbdihs.rbcB;
real VA[4], VB[4];
VA[3] = -0.25 * rbcA[4];
for (int i = 0; i < NR_FOURDIHS; i++)
{
- fprintf(fp, "%sFourA[%d]=%15.8e", i == 0 ? "" : ", ", i, VA[i]);
+ writer->writeStringFormatted("%sFourA[%d]=%15.8e", i == 0 ? "" : ", ", i, VA[i]);
}
- fprintf(fp, "\n");
+ writer->ensureLineBreak();
for (int i = 0; i < NR_FOURDIHS; i++)
{
- fprintf(fp, "%sFourB[%d]=%15.8e", i == 0 ? "" : ", ", i, VB[i]);
+ writer->writeStringFormatted("%sFourB[%d]=%15.8e", i == 0 ? "" : ", ", i, VB[i]);
}
- fprintf(fp, "\n");
+ writer->ensureLineBreak();
break;
}
case F_CONSTR:
case F_CONSTRNC:
- fprintf(fp, "dA=%15.8e, dB=%15.8e\n", iparams->constr.dA, iparams->constr.dB);
+ writer->writeLineFormatted("dA=%15.8e, dB=%15.8e", iparams.constr.dA, iparams.constr.dB);
break;
case F_SETTLE:
- fprintf(fp, "doh=%15.8e, dhh=%15.8e\n", iparams->settle.doh, iparams->settle.dhh);
+ writer->writeLineFormatted("doh=%15.8e, dhh=%15.8e", iparams.settle.doh, iparams.settle.dhh);
break;
- case F_VSITE2: fprintf(fp, "a=%15.8e\n", iparams->vsite.a); break;
+ case F_VSITE2: writer->writeLineFormatted("a=%15.8e", iparams.vsite.a); break;
case F_VSITE3:
case F_VSITE3FD:
case F_VSITE3FAD:
- fprintf(fp, "a=%15.8e, b=%15.8e\n", iparams->vsite.a, iparams->vsite.b);
+ writer->writeLineFormatted("a=%15.8e, b=%15.8e", iparams.vsite.a, iparams.vsite.b);
break;
case F_VSITE3OUT:
case F_VSITE4FD:
case F_VSITE4FDN:
- fprintf(fp, "a=%15.8e, b=%15.8e, c=%15.8e\n", iparams->vsite.a, iparams->vsite.b,
- iparams->vsite.c);
+ writer->writeLineFormatted("a=%15.8e, b=%15.8e, c=%15.8e", iparams.vsite.a,
+ iparams.vsite.b, iparams.vsite.c);
break;
case F_VSITEN:
- fprintf(fp, "n=%2d, a=%15.8e\n", iparams->vsiten.n, iparams->vsiten.a);
+ writer->writeLineFormatted("n=%2d, a=%15.8e", iparams.vsiten.n, iparams.vsiten.a);
break;
case F_GB12_NOLONGERUSED:
case F_GB13_NOLONGERUSED:
// nothing to print.
break;
case F_CMAP:
- fprintf(fp, "cmapA=%1d, cmapB=%1d\n", iparams->cmap.cmapA, iparams->cmap.cmapB);
+ writer->writeLineFormatted("cmapA=%1d, cmapB=%1d", iparams.cmap.cmapA, iparams.cmap.cmapB);
break;
- case F_RESTRANGLES: pr_harm(fp, iparams, "ktheta", "costheta0"); break;
+ case F_RESTRANGLES: printHarmonicInteraction(writer, iparams, "ktheta", "costheta0"); break;
case F_RESTRDIHS:
- fprintf(fp, "phiA=%15.8e, cpA=%15.8e", iparams->pdihs.phiA, iparams->pdihs.cpA);
+ writer->writeLineFormatted("phiA=%15.8e, cpA=%15.8e", iparams.pdihs.phiA, iparams.pdihs.cpA);
break;
case F_CBTDIHS:
- fprintf(fp, "kphi=%15.8e", iparams->cbtdihs.cbtcA[0]);
+ writer->writeLineFormatted("kphi=%15.8e", iparams.cbtdihs.cbtcA[0]);
for (int i = 1; i < NR_CBTDIHS; i++)
{
- fprintf(fp, ", cbtcA[%d]=%15.8e", i - 1, iparams->cbtdihs.cbtcA[i]);
+ writer->writeStringFormatted(", cbtcA[%d]=%15.8e", i - 1, iparams.cbtdihs.cbtcA[i]);
}
- fprintf(fp, "\n");
+ writer->ensureLineBreak();
break;
default:
gmx_fatal(FARGS, "unknown function type %d (%s) in %s line %d", ftype,
if (bShowParameters)
{
fprintf(fp, " ");
- pr_iparams(fp, ftype, &iparams[type]);
+ pr_iparams(fp, ftype, iparams[type]);
}
fprintf(fp, "\n");
i += 1 + interaction_function[ftype].nratoms;
pr_indent(fp, indent + INDENT);
fprintf(fp, "functype[%d]=%s, ", bShowNumbers ? i : -1,
interaction_function[idef->functype[i]].name);
- pr_iparams(fp, idef->functype[i], &idef->iparams[i]);
+ pr_iparams(fp, idef->functype[i], idef->iparams[i]);
}
pr_real(fp, indent, "fudgeQQ", idef->fudgeQQ);
* The state of the sorting of il, values are provided above.
*/
-void pr_iparams(FILE* fp, t_functype ftype, const t_iparams* iparams);
+namespace gmx
+{
+class TextWriter;
+} // namespace gmx
+
+void printInteractionParameters(gmx::TextWriter* writer, t_functype ftype, const t_iparams& iparams);
+void pr_iparams(FILE* fp, t_functype ftype, const t_iparams& iparams);
void pr_ilist(FILE* fp,
int indent,
const char* title,
if (bDiff)
{
fprintf(fp, "%s1: ", s);
- pr_iparams(fp, ft, &ip1);
+ pr_iparams(fp, ft, ip1);
fprintf(fp, "%s2: ", s);
- pr_iparams(fp, ft, &ip2);
+ pr_iparams(fp, ft, ip2);
}
}
if (bDiff)
{
fprintf(fp, "%s: ", s);
- pr_iparams(fp, ft, &ip1);
+ pr_iparams(fp, ft, ip1);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff