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, 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/gmxomp.h"
  54 #include "gromacs/utility/smalloc.h"
  55
  56 t_vcm *init_vcm(FILE *fp, gmx_groups_t *groups, t_inputrec *ir)
  57 {
  58     t_vcm *vcm;
  59     int    g;
  60
  61     snew(vcm, 1);
  62
  63     vcm->mode = (ir->nstcomm > 0) ? ir->comm_mode : ecmNO;
  64     vcm->ndim = ndof_com(ir);
  65
  66     if (vcm->mode == ecmANGULAR && vcm->ndim < 3)
  67     {
  68         gmx_fatal(FARGS, "Can not have angular comm removal with pbc=%s",
  69                   epbc_names[ir->ePBC]);
  70     }
  71
  72     if (vcm->mode != ecmNO)
  73     {
  74         vcm->nr = groups->grps[egcVCM].nr;
  75         /* Allocate one extra for a possible rest group */
  76         vcm->size = vcm->nr + 1;
  77         /* We need vcm->nr+1 elements per thread, but to avoid cache
  78          * invalidation we add 2 elements to get a 152 byte separation.
  79          */
  80         vcm->stride = vcm->nr + 3;
  81         if (vcm->mode == ecmANGULAR)
  82         {
  83             snew(vcm->group_j, vcm->size);
  84             snew(vcm->group_x, vcm->size);
  85             snew(vcm->group_i, vcm->size);
  86             snew(vcm->group_w, vcm->size);
  87         }
  88         snew(vcm->group_p, vcm->size);
  89         snew(vcm->group_v, vcm->size);
  90         snew(vcm->group_mass, vcm->size);
  91         snew(vcm->group_name, vcm->size);
  92         snew(vcm->group_ndf, vcm->size);
  93         for (g = 0; (g < vcm->nr); g++)
  94         {
  95             vcm->group_ndf[g] = ir->opts.nrdf[g];
  96         }
  97
  98         /* Copy pointer to group names and print it. */
  99         if (fp)
 100         {
 101             fprintf(fp, "Center of mass motion removal mode is %s\n",
 102                     ECOM(vcm->mode));
 103             fprintf(fp, "We have the following groups for center of"
 104                     " mass motion removal:\n");
 105         }
 106         for (g = 0; (g < vcm->nr); g++)
 107         {
 108             vcm->group_name[g] = *groups->grpname[groups->grps[egcVCM].nm_ind[g]];
 109             if (fp)
 110             {
 111                 fprintf(fp, "%3d:  %s\n", g, vcm->group_name[g]);
 112             }
 113         }
 114
 115         snew(vcm->thread_vcm, gmx_omp_nthreads_get(emntDefault) * vcm->stride);
 116     }
 117
 118     return vcm;
 119 }
 120
 121 static void update_tensor(rvec x, real m0, tensor I)
 122 {
 123     real xy, xz, yz;
 124
 125     /* Compute inertia tensor contribution due to this atom */
 126     xy         = x[XX]*x[YY]*m0;
 127     xz         = x[XX]*x[ZZ]*m0;
 128     yz         = x[YY]*x[ZZ]*m0;
 129     I[XX][XX] += x[XX]*x[XX]*m0;
 130     I[YY][YY] += x[YY]*x[YY]*m0;
 131     I[ZZ][ZZ] += x[ZZ]*x[ZZ]*m0;
 132     I[XX][YY] += xy;
 133     I[YY][XX] += xy;
 134     I[XX][ZZ] += xz;
 135     I[ZZ][XX] += xz;
 136     I[YY][ZZ] += yz;
 137     I[ZZ][YY] += yz;
 138 }
 139
 140 /* Center of mass code for groups */
 141 void calc_vcm_grp(int start, int homenr, t_mdatoms *md,
 142                   rvec x[], rvec v[], t_vcm *vcm)
 143 {
 144     int nthreads = gmx_omp_nthreads_get(emntDefault);
 145     if (vcm->mode != ecmNO)
 146     {
 147 #pragma omp parallel num_threads(nthreads)
 148         {
 149             int t = gmx_omp_get_thread_num();
 150             for (int g = 0; g < vcm->size; g++)
 151             {
 152                 /* Reset linear momentum */
 153                 t_vcm_thread *vcm_t = &vcm->thread_vcm[t*vcm->stride + g];
 154                 vcm_t->mass = 0;
 155                 clear_rvec(vcm_t->p);
 156                 if (vcm->mode == ecmANGULAR)
 157                 {
 158                     /* Reset angular momentum */
 159                     clear_rvec(vcm_t->j);
 160                     clear_rvec(vcm_t->x);
 161                     clear_mat(vcm_t->i);
 162                 }
 163             }
 164
 165 #pragma omp for schedule(static)
 166             for (int i = start; i < start+homenr; i++)
 167             {
 168                 int  g  = 0;
 169                 real m0 = md->massT[i];
 170                 if (md->cVCM)
 171                 {
 172                     g = md->cVCM[i];
 173                 }
 174                 t_vcm_thread *vcm_t = &vcm->thread_vcm[t*vcm->stride + g];
 175                 /* Calculate linear momentum */
 176                 vcm_t->mass  += m0;
 177                 int m;
 178                 for (m = 0; (m < DIM); m++)
 179                 {
 180                     vcm_t->p[m] += m0*v[i][m];
 181                 }
 182
 183                 if (vcm->mode == ecmANGULAR)
 184                 {
 185                     /* Calculate angular momentum */
 186                     rvec   j0;
 187                     cprod(x[i], v[i], j0);
 188
 189                     for (m = 0; (m < DIM); m++)
 190                     {
 191                         vcm_t->j[m] += m0*j0[m];
 192                         vcm_t->x[m] += m0*x[i][m];
 193                     }
 194                     /* Update inertia tensor */
 195                     update_tensor(x[i], m0, vcm_t->i);
 196                 }
 197             }
 198         }
 199         for (int g = 0; g < vcm->size; g++)
 200         {
 201             /* Reset linear momentum */
 202             vcm->group_mass[g] = 0;
 203             clear_rvec(vcm->group_p[g]);
 204             if (vcm->mode == ecmANGULAR)
 205             {
 206                 /* Reset angular momentum */
 207                 clear_rvec(vcm->group_j[g]);
 208                 clear_rvec(vcm->group_x[g]);
 209                 clear_rvec(vcm->group_w[g]);
 210                 clear_mat(vcm->group_i[g]);
 211             }
 212
 213             for (int t = 0; t < nthreads; t++)
 214             {
 215                 t_vcm_thread *vcm_t = &vcm->thread_vcm[t*vcm->stride + g];
 216                 vcm->group_mass[g] += vcm_t->mass;
 217                 rvec_inc(vcm->group_p[g], vcm_t->p);
 218                 if (vcm->mode == ecmANGULAR)
 219                 {
 220                     rvec_inc(vcm->group_j[g], vcm_t->j);
 221                     rvec_inc(vcm->group_x[g], vcm_t->x);
 222                     m_add(vcm_t->i, vcm->group_i[g], vcm->group_i[g]);
 223                 }
 224             }
 225         }
 226
 227     }
 228 }
 229
 230 void do_stopcm_grp(int start, int homenr, unsigned short *group_id,
 231                    rvec x[], rvec v[], t_vcm *vcm)
 232 {
 233     if (vcm->mode != ecmNO)
 234     {
 235         // cppcheck-suppress unreadVariable
 236         int gmx_unused nth = gmx_omp_nthreads_get(emntDefault);
 237 #pragma omp parallel num_threads(nth)
 238         {
 239             int  i, g = 0;
 240             rvec dv, dx;
 241             /* Subtract linear momentum */
 242             switch (vcm->ndim)
 243             {
 244                 case 1:
 245 #pragma omp for schedule(static)
 246                     for (i = start; i < start+homenr; i++)
 247                     {
 248                         if (group_id)
 249                         {
 250                             g = group_id[i];
 251                         }
 252                         v[i][XX] -= vcm->group_v[g][XX];
 253                     }
 254                     break;
 255                 case 2:
 256 #pragma omp for schedule(static)
 257                     for (i = start; i < start+homenr; i++)
 258                     {
 259                         if (group_id)
 260                         {
 261                             g = group_id[i];
 262                         }
 263                         v[i][XX] -= vcm->group_v[g][XX];
 264                         v[i][YY] -= vcm->group_v[g][YY];
 265                     }
 266                     break;
 267                 case 3:
 268 #pragma omp for schedule(static)
 269                     for (i = start; i < start+homenr; i++)
 270                     {
 271                         if (group_id)
 272                         {
 273                             g = group_id[i];
 274                         }
 275                         rvec_dec(v[i], vcm->group_v[g]);
 276                     }
 277                     break;
 278             }
 279             if (vcm->mode == ecmANGULAR)
 280             {
 281                 /* Subtract angular momentum */
 282 #pragma omp for schedule(static)
 283                 for (i = start; i < start+homenr; i++)
 284                 {
 285                     if (group_id)
 286                     {
 287                         g = group_id[i];
 288                     }
 289                     /* Compute the correction to the velocity for each atom */
 290                     rvec_sub(x[i], vcm->group_x[g], dx);
 291                     cprod(vcm->group_w[g], dx, dv);
 292                     rvec_dec(v[i], dv);
 293                 }
 294             }
 295         }
 296     }
 297 }
 298
 299 static void get_minv(tensor A, tensor B)
 300 {
 301     int    m, n;
 302     double fac, rfac;
 303     tensor tmp;
 304
 305     tmp[XX][XX] =  A[YY][YY] + A[ZZ][ZZ];
 306     tmp[YY][XX] = -A[XX][YY];
 307     tmp[ZZ][XX] = -A[XX][ZZ];
 308     tmp[XX][YY] = -A[XX][YY];
 309     tmp[YY][YY] =  A[XX][XX] + A[ZZ][ZZ];
 310     tmp[ZZ][YY] = -A[YY][ZZ];
 311     tmp[XX][ZZ] = -A[XX][ZZ];
 312     tmp[YY][ZZ] = -A[YY][ZZ];
 313     tmp[ZZ][ZZ] =  A[XX][XX] + A[YY][YY];
 314
 315     /* This is a hack to prevent very large determinants */
 316     rfac  = (tmp[XX][XX]+tmp[YY][YY]+tmp[ZZ][ZZ])/3;
 317     if (rfac == 0.0)
 318     {
 319         gmx_fatal(FARGS, "Can not stop center of mass: maybe 2dimensional system");
 320     }
 321     fac = 1.0/rfac;
 322     for (m = 0; (m < DIM); m++)
 323     {
 324         for (n = 0; (n < DIM); n++)
 325         {
 326             tmp[m][n] *= fac;
 327         }
 328     }
 329     gmx::invertMatrix(tmp, B);
 330     for (m = 0; (m < DIM); m++)
 331     {
 332         for (n = 0; (n < DIM); n++)
 333         {
 334             B[m][n] *= fac;
 335         }
 336     }
 337 }
 338
 339 void check_cm_grp(FILE *fp, t_vcm *vcm, t_inputrec *ir, real Temp_Max)
 340 {
 341     int    m, g;
 342     real   ekcm, ekrot, tm, tm_1, Temp_cm;
 343     rvec   jcm;
 344     tensor Icm;
 345
 346     /* First analyse the total results */
 347     if (vcm->mode != ecmNO)
 348     {
 349         for (g = 0; (g < vcm->nr); g++)
 350         {
 351             tm = vcm->group_mass[g];
 352             if (tm != 0)
 353             {
 354                 tm_1 = 1.0/tm;
 355                 svmul(tm_1, vcm->group_p[g], vcm->group_v[g]);
 356             }
 357             /* Else it's zero anyway! */
 358         }
 359         if (vcm->mode == ecmANGULAR)
 360         {
 361             for (g = 0; (g < vcm->nr); g++)
 362             {
 363                 tm = vcm->group_mass[g];
 364                 if (tm != 0)
 365                 {
 366                     tm_1 = 1.0/tm;
 367
 368                     /* Compute center of mass for this group */
 369                     for (m = 0; (m < DIM); m++)
 370                     {
 371                         vcm->group_x[g][m] *= tm_1;
 372                     }
 373
 374                     /* Subtract the center of mass contribution to the
 375                      * angular momentum
 376                      */
 377                     cprod(vcm->group_x[g], vcm->group_v[g], jcm);
 378                     for (m = 0; (m < DIM); m++)
 379                     {
 380                         vcm->group_j[g][m] -= tm*jcm[m];
 381                     }
 382
 383                     /* Subtract the center of mass contribution from the inertia
 384                      * tensor (this is not as trivial as it seems, but due to
 385                      * some cancellation we can still do it, even in parallel).
 386                      */
 387                     clear_mat(Icm);
 388                     update_tensor(vcm->group_x[g], tm, Icm);
 389                     m_sub(vcm->group_i[g], Icm, vcm->group_i[g]);
 390
 391                     /* Compute angular velocity, using matrix operation
 392                      * Since J = I w
 393                      * we have
 394                      * w = I^-1 J
 395                      */
 396                     get_minv(vcm->group_i[g], Icm);
 397                     mvmul(Icm, vcm->group_j[g], vcm->group_w[g]);
 398                 }
 399                 /* Else it's zero anyway! */
 400             }
 401         }
 402     }
 403     for (g = 0; (g < vcm->nr); g++)
 404     {
 405         ekcm    = 0;
 406         if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0)
 407         {
 408             for (m = 0; m < vcm->ndim; m++)
 409             {
 410                 ekcm += gmx::square(vcm->group_v[g][m]);
 411             }
 412             ekcm   *= 0.5*vcm->group_mass[g];
 413             Temp_cm = 2*ekcm/vcm->group_ndf[g];
 414
 415             if ((Temp_cm > Temp_Max) && fp)
 416             {
 417                 fprintf(fp, "Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
 418                         vcm->group_name[g], vcm->group_v[g][XX],
 419                         vcm->group_v[g][YY], vcm->group_v[g][ZZ], Temp_cm);
 420             }
 421
 422             if (vcm->mode == ecmANGULAR)
 423             {
 424                 ekrot = 0.5*iprod(vcm->group_j[g], vcm->group_w[g]);
 425                 if ((ekrot > 1) && fp && !EI_RANDOM(ir->eI))
 426                 {
 427                     /* if we have an integrator that may not conserve momenta, skip */
 428                     tm    = vcm->group_mass[g];
 429                     fprintf(fp, "Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
 430                             vcm->group_name[g], tm, ekrot, det(vcm->group_i[g]));
 431                     fprintf(fp, "  COM: %12.5f  %12.5f  %12.5f\n",
 432                             vcm->group_x[g][XX], vcm->group_x[g][YY], vcm->group_x[g][ZZ]);
 433                     fprintf(fp, "  P:   %12.5f  %12.5f  %12.5f\n",
 434                             vcm->group_p[g][XX], vcm->group_p[g][YY], vcm->group_p[g][ZZ]);
 435                     fprintf(fp, "  V:   %12.5f  %12.5f  %12.5f\n",
 436                             vcm->group_v[g][XX], vcm->group_v[g][YY], vcm->group_v[g][ZZ]);
 437                     fprintf(fp, "  J:   %12.5f  %12.5f  %12.5f\n",
 438                             vcm->group_j[g][XX], vcm->group_j[g][YY], vcm->group_j[g][ZZ]);
 439                     fprintf(fp, "  w:   %12.5f  %12.5f  %12.5f\n",
 440                             vcm->group_w[g][XX], vcm->group_w[g][YY], vcm->group_w[g][ZZ]);
 441                     pr_rvecs(fp, 0, "Inertia tensor", vcm->group_i[g], DIM);
 442                 }
 443             }
 444         }
 445     }
 446 }