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