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