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 {
  59     mode     = (ir.nstcomm > 0) ? ir.comm_mode : ecmNO;
  60     ndim     = ndof_com(&ir);
  61     timeStep = ir.nstcomm*ir.delta_t;
  62
  63     if (mode == ecmANGULAR && ndim < 3)
  64     {
  65         gmx_fatal(FARGS, "Can not have angular comm removal with pbc=%s",
  66                   epbc_names[ir.ePBC]);
  67     }
  68
  69     if (mode != ecmNO)
  70     {
  71         nr = groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].nr;
  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
  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] = *groups.groupNames[groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].nm_ind[g]];
  97
  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 == ecmANGULAR)
 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 != ecmNO)
 119     {
 120         fprintf(fp, "Center of mass motion removal mode is %s\n",
 121                 ECOM(vcm.mode));
 122         fprintf(fp, "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(int start, int homenr, t_mdatoms *md,
 155                   rvec x[], rvec v[], t_vcm *vcm)
 156 {
 157     int nthreads = gmx_omp_nthreads_get(emntDefault);
 158     if (vcm->mode != ecmNO)
 159     {
 160 #pragma omp parallel num_threads(nthreads)
 161         {
 162             int t = gmx_omp_get_thread_num();
 163             for (int g = 0; g < vcm->size; g++)
 164             {
 165                 /* Reset linear momentum */
 166                 t_vcm_thread *vcm_t = &vcm->thread_vcm[t*vcm->stride + g];
 167                 vcm_t->mass = 0;
 168                 clear_rvec(vcm_t->p);
 169                 if (vcm->mode == ecmANGULAR)
 170                 {
 171                     /* Reset angular momentum */
 172                     clear_rvec(vcm_t->j);
 173                     clear_rvec(vcm_t->x);
 174                     clear_mat(vcm_t->i);
 175                 }
 176             }
 177
 178 #pragma omp for schedule(static)
 179             for (int i = start; i < start+homenr; i++)
 180             {
 181                 int  g  = 0;
 182                 real m0 = md->massT[i];
 183                 if (md->cVCM)
 184                 {
 185                     g = md->cVCM[i];
 186                 }
 187                 t_vcm_thread *vcm_t = &vcm->thread_vcm[t*vcm->stride + g];
 188                 /* Calculate linear momentum */
 189                 vcm_t->mass  += m0;
 190                 int m;
 191                 for (m = 0; (m < DIM); m++)
 192                 {
 193                     vcm_t->p[m] += m0*v[i][m];
 194                 }
 195
 196                 if (vcm->mode == ecmANGULAR)
 197                 {
 198                     /* Calculate angular momentum */
 199                     rvec   j0;
 200                     cprod(x[i], v[i], j0);
 201
 202                     for (m = 0; (m < DIM); m++)
 203                     {
 204                         vcm_t->j[m] += m0*j0[m];
 205                         vcm_t->x[m] += m0*x[i][m];
 206                     }
 207                     /* Update inertia tensor */
 208                     update_tensor(x[i], m0, vcm_t->i);
 209                 }
 210             }
 211         }
 212         for (int g = 0; g < vcm->size; g++)
 213         {
 214             /* Reset linear momentum */
 215             vcm->group_mass[g] = 0;
 216             clear_rvec(vcm->group_p[g]);
 217             if (vcm->mode == ecmANGULAR)
 218             {
 219                 /* Reset angular momentum */
 220                 clear_rvec(vcm->group_j[g]);
 221                 clear_rvec(vcm->group_x[g]);
 222                 clear_rvec(vcm->group_w[g]);
 223                 clear_mat(vcm->group_i[g]);
 224             }
 225
 226             for (int t = 0; t < nthreads; t++)
 227             {
 228                 t_vcm_thread *vcm_t = &vcm->thread_vcm[t*vcm->stride + g];
 229                 vcm->group_mass[g] += vcm_t->mass;
 230                 rvec_inc(vcm->group_p[g], vcm_t->p);
 231                 if (vcm->mode == ecmANGULAR)
 232                 {
 233                     rvec_inc(vcm->group_j[g], vcm_t->j);
 234                     rvec_inc(vcm->group_x[g], vcm_t->x);
 235                     m_add(vcm_t->i, vcm->group_i[g], vcm->group_i[g]);
 236                 }
 237             }
 238         }
 239
 240     }
 241 }
 242
 243 /*! \brief Remove the COM motion velocity from the velocities
 244  *
 245  * \note This routine should be called from within an OpenMP parallel region.
 246  *
 247  * \tparam numDimensions    Correct dimensions 0 to \p numDimensions-1
 248  * \param[in]     mdatoms   The atom property and group information
 249  * \param[in,out] v         The velocities to correct
 250  * \param[in]     vcm       VCM data
 251  */
 252 template<int numDimensions>
 253 static void
 254 doStopComMotionLinear(const t_mdatoms      &mdatoms,
 255                       rvec                 *v,
 256                       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
 317 doStopComMotionAccelerationCorrection(int                   homenr,
 318                                       const unsigned short *group_id,
 319                                       rvec * gmx_restrict   x,
 320                                       rvec * gmx_restrict   v,
 321                                       const t_vcm          &vcm)
 322 {
 323     const real xCorrectionFactor = 0.5*vcm.timeStep;
 324
 325     if (group_id == nullptr)
 326     {
 327 #pragma omp for schedule(static)
 328         for (int i = 0; i < homenr; i++)
 329         {
 330             for (int d = 0; d < numDimensions; d++)
 331             {
 332                 x[i][d] -= vcm.group_v[0][d]*xCorrectionFactor;
 333                 v[i][d] -= vcm.group_v[0][d];
 334             }
 335         }
 336     }
 337     else
 338     {
 339 #pragma omp for schedule(static)
 340         for (int i = 0; i < homenr; i++)
 341         {
 342             const int g = group_id[i];
 343             for (int d = 0; d < numDimensions; d++)
 344             {
 345                 x[i][d] -= vcm.group_v[g][d]*xCorrectionFactor;
 346                 v[i][d] -= vcm.group_v[g][d];
 347             }
 348         }
 349     }
 350 }
 351
 352 void do_stopcm_grp(const t_mdatoms &mdatoms,
 353                    rvec x[], rvec v[], const t_vcm &vcm)
 354 {
 355     if (vcm.mode != ecmNO)
 356     {
 357         const int             homenr   = mdatoms.homenr;
 358         const unsigned short *group_id = mdatoms.cVCM;
 359
 360         int gmx_unused        nth = gmx_omp_nthreads_get(emntDefault);
 361 #pragma omp parallel num_threads(nth)
 362         {
 363             if (vcm.mode == ecmLINEAR ||
 364                 vcm.mode == ecmANGULAR ||
 365                 (vcm.mode == ecmLINEAR_ACCELERATION_CORRECTION && x == nullptr))
 366             {
 367                 /* Subtract linear momentum for v */
 368                 switch (vcm.ndim)
 369                 {
 370                     case 1:
 371                         doStopComMotionLinear<1>(mdatoms, v, vcm);
 372                         break;
 373                     case 2:
 374                         doStopComMotionLinear<2>(mdatoms, v, vcm);
 375                         break;
 376                     case 3:
 377                         doStopComMotionLinear<3>(mdatoms, v, vcm);
 378                         break;
 379                 }
 380             }
 381             else
 382             {
 383                 GMX_ASSERT(vcm.mode == ecmLINEAR_ACCELERATION_CORRECTION, "When the mode is not linear or angular, it should be acceleration correction");
 384                 /* Subtract linear momentum for v and x*/
 385                 switch (vcm.ndim)
 386                 {
 387                     case 1:
 388                         doStopComMotionAccelerationCorrection<1>(homenr, group_id, x, v, vcm);
 389                         break;
 390                     case 2:
 391                         doStopComMotionAccelerationCorrection<2>(homenr, group_id, x, v, vcm);
 392                         break;
 393                     case 3:
 394                         doStopComMotionAccelerationCorrection<3>(homenr, group_id, x, v, vcm);
 395                         break;
 396                 }
 397
 398             }
 399             if (vcm.mode == ecmANGULAR)
 400             {
 401                 /* Subtract angular momentum */
 402                 GMX_ASSERT(x, "Need x to compute angular momentum correction");
 403
 404                 int g = 0;
 405 #pragma omp for schedule(static)
 406                 for (int i = 0; i < homenr; i++)
 407                 {
 408                     if (group_id)
 409                     {
 410                         g = group_id[i];
 411                     }
 412                     /* Compute the correction to the velocity for each atom */
 413                     rvec dv, dx;
 414                     rvec_sub(x[i], vcm.group_x[g], dx);
 415                     cprod(vcm.group_w[g], dx, dv);
 416                     rvec_dec(v[i], dv);
 417                 }
 418             }
 419         }
 420     }
 421 }
 422
 423 static void get_minv(tensor A, tensor B)
 424 {
 425     int    m, n;
 426     double fac, rfac;
 427     tensor tmp;
 428
 429     tmp[XX][XX] =  A[YY][YY] + A[ZZ][ZZ];
 430     tmp[YY][XX] = -A[XX][YY];
 431     tmp[ZZ][XX] = -A[XX][ZZ];
 432     tmp[XX][YY] = -A[XX][YY];
 433     tmp[YY][YY] =  A[XX][XX] + A[ZZ][ZZ];
 434     tmp[ZZ][YY] = -A[YY][ZZ];
 435     tmp[XX][ZZ] = -A[XX][ZZ];
 436     tmp[YY][ZZ] = -A[YY][ZZ];
 437     tmp[ZZ][ZZ] =  A[XX][XX] + A[YY][YY];
 438
 439     /* This is a hack to prevent very large determinants */
 440     rfac  = (tmp[XX][XX]+tmp[YY][YY]+tmp[ZZ][ZZ])/3;
 441     if (rfac == 0.0)
 442     {
 443         gmx_fatal(FARGS, "Can not stop center of mass: maybe 2dimensional system");
 444     }
 445     fac = 1.0/rfac;
 446     for (m = 0; (m < DIM); m++)
 447     {
 448         for (n = 0; (n < DIM); n++)
 449         {
 450             tmp[m][n] *= fac;
 451         }
 452     }
 453     gmx::invertMatrix(tmp, B);
 454     for (m = 0; (m < DIM); m++)
 455     {
 456         for (n = 0; (n < DIM); n++)
 457         {
 458             B[m][n] *= fac;
 459         }
 460     }
 461 }
 462
 463 void check_cm_grp(FILE *fp, t_vcm *vcm, t_inputrec *ir, real Temp_Max)
 464 {
 465     int    m, g;
 466     real   ekcm, ekrot, tm, tm_1, Temp_cm;
 467     rvec   jcm;
 468     tensor Icm;
 469
 470     /* First analyse the total results */
 471     if (vcm->mode != ecmNO)
 472     {
 473         for (g = 0; (g < vcm->nr); g++)
 474         {
 475             tm = vcm->group_mass[g];
 476             if (tm != 0)
 477             {
 478                 tm_1 = 1.0/tm;
 479                 svmul(tm_1, vcm->group_p[g], vcm->group_v[g]);
 480             }
 481             /* Else it's zero anyway! */
 482         }
 483         if (vcm->mode == ecmANGULAR)
 484         {
 485             for (g = 0; (g < vcm->nr); g++)
 486             {
 487                 tm = vcm->group_mass[g];
 488                 if (tm != 0)
 489                 {
 490                     tm_1 = 1.0/tm;
 491
 492                     /* Compute center of mass for this group */
 493                     for (m = 0; (m < DIM); m++)
 494                     {
 495                         vcm->group_x[g][m] *= tm_1;
 496                     }
 497
 498                     /* Subtract the center of mass contribution to the
 499                      * angular momentum
 500                      */
 501                     cprod(vcm->group_x[g], vcm->group_v[g], jcm);
 502                     for (m = 0; (m < DIM); m++)
 503                     {
 504                         vcm->group_j[g][m] -= tm*jcm[m];
 505                     }
 506
 507                     /* Subtract the center of mass contribution from the inertia
 508                      * tensor (this is not as trivial as it seems, but due to
 509                      * some cancellation we can still do it, even in parallel).
 510                      */
 511                     clear_mat(Icm);
 512                     update_tensor(vcm->group_x[g], tm, Icm);
 513                     m_sub(vcm->group_i[g], Icm, vcm->group_i[g]);
 514
 515                     /* Compute angular velocity, using matrix operation
 516                      * Since J = I w
 517                      * we have
 518                      * w = I^-1 J
 519                      */
 520                     get_minv(vcm->group_i[g], Icm);
 521                     mvmul(Icm, vcm->group_j[g], vcm->group_w[g]);
 522                 }
 523                 /* Else it's zero anyway! */
 524             }
 525         }
 526     }
 527     for (g = 0; (g < vcm->nr); g++)
 528     {
 529         ekcm    = 0;
 530         if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0)
 531         {
 532             for (m = 0; m < vcm->ndim; m++)
 533             {
 534                 ekcm += gmx::square(vcm->group_v[g][m]);
 535             }
 536             ekcm   *= 0.5*vcm->group_mass[g];
 537             Temp_cm = 2*ekcm/vcm->group_ndf[g];
 538
 539             if ((Temp_cm > Temp_Max) && fp)
 540             {
 541                 fprintf(fp, "Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
 542                         vcm->group_name[g], vcm->group_v[g][XX],
 543                         vcm->group_v[g][YY], vcm->group_v[g][ZZ], Temp_cm);
 544             }
 545
 546             if (vcm->mode == ecmANGULAR)
 547             {
 548                 ekrot = 0.5*iprod(vcm->group_j[g], vcm->group_w[g]);
 549                 if ((ekrot > 1) && fp && !EI_RANDOM(ir->eI))
 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",
 556                             vcm->group_x[g][XX], vcm->group_x[g][YY], vcm->group_x[g][ZZ]);
 557                     fprintf(fp, "  P:   %12.5f  %12.5f  %12.5f\n",
 558                             vcm->group_p[g][XX], vcm->group_p[g][YY], vcm->group_p[g][ZZ]);
 559                     fprintf(fp, "  V:   %12.5f  %12.5f  %12.5f\n",
 560                             vcm->group_v[g][XX], vcm->group_v[g][YY], vcm->group_v[g][ZZ]);
 561                     fprintf(fp, "  J:   %12.5f  %12.5f  %12.5f\n",
 562                             vcm->group_j[g][XX], vcm->group_j[g][YY], vcm->group_j[g][ZZ]);
 563                     fprintf(fp, "  w:   %12.5f  %12.5f  %12.5f\n",
 564                             vcm->group_w[g][XX], 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 }