Remove velocity from partially frozen atoms in md-vv

md-vv would add some velocity to the frozen dimensions of partially
frozen atoms during constraining. This did not lead to wrong
trajectories, as the frozen dimensions of the positions are kept fixed
during propagation. The non-zero velocities were, however, reported in
trajectories and final configurations. They might also have lead to
slightly wrong kinetic energies, since the reported kinetic energy is
calculated after the velocities are constrained. All effects are
expected to be relatively small, since they did not accumulate, as the
velocities were regularly reset to zero once per step.

Refs #3849

diff --git a/docs/release-notes/2021/2021.1.rst b/docs/release-notes/2021/2021.1.rst
index 0491f63f1578ae143b1fbb08a9e0e9d81d9e40eb..77d39a22e40f25dc23a14a87e6daa171257def11 100644 (file)
--- a/docs/release-notes/2021/2021.1.rst
+++ b/docs/release-notes/2021/2021.1.rst
@@ -44,6 +44,21 @@ properly. This has been disabled in grompp.
 
 :issue:`3946`
 
+Remove velocity from partially frozen atoms in md-vv
+""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+md-vv would add some velocity to the frozen dimensions of partially
+frozen atoms during constraining. This did not lead to wrong
+trajectories, as the frozen dimensions of the positions are kept fixed
+during propagation. The non-zero velocities were, however, reported in
+trajectories and final configurations. They might also have lead to
+slightly wrong kinetic energies, since the reported kinetic energy is
+calculated after the velocities are constrained. All effects are
+expected to be relatively small, since they did not accumulate, as the
+velocities were regularly reset to zero once per step.
+
+:issue:`3849`
+
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp
index 21397f342b3dedc28d3fd458da2e1235255bca42..03d4783380b8566995027e32a774821549138270 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -697,10 +697,19 @@ bool Constraints::Impl::apply(bool                      bLog,
     }
     wallcycle_stop(wcycle, ewcCONSTR);
 
-    if (!v.empty() && cFREEZE_)
+    const bool haveVelocities = (!v.empty() || econq == ConstraintVariable::Velocities);
+    if (haveVelocities && cFREEZE_)
     {
         /* Set the velocities of frozen dimensions to zero */
-        ArrayRef<RVec> vRef = v.unpaddedArrayRef();
+        ArrayRef<RVec> vRef;
+        if (econq == ConstraintVariable::Velocities)
+        {
+            vRef = xprime.unpaddedArrayRef();
+        }
+        else
+        {
+            vRef = v.unpaddedArrayRef();
+        }
 
         int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);