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