src/gromacs/mdlib/vcm.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #include "gmxpre.h"
  40
  41 #include "vcm.h"
  42
  43 #include "gromacs/math/functions.h"
  44 #include "gromacs/math/invertmatrix.h"
  45 #include "gromacs/math/vec.h"
  46 #include "gromacs/math/vecdump.h"
  47 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  48 #include "gromacs/mdtypes/inputrec.h"
  49 #include "gromacs/mdtypes/md_enums.h"
  50 #include "gromacs/mdtypes/mdatom.h"
  51 #include "gromacs/pbcutil/pbc.h"
  52 #include "gromacs/topology/topology.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/gmxassert.h"
  55 #include "gromacs/utility/gmxomp.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 t_vcm::t_vcm(const SimulationGroups& groups, const t_inputrec& ir) :
  59     integratorConservesMomentum(!EI_RANDOM(ir.eI))
  60 {
  61     mode     = (ir.nstcomm > 0) ? ir.comm_mode : ComRemovalAlgorithm::No;
  62     ndim     = ndof_com(&ir);
  63     timeStep = ir.nstcomm * ir.delta_t;
  64
  65     if (mode == ComRemovalAlgorithm::Angular && ndim < 3)
  66     {
  67         gmx_fatal(FARGS, "Can not have angular comm removal with pbc=%s", c_pbcTypeNames[ir.pbcType].c_str());
  68     }
  69
  70     if (mode != ComRemovalAlgorithm::No)
  71     {
  72         nr = groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size();
  73         /* Allocate one extra for a possible rest group */
  74         size = nr + 1;
  75         /* We need vcm->nr+1 elements per thread, but to avoid cache
  76          * invalidation we add 2 elements to get a 152 byte separation.
  77          */
  78         stride = nr + 3;
  79         if (mode == ComRemovalAlgorithm::Angular)
  80         {
  81             snew(group_i, size);
  82
  83             group_j.resize(size);
  84             group_x.resize(size);
  85             group_w.resize(size);
  86         }
  87
  88         group_name.resize(size);
  89         group_p.resize(size);
  90         group_v.resize(size);
  91         group_mass.resize(size);
  92         group_ndf.resize(size);
  93         for (int g = 0; (g < nr); g++)
  94         {
  95             group_ndf[g] = ir.opts.nrdf[g];
  96             group_name[g] =
  97                     *groups.groupNames[groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval][g]];
  98         }
  99
 100         thread_vcm.resize(gmx_omp_nthreads_get(emntDefault) * stride);
 101     }
 102
 103     nFreeze = ir.opts.nFreeze;
 104 }
 105
 106 t_vcm::~t_vcm()
 107 {
 108     if (mode == ComRemovalAlgorithm::Angular)
 109     {
 110         sfree(group_i);
 111     }
 112 }
 113
 114 void reportComRemovalInfo(FILE* fp, const t_vcm& vcm)
 115 {
 116
 117     /* Copy pointer to group names and print it. */
 118     if (fp && vcm.mode != ComRemovalAlgorithm::No)
 119     {
 120         fprintf(fp, "Center of mass motion removal mode is %s\n", enumValueToString(vcm.mode));
 121         fprintf(fp,
 122                 "We have the following groups for center of"
 123                 " mass motion removal:\n");
 124
 125         for (int g = 0; (g < vcm.nr); g++)
 126         {
 127
 128             fprintf(fp, "%3d:  %s\n", g, vcm.group_name[g]);
 129         }
 130     }
 131 }
 132
 133
 134 static void update_tensor(const rvec x, real m0, tensor I)
 135 {
 136     real xy, xz, yz;
 137
 138     /* Compute inertia tensor contribution due to this atom */
 139     xy = x[XX] * x[YY] * m0;
 140     xz = x[XX] * x[ZZ] * m0;
 141     yz = x[YY] * x[ZZ] * m0;
 142     I[XX][XX] += x[XX] * x[XX] * m0;
 143     I[YY][YY] += x[YY] * x[YY] * m0;
 144     I[ZZ][ZZ] += x[ZZ] * x[ZZ] * m0;
 145     I[XX][YY] += xy;
 146     I[YY][XX] += xy;
 147     I[XX][ZZ] += xz;
 148     I[ZZ][XX] += xz;
 149     I[YY][ZZ] += yz;
 150     I[ZZ][YY] += yz;
 151 }
 152
 153 /* Center of mass code for groups */
 154 void calc_vcm_grp(const t_mdatoms&               md,
 155                   gmx::ArrayRef<const gmx::RVec> x,
 156                   gmx::ArrayRef<const gmx::RVec> v,
 157                   t_vcm*                         vcm)
 158 {
 159     if (vcm->mode == ComRemovalAlgorithm::No)
 160     {
 161         return;
 162     }
 163     int nthreads = gmx_omp_nthreads_get(emntDefault);
 164
 165     {
 166 #pragma omp parallel num_threads(nthreads) default(none) shared(x, v, vcm, md)
 167         {
 168             int t = gmx_omp_get_thread_num();
 169             for (int g = 0; g < vcm->size; g++)
 170             {
 171                 /* Reset linear momentum */
 172                 t_vcm_thread* vcm_t = &vcm->thread_vcm[t * vcm->stride + g];
 173                 vcm_t->mass         = 0;
 174                 clear_rvec(vcm_t->p);
 175                 if (vcm->mode == ComRemovalAlgorithm::Angular)
 176                 {
 177                     /* Reset angular momentum */
 178                     clear_rvec(vcm_t->j);
 179                     clear_rvec(vcm_t->x);
 180                     clear_mat(vcm_t->i);
 181                 }
 182             }
 183
 184 #pragma omp for schedule(static)
 185             for (int i = 0; i < md.homenr; i++)
 186             {
 187                 int  g  = 0;
 188                 real m0 = md.massT[i];
 189                 if (md.cVCM)
 190                 {
 191                     g = md.cVCM[i];
 192                 }
 193                 t_vcm_thread* vcm_t = &vcm->thread_vcm[t * vcm->stride + g];
 194                 /* Calculate linear momentum */
 195                 vcm_t->mass += m0;
 196                 int m;
 197                 for (m = 0; (m < DIM); m++)
 198                 {
 199                     vcm_t->p[m] += m0 * v[i][m];
 200                 }
 201
 202                 if (vcm->mode == ComRemovalAlgorithm::Angular)
 203                 {
 204                     /* Calculate angular momentum */
 205                     rvec j0;
 206                     cprod(x[i], v[i], j0);
 207
 208                     for (m = 0; (m < DIM); m++)
 209                     {
 210                         vcm_t->j[m] += m0 * j0[m];
 211                         vcm_t->x[m] += m0 * x[i][m];
 212                     }
 213                     /* Update inertia tensor */
 214                     update_tensor(x[i], m0, vcm_t->i);
 215                 }
 216             }
 217         }
 218         for (int g = 0; g < vcm->size; g++)
 219         {
 220             /* Reset linear momentum */
 221             vcm->group_mass[g] = 0;
 222             clear_rvec(vcm->group_p[g]);
 223             if (vcm->mode == ComRemovalAlgorithm::Angular)
 224             {
 225                 /* Reset angular momentum */
 226                 clear_rvec(vcm->group_j[g]);
 227                 clear_rvec(vcm->group_x[g]);
 228                 clear_rvec(vcm->group_w[g]);
 229                 clear_mat(vcm->group_i[g]);
 230             }
 231
 232             for (int t = 0; t < nthreads; t++)
 233             {
 234                 t_vcm_thread* vcm_t = &vcm->thread_vcm[t * vcm->stride + g];
 235                 vcm->group_mass[g] += vcm_t->mass;
 236                 rvec_inc(vcm->group_p[g], vcm_t->p);
 237                 if (vcm->mode == ComRemovalAlgorithm::Angular)
 238                 {
 239                     rvec_inc(vcm->group_j[g], vcm_t->j);
 240                     rvec_inc(vcm->group_x[g], vcm_t->x);
 241                     m_add(vcm_t->i, vcm->group_i[g], vcm->group_i[g]);
 242                 }
 243             }
 244         }
 245     }
 246 }
 247
 248 /*! \brief Remove the COM motion velocity from the velocities
 249  *
 250  * \note This routine should be called from within an OpenMP parallel region.
 251  *
 252  * \tparam numDimensions    Correct dimensions 0 to \p numDimensions-1
 253  * \param[in]     mdatoms   The atom property and group information
 254  * \param[in,out] v         The velocities to correct
 255  * \param[in]     vcm       VCM data
 256  */
 257 template<int numDimensions>
 258 static void doStopComMotionLinear(const t_mdatoms& mdatoms, gmx::ArrayRef<gmx::RVec> v, const t_vcm& vcm)
 259 {
 260     const int             homenr   = mdatoms.homenr;
 261     const unsigned short* group_id = mdatoms.cVCM;
 262
 263     if (mdatoms.cFREEZE != nullptr)
 264     {
 265         GMX_RELEASE_ASSERT(vcm.nFreeze != nullptr, "Need freeze dimension info with freeze groups");
 266
 267 #pragma omp for schedule(static)
 268         for (int i = 0; i < homenr; i++)
 269         {
 270             unsigned short vcmGroup    = (group_id == nullptr ? 0 : group_id[i]);
 271             unsigned short freezeGroup = mdatoms.cFREEZE[i];
 272             for (int d = 0; d < numDimensions; d++)
 273             {
 274                 if (vcm.nFreeze[freezeGroup][d] == 0)
 275                 {
 276                     v[i][d] -= vcm.group_v[vcmGroup][d];
 277                 }
 278             }
 279         }
 280     }
 281     else if (group_id == nullptr)
 282     { // NOLINT bugprone-branch-clone This is actually a clang-tidy bug
 283 #pragma omp for schedule(static)
 284         for (int i = 0; i < homenr; i++)
 285         {
 286             for (int d = 0; d < numDimensions; d++)
 287             {
 288                 v[i][d] -= vcm.group_v[0][d];
 289             }
 290         }
 291     }
 292     else
 293     {
 294 #pragma omp for schedule(static)
 295         for (int i = 0; i < homenr; i++)
 296         {
 297             const int g = group_id[i];
 298             for (int d = 0; d < numDimensions; d++)
 299             {
 300                 v[i][d] -= vcm.group_v[g][d];
 301             }
 302         }
 303     }
 304 }
 305
 306 /*! \brief Remove the COM motion velocity from the velocities, correct the coordinates assuming constant acceleration
 307  *
 308  * \note This routine should be called from within an OpenMP parallel region.
 309  *
 310  * \tparam numDimensions    Correct dimensions 0 to \p numDimensions-1
 311  * \param[in]     homenr    The number of atoms to correct
 312  * \param[in]     group_id  List of VCM group ids, when nullptr is passed all atoms are assumed to be in group 0
 313  * \param[in,out] x         The coordinates to correct
 314  * \param[in,out] v         The velocities to correct
 315  * \param[in]     vcm       VCM data
 316  */
 317 template<int numDimensions>
 318 static void doStopComMotionAccelerationCorrection(int                      homenr,
 319                                                   const unsigned short*    group_id,
 320                                                   gmx::ArrayRef<gmx::RVec> x,
 321                                                   gmx::ArrayRef<gmx::RVec> v,
 322                                                   const t_vcm&             vcm)
 323 {
 324     const real xCorrectionFactor = 0.5 * vcm.timeStep;
 325
 326     // NOLINTNEXTLINE bugprone-branch-clone This is actually a clang-tidy bug
 327     if (group_id == nullptr)
 328     {
 329 #pragma omp for schedule(static)
 330         for (int i = 0; i < homenr; i++)
 331         {
 332             for (int d = 0; d < numDimensions; d++)
 333             {
 334                 x[i][d] -= vcm.group_v[0][d] * xCorrectionFactor;
 335                 v[i][d] -= vcm.group_v[0][d];
 336             }
 337         }
 338     }
 339     else
 340     {
 341 #pragma omp for schedule(static)
 342         for (int i = 0; i < homenr; i++)
 343         {
 344             const int g = group_id[i];
 345             for (int d = 0; d < numDimensions; d++)
 346             {
 347                 x[i][d] -= vcm.group_v[g][d] * xCorrectionFactor;
 348                 v[i][d] -= vcm.group_v[g][d];
 349             }
 350         }
 351     }
 352 }
 353
 354 static void do_stopcm_grp(const t_mdatoms&         mdatoms,
 355                           gmx::ArrayRef<gmx::RVec> x,
 356                           gmx::ArrayRef<gmx::RVec> v,
 357                           const t_vcm&             vcm)
 358 {
 359     if (vcm.mode == ComRemovalAlgorithm::No)
 360     {
 361         return;
 362     }
 363     {
 364         const int             homenr   = mdatoms.homenr;
 365         const unsigned short* group_id = mdatoms.cVCM;
 366
 367         int gmx_unused nth = gmx_omp_nthreads_get(emntDefault);
 368         // homenr could be shared, but gcc-8 & gcc-9 don't agree how to write that...
 369         // https://www.gnu.org/software/gcc/gcc-9/porting_to.html -> OpenMP data sharing
 370 #pragma omp parallel num_threads(nth) default(none) shared(x, v, vcm, group_id, mdatoms) \
 371         firstprivate(homenr)
 372         {
 373             if (vcm.mode == ComRemovalAlgorithm::Linear || vcm.mode == ComRemovalAlgorithm::Angular
 374                 || (vcm.mode == ComRemovalAlgorithm::LinearAccelerationCorrection && x.empty()))
 375             {
 376                 /* Subtract linear momentum for v */
 377                 switch (vcm.ndim)
 378                 {
 379                     case 1: doStopComMotionLinear<1>(mdatoms, v, vcm); break;
 380                     case 2: doStopComMotionLinear<2>(mdatoms, v, vcm); break;
 381                     case 3: doStopComMotionLinear<3>(mdatoms, v, vcm); break;
 382                 }
 383             }
 384             else
 385             {
 386                 GMX_ASSERT(vcm.mode == ComRemovalAlgorithm::LinearAccelerationCorrection,
 387                            "When the mode is not linear or angular, it should be acceleration "
 388                            "correction");
 389                 /* Subtract linear momentum for v and x*/
 390                 switch (vcm.ndim)
 391                 {
 392                     case 1:
 393                         doStopComMotionAccelerationCorrection<1>(homenr, group_id, x, v, vcm);
 394                         break;
 395                     case 2:
 396                         doStopComMotionAccelerationCorrection<2>(homenr, group_id, x, v, vcm);
 397                         break;
 398                     case 3:
 399                         doStopComMotionAccelerationCorrection<3>(homenr, group_id, x, v, vcm);
 400                         break;
 401                 }
 402             }
 403             if (vcm.mode == ComRemovalAlgorithm::Angular)
 404             {
 405                 /* Subtract angular momentum */
 406                 GMX_ASSERT(!x.empty(), "Need x to compute angular momentum correction");
 407
 408                 int g = 0;
 409 #pragma omp for schedule(static)
 410                 for (int i = 0; i < homenr; i++)
 411                 {
 412                     if (group_id)
 413                     {
 414                         g = group_id[i];
 415                     }
 416                     /* Compute the correction to the velocity for each atom */
 417                     rvec dv, dx;
 418                     rvec_sub(x[i], vcm.group_x[g], dx);
 419                     cprod(vcm.group_w[g], dx, dv);
 420                     rvec_dec(v[i], dv);
 421                 }
 422             }
 423         }
 424     }
 425 }
 426
 427 static void get_minv(tensor A, tensor B)
 428 {
 429     int    m, n;
 430     double fac, rfac;
 431     tensor tmp;
 432
 433     tmp[XX][XX] = A[YY][YY] + A[ZZ][ZZ];
 434     tmp[YY][XX] = -A[XX][YY];
 435     tmp[ZZ][XX] = -A[XX][ZZ];
 436     tmp[XX][YY] = -A[XX][YY];
 437     tmp[YY][YY] = A[XX][XX] + A[ZZ][ZZ];
 438     tmp[ZZ][YY] = -A[YY][ZZ];
 439     tmp[XX][ZZ] = -A[XX][ZZ];
 440     tmp[YY][ZZ] = -A[YY][ZZ];
 441     tmp[ZZ][ZZ] = A[XX][XX] + A[YY][YY];
 442
 443     /* This is a hack to prevent very large determinants */
 444     rfac = (tmp[XX][XX] + tmp[YY][YY] + tmp[ZZ][ZZ]) / 3;
 445     if (rfac == 0.0)
 446     {
 447         gmx_fatal(FARGS, "Can not stop center of mass: maybe 2dimensional system");
 448     }
 449     fac = 1.0 / rfac;
 450     for (m = 0; (m < DIM); m++)
 451     {
 452         for (n = 0; (n < DIM); n++)
 453         {
 454             tmp[m][n] *= fac;
 455         }
 456     }
 457     gmx::invertMatrix(tmp, B);
 458     for (m = 0; (m < DIM); m++)
 459     {
 460         for (n = 0; (n < DIM); n++)
 461         {
 462             B[m][n] *= fac;
 463         }
 464     }
 465 }
 466
 467 /* Processes VCM after reduction over ranks and prints warning with high VMC and fp != nullptr */
 468 static void process_and_check_cm_grp(FILE* fp, t_vcm* vcm, real Temp_Max)
 469 {
 470     int    m, g;
 471     real   ekcm, ekrot, tm, tm_1, Temp_cm;
 472     rvec   jcm;
 473     tensor Icm;
 474
 475     /* First analyse the total results */
 476     if (vcm->mode != ComRemovalAlgorithm::No)
 477     {
 478         for (g = 0; (g < vcm->nr); g++)
 479         {
 480             tm = vcm->group_mass[g];
 481             if (tm != 0)
 482             {
 483                 tm_1 = 1.0 / tm;
 484                 svmul(tm_1, vcm->group_p[g], vcm->group_v[g]);
 485             }
 486             /* Else it's zero anyway! */
 487         }
 488         if (vcm->mode == ComRemovalAlgorithm::Angular)
 489         {
 490             for (g = 0; (g < vcm->nr); g++)
 491             {
 492                 tm = vcm->group_mass[g];
 493                 if (tm != 0)
 494                 {
 495                     tm_1 = 1.0 / tm;
 496
 497                     /* Compute center of mass for this group */
 498                     for (m = 0; (m < DIM); m++)
 499                     {
 500                         vcm->group_x[g][m] *= tm_1;
 501                     }
 502
 503                     /* Subtract the center of mass contribution to the
 504                      * angular momentum
 505                      */
 506                     cprod(vcm->group_x[g], vcm->group_v[g], jcm);
 507                     for (m = 0; (m < DIM); m++)
 508                     {
 509                         vcm->group_j[g][m] -= tm * jcm[m];
 510                     }
 511
 512                     /* Subtract the center of mass contribution from the inertia
 513                      * tensor (this is not as trivial as it seems, but due to
 514                      * some cancellation we can still do it, even in parallel).
 515                      */
 516                     clear_mat(Icm);
 517                     update_tensor(vcm->group_x[g], tm, Icm);
 518                     m_sub(vcm->group_i[g], Icm, vcm->group_i[g]);
 519
 520                     /* Compute angular velocity, using matrix operation
 521                      * Since J = I w
 522                      * we have
 523                      * w = I^-1 J
 524                      */
 525                     get_minv(vcm->group_i[g], Icm);
 526                     mvmul(Icm, vcm->group_j[g], vcm->group_w[g]);
 527                 }
 528                 /* Else it's zero anyway! */
 529             }
 530         }
 531     }
 532     for (g = 0; (g < vcm->nr); g++)
 533     {
 534         ekcm = 0;
 535         if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0)
 536         {
 537             for (m = 0; m < vcm->ndim; m++)
 538             {
 539                 ekcm += gmx::square(vcm->group_v[g][m]);
 540             }
 541             ekcm *= 0.5 * vcm->group_mass[g];
 542             Temp_cm = 2 * ekcm / vcm->group_ndf[g];
 543
 544             if ((Temp_cm > Temp_Max) && fp)
 545             {
 546                 fprintf(fp,
 547                         "Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
 548                         vcm->group_name[g],
 549                         vcm->group_v[g][XX],
 550                         vcm->group_v[g][YY],
 551                         vcm->group_v[g][ZZ],
 552                         Temp_cm);
 553             }
 554
 555             if (vcm->mode == ComRemovalAlgorithm::Angular)
 556             {
 557                 ekrot = 0.5 * iprod(vcm->group_j[g], vcm->group_w[g]);
 558                 // TODO: Change absolute energy comparison to relative
 559                 if ((ekrot > 1) && fp && vcm->integratorConservesMomentum)
 560                 {
 561                     /* if we have an integrator that may not conserve momenta, skip */
 562                     tm = vcm->group_mass[g];
 563                     fprintf(fp,
 564                             "Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
 565                             vcm->group_name[g],
 566                             tm,
 567                             ekrot,
 568                             det(vcm->group_i[g]));
 569                     fprintf(fp,
 570                             "  COM: %12.5f  %12.5f  %12.5f\n",
 571                             vcm->group_x[g][XX],
 572                             vcm->group_x[g][YY],
 573                             vcm->group_x[g][ZZ]);
 574                     fprintf(fp,
 575                             "  P:   %12.5f  %12.5f  %12.5f\n",
 576                             vcm->group_p[g][XX],
 577                             vcm->group_p[g][YY],
 578                             vcm->group_p[g][ZZ]);
 579                     fprintf(fp,
 580                             "  V:   %12.5f  %12.5f  %12.5f\n",
 581                             vcm->group_v[g][XX],
 582                             vcm->group_v[g][YY],
 583                             vcm->group_v[g][ZZ]);
 584                     fprintf(fp,
 585                             "  J:   %12.5f  %12.5f  %12.5f\n",
 586                             vcm->group_j[g][XX],
 587                             vcm->group_j[g][YY],
 588                             vcm->group_j[g][ZZ]);
 589                     fprintf(fp,
 590                             "  w:   %12.5f  %12.5f  %12.5f\n",
 591                             vcm->group_w[g][XX],
 592                             vcm->group_w[g][YY],
 593                             vcm->group_w[g][ZZ]);
 594                     pr_rvecs(fp, 0, "Inertia tensor", vcm->group_i[g], DIM);
 595                 }
 596             }
 597         }
 598     }
 599 }
 600
 601 void process_and_stopcm_grp(FILE*                    fplog,
 602                             t_vcm*                   vcm,
 603                             const t_mdatoms&         mdatoms,
 604                             gmx::ArrayRef<gmx::RVec> x,
 605                             gmx::ArrayRef<gmx::RVec> v)
 606 {
 607     if (vcm->mode != ComRemovalAlgorithm::No)
 608     {
 609         // TODO: Replace fixed temperature of 1 by a system value
 610         process_and_check_cm_grp(fplog, vcm, 1);
 611
 612         do_stopcm_grp(mdatoms, x, v, *vcm);
 613     }
 614 }