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