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