2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
38 /* This file is completely threadsafe - keep it that way! */
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"
58 const char* enumValueToString(ComRemovalAlgorithm enumValue)
60 static constexpr gmx::EnumerationArray<ComRemovalAlgorithm, const char*> comRemovalAlgorithmNames = {
61 "Linear", "Angular", "None", "Linear-acceleration-correction"
63 return comRemovalAlgorithmNames[enumValue];
66 t_vcm::t_vcm(const SimulationGroups& groups, const t_inputrec& ir) :
67 integratorConservesMomentum(!EI_RANDOM(ir.eI))
69 mode = (ir.nstcomm > 0) ? ir.comm_mode : ComRemovalAlgorithm::No;
71 timeStep = ir.nstcomm * ir.delta_t;
73 if (mode == ComRemovalAlgorithm::Angular && ndim < 3)
75 gmx_fatal(FARGS, "Can not have angular comm removal with pbc=%s", c_pbcTypeNames[ir.pbcType].c_str());
78 if (mode != ComRemovalAlgorithm::No)
80 nr = groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size();
81 /* Allocate one extra for a possible rest group */
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.
87 if (mode == ComRemovalAlgorithm::Angular)
96 group_name.resize(size);
99 group_mass.resize(size);
100 group_ndf.resize(size);
101 for (int g = 0; (g < nr); g++)
103 group_ndf[g] = ir.opts.nrdf[g];
105 *groups.groupNames[groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval][g]];
108 thread_vcm.resize(gmx_omp_nthreads_get(ModuleMultiThread::Default) * stride);
111 nFreeze = ir.opts.nFreeze;
116 if (mode == ComRemovalAlgorithm::Angular)
122 void reportComRemovalInfo(FILE* fp, const t_vcm& vcm)
125 /* Copy pointer to group names and print it. */
126 if (fp && vcm.mode != ComRemovalAlgorithm::No)
128 fprintf(fp, "Center of mass motion removal mode is %s\n", enumValueToString(vcm.mode));
130 "We have the following groups for center of"
131 " mass motion removal:\n");
133 for (int g = 0; (g < vcm.nr); g++)
136 fprintf(fp, "%3d: %s\n", g, vcm.group_name[g]);
142 static void update_tensor(const rvec x, real m0, tensor I)
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;
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,
167 if (vcm->mode == ComRemovalAlgorithm::No)
171 int nthreads = gmx_omp_nthreads_get(ModuleMultiThread::Default);
174 #pragma omp parallel num_threads(nthreads) default(none) shared(x, v, vcm, md)
176 int t = gmx_omp_get_thread_num();
177 for (int g = 0; g < vcm->size; g++)
179 /* Reset linear momentum */
180 t_vcm_thread* vcm_t = &vcm->thread_vcm[t * vcm->stride + g];
182 clear_rvec(vcm_t->p);
183 if (vcm->mode == ComRemovalAlgorithm::Angular)
185 /* Reset angular momentum */
186 clear_rvec(vcm_t->j);
187 clear_rvec(vcm_t->x);
192 #pragma omp for schedule(static)
193 for (int i = 0; i < md.homenr; i++)
196 real m0 = md.massT[i];
201 t_vcm_thread* vcm_t = &vcm->thread_vcm[t * vcm->stride + g];
202 /* Calculate linear momentum */
205 for (m = 0; (m < DIM); m++)
207 vcm_t->p[m] += m0 * v[i][m];
210 if (vcm->mode == ComRemovalAlgorithm::Angular)
212 /* Calculate angular momentum */
214 cprod(x[i], v[i], j0);
216 for (m = 0; (m < DIM); m++)
218 vcm_t->j[m] += m0 * j0[m];
219 vcm_t->x[m] += m0 * x[i][m];
221 /* Update inertia tensor */
222 update_tensor(x[i], m0, vcm_t->i);
226 for (int g = 0; g < vcm->size; g++)
228 /* Reset linear momentum */
229 vcm->group_mass[g] = 0;
230 clear_rvec(vcm->group_p[g]);
231 if (vcm->mode == ComRemovalAlgorithm::Angular)
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]);
240 for (int t = 0; t < nthreads; t++)
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)
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]);
256 /*! \brief Remove the COM motion velocity from the velocities
258 * \note This routine should be called from within an OpenMP parallel region.
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
265 template<int numDimensions>
266 static void doStopComMotionLinear(const t_mdatoms& mdatoms, gmx::ArrayRef<gmx::RVec> v, const t_vcm& vcm)
268 const int homenr = mdatoms.homenr;
269 const unsigned short* group_id = mdatoms.cVCM;
271 if (mdatoms.cFREEZE != nullptr)
273 GMX_RELEASE_ASSERT(vcm.nFreeze != nullptr, "Need freeze dimension info with freeze groups");
275 #pragma omp for schedule(static)
276 for (int i = 0; i < homenr; i++)
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++)
282 if (vcm.nFreeze[freezeGroup][d] == 0)
284 v[i][d] -= vcm.group_v[vcmGroup][d];
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++)
294 for (int d = 0; d < numDimensions; d++)
296 v[i][d] -= vcm.group_v[0][d];
302 #pragma omp for schedule(static)
303 for (int i = 0; i < homenr; i++)
305 const int g = group_id[i];
306 for (int d = 0; d < numDimensions; d++)
308 v[i][d] -= vcm.group_v[g][d];
314 /*! \brief Remove the COM motion velocity from the velocities, correct the coordinates assuming constant acceleration
316 * \note This routine should be called from within an OpenMP parallel region.
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
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,
332 const real xCorrectionFactor = 0.5 * vcm.timeStep;
334 // NOLINTNEXTLINE bugprone-branch-clone This is actually a clang-tidy bug
335 if (group_id == nullptr)
337 #pragma omp for schedule(static)
338 for (int i = 0; i < homenr; i++)
340 for (int d = 0; d < numDimensions; d++)
342 x[i][d] -= vcm.group_v[0][d] * xCorrectionFactor;
343 v[i][d] -= vcm.group_v[0][d];
349 #pragma omp for schedule(static)
350 for (int i = 0; i < homenr; i++)
352 const int g = group_id[i];
353 for (int d = 0; d < numDimensions; d++)
355 x[i][d] -= vcm.group_v[g][d] * xCorrectionFactor;
356 v[i][d] -= vcm.group_v[g][d];
362 static void do_stopcm_grp(const t_mdatoms& mdatoms,
363 gmx::ArrayRef<gmx::RVec> x,
364 gmx::ArrayRef<gmx::RVec> v,
367 if (vcm.mode == ComRemovalAlgorithm::No)
372 const int homenr = mdatoms.homenr;
373 const unsigned short* group_id = mdatoms.cVCM;
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) \
381 if (vcm.mode == ComRemovalAlgorithm::Linear || vcm.mode == ComRemovalAlgorithm::Angular
382 || (vcm.mode == ComRemovalAlgorithm::LinearAccelerationCorrection && x.empty()))
384 /* Subtract linear momentum for v */
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;
394 GMX_ASSERT(vcm.mode == ComRemovalAlgorithm::LinearAccelerationCorrection,
395 "When the mode is not linear or angular, it should be acceleration "
397 /* Subtract linear momentum for v and x*/
401 doStopComMotionAccelerationCorrection<1>(homenr, group_id, x, v, vcm);
404 doStopComMotionAccelerationCorrection<2>(homenr, group_id, x, v, vcm);
407 doStopComMotionAccelerationCorrection<3>(homenr, group_id, x, v, vcm);
411 if (vcm.mode == ComRemovalAlgorithm::Angular)
413 /* Subtract angular momentum */
414 GMX_ASSERT(!x.empty(), "Need x to compute angular momentum correction");
417 #pragma omp for schedule(static)
418 for (int i = 0; i < homenr; i++)
424 /* Compute the correction to the velocity for each atom */
426 rvec_sub(x[i], vcm.group_x[g], dx);
427 cprod(vcm.group_w[g], dx, dv);
435 static void get_minv(tensor A, tensor B)
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];
451 /* This is a hack to prevent very large determinants */
452 rfac = (tmp[XX][XX] + tmp[YY][YY] + tmp[ZZ][ZZ]) / 3;
455 gmx_fatal(FARGS, "Can not stop center of mass: maybe 2dimensional system");
458 for (m = 0; (m < DIM); m++)
460 for (n = 0; (n < DIM); n++)
465 gmx::invertMatrix(tmp, B);
466 for (m = 0; (m < DIM); m++)
468 for (n = 0; (n < DIM); n++)
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)
479 real ekcm, ekrot, tm, tm_1, Temp_cm;
483 /* First analyse the total results */
484 if (vcm->mode != ComRemovalAlgorithm::No)
486 for (g = 0; (g < vcm->nr); g++)
488 tm = vcm->group_mass[g];
492 svmul(tm_1, vcm->group_p[g], vcm->group_v[g]);
494 /* Else it's zero anyway! */
496 if (vcm->mode == ComRemovalAlgorithm::Angular)
498 for (g = 0; (g < vcm->nr); g++)
500 tm = vcm->group_mass[g];
505 /* Compute center of mass for this group */
506 for (m = 0; (m < DIM); m++)
508 vcm->group_x[g][m] *= tm_1;
511 /* Subtract the center of mass contribution to the
514 cprod(vcm->group_x[g], vcm->group_v[g], jcm);
515 for (m = 0; (m < DIM); m++)
517 vcm->group_j[g][m] -= tm * jcm[m];
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).
525 update_tensor(vcm->group_x[g], tm, Icm);
526 m_sub(vcm->group_i[g], Icm, vcm->group_i[g]);
528 /* Compute angular velocity, using matrix operation
533 get_minv(vcm->group_i[g], Icm);
534 mvmul(Icm, vcm->group_j[g], vcm->group_w[g]);
536 /* Else it's zero anyway! */
540 for (g = 0; (g < vcm->nr); g++)
543 if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0)
545 for (m = 0; m < vcm->ndim; m++)
547 ekcm += gmx::square(vcm->group_v[g][m]);
549 ekcm *= 0.5 * vcm->group_mass[g];
550 Temp_cm = 2 * ekcm / vcm->group_ndf[g];
552 if ((Temp_cm > Temp_Max) && fp)
555 "Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
563 if (vcm->mode == ComRemovalAlgorithm::Angular)
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)
569 /* if we have an integrator that may not conserve momenta, skip */
570 tm = vcm->group_mass[g];
572 "Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
576 det(vcm->group_i[g]));
578 " COM: %12.5f %12.5f %12.5f\n",
581 vcm->group_x[g][ZZ]);
583 " P: %12.5f %12.5f %12.5f\n",
586 vcm->group_p[g][ZZ]);
588 " V: %12.5f %12.5f %12.5f\n",
591 vcm->group_v[g][ZZ]);
593 " J: %12.5f %12.5f %12.5f\n",
596 vcm->group_j[g][ZZ]);
598 " w: %12.5f %12.5f %12.5f\n",
601 vcm->group_w[g][ZZ]);
602 pr_rvecs(fp, 0, "Inertia tensor", vcm->group_i[g], DIM);
609 void process_and_stopcm_grp(FILE* fplog,
611 const t_mdatoms& mdatoms,
612 gmx::ArrayRef<gmx::RVec> x,
613 gmx::ArrayRef<gmx::RVec> v)
615 if (vcm->mode != ComRemovalAlgorithm::No)
617 // TODO: Replace fixed temperature of 1 by a system value
618 process_and_check_cm_grp(fplog, vcm, 1);
620 do_stopcm_grp(mdatoms, x, v, *vcm);