src/gromacs/mdlib/vcm.c

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