Merge release-2019 branch into master
authorSzilárd Páll <pall.szilard@gmail.com>
Tue, 18 Dec 2018 16:30:28 +0000 (17:30 +0100)
committerSzilárd Páll <pall.szilard@gmail.com>
Tue, 18 Dec 2018 22:48:53 +0000 (23:48 +0100)
Change-Id: Ie08781a0f40b2c92321e605db91d047a4451d2b2

60 files changed:
COPYING
cmake/gmxVersionInfo.cmake
docs/CMakeLists.txt
docs/release-notes/2019/major/deprecated-functionality.rst
docs/release-notes/2020/major/removed-functionality.rst [new file with mode: 0644]
docs/release-notes/bugs-fixed.rst [new file with mode: 0644]
docs/release-notes/deprecated-functionality.rst [new file with mode: 0644]
docs/release-notes/features.rst [new file with mode: 0644]
docs/release-notes/highlights.rst [new file with mode: 0644]
docs/release-notes/index.rst
docs/release-notes/miscellaneous.rst [new file with mode: 0644]
docs/release-notes/performance.rst [new file with mode: 0644]
docs/release-notes/portability.rst [new file with mode: 0644]
docs/release-notes/removed-functionality.rst [new file with mode: 0644]
docs/release-notes/tools.rst [new file with mode: 0644]
src/api/CMakeLists.txt
src/api/cpp/context.cpp
src/gromacs/CMakeLists.txt
src/gromacs/fileio/tngio.cpp
src/gromacs/gmxana/gmx_ana.h
src/gromacs/gmxana/gmx_anadock.cpp [deleted file]
src/gromacs/gmxana/gmx_dyndom.cpp [deleted file]
src/gromacs/gmxana/gmx_morph.cpp [deleted file]
src/gromacs/gmxana/gmx_msd.cpp
src/gromacs/gmxpreprocess/grompp.cpp
src/gromacs/gmxpreprocess/readir.cpp
src/gromacs/gmxpreprocess/tests/solvate.cpp
src/gromacs/listed-forces/orires.cpp
src/gromacs/mdlib/forcerec.cpp
src/gromacs/mdlib/mdatoms.cpp
src/gromacs/mdlib/mdoutf.cpp
src/gromacs/mdlib/ns.cpp
src/gromacs/mdlib/qmmm.cpp
src/gromacs/mdlib/tgroup.cpp
src/gromacs/mdrun/md.cpp
src/gromacs/mdspan/CMakeLists.txt [new file with mode: 0644]
src/gromacs/mdspan/README [new file with mode: 0644]
src/gromacs/mdspan/extents.h [new file with mode: 0644]
src/gromacs/mdspan/tests/CMakeLists.txt [new file with mode: 0644]
src/gromacs/mdspan/tests/extents.cpp [new file with mode: 0644]
src/gromacs/mimic/MimicCommunicator.cpp
src/gromacs/pulling/pull.cpp
src/gromacs/tools/check.cpp
src/gromacs/tools/dump.cpp
src/gromacs/topology/topology.cpp
src/gromacs/topology/topology.h
src/gromacs/utility/coolstuff.cpp
src/gromacs/utility/pleasecite.cpp
src/programs/legacymodules.cpp
src/programs/mdrun/tests/CMakeLists.txt
src/programs/mdrun/tests/moduletest.cpp
src/programs/mdrun/tests/moduletest.h
src/programs/mdrun/tests/normalmodes.cpp [new file with mode: 0644]
src/programs/mdrun/tests/refdata/NormalModesWorks_NormalModesTest_WithinTolerances_0.xml [new file with mode: 0644]
src/programs/mdrun/tests/simulationdatabase.cpp
src/testutils/cmdlinetest.cpp
src/testutils/cmdlinetest.h
src/testutils/simulationdatabase/scaled-water.g96 [new file with mode: 0644]
src/testutils/simulationdatabase/scaled-water.ndx [new file with mode: 0644]
src/testutils/simulationdatabase/scaled-water.top [new file with mode: 0644]

diff --git a/COPYING b/COPYING
index 1d94bc08a950878a1cd095de1770a48aa4aa1103..7203cc4edfa9082eb165e2bb5e197f65675fe41d 100644 (file)
--- a/COPYING
+++ b/COPYING
@@ -25,6 +25,7 @@ This file contains the licenses for the following bodies of code:
 14. lmfit
 15. clFFT
 16. Guidelines Support Library (GSL)
+17. P0009 reference implementation (Sandia Corporation)
 
 Our chosen method for packaging distributions (CPack) only permits a
 package to have a single license file, so we are unfortunately forced
@@ -1322,3 +1323,45 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 --
+
+17. P0009 reference implementation (Sandia Corporation)
+=======================================================
+
+   Files: src/gromacs/mdspan/*
+
+                         Kokkos v. 2.0
+              Copyright (2014) Sandia Corporation
+
+Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+the U.S. Government retains certain rights in this software.
+
+Kokkos is licensed under 3-clause BSD terms of use:
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the Corporation nor the names of the
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+Questions? Contact Christian R. Trott (crtrott@sandia.gov)
index f6c8c18e1c92bbae047ac6e825bb88978b89cac7..00c28c65a07f50895a41fb3b903a255deb55362b 100644 (file)
@@ -57,6 +57,7 @@
 #         GROMACS     2016   2
 #         GROMACS     2018   3
 #         GROMACS     2019   4
+#         GROMACS     2020   5
 #   LIBRARY_SOVERSION_MINOR so minor version for the built libraries.
 #       Should be increased for each release that changes only the implementation.
 #       In GROMACS, the typical policy is to increase it for each patch version
 
 # The GROMACS convention is that these are the version number of the next
 # release that is going to be made from this branch.
-set(GMX_VERSION_MAJOR 2019)
+set(GMX_VERSION_MAJOR 2020)
 set(GMX_VERSION_PATCH 0)
 # The suffix, on the other hand, is used mainly for betas and release
 # candidates, where it signifies the most recent such release from
 # this branch; it will be empty before the first such release, as well
 # as after the final release is out.
-set(GMX_VERSION_SUFFIX "-rc1")
+set(GMX_VERSION_SUFFIX "")
 
 # Conventionally with libtool, any ABI change must change the major
 # version number, the minor version number should change if it's just
@@ -212,7 +213,7 @@ set(GMX_VERSION_SUFFIX "-rc1")
 # here. The important thing is to minimize the chance of third-party
 # code being able to dynamically link with a version of libgromacs
 # that might not work.
-set(LIBRARY_SOVERSION_MAJOR 4)
+set(LIBRARY_SOVERSION_MAJOR 5)
 set(LIBRARY_SOVERSION_MINOR 0)
 set(LIBRARY_VERSION ${LIBRARY_SOVERSION_MAJOR}.${LIBRARY_SOVERSION_MINOR}.0)
 
@@ -235,13 +236,13 @@ if (NOT SOURCE_IS_SOURCE_DISTRIBUTION AND
 endif()
 
 set(REGRESSIONTEST_VERSION "${GMX_VERSION_STRING}")
-set(REGRESSIONTEST_BRANCH "refs/heads/release-2019")
+set(REGRESSIONTEST_BRANCH "refs/heads/master")
 # Run the regressiontests packaging job with the correct pakage
 # version string, and the release box checked, in order to have it
 # 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 "69da7f6a5bfcdb38b202eb9f9df69d01" CACHE INTERNAL "MD5 sum of the regressiontests tarball for this GROMACS version")
+set(REGRESSIONTEST_MD5SUM "3d06d41e07f523d70ae575b9ad75c670" CACHE INTERNAL "MD5 sum of the regressiontests tarball for this GROMACS version")
 
 math(EXPR GMX_VERSION_NUMERIC
      "${GMX_VERSION_MAJOR}*10000 + ${GMX_VERSION_PATCH}")
index 36e184bc76c48b870a2de2cd7cc49b61bece5a7d..5c4fc632c7c67abc8838564c64c7f4b6b86324b2 100644 (file)
@@ -367,6 +367,15 @@ if (SPHINX_FOUND)
         how-to/visualize.rst
         install-guide/index.rst
         release-notes/index.rst
+        release-notes/highlights.rst
+        release-notes/features.rst
+        release-notes/performance.rst
+        release-notes/tools.rst
+        release-notes/bugs-fixed.rst
+        release-notes/removed-functionality.rst
+        release-notes/deprecated-functionality.rst
+        release-notes/portability.rst
+        release-notes/miscellaneous.rst
         release-notes/2019/major/highlights.rst
         release-notes/2019/major/features.rst
         release-notes/2019/major/performance.rst
index bcb0a362ce6c9479ea8ac8c61c1d83c602f92bfb..182ed840477078998bbc688640fdee7da8a6a466 100644 (file)
@@ -1,5 +1,3 @@
-.. _anticipated-changes:
-
 Changes anticipated to GROMACS 2019 functionality
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/docs/release-notes/2020/major/removed-functionality.rst b/docs/release-notes/2020/major/removed-functionality.rst
new file mode 100644 (file)
index 0000000..72f7a86
--- /dev/null
@@ -0,0 +1,17 @@
+Removed functionality
+^^^^^^^^^^^^^^^^^^^^^
+
+gmx anadock
+"""""""""""
+The gmx anadock tool was removed since it does not belong in gromacs
+(it analyzes AutoDock outputs).
+
+gmx dyndom
+""""""""""
+The gmx dyndom tool was removed since it does not belong in gromacs
+(it analyzes DynDom outputs).
+
+gmx morph
+"""""""""
+The gmx morph tool was removed since it yields non-physical structures
+that can easily be done by a script.
diff --git a/docs/release-notes/bugs-fixed.rst b/docs/release-notes/bugs-fixed.rst
new file mode 100644 (file)
index 0000000..81b559a
--- /dev/null
@@ -0,0 +1,2 @@
+Bugs fixed
+^^^^^^^^^^
diff --git a/docs/release-notes/deprecated-functionality.rst b/docs/release-notes/deprecated-functionality.rst
new file mode 100644 (file)
index 0000000..90218f4
--- /dev/null
@@ -0,0 +1,52 @@
+.. _anticipated-changes:
+
+Changes anticipated to |Gromacs| NEXT functionality
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Functionality deprecated in |Gromacs| NEXT
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Functionality deprecated in |Gromacs| 2019
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Generation of virtual sites to replace aromatic rings in standard residues
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+These are thought to produce artefacts under some circumstances
+(unpublished results), were never well tested, are not widely used,
+and we need to simplify pdb2gmx.
+
+``gmx mdrun -gcom``
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+This feature sometimes overrides the effects of various .mdp settings
+in a way that is difficult to understand and report. A user who wants
+to do communication between PP ranks less often should choose their
+``nst*`` mdp options accordingly.
+
+Benchmarking options only available with ``gmx benchmark``
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+Options such as ``-confout``, ``-resethway``, ``-resetstep`` are not
+intended for use by regular mdrun users, so making them only available
+with a dedicated tool is more clear. Also, this permits us to customize
+defaults for e.g. writing files at the end of a simulation part in ways
+that suit the respective mdrun and benchmark use cases, so ``-confout``
+will no longer be required.
+
+``gmx mdrun -nsteps``
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+The number of simulation steps described by the .tpr file can be
+changed with ``gmx convert-tpr``, or altered in .mdp file before the
+call to ``gmx grompp``. The convenience of this mdrun option was
+outweighted by the doubtful quality of its implementation, no clear
+record in the log file, and lack of maintenance.
+
+Functionality deprecated before |Gromacs| 2019
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+This functionality has been declared as deprecated in previous versions
+of |Gromacs|, but has not yet been removed.
+
+The group cutoff scheme
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+All remaining aspects of the group cutoff scheme will be removed, once
+a few remaining features (e.g. tabulated interactions, energy-group
+exclusions, and vacuum simulations) are available with the Verlet
+scheme. Deprecated in GROMACS 5.1
diff --git a/docs/release-notes/features.rst b/docs/release-notes/features.rst
new file mode 100644 (file)
index 0000000..5aa11eb
--- /dev/null
@@ -0,0 +1,2 @@
+New and improved features
+^^^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/docs/release-notes/highlights.rst b/docs/release-notes/highlights.rst
new file mode 100644 (file)
index 0000000..7026520
--- /dev/null
@@ -0,0 +1,14 @@
+Highlights
+^^^^^^^^^^
+
+|Gromacs| NEXT was released on INSERT DATE HERE. Patch releases may
+have been made since then, please use the updated versions!  Here are
+some highlights of what you can expect, along with more detail in the
+links below!
+
+As always, we've got several useful performance improvements, with or
+without GPUs, all enabled and automated by default. We are extremely
+interested in your feedback on how well this worked on your
+simulations and hardware. They are:
+
+* Cool quotes autogenerator 
index 014de1170b01b9638e5f3f74e62c6396824f216f..efa62f1b555acbc1d64ea676c3d93defe2cd0801 100644 (file)
@@ -8,18 +8,35 @@ releases of |Gromacs|. Major releases contain changes to the
 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 2019
-series and the 2018 series. In the latter, only highly conservative
+Two versions of |Gromacs| are under active maintenance, the NEXT
+series and the 2019 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 2018 ends, but we keep 2018 in the name so users understand how
+year 2019 ends, but we keep 2018 in the name so users understand how
 up to date their version is. Such fixes will also be incorporated into
-the 2019 release series, as appropriate. Around the time the 2020
-release is made, the 2018 series will no longer be maintained.
+the NEXT release series, as appropriate. Around the time the NEXT+1
+release is made, the 2019 series will no longer be maintained.
 
 Where issue numbers are reported in these release notes, more details
 can be found at https://redmine.gromacs.org at that issue number.
 
+|Gromacs| NEXT series
+---------------------
+
+.. toctree::
+   :maxdepth: 1
+
+   highlights
+   features
+   performance
+   tools
+   bugs-fixed
+   deprecated-functionality
+   removed-functionality
+   portability
+   miscellaneous
+
+
 |Gromacs| 2019 series
 ---------------------
 
diff --git a/docs/release-notes/miscellaneous.rst b/docs/release-notes/miscellaneous.rst
new file mode 100644 (file)
index 0000000..27789e4
--- /dev/null
@@ -0,0 +1,2 @@
+Miscellaneous
+^^^^^^^^^^^^^
diff --git a/docs/release-notes/performance.rst b/docs/release-notes/performance.rst
new file mode 100644 (file)
index 0000000..c9d7e7c
--- /dev/null
@@ -0,0 +1,2 @@
+Performance improvements
+^^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/docs/release-notes/portability.rst b/docs/release-notes/portability.rst
new file mode 100644 (file)
index 0000000..bc64842
--- /dev/null
@@ -0,0 +1,2 @@
+Portability
+^^^^^^^^^^^
diff --git a/docs/release-notes/removed-functionality.rst b/docs/release-notes/removed-functionality.rst
new file mode 100644 (file)
index 0000000..8d99815
--- /dev/null
@@ -0,0 +1,2 @@
+Removed functionality
+^^^^^^^^^^^^^^^^^^^^^
diff --git a/docs/release-notes/tools.rst b/docs/release-notes/tools.rst
new file mode 100644 (file)
index 0000000..fb2512b
--- /dev/null
@@ -0,0 +1,2 @@
+Improvements to |Gromacs| tools
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
index 67d428cbb10b3758fb168c19bcb31e157635a5cd..8f497476661bcde34c71c9e7fec19c6593d15863 100644 (file)
@@ -46,7 +46,7 @@
 # release string. For roadmap details, see https://redmine.gromacs.org/issues/2585
 set(GMXAPI_MAJOR 0)
 set(GMXAPI_MINOR 0)
-set(GMXAPI_PATCH 7)
+set(GMXAPI_PATCH 8)
 set(GMXAPI_RELEASE ${GMXAPI_MAJOR}.${GMXAPI_MINOR}.${GMXAPI_PATCH})
 
 ###############################
index 8b455db34d21074c2435b7d384df0958fe33a274..bea12dbf0f65ea0b57fcee99840275f772d3b693 100644 (file)
@@ -120,23 +120,22 @@ std::shared_ptr<Session> ContextImpl::launch(const Workflow &work)
          * a microstate for gmxapi interfaces.
          */
 
-        // Note: these output options normalize the file names, but not their
-        // paths. gmxapi 0.0.7 changes working directory for each session, so the
-        // relative paths are appropriate, but a near-future version will avoid
-        // changing directories once the process starts and manage file paths explicitly.
-        using gmxapi::c_majorVersion;
-        using gmxapi::c_minorVersion;
-        using gmxapi::c_patchVersion;
-        static_assert(!(c_majorVersion != 0 || c_minorVersion != 0 || c_patchVersion > 7),
-                      "Developer notice: check assumptions about working directory and relative file paths for this "
-                      "software version.");
-
         // Set input TPR name
         mdArgs_.emplace_back("-s");
         mdArgs_.emplace_back(filename);
+
         // Set checkpoint file name
         mdArgs_.emplace_back("-cpi");
         mdArgs_.emplace_back("state.cpt");
+        /* Note: we normalize the checkpoint file name, but not its full path.
+         * Through version 0.0.8, gmxapi clients change working directory
+         * for each session, so relative path(s) below are appropriate.
+         * A future gmxapi version should avoid changing directories once the
+         * process starts and instead manage files (paths) in an absolute and
+         * immutable way, with abstraction provided through the Context chain-of-responsibility.
+         * TODO: API abstractions for initializing simulations that may be new or partially complete.
+         * Reference gmxapi milestone 13 at https://redmine.gromacs.org/issues/2585
+         */
 
         // Create a mock argv. Note that argv[0] is expected to hold the program name.
         const int  offset = 1;
index 3430faf749065e9a97f7e985d0bf94ce87c6ea62..42f54f1fb770ac25e9c4cc31880c836c8349ad59 100644 (file)
@@ -107,6 +107,7 @@ add_subdirectory(linearalgebra)
 add_subdirectory(math)
 add_subdirectory(mdrun)
 add_subdirectory(mdrunutility)
+add_subdirectory(mdspan)
 add_subdirectory(mdtypes)
 add_subdirectory(onlinehelp)
 add_subdirectory(options)
index b69f38b7fb797536f2e54c2ffc2cd2db0b40b3ed..ee1b2d2e7aa74c0730370b9696b53f512ab33e13 100644 (file)
@@ -588,7 +588,7 @@ static gmx_bool all_atoms_selected(const gmx_mtop_t *mtop,
         {
             for (int atomIndex = 0; atomIndex < atoms.nr; atomIndex++, i++)
             {
-                if (getGroupType(&mtop->groups, gtype, i) != 0)
+                if (getGroupType(mtop->groups, gtype, i) != 0)
                 {
                     return FALSE;
                 }
@@ -664,7 +664,7 @@ static void add_selection_groups(gmx_tng_trajectory_t  gmx_tng,
                 char *res_name;
                 int   res_id;
 
-                if (getGroupType(&mtop->groups, egcCompressedX, i) != 0)
+                if (getGroupType(mtop->groups, egcCompressedX, i) != 0)
                 {
                     continue;
                 }
@@ -707,8 +707,8 @@ static void add_selection_groups(gmx_tng_trajectory_t  gmx_tng,
                             int atom1, atom2;
                             atom1 = ilist.iatoms[l] + atom_offset;
                             atom2 = ilist.iatoms[l + 1] + atom_offset;
-                            if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0 &&
-                                getGroupType(&mtop->groups, egcCompressedX, atom2) == 0)
+                            if (getGroupType(mtop->groups, egcCompressedX, atom1) == 0 &&
+                                getGroupType(mtop->groups, egcCompressedX, atom2) == 0)
                             {
                                 tng_molecule_bond_add(tng, mol, ilist.iatoms[l],
                                                       ilist.iatoms[l + 1], &tngBond);
@@ -724,14 +724,14 @@ static void add_selection_groups(gmx_tng_trajectory_t  gmx_tng,
                     atom1 = ilist.iatoms[l] + atom_offset;
                     atom2 = ilist.iatoms[l + 1] + atom_offset;
                     atom3 = ilist.iatoms[l + 2] + atom_offset;
-                    if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0)
+                    if (getGroupType(mtop->groups, egcCompressedX, atom1) == 0)
                     {
-                        if (getGroupType(&mtop->groups, egcCompressedX, atom2) == 0)
+                        if (getGroupType(mtop->groups, egcCompressedX, atom2) == 0)
                         {
                             tng_molecule_bond_add(tng, mol, atom1,
                                                   atom2, &tngBond);
                         }
-                        if (getGroupType(&mtop->groups, egcCompressedX, atom3) == 0)
+                        if (getGroupType(mtop->groups, egcCompressedX, atom3) == 0)
                         {
                             tng_molecule_bond_add(tng, mol, atom1,
                                                   atom3, &tngBond);
index 0517719bcc2da62174af8f29f87ad73ffedea207..0924f04f19940a4167b967195aec831925f51299 100644 (file)
@@ -38,9 +38,6 @@
 #ifndef _gmx_ana_h
 #define _gmx_ana_h
 
-int
-gmx_anadock(int argc, char *argv[]);
-
 int
 gmx_analyze(int argc, char *argv[]);
 
@@ -101,9 +98,6 @@ gmx_dos(int argc, char *argv[]);
 int
 gmx_dyecoupl(int argc, char *argv[]);
 
-int
-gmx_dyndom(int argc, char *argv[]);
-
 int
 gmx_editconf(int argc, char *argv[]);
 
@@ -161,9 +155,6 @@ gmx_mk_angndx(int argc, char *argv[]);
 int
 gmx_msd(int argc, char *argv[]);
 
-int
-gmx_morph(int argc, char *argv[]);
-
 int
 gmx_nmeig(int argc, char *argv[]);
 
diff --git a/src/gromacs/gmxana/gmx_anadock.cpp b/src/gromacs/gmxana/gmx_anadock.cpp
deleted file mode 100644 (file)
index b8ca3f8..0000000
+++ /dev/null
@@ -1,391 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2013,2014,2015,2017,2018, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include <cmath>
-#include <cstdlib>
-#include <cstring>
-
-#include <string>
-
-#include "gromacs/commandline/pargs.h"
-#include "gromacs/fileio/confio.h"
-#include "gromacs/fileio/xvgr.h"
-#include "gromacs/gmxana/gmx_ana.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/statistics/statistics.h"
-#include "gromacs/topology/topology.h"
-#include "gromacs/utility/arraysize.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/futil.h"
-#include "gromacs/utility/path.h"
-#include "gromacs/utility/pleasecite.h"
-#include "gromacs/utility/smalloc.h"
-
-static const char *etitles[] = { "E-docked", "Free Energy" };
-
-typedef struct {
-    real     edocked, efree;
-    int      index, cluster_id;
-    t_atoms  atoms;
-    rvec    *x;
-    int      ePBC;
-    matrix   box;
-} t_pdbfile;
-
-static t_pdbfile *read_pdbf(const char *fn)
-{
-    t_pdbfile *pdbf;
-    double     e;
-    FILE      *fp;
-
-    snew(pdbf, 1);
-    t_topology top;
-    read_tps_conf(fn, &top, &pdbf->ePBC, &pdbf->x, nullptr, pdbf->box, FALSE);
-    pdbf->atoms = top.atoms;
-    fp          = gmx_ffopen(fn, "r");
-    char       buf[256], *ptr;
-    while ((ptr = fgets2(buf, 255, fp)) != nullptr)
-    {
-        if (std::strstr(buf, "Intermolecular") != nullptr)
-        {
-            ptr = std::strchr(buf, '=');
-            sscanf(ptr+1, "%lf", &e);
-            pdbf->edocked = e;
-        }
-        else if (std::strstr(buf, "Estimated Free") != nullptr)
-        {
-            ptr = std::strchr(buf, '=');
-            sscanf(ptr+1, "%lf", &e);
-            pdbf->efree = e;
-        }
-    }
-    gmx_ffclose(fp);
-
-    return pdbf;
-}
-
-static t_pdbfile **read_em_all(const std::string &fn, int *npdbf)
-{
-    t_pdbfile **pdbf = nullptr;
-    int         i, maxpdbf;
-    gmx_bool    bExist;
-
-    maxpdbf                 = 100;
-    snew(pdbf, maxpdbf);
-    i = 0;
-    do
-    {
-        auto newFilename = gmx::Path::concatenateBeforeExtension(fn, gmx::formatString("_%d", i+1));
-        if ((bExist = gmx_fexist(newFilename)))
-        {
-            pdbf[i]        = read_pdbf(newFilename.c_str());
-            pdbf[i]->index = i+1;
-            i++;
-            if (i >= maxpdbf)
-            {
-                maxpdbf += 100;
-                srenew(pdbf, maxpdbf);
-            }
-        }
-    }
-    while (bExist);
-
-    *npdbf = i;
-
-    printf("Found %d pdb files\n", i);
-
-    return pdbf;
-}
-
-static gmx_bool bFreeSort = FALSE;
-
-static int pdbf_comp(const void *a, const void *b)
-{
-    const t_pdbfile *pa, *pb;
-    real             x;
-    int              dc;
-
-    pa = *static_cast<t_pdbfile *const*>(a);
-    pb = *static_cast<t_pdbfile *const*>(b);
-
-    dc = pa->cluster_id - pb->cluster_id;
-
-    if (dc == 0)
-    {
-        if (bFreeSort)
-        {
-            x = pa->efree   - pb->efree;
-        }
-        else
-        {
-            x = pa->edocked - pb->edocked;
-        }
-
-        if (x < 0)
-        {
-            return -1;
-        }
-        else if (x > 0)
-        {
-            return 1;
-        }
-        else
-        {
-            return 0;
-        }
-    }
-    else
-    {
-        return dc;
-    }
-}
-
-static void analyse_em_all(int npdb, t_pdbfile *pdbf[], const char *edocked,
-                           const char *efree, const gmx_output_env_t *oenv)
-{
-    FILE *fp;
-    int   i;
-
-    for (int freeSort = 0; freeSort < 2; freeSort++)
-    {
-        bFreeSort = (freeSort == 1);
-        qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
-        fp = xvgropen(bFreeSort ? efree : edocked,
-                      etitles[bFreeSort], "()", "E (kJ/mol)", oenv);
-        for (i = 0; (i < npdb); i++)
-        {
-            fprintf(fp, "%12lf\n", bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
-        }
-        xvgrclose(fp);
-    }
-}
-
-static void clust_stat(FILE *fp, int start, int end, t_pdbfile *pdbf[])
-{
-    int         i;
-    gmx_stats_t ed, ef;
-    real        aver, sigma;
-
-    ed = gmx_stats_init();
-    ef = gmx_stats_init();
-    for (i = start; (i < end); i++)
-    {
-        gmx_stats_add_point(ed, i-start, pdbf[i]->edocked, 0, 0);
-        gmx_stats_add_point(ef, i-start, pdbf[i]->efree, 0, 0);
-    }
-    gmx_stats_get_ase(ed, &aver, &sigma, nullptr);
-    fprintf(fp, "  <%12s> = %8.3f (+/- %6.3f)\n", etitles[FALSE], aver, sigma);
-    gmx_stats_get_ase(ef, &aver, &sigma, nullptr);
-    fprintf(fp, "  <%12s> = %8.3f (+/- %6.3f)\n", etitles[TRUE], aver, sigma);
-    gmx_stats_free(ed);
-    gmx_stats_free(ef);
-}
-
-static real rmsd_dist(t_pdbfile *pa, t_pdbfile *pb, gmx_bool bRMSD)
-{
-    int  i;
-    real rmsd;
-    rvec acm, bcm, dx;
-
-    if (bRMSD)
-    {
-        rmsd = 0;
-        for (i = 0; (i < pa->atoms.nr); i++)
-        {
-            rvec_sub(pa->x[i], pb->x[i], dx);
-            rmsd += iprod(dx, dx);
-        }
-        rmsd = std::sqrt(rmsd/pa->atoms.nr);
-    }
-    else
-    {
-        clear_rvec(acm);
-        clear_rvec(bcm);
-        for (i = 0; (i < pa->atoms.nr); i++)
-        {
-            rvec_inc(acm, pa->x[i]);
-            rvec_inc(bcm, pb->x[i]);
-        }
-        rvec_sub(acm, bcm, dx);
-        for (i = 0; (i < DIM); i++)
-        {
-            dx[i] /= pa->atoms.nr;
-        }
-        rmsd = norm(dx);
-    }
-    return rmsd;
-}
-
-static void line(FILE *fp)
-{
-    fprintf(fp, "                   - - - - - - - - - -\n");
-}
-
-static void cluster_em_all(FILE *fp, int npdb, t_pdbfile *pdbf[],
-                           gmx_bool bFree, gmx_bool bRMSD, real cutoff)
-{
-    int   i, j, k;
-    int  *cndx, ncluster;
-    real  rmsd;
-
-    bFreeSort = bFree;
-    qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
-
-    fprintf(fp, "Statistics over all structures:\n");
-    clust_stat(fp, 0, npdb, pdbf);
-    line(fp);
-
-    /* Index to first structure in a cluster */
-    snew(cndx, npdb);
-    ncluster = 1;
-
-    for (i = 1; (i < npdb); i++)
-    {
-        for (j = 0; (j < ncluster); j++)
-        {
-            rmsd = rmsd_dist(pdbf[cndx[j]], pdbf[i], bRMSD);
-            if (rmsd <= cutoff)
-            {
-                /* Structure i is in cluster j */
-                pdbf[i]->cluster_id = pdbf[cndx[j]]->cluster_id;
-                break;
-            }
-        }
-        if (j == ncluster)
-        {
-            /* New cluster! Cool! */
-            cndx[ncluster]      = i;
-            pdbf[i]->cluster_id = ncluster;
-            ncluster++;
-        }
-    }
-    cndx[ncluster] = npdb;
-    qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
-
-    j         = 0;
-    cndx[j++] = 0;
-    for (i = 1; (i < npdb); i++)
-    {
-        if (pdbf[i]->cluster_id != pdbf[i-1]->cluster_id)
-        {
-            cndx[j++] = i;
-        }
-    }
-    cndx[j] = npdb;
-    if (j != ncluster)
-    {
-        gmx_fatal(FARGS, "Consistency error: j = %d, ncluster = %d", j, ncluster);
-    }
-
-    fprintf(fp, "I found %d clusters based on %s and %s with a %.3f nm cut-off\n",
-            ncluster, etitles[bFree], bRMSD ? "RMSD" : "distance", cutoff);
-    line(fp);
-    for (j = 0; (j < ncluster); j++)
-    {
-        fprintf(fp, "Cluster: %3d  %s: %10.5f kJ/mol %3d elements\n",
-                j, etitles[bFree],
-                bFree ? pdbf[cndx[j]]->efree : pdbf[cndx[j]]->edocked,
-                cndx[j+1]-cndx[j]);
-        clust_stat(fp, cndx[j], cndx[j+1], pdbf);
-        for (k = cndx[j]; (k < cndx[j+1]); k++)
-        {
-            fprintf(fp, "  %3d", pdbf[k]->index);
-        }
-        fprintf(fp, "\n");
-        line(fp);
-    }
-    sfree(cndx);
-}
-
-int gmx_anadock(int argc, char *argv[])
-{
-    const char       *desc[] = {
-        "[THISMODULE] analyses the results of an Autodock run and clusters the",
-        "structures together, based on distance or RMSD. The docked energy",
-        "and free energy estimates are analysed, and for each cluster the",
-        "energy statistics are printed.[PAR]",
-        "An alternative approach to this is to cluster the structures first",
-        "using [gmx-cluster] and then sort the clusters on either lowest",
-        "energy or average energy."
-    };
-    t_filenm          fnm[] = {
-        { efPDB, "-f", nullptr,       ffREAD  },
-        { efXVG, "-od", "edocked", ffWRITE },
-        { efXVG, "-of", "efree",   ffWRITE },
-        { efLOG, "-g",  "anadock", ffWRITE }
-    };
-    gmx_output_env_t *oenv;
-#define NFILE asize(fnm)
-    static gmx_bool   bFree  = FALSE, bRMS = TRUE;
-    static real       cutoff = 0.2;
-    t_pargs           pa[]   = {
-        { "-free",   FALSE, etBOOL, {&bFree},
-          "Use Free energy estimate from autodock for sorting the classes" },
-        { "-rms",    FALSE, etBOOL, {&bRMS},
-          "Cluster on RMS or distance" },
-        { "-cutoff", FALSE, etREAL, {&cutoff},
-          "Maximum RMSD/distance for belonging to the same cluster" }
-    };
-#define NPA asize(pa)
-
-    FILE       *fp;
-    t_pdbfile **pdbf = nullptr;
-    int         npdbf;
-
-    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0,
-                           nullptr, &oenv))
-    {
-        return 0;
-    }
-
-    fp = gmx_ffopen(opt2fn("-g", NFILE, fnm), "w");
-    please_cite(stdout, "Hetenyi2002b");
-    please_cite(fp, "Hetenyi2002b");
-
-    pdbf = read_em_all(opt2fn("-f", NFILE, fnm), &npdbf);
-
-    analyse_em_all(npdbf, pdbf, opt2fn("-od", NFILE, fnm), opt2fn("-of", NFILE, fnm),
-                   oenv);
-
-    cluster_em_all(fp, npdbf, pdbf, bFree, bRMS, cutoff);
-
-    gmx_ffclose(fp);
-
-    return 0;
-}
diff --git a/src/gromacs/gmxana/gmx_dyndom.cpp b/src/gromacs/gmxana/gmx_dyndom.cpp
deleted file mode 100644 (file)
index 55219d8..0000000
+++ /dev/null
@@ -1,266 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * 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, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include <cmath>
-
-#include "gromacs/commandline/pargs.h"
-#include "gromacs/fileio/confio.h"
-#include "gromacs/fileio/trxio.h"
-#include "gromacs/gmxana/gmx_ana.h"
-#include "gromacs/math/3dtransforms.h"
-#include "gromacs/math/units.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/topology/atoms.h"
-#include "gromacs/topology/index.h"
-#include "gromacs/topology/topology.h"
-#include "gromacs/utility/arraysize.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-static void rot_conf(t_atoms *atoms, const rvec x[], const rvec v[], real trans, real angle,
-                     rvec head, rvec tail, int isize, const int index[],
-                     rvec xout[], rvec vout[])
-{
-    rvec     arrow, xcm;
-    real     theta, phi, arrow_len;
-    mat4     Rx, Ry, Rz, Rinvy, Rinvz, Mtot;
-    mat4     temp1, temp2, temp3;
-    vec4     xv;
-    int      i, ai;
-
-    rvec_sub(tail, head, arrow);
-    arrow_len = norm(arrow);
-    if (debug)
-    {
-        fprintf(debug, "Arrow vector:   %10.4f  %10.4f  %10.4f\n",
-                arrow[XX], arrow[YY], arrow[ZZ]);
-        fprintf(debug, "Effective translation %g nm\n", trans);
-    }
-    if (arrow_len == 0.0)
-    {
-        gmx_fatal(FARGS, "Arrow vector not given");
-    }
-
-    /* Copy all aoms to output */
-    for (i = 0; (i < atoms->nr); i++)
-    {
-        copy_rvec(x[i], xout[i]);
-        copy_rvec(v[i], vout[i]);
-    }
-
-    /* Compute center of mass and move atoms there */
-    clear_rvec(xcm);
-    for (i = 0; (i < isize); i++)
-    {
-        rvec_inc(xcm, x[index[i]]);
-    }
-    for (i = 0; (i < DIM); i++)
-    {
-        xcm[i] /= isize;
-    }
-    if (debug)
-    {
-        fprintf(debug, "Center of mass: %10.4f  %10.4f  %10.4f\n",
-                xcm[XX], xcm[YY], xcm[ZZ]);
-    }
-    for (i = 0; (i < isize); i++)
-    {
-        rvec_sub(x[index[i]], xcm, xout[index[i]]);
-    }
-
-    /* Compute theta and phi that describe the arrow */
-    theta = std::acos(arrow[ZZ]/arrow_len);
-    phi   = std::atan2(arrow[YY]/arrow_len, arrow[XX]/arrow_len);
-    if (debug)
-    {
-        fprintf(debug, "Phi = %.1f, Theta = %.1f\n", RAD2DEG*phi, RAD2DEG*theta);
-    }
-
-    /* Now the total rotation matrix: */
-    /* Rotate a couple of times */
-    gmx_mat4_init_rotation(ZZ, -phi, Rz);
-    gmx_mat4_init_rotation(YY, M_PI/2-theta, Ry);
-    gmx_mat4_init_rotation(XX, angle*DEG2RAD, Rx);
-    Rx[WW][XX] = trans;
-    gmx_mat4_init_rotation(YY, theta-M_PI/2, Rinvy);
-    gmx_mat4_init_rotation(ZZ, phi, Rinvz);
-
-    gmx_mat4_mmul(temp1, Ry, Rz);
-    gmx_mat4_mmul(temp2, Rinvy, Rx);
-    gmx_mat4_mmul(temp3, temp2, temp1);
-    gmx_mat4_mmul(Mtot, Rinvz, temp3);
-
-    if (debug)
-    {
-        gmx_mat4_print(debug, "Rz", Rz);
-        gmx_mat4_print(debug, "Ry", Ry);
-        gmx_mat4_print(debug, "Rx", Rx);
-        gmx_mat4_print(debug, "Rinvy", Rinvy);
-        gmx_mat4_print(debug, "Rinvz", Rinvz);
-        gmx_mat4_print(debug, "Mtot", Mtot);
-    }
-
-    for (i = 0; (i < isize); i++)
-    {
-        ai = index[i];
-        gmx_mat4_transform_point(Mtot, xout[ai], xv);
-        rvec_add(xv, xcm, xout[ai]);
-        gmx_mat4_transform_point(Mtot, v[ai], xv);
-        copy_rvec(xv, vout[ai]);
-    }
-}
-
-int gmx_dyndom(int argc, char *argv[])
-{
-    const char       *desc[] = {
-        "[THISMODULE] reads a [REF].pdb[ref] file output from DynDom",
-        "(http://www.cmp.uea.ac.uk/dyndom/).",
-        "It reads the coordinates, the coordinates of the rotation axis,",
-        "and an index file containing the domains.",
-        "Furthermore, it takes the first and last atom of the arrow file",
-        "as command line arguments (head and tail) and",
-        "finally it takes the translation vector (given in DynDom info file)",
-        "and the angle of rotation (also as command line arguments). If the angle",
-        "determined by DynDom is given, one should be able to recover the",
-        "second structure used for generating the DynDom output.",
-        "Because of limited numerical accuracy this should be verified by",
-        "computing an all-atom RMSD (using [gmx-confrms]) rather than by file",
-        "comparison (using diff).[PAR]",
-        "The purpose of this program is to interpolate and extrapolate the",
-        "rotation as found by DynDom. As a result unphysical structures with",
-        "long or short bonds, or overlapping atoms may be produced. Visual",
-        "inspection, and energy minimization may be necessary to",
-        "validate the structure."
-    };
-    static real       trans0 = 0;
-    static rvec       head   = { 0, 0, 0 };
-    static rvec       tail   = { 0, 0, 0 };
-    static real       angle0 = 0, angle1 = 0, maxangle = 0;
-    static int        label  = 0, nframes = 11;
-    t_pargs           pa[]   = {
-        { "-firstangle",    FALSE, etREAL, {&angle0},
-          "Angle of rotation about rotation vector" },
-        { "-lastangle",    FALSE, etREAL, {&angle1},
-          "Angle of rotation about rotation vector" },
-        { "-nframe",   FALSE, etINT,  {&nframes},
-          "Number of steps on the pathway" },
-        { "-maxangle", FALSE, etREAL, {&maxangle},
-          "DymDom dtermined angle of rotation about rotation vector" },
-        { "-trans",    FALSE, etREAL, {&trans0},
-          "Translation (Angstrom) along rotation vector (see DynDom info file)" },
-        { "-head",     FALSE, etRVEC, {head},
-          "First atom of the arrow vector" },
-        { "-tail",     FALSE, etRVEC, {tail},
-          "Last atom of the arrow vector" }
-    };
-    int               i, j, natoms, isize;
-    t_trxstatus      *status;
-    int              *index = nullptr, *index_all;
-    char             *grpname;
-    real              angle, trans;
-    rvec             *x, *v, *xout, *vout;
-    matrix            box;
-    gmx_output_env_t *oenv;
-
-    t_filenm          fnm[] = {
-        { efPDB, "-f", "dyndom",  ffREAD },
-        { efTRO, "-o", "rotated", ffWRITE },
-        { efNDX, "-n", "domains", ffREAD }
-    };
-#define NFILE asize(fnm)
-
-    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
-                           asize(desc), desc, 0, nullptr, &oenv))
-    {
-        return 0;
-    }
-
-    if (maxangle == 0)
-    {
-        gmx_fatal(FARGS, "maxangle not given");
-    }
-
-    t_topology *top;
-    snew(top, 1);
-    read_tps_conf(opt2fn("-f", NFILE, fnm), top, nullptr, &x, &v, box, FALSE);
-    t_atoms  &atoms = top->atoms;
-    if (atoms.pdbinfo == nullptr)
-    {
-        snew(atoms.pdbinfo, atoms.nr);
-    }
-    atoms.havePdbInfo = TRUE;
-    natoms            = atoms.nr;
-    snew(xout, natoms);
-    snew(vout, natoms);
-
-    printf("Select group to rotate:\n");
-    rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
-    printf("Going to rotate %s containing %d atoms\n", grpname, isize);
-
-    snew(index_all, atoms.nr);
-    for (i = 0; (i < atoms.nr); i++)
-    {
-        index_all[i] = i;
-    }
-
-    status = open_trx(opt2fn("-o", NFILE, fnm), "w");
-
-    label = 'A';
-    for (i = 0; (i < nframes); i++, label++)
-    {
-        angle = angle0 + (i*(angle1-angle0))/(nframes-1);
-        trans = trans0*0.1*angle/maxangle;
-        printf("Frame: %2d (label %c), angle: %8.3f deg., trans: %8.3f nm\n",
-               i, label, angle, trans);
-        rot_conf(&atoms, x, v, trans, angle, head, tail, isize, index, xout, vout);
-
-        if (label > 'Z')
-        {
-            label -= 26;
-        }
-        for (j = 0; (j < atoms.nr); j++)
-        {
-            atoms.resinfo[atoms.atom[j].resind].chainid = label;
-        }
-
-        write_trx(status, atoms.nr, index_all, &atoms, i, angle, box, xout, vout, nullptr);
-    }
-    close_trx(status);
-
-    return 0;
-}
diff --git a/src/gromacs/gmxana/gmx_morph.cpp b/src/gromacs/gmxana/gmx_morph.cpp
deleted file mode 100644 (file)
index 20618e0..0000000
+++ /dev/null
@@ -1,204 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "gromacs/commandline/pargs.h"
-#include "gromacs/commandline/viewit.h"
-#include "gromacs/fileio/confio.h"
-#include "gromacs/fileio/trxio.h"
-#include "gromacs/fileio/xvgr.h"
-#include "gromacs/gmxana/gmx_ana.h"
-#include "gromacs/math/do_fit.h"
-#include "gromacs/topology/atoms.h"
-#include "gromacs/topology/index.h"
-#include "gromacs/topology/topology.h"
-#include "gromacs/utility/arraysize.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-static real dointerp(int n, rvec x1[], rvec x2[], rvec xx[],
-                     int I, int N, real first, real last)
-{
-    int    i, j;
-    double fac, fac0, fac1;
-
-    fac  = first + (I*(last-first))/(N-1);
-    fac0 = 1-fac;
-    fac1 = fac;
-    for (i = 0; (i < n); i++)
-    {
-        for (j = 0; (j < DIM); j++)
-        {
-            xx[i][j] = fac0*x1[i][j] + fac1*x2[i][j];
-        }
-    }
-
-    return fac;
-}
-
-int gmx_morph(int argc, char *argv[])
-{
-    const char       *desc[] = {
-        "[THISMODULE] does a linear interpolation of conformations in order to",
-        "create intermediates. Of course these are completely unphysical, but",
-        "that you may try to justify yourself. Output is in the form of a ",
-        "generic trajectory. The number of intermediates can be controlled with",
-        "the [TT]-ninterm[tt] flag. The first and last flag correspond to the way of",
-        "interpolating: 0 corresponds to input structure 1 while",
-        "1 corresponds to input structure 2.",
-        "If you specify [TT]-first[tt] < 0 or [TT]-last[tt] > 1 extrapolation will be",
-        "on the path from input structure x[SUB]1[sub] to x[SUB]2[sub]. In general, the coordinates",
-        "of the intermediate x(i) out of N total intermediates correspond to:[PAR]",
-        "x(i) = x[SUB]1[sub] + (first+(i/(N-1))*(last-first))*(x[SUB]2[sub]-x[SUB]1[sub])[PAR]",
-        "Finally the RMSD with respect to both input structures can be computed",
-        "if explicitly selected ([TT]-or[tt] option). In that case, an index file may be",
-        "read to select the group from which the RMS is computed."
-    };
-    t_filenm          fnm[] = {
-        { efSTX, "-f1", "conf1",  ffREAD },
-        { efSTX, "-f2", "conf2",  ffREAD },
-        { efTRX, "-o",  "interm", ffWRITE },
-        { efXVG, "-or", "rms-interm", ffOPTWR },
-        { efNDX, "-n",  "index",  ffOPTRD }
-    };
-#define NFILE asize(fnm)
-    static  int       ninterm = 11;
-    static  real      first   = 0.0;
-    static  real      last    = 1.0;
-    static  gmx_bool  bFit    = TRUE;
-    t_pargs           pa []   = {
-        { "-ninterm", FALSE, etINT,  {&ninterm},
-          "Number of intermediates" },
-        { "-first",   FALSE, etREAL, {&first},
-          "Corresponds to first generated structure (0 is input x[SUB]1[sub], see above)" },
-        { "-last",    FALSE, etREAL, {&last},
-          "Corresponds to last generated structure (1 is input x[SUB]2[sub], see above)" },
-        { "-fit",     FALSE, etBOOL, {&bFit},
-          "Do a least squares fit of the second to the first structure before interpolating" }
-    };
-    const char       *leg[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
-    FILE             *fp    = nullptr;
-    int               i, isize, is_lsq, nat1, nat2;
-    t_trxstatus      *status;
-    int              *index, *index_lsq, *index_all, *dummy;
-    rvec             *x1, *x2, *xx;
-    matrix            box;
-    real              rms1, rms2, fac, *mass;
-    char             *grpname;
-    gmx_bool          bRMS;
-    gmx_output_env_t *oenv;
-
-    if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
-                           NFILE, fnm, asize(pa), pa, asize(desc), desc,
-                           0, nullptr, &oenv))
-    {
-        return 0;
-    }
-
-    t_topology *top;
-    snew(top, 1);
-    read_tps_conf(opt2fn("-f1", NFILE, fnm), top, nullptr, &x1, nullptr, box, FALSE);
-    nat1 = top->atoms.nr;
-    read_tps_conf(opt2fn("-f2", NFILE, fnm), top, nullptr, &x2, nullptr, box, FALSE);
-    nat2 = top->atoms.nr;
-    if (nat1 != nat2)
-    {
-        gmx_fatal(FARGS, "Number of atoms in first structure is %d, in second %d",
-                  nat1, nat2);
-    }
-    snew(xx, nat1);
-    t_atoms  &atoms = top->atoms;
-
-    snew(mass, nat1);
-    snew(index_all, nat1);
-    for (i = 0; (i < nat1); i++)
-    {
-        mass[i]      = 1;
-        index_all[i] = i;
-    }
-    if (bFit)
-    {
-        printf("Select group for LSQ superposition:\n");
-        get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &is_lsq, &index_lsq,
-                  &grpname);
-        reset_x(is_lsq, index_lsq, nat1, index_all, x1, mass);
-        reset_x(is_lsq, index_lsq, nat1, index_all, x2, mass);
-        do_fit(nat1, mass, x1, x2);
-    }
-
-    bRMS = opt2bSet("-or", NFILE, fnm);
-    if (bRMS)
-    {
-        fp = xvgropen(opt2fn("-or", NFILE, fnm), "RMSD", "Conf", "(nm)", oenv);
-        xvgr_legend(fp, asize(leg), leg, oenv);
-        printf("Select group for RMSD calculation:\n");
-        get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &isize, &index, &grpname);
-        printf("You selected group %s, containing %d atoms\n", grpname, isize);
-        rms1 = rmsdev_ind(isize, index, mass, x1, x2);
-        fprintf(stderr, "RMSD between input conformations is %g nm\n", rms1);
-    }
-
-    snew(dummy, nat1);
-    for (i = 0; (i < nat1); i++)
-    {
-        dummy[i] = i;
-    }
-    status = open_trx(ftp2fn(efTRX, NFILE, fnm), "w");
-
-    for (i = 0; (i < ninterm); i++)
-    {
-        fac = dointerp(nat1, x1, x2, xx, i, ninterm, first, last);
-        write_trx(status, nat1, dummy, &atoms, i, fac, box, xx, nullptr, nullptr);
-        if (bRMS)
-        {
-            rms1 = rmsdev_ind(isize, index, mass, x1, xx);
-            rms2 = rmsdev_ind(isize, index, mass, x2, xx);
-            fprintf(fp, "%10g  %10g  %10g\n", fac, rms1, rms2);
-        }
-    }
-
-    close_trx(status);
-
-    if (bRMS)
-    {
-        xvgrclose(fp);
-        do_view(oenv, opt2fn("-or", NFILE, fnm), "-nxy");
-    }
-
-    return 0;
-}
index ec7feb678f3c5790602819a41d4160b14025df60..0ffe400b984f6491d1ea6bfa4c037fc12b59e09b 100644 (file)
@@ -41,6 +41,7 @@
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/commandline/viewit.h"
+#include "gromacs/compat/make_unique.h"
 #include "gromacs/fileio/confio.h"
 #include "gromacs/fileio/trxio.h"
 #include "gromacs/fileio/xvgr.h"
@@ -68,7 +69,7 @@ typedef enum {
     NOT_USED, NORMAL, X, Y, Z, LATERAL
 } msd_type;
 
-typedef struct {
+struct t_corr {
     real          t0;         /* start time and time increment between  */
     real          delta_t;    /* time between restart points */
     real          beginfit,   /* the begin/end time for fits as reals between */
@@ -94,82 +95,81 @@ typedef struct {
     int          *n_offs;
     int         **ndata;      /* the number of msds (particles/mols) per data
                                  point. */
-} t_corr;
-
-typedef real t_calc_func (t_corr *curr, int nx, const int index[], int nx0, rvec xc[],
-                          const rvec dcom, gmx_bool bTen, matrix mat);
-
-static real thistime(t_corr *curr)
-{
-    return curr->time[curr->nframes];
-}
-
-static int in_data(t_corr *curr, int nx00)
-{
-    return curr->nframes-curr->n_offs[nx00];
-}
-
-static t_corr *init_corr(int nrgrp, int type, int axis, real dim_factor,
-                         int nmol, gmx_bool bTen, gmx_bool bMass, real dt, const t_topology *top,
-                         real beginfit, real endfit)
-{
-    t_corr  *curr;
-    int      i;
-
-    snew(curr, 1);
-    curr->type       = static_cast<msd_type>(type);
-    curr->axis       = axis;
-    curr->ngrp       = nrgrp;
-    curr->nrestart   = 0;
-    curr->delta_t    = dt;
-    curr->beginfit   = (1 - 2*GMX_REAL_EPS)*beginfit;
-    curr->endfit     = (1 + 2*GMX_REAL_EPS)*endfit;
-    curr->x0         = nullptr;
-    curr->n_offs     = nullptr;
-    curr->nframes    = 0;
-    curr->nlast      = 0;
-    curr->dim_factor = dim_factor;
-
-    snew(curr->ndata, nrgrp);
-    snew(curr->data, nrgrp);
-    if (bTen)
+    t_corr(int nrgrp, int type, int axis, real dim_factor, int nmol,
+           gmx_bool bTen, gmx_bool bMass, real dt, const t_topology *top,
+           real beginfit, real endfit) :
+        t0(0),
+        delta_t(dt),
+        beginfit((1 - 2*GMX_REAL_EPS)*beginfit),
+        endfit((1 + 2*GMX_REAL_EPS)*endfit),
+        dim_factor(dim_factor),
+        data(nullptr),
+        time(nullptr),
+        mass(nullptr),
+        datam(nullptr),
+        x0(nullptr),
+        com(nullptr),
+        lsq(nullptr),
+        type(static_cast<msd_type>(type)),
+        axis(axis),
+        ncoords(0),
+        nrestart(0),
+        nmol(nmol),
+        nframes(0),
+        nlast(0),
+        ngrp(nrgrp),
+        n_offs(nullptr),
+        ndata(nullptr)
     {
-        snew(curr->datam, nrgrp);
-    }
-    for (i = 0; (i < nrgrp); i++)
-    {
-        curr->ndata[i] = nullptr;
-        curr->data[i]  = nullptr;
+        snew(ndata, nrgrp);
+        snew(data, nrgrp);
         if (bTen)
         {
-            curr->datam[i] = nullptr;
+            snew(datam, nrgrp);
         }
-    }
-    curr->time = nullptr;
-    curr->lsq  = nullptr;
-    curr->nmol = nmol;
-    if (curr->nmol > 0)
-    {
-        snew(curr->mass, curr->nmol);
-        for (i = 0; i < curr->nmol; i++)
+        for (int i = 0; (i < nrgrp); i++)
+        {
+            ndata[i] = nullptr;
+            data[i]  = nullptr;
+            if (bTen)
+            {
+                datam[i] = nullptr;
+            }
+        }
+        if (nmol > 0)
         {
-            curr->mass[i] = 1;
+            snew(mass, nmol);
+            for (int i = 0; i < nmol; i++)
+            {
+                mass[i] = 1;
+            }
         }
-    }
-    else
-    {
-        if (bMass)
+        else
         {
-            const t_atoms *atoms = &top->atoms;
-            snew(curr->mass, atoms->nr);
-            for (i = 0; (i < atoms->nr); i++)
+            if (bMass)
             {
-                curr->mass[i] = atoms->atom[i].m;
+                const t_atoms *atoms = &top->atoms;
+                snew(mass, atoms->nr);
+                for (int i = 0; (i < atoms->nr); i++)
+                {
+                    mass[i] = atoms->atom[i].m;
+                }
             }
         }
     }
+};
 
-    return curr;
+typedef real t_calc_func (t_corr *curr, int nx, const int index[], int nx0, rvec xc[],
+                          const rvec dcom, gmx_bool bTen, matrix mat);
+
+static real thistime(t_corr *curr)
+{
+    return curr->time[curr->nframes];
+}
+
+static int in_data(t_corr *curr, int nx00)
+{
+    return curr->nframes-curr->n_offs[nx00];
 }
 
 static void corr_print(t_corr *curr, gmx_bool bTen, const char *fn, const char *title,
@@ -878,17 +878,17 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_
                     int type, real dim_factor, int axis,
                     real dt, real beginfit, real endfit, const gmx_output_env_t *oenv)
 {
-    t_corr        *msd;
-    int           *gnx;   /* the selected groups' sizes */
-    int          **index; /* selected groups' indices */
-    char         **grpname;
-    int            i, i0, i1, j, N, nat_trx;
-    real          *DD, *SigmaD, a, a2, b, r, chi2;
-    rvec          *x = nullptr;
-    matrix         box;
-    int           *gnx_com     = nullptr; /* the COM removal group size  */
-    int          **index_com   = nullptr; /* the COM removal group atom indices */
-    char         **grpname_com = nullptr; /* the COM removal group name */
+    std::unique_ptr<t_corr> msd;
+    int                    *gnx;   /* the selected groups' sizes */
+    int                   **index; /* selected groups' indices */
+    char                  **grpname;
+    int                     i, i0, i1, j, N, nat_trx;
+    real                   *DD, *SigmaD, a, a2, b, r, chi2;
+    rvec                   *x = nullptr;
+    matrix                  box;
+    int                    *gnx_com     = nullptr; /* the COM removal group size  */
+    int                   **index_com   = nullptr; /* the COM removal group atom indices */
+    char                  **grpname_com = nullptr; /* the COM removal group name */
 
     snew(gnx, nrgrp);
     snew(index, nrgrp);
@@ -912,12 +912,12 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_
         index_atom2mol(&gnx[0], index[0], &top->mols);
     }
 
-    msd = init_corr(nrgrp, type, axis, dim_factor,
-                    mol_file == nullptr ? 0 : gnx[0], bTen, bMW, dt, top,
-                    beginfit, endfit);
+    msd = gmx::compat::make_unique<t_corr>(nrgrp, type, axis, dim_factor,
+                                           mol_file == nullptr ? 0 : gnx[0],
+                                           bTen, bMW, dt, top, beginfit, endfit);
 
     nat_trx =
-        corr_loop(msd, trx_file, top, ePBC, mol_file ? gnx[0] != 0 : false, gnx, index,
+        corr_loop(msd.get(), trx_file, top, ePBC, mol_file ? gnx[0] != 0 : false, gnx, index,
                   (mol_file != nullptr) ? calc1_mol : (bMW ? calc1_mw : calc1_norm),
                   bTen, gnx_com, index_com, dt, t_pdb,
                   pdb_file ? &x : nullptr, box, oenv);
@@ -948,7 +948,7 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_
         {
             snew(top->atoms.pdbinfo, top->atoms.nr);
         }
-        printmol(msd, mol_file, pdb_file, index[0], top, x, ePBC, box, oenv);
+        printmol(msd.get(), mol_file, pdb_file, index[0], top, x, ePBC, box, oenv);
         top->atoms.nr = i;
     }
 
@@ -1021,7 +1021,7 @@ static void do_corr(const char *trx_file, const char *ndx_file, const char *msd_
         }
     }
     /* Print mean square displacement */
-    corr_print(msd, bTen, msd_file,
+    corr_print(msd.get(), bTen, msd_file,
                "Mean Square Displacement",
                "MSD (nm\\S2\\N)",
                msd->time[msd->nframes-1], beginfit, endfit, DD, SigmaD, grpname, oenv);
index 2485d5a2e222325751ade6e26cf086dd1c903ca9..970ae5c1402faae5ddc0f521e3a7827c82cccb86 100644 (file)
@@ -180,10 +180,10 @@ static void check_eg_vs_cg(gmx_mtop_t *mtop)
             {
                 /* Get the energy group of the first atom in this charge group */
                 firstj  = astart + molt->cgs.index[cg];
-                firsteg = getGroupType(&mtop->groups, egcENER, firstj);
+                firsteg = getGroupType(mtop->groups, egcENER, firstj);
                 for (j = molt->cgs.index[cg]+1; j < molt->cgs.index[cg+1]; j++)
                 {
-                    eg = getGroupType(&mtop->groups, egcENER, astart+j);
+                    eg = getGroupType(mtop->groups, egcENER, astart+j);
                     if (eg != firsteg)
                     {
                         gmx_fatal(FARGS, "atoms %d and %d in charge group %d of molecule type '%s' are in different energy groups",
index d0df5b43234ab3b5b19f0c7543bb8c52f6cd5be5..9d59e2f4d94758b1701b50b7e17bb101182273c5 100644 (file)
@@ -2760,7 +2760,6 @@ static bool do_numbering(int natoms, gmx_groups_t *groups,
 static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
 {
     t_grpopts              *opts;
-    const gmx_groups_t     *groups;
     pull_params_t          *pull;
     int                     natoms, ai, aj, i, j, d, g, imin, jmin;
     int                    *nrdf2, *na_vcm, na_tot;
@@ -2778,24 +2777,24 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
 
     opts = &ir->opts;
 
-    groups = &mtop->groups;
+    const gmx_groups_t &groups = mtop->groups;
     natoms = mtop->natoms;
 
     /* Allocate one more for a possible rest group */
     /* We need to sum degrees of freedom into doubles,
      * since floats give too low nrdf's above 3 million atoms.
      */
-    snew(nrdf_tc, groups->grps[egcTC].nr+1);
-    snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
-    snew(dof_vcm, groups->grps[egcVCM].nr+1);
-    snew(na_vcm, groups->grps[egcVCM].nr+1);
-    snew(nrdf_vcm_sub, groups->grps[egcVCM].nr+1);
+    snew(nrdf_tc, groups.grps[egcTC].nr+1);
+    snew(nrdf_vcm, groups.grps[egcVCM].nr+1);
+    snew(dof_vcm, groups.grps[egcVCM].nr+1);
+    snew(na_vcm, groups.grps[egcVCM].nr+1);
+    snew(nrdf_vcm_sub, groups.grps[egcVCM].nr+1);
 
-    for (i = 0; i < groups->grps[egcTC].nr; i++)
+    for (i = 0; i < groups.grps[egcTC].nr; i++)
     {
         nrdf_tc[i] = 0;
     }
-    for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
+    for (i = 0; i < groups.grps[egcVCM].nr+1; i++)
     {
         nrdf_vcm[i]     = 0;
         clear_ivec(dof_vcm[i]);
@@ -2932,7 +2931,7 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
                     nrdf_vcm[getGroupType(groups, egcVCM, ai)] -= 0.5*imin;
                     if (nrdf_tc[getGroupType(groups, egcTC, ai)] < 0)
                     {
-                        gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[getGroupType(groups, egcTC, ai)]]);
+                        gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups.grps[egcTC].nm_ind[getGroupType(groups, egcTC, ai)]]);
                     }
                 }
                 else
@@ -2955,7 +2954,7 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
          * the number of degrees of freedom in each vcm group when COM
          * translation is removed and 6 when rotation is removed as well.
          */
-        for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
+        for (j = 0; j < groups.grps[egcVCM].nr+1; j++)
         {
             switch (ir->comm_mode)
             {
@@ -2978,10 +2977,10 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
             }
         }
 
-        for (i = 0; i < groups->grps[egcTC].nr; i++)
+        for (i = 0; i < groups.grps[egcTC].nr; i++)
         {
             /* Count the number of atoms of TC group i for every VCM group */
-            for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
+            for (j = 0; j < groups.grps[egcVCM].nr+1; j++)
             {
                 na_vcm[j] = 0;
             }
@@ -2999,7 +2998,7 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
              */
             nrdf_uc    = nrdf_tc[i];
             nrdf_tc[i] = 0;
-            for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
+            for (j = 0; j < groups.grps[egcVCM].nr+1; j++)
             {
                 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
                 {
@@ -3009,7 +3008,7 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
             }
         }
     }
-    for (i = 0; (i < groups->grps[egcTC].nr); i++)
+    for (i = 0; (i < groups.grps[egcTC].nr); i++)
     {
         opts->nrdf[i] = nrdf_tc[i];
         if (opts->nrdf[i] < 0)
@@ -3018,7 +3017,7 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
         }
         fprintf(stderr,
                 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
-                gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
+                gnames[groups.grps[egcTC].nm_ind[i]], opts->nrdf[i]);
     }
 
     sfree(nrdf2);
@@ -4147,7 +4146,7 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
         const t_atom *atom;
         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
         {
-            mgrp[getGroupType(&sys->groups, egcACC, i)] += atom->m;
+            mgrp[getGroupType(sys->groups, egcACC, i)] += atom->m;
         }
         mt = 0.0;
         for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
index 8b8e6594593cded197500e4a7449b8013eb2a631..c4e5f36e9e18f0374ec32f333fc6dcb99d486266 100644 (file)
@@ -104,12 +104,7 @@ TEST_F(SolvateTest, cs_cp_p_Works)
         "solvate", "-cs"
     };
     setInputFile("-cp", "spc-and-methanol.gro");
-
-    // TODO: Consider adding a convenience method for this.
-    std::string topFileName           = gmx::test::TestFileManager::getInputFilePath("spc-and-methanol.top");
-    std::string modifiableTopFileName = fileManager().getTemporaryFilePath(".top");
-    gmx_file_copy(topFileName.c_str(), modifiableTopFileName.c_str(), true);
-    commandLine().addOption("-p", modifiableTopFileName);
+    setModifiableInputFile("-p", "spc-and-methanol.top");
 
     runTest(CommandLine(cmdline));
 }
@@ -134,14 +129,7 @@ TEST_F(SolvateTest, update_Topology_Works)
     };
     setInputFile("-cs", "mixed_solvent.gro");
     setInputFile("-cp", "simple.gro");
-
-    // TODO: Consider adding a convenience method for this.
-    // Copies topology file to where it would be found as an output file, so the copied
-    // .top file is used as both input and output
-    std::string topFileName           = gmx::test::TestFileManager::getInputFilePath("simple.top");
-    std::string modifiableTopFileName = fileManager().getTemporaryFilePath("simple.top");
-    gmx_file_copy(topFileName.c_str(), modifiableTopFileName.c_str(), true);
-    setOutputFile("-p", "simple.top", ExactTextMatch());
+    setInputAndOutputFile("-p", "simple.top", ExactTextMatch());
 
     runTest(CommandLine(cmdline));
 }
index 93d1c5772977dcec939fb6381976f4bfaa2f1b44..786cdab1b913d4a724561c282857ed58ec925d76 100644 (file)
@@ -202,7 +202,7 @@ void init_orires(FILE                 *fplog,
     od->nref = 0;
     for (int i = 0; i < mtop->natoms; i++)
     {
-        if (getGroupType(&mtop->groups, egcORFIT, i) == 0)
+        if (getGroupType(mtop->groups, egcORFIT, i) == 0)
         {
             od->nref++;
         }
index 6475581325ad4de78ec3e9c3d06f211b16b4e038..0c9cefce72cbefea53e08f11668621bacf178185 100644 (file)
@@ -675,8 +675,8 @@ static cginfo_mb_t *init_cginfo_mb(FILE *fplog, const gmx_mtop_t *mtop,
             {
                 a0 = cgs->index[cg];
                 a1 = cgs->index[cg+1];
-                if (getGroupType(&mtop->groups, egcENER, a_offset+am+a0) !=
-                    getGroupType(&mtop->groups, egcENER, a_offset   +a0))
+                if (getGroupType(mtop->groups, egcENER, a_offset+am+a0) !=
+                    getGroupType(mtop->groups, egcENER, a_offset   +a0))
                 {
                     bId = FALSE;
                 }
@@ -732,7 +732,7 @@ static cginfo_mb_t *init_cginfo_mb(FILE *fplog, const gmx_mtop_t *mtop,
                 a1 = cgs->index[cg+1];
 
                 /* Store the energy group in cginfo */
-                gid = getGroupType(&mtop->groups, egcENER, a_offset+am+a0);
+                gid = getGroupType(mtop->groups, egcENER, a_offset+am+a0);
                 SET_CGINFO_GID(cginfo[cgm+cg], gid);
 
                 /* Check the intra/inter charge group exclusions */
index 10a6359521792c896ad79016d812f400cd7777e6..6c015cc6eaae22cf7e20491c93e50f6be6339ee3 100644 (file)
@@ -129,7 +129,7 @@ makeMDAtoms(FILE *fp, const gmx_mtop_t &mtop, const t_inputrec &ir,
     md->bVCMgrps = FALSE;
     for (int i = 0; i < mtop.natoms; i++)
     {
-        if (getGroupType(&mtop.groups, egcVCM, i) > 0)
+        if (getGroupType(mtop.groups, egcVCM, i) > 0)
         {
             md->bVCMgrps = TRUE;
         }
@@ -207,16 +207,15 @@ void atoms2md(const gmx_mtop_t *mtop, const t_inputrec *ir,
 {
     gmx_bool              bLJPME;
     const t_grpopts      *opts;
-    const gmx_groups_t   *groups;
     int                   nthreads gmx_unused;
 
     bLJPME = EVDW_PME(ir->vdwtype);
 
     opts = &ir->opts;
 
-    groups = &mtop->groups;
+    const gmx_groups_t &groups = mtop->groups;
 
-    auto md = mdAtoms->mdatoms();
+    auto                md = mdAtoms->mdatoms();
     /* nindex>=0 indicates DD where we use an index */
     if (nindex >= 0)
     {
@@ -369,7 +368,7 @@ void atoms2md(const gmx_mtop_t *mtop, const t_inputrec *ir,
                 else
                 {
                     /* The friction coefficient is mass/tau_t */
-                    fac = ir->delta_t/opts->tau_t[md->cTC ? groups->grpnr[egcTC][ag] : 0];
+                    fac = ir->delta_t/opts->tau_t[md->cTC ? groups.grpnr[egcTC][ag] : 0];
                     mA  = 0.5*atom.m*fac;
                     mB  = 0.5*atom.mB*fac;
                 }
@@ -466,16 +465,16 @@ void atoms2md(const gmx_mtop_t *mtop, const t_inputrec *ir,
             md->ptype[i]    = atom.ptype;
             if (md->cTC)
             {
-                md->cTC[i]    = groups->grpnr[egcTC][ag];
+                md->cTC[i]    = groups.grpnr[egcTC][ag];
             }
             md->cENER[i]    = getGroupType(groups, egcENER, ag);
             if (md->cACC)
             {
-                md->cACC[i]   = groups->grpnr[egcACC][ag];
+                md->cACC[i]   = groups.grpnr[egcACC][ag];
             }
             if (md->cVCM)
             {
-                md->cVCM[i]       = groups->grpnr[egcVCM][ag];
+                md->cVCM[i]       = groups.grpnr[egcVCM][ag];
             }
             if (md->cORF)
             {
@@ -484,17 +483,17 @@ void atoms2md(const gmx_mtop_t *mtop, const t_inputrec *ir,
 
             if (md->cU1)
             {
-                md->cU1[i]        = groups->grpnr[egcUser1][ag];
+                md->cU1[i]        = groups.grpnr[egcUser1][ag];
             }
             if (md->cU2)
             {
-                md->cU2[i]        = groups->grpnr[egcUser2][ag];
+                md->cU2[i]        = groups.grpnr[egcUser2][ag];
             }
 
             if (ir->bQMMM)
             {
-                if (groups->grpnr[egcQMMM] == nullptr ||
-                    groups->grpnr[egcQMMM][ag] < groups->grps[egcQMMM].nr-1)
+                if (groups.grpnr[egcQMMM] == nullptr ||
+                    groups.grpnr[egcQMMM][ag] < groups.grps[egcQMMM].nr-1)
                 {
                     md->bQM[i]      = TRUE;
                 }
index 2e4286624a8bc9fb98928c5b9be4919d63bc3499..6ecf24a2e260208fc36ca046738dc693d346eaeb 100644 (file)
@@ -209,7 +209,7 @@ gmx_mdoutf_t init_mdoutf(FILE *fplog, int nfile, const t_filenm fnm[],
         of->natoms_x_compressed = 0;
         for (i = 0; (i < top_global->natoms); i++)
         {
-            if (getGroupType(of->groups, egcCompressedX, i) == 0)
+            if (getGroupType(*of->groups, egcCompressedX, i) == 0)
             {
                 of->natoms_x_compressed++;
             }
@@ -361,7 +361,7 @@ void mdoutf_write_to_trajectory_files(FILE *fplog, const t_commrec *cr,
                 auto x = makeArrayRef(state_global->x);
                 for (i = 0, j = 0; (i < of->natoms_global); i++)
                 {
-                    if (getGroupType(of->groups, egcCompressedX, i) == 0)
+                    if (getGroupType(*of->groups, egcCompressedX, i) == 0)
                     {
                         copy_rvec(x[i], xxtc[j++]);
                     }
index 6ad71ef9f906d3bdf0742e5408739763358cda4f..51fc6cfbcc89fe03a0fe807a0f3981533970d1c1 100644 (file)
@@ -2194,6 +2194,10 @@ void init_ns(FILE *fplog, const t_commrec *cr,
 
 void done_ns(gmx_ns_t *ns, int numEnergyGroups)
 {
+    if (ns->bexcl != nullptr)
+    {
+        sfree(ns->bexcl);
+    }
     sfree(ns->bExcludeAlleg);
     if (ns->ns_buf)
     {
@@ -2203,8 +2207,20 @@ void done_ns(gmx_ns_t *ns, int numEnergyGroups)
         }
         sfree(ns->ns_buf);
     }
+    if (ns->nl_sr != nullptr)
+    {
+        for (int i = 0; i < numEnergyGroups; i++)
+        {
+            if (ns->nl_sr[i] != nullptr)
+            {
+                sfree(ns->nl_sr[i]);
+            }
+        }
+        sfree(ns->nl_sr);
+    }
     sfree(ns->simple_aaj);
     sfree(ns->bHaveVdW);
+    sfree(ns->nsr);
     done_grid(ns->grid);
     sfree(ns);
 }
index a0081e6fe2b28e6a313b3b97807f62608cfc76b9..0b751cc86ca8e52eb89690946e73733812a884dd 100644 (file)
@@ -420,7 +420,7 @@ void init_QMMMrec(const t_commrec  *cr,
         qr->nrQMlayers = 1;
     }
 
-    const gmx_groups_t *groups = &mtop->groups;
+    const gmx_groups_t &groups = mtop->groups;
 
     /* there are jmax groups of QM atoms. In case of multiple QM groups
      * I assume that the users wants to do ONIOM. However, maybe it
index e52e74c137518f39241d8f138c7b27d33c2d39ab..fdae985a973050a5fc05da4a4d14a423fba75bef 100644 (file)
@@ -71,14 +71,13 @@ static void init_grptcstat(int ngtc, t_grp_tcstat tcstat[])
 
 static void init_grpstat(const gmx_mtop_t *mtop, int ngacc, t_grp_acc gstat[])
 {
-    const gmx_groups_t     *groups;
     gmx_mtop_atomloop_all_t aloop;
     int                     i, grp;
     const t_atom           *atom;
 
     if (ngacc > 0)
     {
-        groups = &mtop->groups;
+        const gmx_groups_t &groups = mtop->groups;
         aloop  = gmx_mtop_atomloop_all_init(mtop);
         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
         {
index ee93053af17439e8af65a94aa54a298cdcff8e79..765ffd18a627c57c52c7687a4235ee1284a6a9f2 100644 (file)
@@ -1526,6 +1526,13 @@ void gmx::Integrator::do_md()
     walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
 
     destroy_enerdata(enerd);
+
     sfree(enerd);
+
+    /* Clean up topology. top->atomtypes has an allocated pointer if no domain decomposition*/
+    if (!DOMAINDECOMP(cr))
+    {
+        done_atomtypes(&top->atomtypes);
+    }
     sfree(top);
 }
diff --git a/src/gromacs/mdspan/CMakeLists.txt b/src/gromacs/mdspan/CMakeLists.txt
new file mode 100644 (file)
index 0000000..7af2f8f
--- /dev/null
@@ -0,0 +1,37 @@
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2018, by the GROMACS development team, led by
+# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+# and including many others, as listed in the AUTHORS file in the
+# top-level source directory and at http://www.gromacs.org.
+#
+# GROMACS is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public License
+# as published by the Free Software Foundation; either version 2.1
+# of the License, or (at your option) any later version.
+#
+# GROMACS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with GROMACS; if not, see
+# http://www.gnu.org/licenses, or write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+#
+# If you want to redistribute modifications to GROMACS, please
+# consider that scientific software is very special. Version
+# control is crucial - bugs must be traceable. We will be happy to
+# consider code for inclusion in the official distribution, but
+# derived work must not be called official GROMACS. Details are found
+# in the README & COPYING files - if they are missing, get the
+# official version at http://www.gromacs.org.
+#
+# To help us fund GROMACS development, we humbly ask that you cite
+# the research papers on the package. Check out http://www.gromacs.org.
+
+if (BUILD_TESTING)
+    add_subdirectory(tests)
+endif()
diff --git a/src/gromacs/mdspan/README b/src/gromacs/mdspan/README
new file mode 100644 (file)
index 0000000..56607fb
--- /dev/null
@@ -0,0 +1,20 @@
+This directory contains a stripped-down and modified-for-GROMACS version
+of a reference implementation of the LEWG P0009 proposal, revision 8:
+"mdspan: A Non-Owning Multidimensional Array Reference"
+
+The full version is available from
+https://github.com/ORNL/cpp-proposals-pub/tree/master/P0009
+SHA1 of the commit used is f433984ba897298069b510aa1311145040dd860c.
+
+This work orginated from "Kokkos" and is distributed under the BSD license as
+stated in COPYING.
+
+The following modifications were carried out:
+* reformatted code with uncrustify
+* added GROMACS header to all files
+* added doxygen comments to all files
+* temporarily removed constexpr from functions, so they compile with C++11
+* put all headers from bits/ into this directory
+* use minimial includes in all header files instead of just including mdspan
+* move tests into a subdirectory tests and rename test_foo.cpp to foo.cpp
+* renamed tests and test fixture class to match GROMACS conventions
diff --git a/src/gromacs/mdspan/extents.h b/src/gromacs/mdspan/extents.h
new file mode 100644 (file)
index 0000000..a7a99b3
--- /dev/null
@@ -0,0 +1,430 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*
+ * This file is a modified version of original work of Sandia Corporation.
+ * In the spirit of the original code, this particular file can be distributed
+ * on the terms of Sandia Corporation.
+ */
+/*
+ *                          Kokkos v. 2.0
+ *               Copyright (2014) Sandia Corporation
+ *
+ * Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+ * the U.S. Government retains certain rights in this software.
+ *
+ * Kokkos is licensed under 3-clause BSD terms of use:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * 3. Neither the name of the Corporation nor the names of the
+ * contributors may be used to endorse or promote products derived from
+ * this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
+ * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * Questions? Contact Christian R. Trott (crtrott@sandia.gov)
+ */
+/*! \libinternal \file
+ * \brief Declares gmx::extents for mdspan.
+ *
+ * \author Christian Trott <crtrott@sandia.gov>
+ * \author Ronan Keryell <ronan.keryell@xilinx.com>
+ * \author Carter Edwards <hedwards@nvidia.com>
+ * \author David Hollman <dshollm@sandia.gov>
+ * \author Christian Blau <cblau@gwdg.de>
+ * \ingroup mdspan
+ */
+#ifndef MDSPAN_EXTENTS_H
+#define MDSPAN_EXTENTS_H
+
+#include <cstddef>
+
+#include <array>
+
+namespace gmx
+{
+
+/*! \brief Define constant that signals dynamic extent.
+ */
+enum : std::ptrdiff_t {
+    dynamic_extent = -1
+};
+
+template< std::ptrdiff_t ... StaticExtents >
+class extents;
+
+template<std::ptrdiff_t... LHS, std::ptrdiff_t... RHS>
+bool operator==(const extents<LHS...> &lhs,
+                const extents<RHS...> &rhs) noexcept;
+
+template<std::ptrdiff_t... LHS, std::ptrdiff_t... RHS>
+constexpr bool operator!=(const extents<LHS...> &lhs,
+                          const extents<RHS...> &rhs) noexcept;
+
+namespace detail
+{
+
+template< int R, std::ptrdiff_t ... StaticExtents >
+struct extents_analyse;
+
+/*! \libinternal \brief Enable querying extent of specific rank by splitting
+ * a static extents off the variadic template arguments.
+ *
+ */
+template< int R, std::ptrdiff_t E0, std::ptrdiff_t ... StaticExtents >
+struct extents_analyse<R, E0, StaticExtents...> {
+
+    //! The extent analysis of the next lower rank.
+    using next_extents_analyse = extents_analyse<R-1, StaticExtents...>;
+
+    /*! \brief Accumulate the total rank from all extents.
+     * \returns incremented rank of the next extent
+     */
+    static constexpr std::size_t rank() noexcept { return next_extents_analyse::rank()+1; }
+    /*! \brief Accumulate the dynamic rank from all extents.
+     * This extent is static, so hand down query to the next extent analysis.
+     * \returns the dynamic rank of the next extent analysis.
+     */
+    static constexpr std::size_t rank_dynamic() noexcept { return next_extents_analyse::rank_dynamic(); }
+
+    //! Store analysis of the next extent of next lower rank.
+    next_extents_analyse next;
+
+    //! Trivial constructor.
+    extents_analyse() : next() {}
+
+    /*! \brief Construction from dynamic extents hands the extents down
+     * to the next extents analysis of lower rank.
+     * \param[in] de dynamic extents
+     */
+    template<class ... DynamicExtents>
+    extents_analyse(DynamicExtents... de) : next(de ...) {}
+
+    /*! \brief Construct from an array of dynamic extentes and rank.
+     * Hand down the dynamic rank parameters to the next extents analysis rank
+     * \param[in] de dynamic extents
+     * \param[in] r rank to read from the dynamic extent
+     */
+    template<std::size_t Rank>
+    extents_analyse(const std::array<std::ptrdiff_t, Rank> &de, const std::size_t r) : next(de, r) {}
+
+    //! Copy constructor.
+    template<std::ptrdiff_t... OtherStaticExtents>
+    extents_analyse(extents_analyse<R, OtherStaticExtents...> rhs) : next(rhs.next) {}
+
+    //! Assignment operator.
+    template<std::ptrdiff_t... OtherStaticExtents>
+    extents_analyse &operator= (extents_analyse<R, OtherStaticExtents...> rhs)
+    {
+        next = rhs.next;
+        return *this;
+    }
+
+    /*! \brief Report extent of dimension r.
+     * \param[in] r the dimension to query
+     * \returns the extent in dimension r.
+     */
+    constexpr std::ptrdiff_t extent(const std::size_t r) const noexcept
+    {
+        // TODO use early return and if instead of ternary operator
+        // after bumping requirements to C++14
+        return (r == R) ? E0 : next.extent(r);
+    }
+    /*! \brief Report the static extent of dimension r.
+     * \param[in] r the dimension to query
+     * \returns the static extent in dimension r.
+     */
+    static constexpr std::ptrdiff_t static_extent(const std::size_t r) noexcept
+    {
+        // TODO use early return and if instead of ternary operator
+        // after bumping requirements to C++14
+        return (r == R) ? E0 : next_extents_analyse::static_extent(r);
+    }
+};
+
+/*! \libinternal \brief Enable querying extent of specific rank by splitting
+ * a static extent off the variadic template arguments.
+ */
+template< int R, std::ptrdiff_t ... StaticExtents >
+struct extents_analyse<R, dynamic_extent, StaticExtents...> {
+    //! The extent analysis of the next lower rank.
+    using next_extents_analyse = extents_analyse<R-1, StaticExtents...>;
+    /*! \brief Accumulate the total rank from all extents.
+     * \returns incremented rank of the next extent
+     */
+    static constexpr std::size_t rank() noexcept { return next_extents_analyse::rank()+1; }
+    /*! \brief Accumulate the dynamic rank from all extents.
+     * \returns the dynamic rank of the next extent analysis.
+     */
+    static constexpr std::size_t rank_dynamic() noexcept { return next_extents_analyse::rank_dynamic()+1; }
+
+    //! Store analysis of the next extent of next lower rank.
+    next_extents_analyse next;
+    //! The dynamic extent of this rank
+    std::ptrdiff_t       this_extent;
+
+    //! Trivial constructor.
+    extents_analyse() : next(), this_extent(0) {}
+
+    /*! \brief Construction from dynamic extents hands the extents down
+     * to the next extents analysis of lower rank.
+     * \param[in] E the dynamic extent of this rank.
+     * \param[in] de dynamic extents
+     */
+    template<class ... DynamicExtents>
+    extents_analyse(std::ptrdiff_t E, DynamicExtents... de) : next(de ...), this_extent(E) {}
+
+    /*! \brief Construct from an array of dynamic extentes and rank.
+     * Hand down the dynamic rank parameters to the next extents analysis rank
+     * \param[in] de dynamic extents
+     * \param[in] r rank to read from the dynamic extent
+     */
+    template<std::size_t Rank>
+    extents_analyse(const std::array<std::ptrdiff_t, Rank> &de, const std::size_t r) : next(de, r+1), this_extent(de[r]) {}
+
+    //! Copy constructor.
+    template<std::ptrdiff_t... OtherStaticExtents>
+    extents_analyse(extents_analyse<R, OtherStaticExtents...> rhs) : next(rhs.next), this_extent(rhs.extent(R)) {}
+
+    //! Assignment operator.
+    template<std::ptrdiff_t... OtherStaticExtents>
+    extents_analyse &operator= (extents_analyse<R, OtherStaticExtents...> rhs)
+    {
+        next        = rhs.next;
+        this_extent = rhs.extent(R);
+        return *this;
+    }
+
+    /*! \brief Report extent of dimension r.
+     * \param[in] r the dimension to query
+     * \returns the extent in dimension r.
+     */
+    constexpr std::ptrdiff_t extent(const std::size_t r) const noexcept
+    {
+        return (r == R) ? this_extent : next.extent(r); // TODO use early return and if instead of ternary operator after bumping requirements to C++14
+    }
+    /*! \brief Report the static extent of dimension r.
+     * \param[in] r the dimension to query
+     * \returns the static extent in dimension r.
+     */
+    static constexpr std::ptrdiff_t static_extent(const std::size_t r) noexcept
+    {
+        return (r == R) ? dynamic_extent : next_extents_analyse::static_extent(r); // TODO use early return and if instead of ternary operator after bumping requirements to C++14
+    }
+};
+
+/*! \libinternal \brief Specialisation for rank 0 extents analysis.
+ * Ends recursive rank analysis.
+ */
+template<>
+struct extents_analyse<0> {
+    /*! \brief Rank of extent of rank 0.
+     * \returns 0
+     */
+    static constexpr std::size_t rank() noexcept { return 0; }
+    /*! \brief Dynamic rank of extent of rank 0.
+     * \returns 0
+     */
+    static constexpr std::size_t rank_dynamic() noexcept { return 0; }
+
+    //! Trivial constructor.
+    extents_analyse() {}
+
+    //! Construct from array and rank, doing nothing.
+    template<std::size_t Rank>
+    extents_analyse(const std::array<std::ptrdiff_t, Rank> & /*de*/, const std::size_t /*r*/) {}
+
+    //extents_analyse & operator=(extents_analyse) = default;
+
+    /*! \brief Extent of rank 0 is 1, ensuring that product of extents yields required size and not zero.
+     * NOTE changed from ORNL reference implementation in making this static constexpr instead of constexpr .. const
+     */
+    static constexpr std::ptrdiff_t extent(const std::size_t /*r*/) noexcept
+    {
+        return 1;
+    }
+
+    //! Static extent of rank 0 is 1, ensuring that product of extents yields required size and not zero.
+    static constexpr std::ptrdiff_t static_extent(const std::size_t /*r*/) noexcept
+    {
+        return 1;
+    }
+
+};
+}   // namespace detail
+
+
+/*! \libinternal \brief Multidimensional extents with static and dynamic dimensions.
+ *
+ * Describes a multidimensional index space of rank R.
+ * This is equivalent to the Cartesian product space of integer intervals
+ * [0, N_0) x [0, N_1) x ... x [0,N_{R-1} )
+ *
+ * Confer to P0009r8 of the Library Evolution Working Group and mdspan.extents
+ *
+ * \tparam StaticExtents rank number of extents, where the dynamic_extent
+ * constant for static extent is used to signal a dynamic extent.
+ */
+template< std::ptrdiff_t ... StaticExtents >
+class extents
+{
+    private:
+        using extents_analyse_t = detail::extents_analyse<sizeof ... (StaticExtents), StaticExtents...>;
+
+    public:
+        //! Type used to index elements.
+        using index_type = std::ptrdiff_t;
+        //! Trivial constructor
+        constexpr extents() noexcept {}
+        //! Move constructor
+        constexpr extents( extents && ) noexcept = default;
+        //! Copy constructor.
+        constexpr extents( const extents & ) noexcept = default;
+        /*! \brief Construct with dynamic extents.
+         *
+         * Allows for extents(u,v,w..) syntax when setting dynamic extents
+         *
+         * \tparam IndexType type of index
+         * \param[in] dn first dynamic index
+         * \param[in] DynamicExtents parameter pack
+         */
+        template< class ... IndexType >
+        constexpr extents( std::ptrdiff_t dn,
+                           IndexType ...  DynamicExtents ) noexcept
+            : impl( dn, DynamicExtents ... )
+        { static_assert( 1+sizeof ... (DynamicExtents) == rank_dynamic(), "" ); }
+
+        /*! \brief Construct from array of dynamic extents.
+         *
+         * Allows for extents({{u,v,w..}}) syntax when setting dynamic extents
+         *
+         * \param[in] dynamic_extents array of dynamic rank size containing extents
+         */
+        constexpr extents( const std::array<std::ptrdiff_t, extents_analyse_t::rank_dynamic()> dynamic_extents) noexcept
+            : impl(dynamic_extents, 0) {}
+
+        //! Copy constructor
+        template<std::ptrdiff_t... OtherStaticExtents>
+        extents( const extents<OtherStaticExtents...> &other )
+            : impl( other.impl ) {}
+
+        //! Default move assignment
+        extents &operator= ( extents && ) noexcept = default;
+        //! Default copy assignment
+        extents &operator= ( const extents & ) noexcept = default;
+        //! Copy assignment
+        template<std::ptrdiff_t... OtherStaticExtents>
+        extents &operator= ( const extents<OtherStaticExtents...> &other )
+        { impl = other.impl; return *this; }
+        //! Default destructor
+        ~extents() = default;
+
+        // [mdspan.extents.obs]
+        /*! \brief The rank of the extent.
+         * \returns the rank all extents together
+         */
+        static constexpr std::size_t rank() noexcept
+        { return sizeof ... (StaticExtents); }
+        /*! \brief The rank of the dynamic extents.
+         * \returns Only the dynamic extents.
+         */
+        static constexpr std::size_t rank_dynamic() noexcept
+        { return extents_analyse_t::rank_dynamic(); }
+        /*! \brief The rank of the static extents.
+         * \returns Only the static extents.
+         */
+        static constexpr index_type static_extent(std::size_t k) noexcept
+        { return extents_analyse_t::static_extent(rank()-k); }
+        /*! \brief The extent along a specific dimension.
+         * \param[in] k the dimension
+         * \returns the extent along that dimension
+         */
+        constexpr index_type extent(std::size_t k) const noexcept
+        { return impl.extent(rank()-k); }
+    private:
+        //! For copy assignment, extents are friends of extents.
+        template< std::ptrdiff_t... > friend class extents;
+        //! The implementation class.
+        extents_analyse_t impl;
+};
+
+
+/*! \brief Comparison operator.
+ * \returns true if extents are equal
+ * TODO add constexpr when C++14 is required
+ */
+template<std::ptrdiff_t... LHS, std::ptrdiff_t... RHS>
+bool operator==(const extents<LHS...> &lhs,
+                const extents<RHS...> &rhs) noexcept
+{
+    bool equal = lhs.rank() == rhs.rank();
+    for (std::size_t r = 0; r < lhs.rank(); r++)
+    {
+        equal = equal && ( lhs.extent(r) == rhs.extent(r) );
+    }
+    return equal;
+}
+
+/*! \brief Check for non-equality.
+ * \returns true if extents are unequal
+ */
+template<std::ptrdiff_t... LHS, std::ptrdiff_t... RHS>
+constexpr bool operator!=(const extents<LHS...> &lhs,
+                          const extents<RHS...> &rhs) noexcept
+{
+    return !(lhs == rhs);
+}
+
+}      // namespace gmx
+#endif /* end of include guard: MDSPAN_EXTENTS_H */
diff --git a/src/gromacs/mdspan/tests/CMakeLists.txt b/src/gromacs/mdspan/tests/CMakeLists.txt
new file mode 100644 (file)
index 0000000..2ac6063
--- /dev/null
@@ -0,0 +1,37 @@
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2018, by the GROMACS development team, led by
+# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+# and including many others, as listed in the AUTHORS file in the
+# top-level source directory and at http://www.gromacs.org.
+#
+# GROMACS is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public License
+# as published by the Free Software Foundation; either version 2.1
+# of the License, or (at your option) any later version.
+#
+# GROMACS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with GROMACS; if not, see
+# http://www.gnu.org/licenses, or write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+#
+# If you want to redistribute modifications to GROMACS, please
+# consider that scientific software is very special. Version
+# control is crucial - bugs must be traceable. We will be happy to
+# consider code for inclusion in the official distribution, but
+# derived work must not be called official GROMACS. Details are found
+# in the README & COPYING files - if they are missing, get the
+# official version at http://www.gromacs.org.
+#
+# To help us fund GROMACS development, we humbly ask that you cite
+# the research papers on the package. Check out http://www.gromacs.org.
+
+gmx_add_unit_test(MDSpanTests mdspan-test
+                  extents.cpp
+                 )
diff --git a/src/gromacs/mdspan/tests/extents.cpp b/src/gromacs/mdspan/tests/extents.cpp
new file mode 100644 (file)
index 0000000..b7ccb34
--- /dev/null
@@ -0,0 +1,203 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*
+ * This file is a modified version of original work of Sandia Corporation.
+ * In the spirit of the original code, this particular file can be distributed
+ * on the terms of Sandia Corporation.
+ */
+/*
+ *                         Kokkos v. 2.0
+ *               Copyright (2014) Sandia Corporation
+ *
+ * Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+ * the U.S. Government retains certain rights in this software.
+ *
+ * Kokkos is licensed under 3-clause BSD terms of use:
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are
+ * met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * 3. Neither the name of the Corporation nor the names of the
+ * contributors may be used to endorse or promote products derived from
+ * this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
+ * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * Questions? Contact Christian R. Trott (crtrott@sandia.gov)
+ */
+/*! \internal \file
+ * \brief Testing gmx::extents.
+ *
+ * \author Christian Trott <crtrott@sandia.gov>
+ * \author Carter Edwards <hedwards@nvidia.com>
+ * \author David Hollman <dshollm@sandia.gov>
+ * \author Christian Blau <cblau@gwdg.de>
+ */
+#include "gmxpre.h"
+
+#include "gromacs/mdspan/extents.h"
+
+#include <gtest/gtest.h>
+
+namespace gmx
+{
+
+template<ptrdiff_t ... E_STATIC>
+class ExtentsTest
+{
+    public:
+        using extents_type  = gmx::extents<E_STATIC...>;
+
+        extents_type my_extents_explicit;
+        extents_type my_extents_array;
+        extents_type my_extents_copy;
+
+        ExtentsTest()
+        {
+            my_extents_explicit = extents<E_STATIC...>();
+            my_extents_array    = extents<E_STATIC...>(std::array<ptrdiff_t, 0>());
+            my_extents_copy     = extents<E_STATIC...>(my_extents_explicit);
+        }
+
+        template<class ... E>
+        ExtentsTest(E ... e)
+        {
+            my_extents_explicit = extents<E_STATIC...>(e ...);
+            my_extents_array    = extents<E_STATIC...>(std::array<ptrdiff_t, 2>({{e ...}}));
+            my_extents_copy     = extents<E_STATIC...>(my_extents_explicit);
+        }
+
+        void check_rank(ptrdiff_t r)
+        {
+            EXPECT_EQ(my_extents_explicit.rank(), r);
+            EXPECT_EQ(my_extents_array.rank(), r);
+            EXPECT_EQ(my_extents_copy.rank(), r);
+        }
+        void check_rank_dynamic(ptrdiff_t r)
+        {
+            EXPECT_EQ(my_extents_explicit.rank_dynamic(), r);
+            EXPECT_EQ(my_extents_array.rank_dynamic(), r);
+            EXPECT_EQ(my_extents_copy.rank_dynamic(), r);
+        }
+        template<class ... E>
+        void check_extents(E ... e)
+        {
+            std::array<ptrdiff_t, extents_type::rank()> s = {{E_STATIC ...}};
+            std::array<ptrdiff_t, extents_type::rank()> a = {{e ...}};
+            for (size_t r = 0; r < extents_type::rank(); r++)
+            {
+                EXPECT_EQ(my_extents_explicit.static_extent(r), s[r]);
+                EXPECT_EQ(my_extents_explicit.extent(r), a[r]);
+
+                EXPECT_EQ(my_extents_array.static_extent(r), s[r]);
+                EXPECT_EQ(my_extents_array.extent(r), a[r]);
+
+                EXPECT_EQ(my_extents_copy.static_extent(r), s[r]);
+                EXPECT_EQ(my_extents_copy.extent(r), a[r]);
+            }
+            EXPECT_EQ(my_extents_explicit.static_extent(extents_type::rank()+1), 1);
+            EXPECT_EQ(my_extents_explicit.extent(extents_type::rank()+1), 1);
+
+            EXPECT_EQ(my_extents_array.static_extent(extents_type::rank()+1), 1);
+            EXPECT_EQ(my_extents_array.extent(extents_type::rank()+1), 1);
+
+            EXPECT_EQ(my_extents_copy.static_extent(extents_type::rank()+1), 1);
+            EXPECT_EQ(my_extents_copy.extent(extents_type::rank()+1), 1);
+        }
+
+};
+
+TEST(ExtentsTest, Construction) {
+
+    // setting two dynamic extents
+    ExtentsTest<5, dynamic_extent, 3, dynamic_extent, 1> test(4, 2);
+
+    test.check_rank(5);
+    test.check_rank_dynamic(2);
+    test.check_extents(5, 4, 3, 2, 1);
+
+}
+
+TEST(ExtentsTest, PurelyStatic) {
+    ExtentsTest<5, 4, 3> test;
+    test.check_rank(3);
+    test.check_rank_dynamic(0);
+    test.check_extents(5, 4, 3);
+}
+
+TEST(ExtentsTest, RankNought) {
+    // Can construct extents of rank nought
+    ExtentsTest<> test;
+    test.check_rank(0);
+    test.check_rank_dynamic(0);
+}
+TEST(ExtentsTest, Assignment) {
+    extents<5, dynamic_extent, 3, dynamic_extent, 1> e1(4, 2);
+    extents<5, 4, 3, 2, 1> e2;
+    e2 = e1;
+    for (size_t r = 0; r < 5; r++)
+    {
+        EXPECT_EQ(e2.extent(r), e1.extent(r));
+    }
+    extents<dynamic_extent, dynamic_extent, dynamic_extent, dynamic_extent, dynamic_extent> e3(9, 8, 7, 6, 5);
+    for (int r = 0; r < 5; r++)
+    {
+        EXPECT_EQ(e3.extent(r), 9-r);
+    }
+    e3 = e1;
+    for (int r = 0; r < 5; r++)
+    {
+        EXPECT_EQ(e3.extent(r), e1.extent(r));
+    }
+}
+} // namespace gmx
index cb4bca655595f7ae0699006a45f36945144a6491..19ba2b84ec10deaca1cd152ae7ec116ae08b887e 100644 (file)
@@ -62,7 +62,7 @@ constexpr int TYPE_INT = 0, TYPE_DOUBLE = 0;
  */
 static void MCL_init_client(const char *) // NOLINT(readability-named-parameter)
 {
-    gmx_fatal(FARGS, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
+    GMX_RELEASE_ASSERT(GMX_MIMIC, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
 }
 
 /*! \brief Stub communication library function to call in case if
@@ -70,7 +70,7 @@ static void MCL_init_client(const char *) // NOLINT(readability-named-parameter)
  */
 static void MCL_send(void *, int, int, int) // NOLINT(readability-named-parameter)
 {
-    gmx_fatal(FARGS, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
+    GMX_RELEASE_ASSERT(GMX_MIMIC, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
 }
 
 /*! \brief Stub communication library function to call in case if
@@ -78,15 +78,15 @@ static void MCL_send(void *, int, int, int) // NOLINT(readability-named-paramete
  */
 static void MCL_receive(void *, int, int, int) // NOLINT(readability-named-parameter)
 {
-    gmx_fatal(FARGS, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
+    GMX_RELEASE_ASSERT(GMX_MIMIC, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
 }
 
 /*! \brief Stub communication library function to call in case if
  * GROMACS is compiled without MiMiC. Calling causes GROMACS to exit!
  */
-static void MCL_destroy() // NOLINT(readability-named-parameter)
+static void MCL_destroy()
 {
-    gmx_fatal(FARGS, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
+    GMX_RELEASE_ASSERT(GMX_MIMIC, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
 }
 #endif
 
index f9ae38483a80b526a95aa065f69a861044e25bc2..2b6017d56efc99399746ccf8d5847970a76738a5 100644 (file)
@@ -1731,7 +1731,7 @@ static void init_pull_group_index(FILE *fplog, const t_commrec *cr,
     /* In parallel, store we need to extract localWeights from weights at DD time */
     std::vector<real>  &weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
 
-    const gmx_groups_t *groups  = &mtop->groups;
+    const gmx_groups_t &groups  = mtop->groups;
 
     /* Count frozen dimensions and (weighted) mass */
     int    nfrozen = 0;
@@ -1787,13 +1787,13 @@ static void init_pull_group_index(FILE *fplog, const t_commrec *cr,
             }
             else
             {
-                if (groups->grpnr[egcTC] == nullptr)
+                if (groups.grpnr[egcTC] == nullptr)
                 {
                     mbd = ir->delta_t/ir->opts.tau_t[0];
                 }
                 else
                 {
-                    mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
+                    mbd = ir->delta_t/ir->opts.tau_t[groups.grpnr[egcTC][ii]];
                 }
             }
             w                   *= m/mbd;
index 67b273e4188d5c48f6a9a39105cb5e26c3f16237..47ff776549039294f7df7f87ccc08ddaa462a6cf 100644 (file)
@@ -118,8 +118,8 @@ static void comp_tpx(const char *fn1, const char *fn2,
         top[0] = gmx_mtop_t_to_t_topology(&mtop[0], false);
         top[1] = gmx_mtop_t_to_t_topology(&mtop[1], false);
         cmp_top(stdout, &top[0], &top[1], ftol, abstol);
-        cmp_groups(stdout, &mtop[0].groups, &mtop[1].groups,
-                   mtop[0].natoms, mtop[1].natoms);
+        compareAtomGroups(stdout, mtop[0].groups, mtop[1].groups,
+                          mtop[0].natoms, mtop[1].natoms);
         comp_state(&state[0], &state[1], bRMSD, ftol, abstol);
     }
     else
index 24a67a3d6a6058cd7288be9d5a4da14088892c7a..117d094a989775510f5ca106a2262bd31e5fdcc0 100644 (file)
@@ -85,7 +85,6 @@ static void list_tpx(const char *fn,
     t_state       state;
     t_tpxheader   tpx;
     gmx_mtop_t    mtop;
-    gmx_groups_t *groups;
     t_topology    top;
 
     read_tpxheader(fn, &tpx, TRUE);
@@ -144,12 +143,12 @@ static void list_tpx(const char *fn,
             pr_rvecs(stdout, indent, "v", tpx.bV ? state.v.rvec_array() : nullptr, state.natoms);
         }
 
-        groups = &mtop.groups;
+        const gmx_groups_t &groups = mtop.groups;
 
         snew(gcount, egcNR);
         for (i = 0; (i < egcNR); i++)
         {
-            snew(gcount[i], groups->grps[i].nr);
+            snew(gcount[i], groups.grps[i].nr);
         }
 
         for (i = 0; (i < mtop.natoms); i++)
@@ -164,7 +163,7 @@ static void list_tpx(const char *fn,
         {
             atot = 0;
             printf("%-12s: ", gtypes[i]);
-            for (j = 0; (j < groups->grps[i].nr); j++)
+            for (j = 0; (j < groups.grps[i].nr); j++)
             {
                 printf("  %5d", gcount[i][j]);
                 atot += gcount[i][j];
index 2dc6509752b9d215b46567d08bab98f232a5657d..9e2caf22daa9e6a5f49b932ca3853c617c343798 100644 (file)
@@ -599,30 +599,28 @@ void cmp_top(FILE *fp, const t_topology *t1, const t_topology *t2, real ftol, re
     }
 }
 
-void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1,
-                int natoms0, int natoms1)
+void compareAtomGroups(FILE *fp, const gmx_groups_t &g0, const gmx_groups_t &g1,
+                       int natoms0, int natoms1)
 {
-    char buf[32];
-
     fprintf(fp, "comparing groups\n");
 
     for (int i = 0; i < egcNR; i++)
     {
-        sprintf(buf, "grps[%d].nr", i);
-        cmp_int(fp, buf, -1, g0->grps[i].nr, g1->grps[i].nr);
-        if (g0->grps[i].nr == g1->grps[i].nr)
+        std::string buf = gmx::formatString("grps[%d].nr", i);
+        cmp_int(fp, buf.c_str(), -1, g0.grps[i].nr, g1.grps[i].nr);
+        if (g0.grps[i].nr == g1.grps[i].nr)
         {
-            for (int j = 0; j < g0->grps[i].nr; j++)
+            for (int j = 0; j < g0.grps[i].nr; j++)
             {
-                sprintf(buf, "grps[%d].name[%d]", i, j);
-                cmp_str(fp, buf, -1,
-                        *g0->grpname[g0->grps[i].nm_ind[j]],
-                        *g1->grpname[g1->grps[i].nm_ind[j]]);
+                buf = gmx::formatString("grps[%d].name[%d]", i, j);
+                cmp_str(fp, buf.c_str(), -1,
+                        *g0.grpname[g0.grps[i].nm_ind[j]],
+                        *g1.grpname[g1.grps[i].nm_ind[j]]);
             }
         }
-        cmp_int(fp, "ngrpnr", i, g0->ngrpnr[i], g1->ngrpnr[i]);
-        if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
-            (g0->grpnr[i] != nullptr || g1->grpnr[i] != nullptr))
+        cmp_int(fp, "ngrpnr", i, g0.ngrpnr[i], g1.ngrpnr[i]);
+        if (g0.ngrpnr[i] == g1.ngrpnr[i] && natoms0 == natoms1 &&
+            (g0.grpnr[i] != nullptr || g1.grpnr[i] != nullptr))
         {
             for (int j = 0; j < natoms0; j++)
             {
@@ -635,9 +633,9 @@ void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1,
      */
 }
 
-int getGroupType(const gmx_groups_t *group, int type, int atom)
+int getGroupType(const gmx_groups_t &group, int type, int atom)
 {
-    return (group->grpnr[type] ? group->grpnr[type][atom] : 0);
+    return (group.grpnr[type] ? group.grpnr[type][atom] : 0);
 }
 
 void copy_moltype(const gmx_moltype_t *src, gmx_moltype_t *dst)
index d958ffe750cbbd38a812e58f33e3e548c81dd804..8d46a12e9d80b107b8e5f95d96454343333aadb3 100644 (file)
@@ -118,7 +118,7 @@ typedef struct gmx_groups_t
  * \param[in] type  Type of group to check.
  * \param[in] atom  Atom to check if it has an entry.
  */
-int getGroupType (const gmx_groups_t *group, int type, int atom);
+int getGroupType (const gmx_groups_t &group, int type, int atom);
 
 /* The global, complete system topology struct, based on molecule types.
  * This structure should contain no data that is O(natoms) in memory.
@@ -226,8 +226,17 @@ void pr_top(FILE *fp, int indent, const char *title, const t_topology *top,
             gmx_bool bShowNumbers, gmx_bool bShowParameters);
 
 void cmp_top(FILE *fp, const t_topology *t1, const t_topology *t2, real ftol, real abstol);
-void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1,
-                int natoms0, int natoms1);
+
+/*! \brief Compare groups.
+ *
+ * \param[in] fp File pointer to write to.
+ * \param[in] g0 First group for comparison.
+ * \param[in] g1 Second group for comparison.
+ * \param[in] natoms0 Number of atoms for first group.
+ * \param[in] natoms1 Number of atoms for second group.
+ */
+void compareAtomGroups(FILE *fp, const gmx_groups_t &g0, const gmx_groups_t &g1,
+                       int natoms0, int natoms1);
 
 //! Deleter for gmx_localtop_t, needed until it has a proper destructor.
 using ExpandedTopologyPtr = gmx::unique_cptr<gmx_localtop_t, done_and_sfree_localtop>;
index e331c0f8811e50f543d5e542c1240c57a8116a45..21fc96e7ac9c1779bcc80d67a61e4599629278f0 100644 (file)
@@ -814,6 +814,7 @@ std::string getCoolQuote()
         { "Why would the backup server database get corrupted anyway?", "Stefan Fleischmann -- system administrator, physicist, optimist." },
         { "Teaching quantum computing is like teaching computer science at Hogwarts.", "Thomas Sterling, ISC2018 keynote" },
         { "It is unfortunate that the authors did not make better use of all the electric power energy that went into these massive computations.", "An anonymous referee" },
+        { "Doctor, doctor, it hurts when I hit myself in the head with the hammer! - So don't do it!", "Bjarne Stroustrup at CppCon2015" },
         { "This is extremely unlikely.", "Berk Hess" },
         { "Nothing is more anarchic than power.", "Pier Paolo Pasolini" },
     };
index 232c4c44a0d489f23ccdf8d612f448226a9b5b6b..d2f11a6b9d7bb7891e5bfa7723c7bf633b07b8da 100644 (file)
@@ -212,11 +212,6 @@ void please_cite(FILE *fp, const char *key)
           "Flexible constraints: an adiabatic treatment of quantum degrees of freedom, with application to the flexible and polarizable MCDHO model for water",
           "J. Chem. Phys.",
           116, 2002, "9602-9610" },
-        { "Hetenyi2002b",
-          "Csaba Hetenyi and David van der Spoel",
-          "Efficient docking of peptides to proteins without prior knowledge of the binding site.",
-          "Prot. Sci.",
-          11, 2002, "1729-1737" },
         { "Hess2003",
           "B. Hess and R.M. Scheek",
           "Orientation restraints in molecular dynamics simulations using time and ensemble averaging",
index ead9b6d4e511d7bf2cb01c8893e110269821199b..24606ce58a52708647b03b7b943ea142259d3889 100644 (file)
@@ -231,8 +231,6 @@ void registerLegacyModules(gmx::CommandLineModuleManager *manager)
     registerModule(manager, &gmx_xpm2ps, "xpm2ps",
                    "Convert XPM (XPixelMap) matrices to postscript or XPM");
 
-    registerModule(manager, &gmx_anadock, "anadock",
-                   "Cluster structures from Autodock runs");
     registerModule(manager, &gmx_anaeig, "anaeig",
                    "Analyze eigenvectors/normal modes");
     registerModule(manager, &gmx_analyze, "analyze",
@@ -278,8 +276,6 @@ void registerLegacyModules(gmx::CommandLineModuleManager *manager)
                    "Analyze density of states and properties based on that");
     registerModule(manager, &gmx_dyecoupl, "dyecoupl",
                    "Extract dye dynamics from trajectories");
-    registerModule(manager, &gmx_dyndom, "dyndom",
-                   "Interpolate and extrapolate structure rotations");
     registerModule(manager, &gmx_enemat, "enemat",
                    "Extract an energy matrix from an energy file");
     registerModule(manager, &gmx_energy, "energy",
@@ -304,8 +300,6 @@ void registerLegacyModules(gmx::CommandLineModuleManager *manager)
                    "Calculate residue contact maps");
     registerModule(manager, &gmx_mindist, "mindist",
                    "Calculate the minimum distance between two groups");
-    registerModule(manager, &gmx_morph, "morph",
-                   "Interpolate linearly between conformations");
     registerModule(manager, &gmx_msd, "msd",
                    "Calculates mean square displacements");
     registerModule(manager, &gmx_nmeig, "nmeig",
@@ -418,10 +412,8 @@ void registerLegacyModules(gmx::CommandLineModuleManager *manager)
             manager->addModuleGroup("Tools");
         group.addModule("analyze");
         group.addModule("awh");
-        group.addModule("dyndom");
         group.addModule("filter");
         group.addModule("lie");
-        group.addModule("morph");
         group.addModule("pme_error");
         group.addModule("sham");
         group.addModule("spatial");
@@ -475,7 +467,6 @@ void registerLegacyModules(gmx::CommandLineModuleManager *manager)
     {
         gmx::CommandLineModuleGroup group =
             manager->addModuleGroup("Structural properties");
-        group.addModule("anadock");
         group.addModule("bundle");
         group.addModule("clustsize");
         group.addModule("disre");
index 37b0c0e4620577e0e64625236635011fb27ec3d3..e56de915683660fe3e705902073c2e06c1e1f131 100644 (file)
@@ -77,6 +77,7 @@ gmx_add_gtest_executable(
     ${exename}
     # files with code for tests
     minimize.cpp
+    normalmodes.cpp
     rerun.cpp
     # pseudo-library for code for testing mdrun
     $<TARGET_OBJECTS:mdrun_test_objlib>
index 88e64ffd8db114b7ea3b6698e4ea76b3699d2e25..a3034dd19baccf5b73ab364cbf6092094c94165c 100644 (file)
@@ -47,6 +47,7 @@
 
 #include <cstdio>
 
+#include "gromacs/gmxana/gmx_ana.h"
 #include "gromacs/gmxpreprocess/grompp.h"
 #include "gromacs/hardware/detecthardware.h"
 #include "gromacs/options/basicoptions.h"
@@ -96,6 +97,7 @@ SimulationRunner::SimulationRunner(TestFileManager *fileManager) :
     tprFileName_(fileManager->getTemporaryFilePath(".tpr")),
     logFileName_(fileManager->getTemporaryFilePath(".log")),
     edrFileName_(fileManager->getTemporaryFilePath(".edr")),
+    mtxFileName_(fileManager->getTemporaryFilePath(".mtx")),
     nsteps_(-2),
     fileManager_(*fileManager)
 {
@@ -136,6 +138,14 @@ SimulationRunner::useStringAsNdxFile(const char *ndxString)
     gmx::TextWriter::writeFileFromString(ndxFileName_, ndxString);
 }
 
+void
+SimulationRunner::useTopG96AndNdxFromDatabase(const std::string &name)
+{
+    topFileName_ = gmx::test::TestFileManager::getInputFilePath(name + ".top");
+    groFileName_ = gmx::test::TestFileManager::getInputFilePath(name + ".g96");
+    ndxFileName_ = gmx::test::TestFileManager::getInputFilePath(name + ".ndx");
+}
+
 void
 SimulationRunner::useTopGroAndNdxFromDatabase(const std::string &name)
 {
@@ -208,6 +218,28 @@ SimulationRunner::callGrompp()
     return callGrompp(CommandLine());
 }
 
+int
+SimulationRunner::callNmeig()
+{
+    /* Conforming to style guide by not passing a non-const reference
+       to this function. Passing a non-const reference might make it
+       easier to write code that incorrectly re-uses callerRef after
+       the call to this function. */
+
+    CommandLine caller;
+    caller.append("nmeig");
+    caller.addOption("-s", tprFileName_);
+    caller.addOption("-f", mtxFileName_);
+    // Ignore the overall translation and rotation in the
+    // first six eigenvectors.
+    caller.addOption("-first", "7");
+    // No need to check more than a number of output values.
+    caller.addOption("-last", "50");
+    caller.addOption("-xvg", "none");
+
+    return gmx_nmeig(caller.argc(), caller.argv());
+}
+
 int
 SimulationRunner::callMdrun(const CommandLine &callerRef)
 {
@@ -223,6 +255,7 @@ SimulationRunner::callMdrun(const CommandLine &callerRef)
 
     caller.addOption("-g", logFileName_);
     caller.addOption("-e", edrFileName_);
+    caller.addOption("-mtx", mtxFileName_);
     caller.addOption("-o", fullPrecisionTrajectoryFileName_);
     caller.addOption("-x", reducedPrecisionTrajectoryFileName_);
 
index edba08915b5e30da553d19beda86fbbace90a4f3..450decba483299a58d9271bcfe69c594869ecb6e 100644 (file)
@@ -98,6 +98,8 @@ class SimulationRunner
         void useStringAsMdpFile(const std::string &mdpString);
         //! Use a string as -n input to grompp
         void useStringAsNdxFile(const char *ndxString);
+        //! Use a standard .top and .g96 file as input to grompp
+        void useTopG96AndNdxFromDatabase(const std::string &name);
         //! Use a standard .top and .gro file as input to grompp
         void useTopGroAndNdxFromDatabase(const std::string &name);
         //! Use a standard .gro file as input to grompp
@@ -110,6 +112,8 @@ class SimulationRunner
         int callGromppOnThisRank(const CommandLine &callerRef);
         //! Convenience wrapper for a default call to \c callGromppOnThisRank
         int callGromppOnThisRank();
+        //! Calls nmeig for testing
+        int callNmeig();
         //! Calls mdrun for testing with a customized command line
         int callMdrun(const CommandLine &callerRef);
         /*! \brief Convenience wrapper for calling mdrun for testing
@@ -137,6 +141,7 @@ class SimulationRunner
         std::string tprFileName_;
         std::string logFileName_;
         std::string edrFileName_;
+        std::string mtxFileName_;
         std::string cptFileName_;
         std::string swapFileName_;
         int         nsteps_;
diff --git a/src/programs/mdrun/tests/normalmodes.cpp b/src/programs/mdrun/tests/normalmodes.cpp
new file mode 100644 (file)
index 0000000..4535c62
--- /dev/null
@@ -0,0 +1,170 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Tests for the normal modes functionality.
+ *
+ * \author David van der Spoel <david.vanderspoel@icm.uu.se>
+ * \ingroup module_mdrun_integration_tests
+ */
+#include "gmxpre.h"
+
+#include <map>
+#include <memory>
+#include <string>
+#include <tuple>
+#include <vector>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/compat/make_unique.h"
+#include "gromacs/options/filenameoption.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/ifunc.h"
+#include "gromacs/trajectory/energyframe.h"
+#include "gromacs/trajectory/trajectoryframe.h"
+#include "gromacs/utility/basenetwork.h"
+#include "gromacs/utility/filestream.h"
+#include "gromacs/utility/stringutil.h"
+
+#include "testutils/mpitest.h"
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+#include "testutils/xvgtest.h"
+
+#include "energycomparison.h"
+#include "energyreader.h"
+#include "moduletest.h"
+#include "simulationdatabase.h"
+
+namespace gmx
+{
+namespace test
+{
+namespace
+{
+
+//! Helper type
+using MdpField = MdpFieldValues::value_type;
+
+/*! \brief Test fixture base for normal mode analysis
+ *
+ * This test ensures mdrun can run a normal mode analys, reaching
+ * a reproducible eigenvalues following diagonalization.
+ *
+ * The choices for tolerance are arbitrary but sufficient. */
+class NormalModesTest : public MdrunTestFixture,
+                        public ::testing::WithParamInterface <
+                        std::tuple < std::string, std::string>>
+{
+};
+
+TEST_P(NormalModesTest, WithinTolerances)
+{
+    auto params         = GetParam();
+    auto simulationName = std::get<0>(params);
+    auto integrator     = std::get<1>(params);
+    SCOPED_TRACE(formatString("Comparing normal modes for '%s'",
+                              simulationName.c_str()));
+
+    // TODO At some point we should also test PME-only ranks.
+    int numRanksAvailable = getNumberOfTestMpiRanks();
+    if (!isNumberOfPpRanksSupported(simulationName, numRanksAvailable))
+    {
+        fprintf(stdout, "Test system '%s' cannot run with %d ranks.\n"
+                "The supported numbers are: %s\n",
+                simulationName.c_str(), numRanksAvailable,
+                reportNumbersOfPpRanksSupported(simulationName).c_str());
+        return;
+    }
+    auto mdpFieldValues = prepareMdpFieldValues(simulationName.c_str(),
+                                                integrator.c_str(),
+                                                "no", "no");
+    mdpFieldValues["nsteps"]      = "1";
+    mdpFieldValues["rcoulomb"]    = "5.6";
+    mdpFieldValues["rlist"]       = "5.6";
+    mdpFieldValues["rvdw"]        = "5.6";
+    mdpFieldValues["constraints"] = "none";
+    mdpFieldValues.insert(MdpField("coulombtype", "Cut-off"));
+    mdpFieldValues.insert(MdpField("vdwtype", "Cut-off"));
+
+    // prepare the .tpr file
+    {
+        CommandLine caller;
+        caller.append("grompp");
+        runner_.useTopG96AndNdxFromDatabase(simulationName);
+        runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
+        EXPECT_EQ(0, runner_.callGrompp(caller));
+    }
+    // Do mdrun, preparing to check the normal modes later
+    {
+        CommandLine mdrunCaller;
+        mdrunCaller.append("mdrun");
+        ASSERT_EQ(0, runner_.callMdrun(mdrunCaller));
+    }
+    // Now run gmx nmeig and check the output
+    {
+        ASSERT_EQ(0, runner_.callNmeig());
+        TestReferenceData refData;
+        auto              checker  = refData.rootChecker()
+                .checkCompound("System", simulationName)
+                .checkCompound("Integrator", integrator);
+        auto          settings = XvgMatchSettings();
+        TextInputFile input("eigenval.xvg");
+        checkXvgFile(&input, &checker, settings);
+    }
+}
+
+//! Containers of systems and integrators to test.
+//! \{
+std::vector<std::string> systemsToTest_g     = { "scaled-water" };
+std::vector<std::string> integratorsToTest_g = { "nm" };
+
+//! \}
+
+// The time for OpenCL kernel compilation means these tests might time
+// out. If that proves to be a problem, these can be disabled for
+// OpenCL builds. However, once that compilation is cached for the
+// lifetime of the whole test binary process, these tests should run in
+// such configurations.
+#if GMX_DOUBLE
+INSTANTIATE_TEST_CASE_P(NormalModesWorks, NormalModesTest,
+                            ::testing::Combine(::testing::ValuesIn(systemsToTest_g),
+                                                   ::testing::ValuesIn(integratorsToTest_g)));
+#endif
+} // namespace
+} // namespace test
+} // namespace gmx
diff --git a/src/programs/mdrun/tests/refdata/NormalModesWorks_NormalModesTest_WithinTolerances_0.xml b/src/programs/mdrun/tests/refdata/NormalModesWorks_NormalModesTest_WithinTolerances_0.xml
new file mode 100644 (file)
index 0000000..11948dd
--- /dev/null
@@ -0,0 +1,74 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <System Name="scaled-water">
+    <Integrator Name="nm">
+      <XvgLegend Name="Legend">
+        <String Name="XvgLegend"><![CDATA[
+]]></String>
+      </XvgLegend>
+      <XvgData Name="Data">
+        <Sequence Name="Row0">
+          <Int Name="Length">2</Int>
+          <Real>7</Real>
+          <Real>0.068698</Real>
+        </Sequence>
+        <Sequence Name="Row1">
+          <Int Name="Length">2</Int>
+          <Real>8</Real>
+          <Real>0.311942</Real>
+        </Sequence>
+        <Sequence Name="Row2">
+          <Int Name="Length">2</Int>
+          <Real>9</Real>
+          <Real>0.318822</Real>
+        </Sequence>
+        <Sequence Name="Row3">
+          <Int Name="Length">2</Int>
+          <Real>10</Real>
+          <Real>0.916834</Real>
+        </Sequence>
+        <Sequence Name="Row4">
+          <Int Name="Length">2</Int>
+          <Real>11</Real>
+          <Real>3.24034</Real>
+        </Sequence>
+        <Sequence Name="Row5">
+          <Int Name="Length">2</Int>
+          <Real>12</Real>
+          <Real>13.6069</Real>
+        </Sequence>
+        <Sequence Name="Row6">
+          <Int Name="Length">2</Int>
+          <Real>13</Real>
+          <Real>14.2014</Real>
+        </Sequence>
+        <Sequence Name="Row7">
+          <Int Name="Length">2</Int>
+          <Real>14</Real>
+          <Real>39.3613</Real>
+        </Sequence>
+        <Sequence Name="Row8">
+          <Int Name="Length">2</Int>
+          <Real>15</Real>
+          <Real>5257.53</Real>
+        </Sequence>
+        <Sequence Name="Row9">
+          <Int Name="Length">2</Int>
+          <Real>16</Real>
+          <Real>5280.78</Real>
+        </Sequence>
+        <Sequence Name="Row10">
+          <Int Name="Length">2</Int>
+          <Real>17</Real>
+          <Real>5315.76</Real>
+        </Sequence>
+        <Sequence Name="Row11">
+          <Int Name="Length">2</Int>
+          <Real>18</Real>
+          <Real>5340.06</Real>
+        </Sequence>
+      </XvgData>
+    </Integrator>
+  </System>
+</ReferenceData>
index bc25a861399e0f1c38fbf8ffcf8817b5a218cd62..edbaa8f8a208958e3194a292e41b093f012a41e8 100644 (file)
@@ -168,6 +168,13 @@ const MdpFileValues mdpFileValueDatabase_g
                                               1, 2, 3, 4, 5, 6, 7, 8, 9
                                           } }
     },
+    // Scaled water for NMA
+    {
+        "scaled-water", { { },
+                          {
+                              1, 2, 3, 4, 5, 6
+                          } }
+    },
     // Nonanol molecule in vacuo, topology suitable for testing FEP
     // on KE, angles, dihedral restraints, coulomb and vdw
     {
@@ -229,6 +236,8 @@ MdpFieldValues prepareDefaultMdpFieldValues(const char *simulationName)
     mdpFieldValues.insert(MdpField("compressibility", "5e-5"));
     mdpFieldValues.insert(MdpField("constraints", "none"));
     mdpFieldValues.insert(MdpField("other", ""));
+    mdpFieldValues.insert(MdpField("rcoulomb", "0.7"));
+    mdpFieldValues.insert(MdpField("rvdw", "0.7"));
 
     return mdpFieldValues;
 }
@@ -289,8 +298,8 @@ prepareMdpFileContents(const MdpFieldValues &mdpFieldValues)
      * energies were not computed with those from rerun on the same
      * coordinates.
      */
-    return formatString(R"(rcoulomb                = 0.7
-                           rvdw                    = 0.7
+    return formatString(R"(rcoulomb                = %s
+                           rvdw                    = %s
                            rlist                   = -1
                            bd-fric                 = 1000
                            verlet-buffer-tolerance = 0.000001
@@ -316,6 +325,8 @@ prepareMdpFileContents(const MdpFieldValues &mdpFieldValues)
                            lincs-order             = 2
                            lincs-iter              = 5
                            %s)",
+                        mdpFieldValues.at("rcoulomb").c_str(),
+                        mdpFieldValues.at("rvdw").c_str(),
                         mdpFieldValues.at("nsteps").c_str(),
                         mdpFieldValues.at("integrator").c_str(),
                         mdpFieldValues.at("tcoupl").c_str(),
index c2cdc40520734e132782066db68231deefa73f82..6d6a3b5efeadee714ac924ac0e2a1885d3817e82 100644 (file)
@@ -55,6 +55,7 @@
 #include "gromacs/commandline/cmdlineoptionsmodule.h"
 #include "gromacs/commandline/cmdlineprogramcontext.h"
 #include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/futil.h"
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/strconvert.h"
 #include "gromacs/utility/stringstream.h"
@@ -441,6 +442,21 @@ void CommandLineTestBase::setInputFile(
     setInputFile(option, filename.c_str());
 }
 
+void CommandLineTestBase::setModifiableInputFile(
+        const char *option, const std::string &filename)
+{
+    setModifiableInputFile(option, filename.c_str());
+}
+
+void CommandLineTestBase::setModifiableInputFile(
+        const char *option, const char *filename)
+{
+    std::string originalFileName   = gmx::test::TestFileManager::getInputFilePath(filename);
+    std::string modifiableFileName = fileManager().getTemporaryFilePath(filename);
+    gmx_file_copy(originalFileName.c_str(), modifiableFileName.c_str(), true);
+    impl_->cmdline_.addOption(option, modifiableFileName);
+}
+
 void CommandLineTestBase::setInputFileContents(
         const char *option, const char *extension, const std::string &contents)
 {
@@ -470,6 +486,26 @@ void CommandLineTestBase::setOutputFile(
     impl_->helper_.setOutputFile(&impl_->cmdline_, option, filename, matcher);
 }
 
+void CommandLineTestBase::setInputAndOutputFile(
+        const char *option, const char *filename,
+        const ITextBlockMatcherSettings &matcher)
+{
+    std::string originalFileName   = gmx::test::TestFileManager::getInputFilePath(filename);
+    std::string modifiableFileName = fileManager().getTemporaryFilePath(filename);
+    gmx_file_copy(originalFileName.c_str(), modifiableFileName.c_str(), true);
+    impl_->helper_.setOutputFile(&impl_->cmdline_, option, filename, matcher);
+}
+
+void CommandLineTestBase::setInputAndOutputFile(
+        const char *option, const char *filename,
+        const IFileMatcherSettings &matcher)
+{
+    std::string originalFileName   = gmx::test::TestFileManager::getInputFilePath(filename);
+    std::string modifiableFileName = fileManager().getTemporaryFilePath(filename);
+    gmx_file_copy(originalFileName.c_str(), modifiableFileName.c_str(), true);
+    impl_->helper_.setOutputFile(&impl_->cmdline_, option, filename, matcher);
+}
+
 CommandLine &CommandLineTestBase::commandLine()
 {
     return impl_->cmdline_;
index 1065022fd40dd12b252c9c1389f72caa5332cf56..4febb0d70e645ff95f8a431d37c66067be9c7577 100644 (file)
@@ -389,6 +389,17 @@ class CommandLineTestBase : public ::testing::Test
         void setInputFile(const char *option, const char *filename);
         //! \copydoc setInputFile(const char *, const char *);
         void setInputFile(const char *option, const std::string &filename);
+        /*! \brief
+         * Sets an input file that may be modified. The file is copied to a
+         * temporary file, which is used as the test input
+         *
+         * \param[in]     option    Option to set.
+         * \param[in]     filename  Name of the input file.
+         *
+         */
+        void setModifiableInputFile(const char *option, const char *filename);
+        //! \copydoc setModifiableInputFile(const char *, const char *);
+        void setModifiableInputFile(const char *option, const std::string &filename);
         /*! \brief
          * Generates and sets an input file.
          *
@@ -419,6 +430,15 @@ class CommandLineTestBase : public ::testing::Test
          */
         void setOutputFile(const char *option, const char *filename,
                            const IFileMatcherSettings &matcher);
+        /*! \brief
+         * Sets a file parameter that is used for input and modified as output. The input file
+         * is copied to a temporary file that is used as input and can be modified.
+         */
+        void setInputAndOutputFile(const char *option, const char *filename,
+                                   const ITextBlockMatcherSettings &matcher);
+        //! \copydoc setInputAndOutputFile(const char *, const char *, const ITextBlockMatcherSettings&);
+        void setInputAndOutputFile(const char *option, const char *filename,
+                                   const IFileMatcherSettings &matcher);
 
         /*! \brief
          * Returns the internal CommandLine object used to construct the
diff --git a/src/testutils/simulationdatabase/scaled-water.g96 b/src/testutils/simulationdatabase/scaled-water.g96
new file mode 100644 (file)
index 0000000..4aee8f4
--- /dev/null
@@ -0,0 +1,14 @@
+TITLE
+2 scaled waters
+END
+POSITION
+    1 SOL   OW         1    0.914818904   17.114674936   24.974532564
+    1 SOL   HW1        2    1.207239659   16.434171182   24.367718701
+    1 SOL   HW2        3    1.794362088   17.256320895    0.323592571
+    2 SOL   OW         4    2.688663748   14.887261002   23.182636745
+    2 SOL   HW1        5    2.999853477   14.939893657   22.278872288
+    2 SOL   HW2        6    2.760097525   13.938314213   23.286543655
+END
+BOX
+   25.001000000   25.001000000   25.001000000
+END
diff --git a/src/testutils/simulationdatabase/scaled-water.ndx b/src/testutils/simulationdatabase/scaled-water.ndx
new file mode 100644 (file)
index 0000000..03ad1d8
--- /dev/null
@@ -0,0 +1,6 @@
+[ System ]
+   1    2    3    4    5    6
+[ Water ]
+   1    2    3    4    5    6
+[ SOL ]
+   1    2    3    4    5    6
diff --git a/src/testutils/simulationdatabase/scaled-water.top b/src/testutils/simulationdatabase/scaled-water.top
new file mode 100644 (file)
index 0000000..e7add3c
--- /dev/null
@@ -0,0 +1,34 @@
+[ defaults ]
+; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
+1               3               yes             0.5     0.5
+
+[ atomtypes ]
+; full atom descriptions are available in ffoplsaa.atp
+; name  bond_type    mass    charge   ptype   sigma      epsilon
+ OW   O  8      15.99940     0.000       A    3          50
+ HW   H  1       1.00800     0.000       A    0.00000e+00  0.00000e+00
+
+[ moleculetype ]
+; molname       nrexcl
+SOL             2
+
+[ atoms ]
+; id    at type res nr  residu name     at name         cg nr   charge
+1       OW     1       SOL              OW             1      -0.6
+2       HW     1       SOL             HW1             1       0.3
+3       HW     1       SOL             HW2             1       0.3
+
+[ bonds ]
+; i     j       funct   length  force.c.
+1       2       1       0.9572 5024.0
+1       3       1       0.9572 5024.0
+
+[ angles ]
+; i     j       k       funct   angle   force.c.
+2       1       3       1       100.0  6.0
+
+[ system ]
+2 scaled waters
+
+[ molecules ]
+SOL 2