src/gromacs/selection/centerofmass.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009-2018, The GROMACS development team.
   5  * Copyright (c) 2019,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements functions in centerofmass.h.
  39  *
  40  * \author Teemu Murtola <teemu.murtola@gmail.com>
  41  * \ingroup module_selection
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "centerofmass.h"
  46
  47 #include <cmath>
  48
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/topology/block.h"
  52 #include "gromacs/topology/mtop_lookup.h"
  53 #include "gromacs/topology/topology.h"
  54 #include "gromacs/utility/gmxassert.h"
  55
  56 void gmx_calc_cog(const gmx_mtop_t* /* top */, rvec x[], int nrefat, const int index[], rvec xout)
  57 {
  58     int m, ai;
  59
  60     clear_rvec(xout);
  61     for (m = 0; m < nrefat; ++m)
  62     {
  63         ai = index[m];
  64         rvec_inc(xout, x[ai]);
  65     }
  66     svmul(1.0 / nrefat, xout, xout);
  67 }
  68
  69 /*!
  70  * \param[in]  top    Topology structure with masses.
  71  * \param[in]  x      Position vectors of all atoms.
  72  * \param[in]  nrefat Number of atoms in the index.
  73  * \param[in]  index  Indices of atoms.
  74  * \param[out] xout   COM position for the indexed atoms.
  75  *
  76  * Works exactly as gmx_calc_cog() with the exception that a center of
  77  * mass are calculated, and hence a topology with masses is required.
  78  */
  79 void gmx_calc_com(const gmx_mtop_t* top, rvec x[], int nrefat, const int index[], rvec xout)
  80 {
  81     GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
  82                        "No masses available while mass weighting was requested");
  83     clear_rvec(xout);
  84     real mtot = 0;
  85     int  molb = 0;
  86     for (int m = 0; m < nrefat; ++m)
  87     {
  88         const int  ai   = index[m];
  89         const real mass = mtopGetAtomMass(*top, ai, &molb);
  90         for (int j = 0; j < DIM; ++j)
  91         {
  92             xout[j] += mass * x[ai][j];
  93         }
  94         mtot += mass;
  95     }
  96     svmul(1.0 / mtot, xout, xout);
  97 }
  98
  99 /*!
 100  * \param[in]  top    Topology structure with masses.
 101  * \param[in]  f      Forces on all atoms.
 102  * \param[in]  nrefat Number of atoms in the index.
 103  * \param[in]  index  Indices of atoms.
 104  * \param[out] fout   Force on the COG position for the indexed atoms.
 105  */
 106 void gmx_calc_cog_f(const gmx_mtop_t* top, rvec f[], int nrefat, const int index[], rvec fout)
 107 {
 108     GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
 109                        "No masses available while mass weighting was requested");
 110     clear_rvec(fout);
 111     real mtot = 0;
 112     int  molb = 0;
 113     for (int m = 0; m < nrefat; ++m)
 114     {
 115         const int  ai   = index[m];
 116         const real mass = mtopGetAtomMass(*top, ai, &molb);
 117         for (int j = 0; j < DIM; ++j)
 118         {
 119             fout[j] += f[ai][j] / mass;
 120         }
 121         mtot += mass;
 122     }
 123     svmul(mtot / nrefat, fout, fout);
 124 }
 125
 126 void gmx_calc_com_f(const gmx_mtop_t* /* top */, rvec f[], int nrefat, const int index[], rvec fout)
 127 {
 128     clear_rvec(fout);
 129     for (int m = 0; m < nrefat; ++m)
 130     {
 131         const int ai = index[m];
 132         rvec_inc(fout, f[ai]);
 133     }
 134 }
 135
 136 /*!
 137  * \param[in]  top   Topology structure with masses
 138  *   (can be NULL if \p bMASS==false).
 139  * \param[in]  x     Position vectors of all atoms.
 140  * \param[in]  nrefat Number of atoms in the index.
 141  * \param[in]  index Indices of atoms.
 142  * \param[in]  bMass If true, mass weighting is used.
 143  * \param[out] xout  COM/COG position for the indexed atoms.
 144  *
 145  * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
 146  * \p bMass.
 147  * Other parameters are passed unmodified to these functions.
 148  */
 149 void gmx_calc_comg(const gmx_mtop_t* top, rvec x[], int nrefat, const int index[], bool bMass, rvec xout)
 150 {
 151     if (bMass)
 152     {
 153         gmx_calc_com(top, x, nrefat, index, xout);
 154     }
 155     else
 156     {
 157         gmx_calc_cog(top, x, nrefat, index, xout);
 158     }
 159 }
 160
 161 /*!
 162  * \param[in]  top   Topology structure with masses
 163  *   (can be NULL if \p bMASS==true).
 164  * \param[in]  f     Forces on all atoms.
 165  * \param[in]  nrefat Number of atoms in the index.
 166  * \param[in]  index Indices of atoms.
 167  * \param[in]  bMass If true, force on COM is calculated.
 168  * \param[out] fout  Force on the COM/COG position for the indexed atoms.
 169  *
 170  * Calls either gmx_calc_cog_f() or gmx_calc_com_f() depending on the value of
 171  * \p bMass.
 172  * Other parameters are passed unmodified to these functions.
 173  */
 174 void gmx_calc_comg_f(const gmx_mtop_t* top, rvec f[], int nrefat, const int index[], bool bMass, rvec fout)
 175 {
 176     if (bMass)
 177     {
 178         gmx_calc_com_f(top, f, nrefat, index, fout);
 179     }
 180     else
 181     {
 182         gmx_calc_cog_f(top, f, nrefat, index, fout);
 183     }
 184 }
 185
 186
 187 /*!
 188  * \param[in]  top    Topology structure (unused, can be NULL).
 189  * \param[in]  x      Position vectors of all atoms.
 190  * \param[in]  pbc    Periodic boundary conditions structure.
 191  * \param[in]  nrefat Number of atoms in the index.
 192  * \param[in]  index  Indices of atoms.
 193  * \param[out] xout   COG position for the indexed atoms.
 194  *
 195  * Works exactly as gmx_calc_com_pbc(), but calculates the center of geometry.
 196  */
 197 void gmx_calc_cog_pbc(const gmx_mtop_t* top, rvec x[], const t_pbc* pbc, int nrefat, const int index[], rvec xout)
 198 {
 199     const real tol = 1e-4;
 200     bool       bChanged;
 201     int        m, j, ai, iter;
 202     rvec       dx, xtest;
 203
 204     /* First simple calculation */
 205     gmx_calc_cog(top, x, nrefat, index, xout);
 206     /* Now check if any atom is more than half the box from the COM */
 207     if (pbc)
 208     {
 209         iter = 0;
 210         do
 211         {
 212             bChanged = false;
 213             for (m = 0; m < nrefat; ++m)
 214             {
 215                 ai = index[m];
 216                 pbc_dx(pbc, x[ai], xout, dx);
 217                 rvec_add(xout, dx, xtest);
 218                 for (j = 0; j < DIM; ++j)
 219                 {
 220                     if (std::fabs(xtest[j] - x[ai][j]) > tol)
 221                     {
 222                         /* Here we have used the wrong image for contributing to the COM */
 223                         xout[j] += (xtest[j] - x[ai][j]) / nrefat;
 224                         x[ai][j] = xtest[j];
 225                         bChanged = true;
 226                     }
 227                 }
 228             }
 229             iter++;
 230         } while (bChanged);
 231     }
 232 }
 233
 234 /*!
 235  * \param[in]  top    Topology structure with masses.
 236  * \param[in]  x      Position vectors of all atoms.
 237  * \param[in]  pbc    Periodic boundary conditions structure.
 238  * \param[in]  nrefat Number of atoms in the index.
 239  * \param[in]  index  Indices of atoms.
 240  * \param[out] xout   COM position for the indexed atoms.
 241  *
 242  * Works as gmx_calc_com(), but takes into account periodic boundary
 243  * conditions: If any atom is more than half the box from the COM,
 244  * it is wrapped around and a new COM is calculated. This is repeated
 245  * until no atoms violate the condition.
 246  *
 247  * Modified from src/tools/gmx_sorient.c in Gromacs distribution.
 248  */
 249 void gmx_calc_com_pbc(const gmx_mtop_t* top, rvec x[], const t_pbc* pbc, int nrefat, const int index[], rvec xout)
 250 {
 251     GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
 252                        "No masses available while mass weighting was requested");
 253     /* First simple calculation */
 254     clear_rvec(xout);
 255     real mtot = 0;
 256     int  molb = 0;
 257     for (int m = 0; m < nrefat; ++m)
 258     {
 259         const int  ai   = index[m];
 260         const real mass = mtopGetAtomMass(*top, ai, &molb);
 261         for (int j = 0; j < DIM; ++j)
 262         {
 263             xout[j] += mass * x[ai][j];
 264         }
 265         mtot += mass;
 266     }
 267     svmul(1.0 / mtot, xout, xout);
 268     /* Now check if any atom is more than half the box from the COM */
 269     if (pbc)
 270     {
 271         const real tol = 1e-4;
 272         bool       bChanged;
 273         do
 274         {
 275             bChanged = false;
 276             molb     = 0;
 277             for (int m = 0; m < nrefat; ++m)
 278             {
 279                 rvec       dx, xtest;
 280                 const int  ai   = index[m];
 281                 const real mass = mtopGetAtomMass(*top, ai, &molb) / mtot;
 282                 pbc_dx(pbc, x[ai], xout, dx);
 283                 rvec_add(xout, dx, xtest);
 284                 for (int j = 0; j < DIM; ++j)
 285                 {
 286                     if (std::fabs(xtest[j] - x[ai][j]) > tol)
 287                     {
 288                         /* Here we have used the wrong image for contributing to the COM */
 289                         xout[j] += mass * (xtest[j] - x[ai][j]);
 290                         x[ai][j] = xtest[j];
 291                         bChanged = true;
 292                     }
 293                 }
 294             }
 295         } while (bChanged);
 296     }
 297 }
 298
 299 /*!
 300  * \param[in]  top   Topology structure with masses
 301  *   (can be NULL if \p bMASS==false).
 302  * \param[in]  x     Position vectors of all atoms.
 303  * \param[in]  pbc    Periodic boundary conditions structure.
 304  * \param[in]  nrefat Number of atoms in the index.
 305  * \param[in]  index Indices of atoms.
 306  * \param[in]  bMass If true, mass weighting is used.
 307  * \param[out] xout  COM/COG position for the indexed atoms.
 308  *
 309  * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
 310  * \p bMass.
 311  * Other parameters are passed unmodified to these functions.
 312  */
 313 void gmx_calc_comg_pbc(const gmx_mtop_t* top, rvec x[], const t_pbc* pbc, int nrefat, const int index[], bool bMass, rvec xout)
 314 {
 315     if (bMass)
 316     {
 317         gmx_calc_com_pbc(top, x, pbc, nrefat, index, xout);
 318     }
 319     else
 320     {
 321         gmx_calc_cog_pbc(top, x, pbc, nrefat, index, xout);
 322     }
 323 }
 324
 325
 326 void gmx_calc_cog_block(const gmx_mtop_t* /* top */, rvec x[], const t_block* block, const int index[], rvec xout[])
 327 {
 328     int  b, i, ai;
 329     rvec xb;
 330
 331     for (b = 0; b < block->nr; ++b)
 332     {
 333         clear_rvec(xb);
 334         for (i = block->index[b]; i < block->index[b + 1]; ++i)
 335         {
 336             ai = index[i];
 337             rvec_inc(xb, x[ai]);
 338         }
 339         svmul(1.0 / (block->index[b + 1] - block->index[b]), xb, xout[b]);
 340     }
 341 }
 342
 343 /*!
 344  * \param[in]  top   Topology structure with masses.
 345  * \param[in]  x     Position vectors of all atoms.
 346  * \param[in]  block t_block structure that divides \p index into blocks.
 347  * \param[in]  index Indices of atoms.
 348  * \param[out] xout  \p block->nr COM positions.
 349  *
 350  * Works exactly as gmx_calc_cog_block() with the exception that centers of
 351  * mass are calculated, and hence a topology with masses is required.
 352  */
 353 void gmx_calc_com_block(const gmx_mtop_t* top, rvec x[], const t_block* block, const int index[], rvec xout[])
 354 {
 355     GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
 356                        "No masses available while mass weighting was requested");
 357     int molb = 0;
 358     for (int b = 0; b < block->nr; ++b)
 359     {
 360         rvec xb;
 361         clear_rvec(xb);
 362         real mtot = 0;
 363         for (int i = block->index[b]; i < block->index[b + 1]; ++i)
 364         {
 365             const int  ai   = index[i];
 366             const real mass = mtopGetAtomMass(*top, ai, &molb);
 367             for (int d = 0; d < DIM; ++d)
 368             {
 369                 xb[d] += mass * x[ai][d];
 370             }
 371             mtot += mass;
 372         }
 373         svmul(1.0 / mtot, xb, xout[b]);
 374     }
 375 }
 376
 377 /*!
 378  * \param[in]  top   Topology structure with masses.
 379  * \param[in]  f     Forces on all atoms.
 380  * \param[in]  block t_block structure that divides \p index into blocks.
 381  * \param[in]  index Indices of atoms.
 382  * \param[out] fout  \p block->nr Forces on COG positions.
 383  */
 384 void gmx_calc_cog_f_block(const gmx_mtop_t* top, rvec f[], const t_block* block, const int index[], rvec fout[])
 385 {
 386     GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
 387                        "No masses available while mass weighting was requested");
 388     int molb = 0;
 389     for (int b = 0; b < block->nr; ++b)
 390     {
 391         rvec fb;
 392         clear_rvec(fb);
 393         real mtot = 0;
 394         for (int i = block->index[b]; i < block->index[b + 1]; ++i)
 395         {
 396             const int  ai   = index[i];
 397             const real mass = mtopGetAtomMass(*top, ai, &molb);
 398             for (int d = 0; d < DIM; ++d)
 399             {
 400                 fb[d] += f[ai][d] / mass;
 401             }
 402             mtot += mass;
 403         }
 404         svmul(mtot / (block->index[b + 1] - block->index[b]), fb, fout[b]);
 405     }
 406 }
 407
 408 void gmx_calc_com_f_block(const gmx_mtop_t* /* top */,
 409                           rvec           f[],
 410                           const t_block* block,
 411                           const int      index[],
 412                           rvec           fout[])
 413 {
 414     for (int b = 0; b < block->nr; ++b)
 415     {
 416         rvec fb;
 417         clear_rvec(fb);
 418         for (int i = block->index[b]; i < block->index[b + 1]; ++i)
 419         {
 420             const int ai = index[i];
 421             rvec_inc(fb, f[ai]);
 422         }
 423         copy_rvec(fb, fout[b]);
 424     }
 425 }
 426
 427 /*!
 428  * \param[in]  top   Topology structure with masses
 429  *   (can be NULL if \p bMASS==false).
 430  * \param[in]  x     Position vectors of all atoms.
 431  * \param[in]  block t_block structure that divides \p index into blocks.
 432  * \param[in]  index Indices of atoms.
 433  * \param[in]  bMass If true, mass weighting is used.
 434  * \param[out] xout  \p block->nr COM/COG positions.
 435  *
 436  * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
 437  * value of \p bMass.
 438  * Other parameters are passed unmodified to these functions.
 439  */
 440 void gmx_calc_comg_block(const gmx_mtop_t* top,
 441                          rvec              x[],
 442                          const t_block*    block,
 443                          const int         index[],
 444                          bool              bMass,
 445                          rvec              xout[])
 446 {
 447     if (bMass)
 448     {
 449         gmx_calc_com_block(top, x, block, index, xout);
 450     }
 451     else
 452     {
 453         gmx_calc_cog_block(top, x, block, index, xout);
 454     }
 455 }
 456
 457 /*!
 458  * \param[in]  top   Topology structure with masses
 459  *   (can be NULL if \p bMASS==false).
 460  * \param[in]  f     Forces on all atoms.
 461  * \param[in]  block t_block structure that divides \p index into blocks.
 462  * \param[in]  index Indices of atoms.
 463  * \param[in]  bMass If true, force on COM is calculated.
 464  * \param[out] fout  \p block->nr forces on the COM/COG positions.
 465  *
 466  * Calls either gmx_calc_com_f_block() or gmx_calc_cog_f_block() depending on
 467  * the value of \p bMass.
 468  * Other parameters are passed unmodified to these functions.
 469  */
 470 void gmx_calc_comg_f_block(const gmx_mtop_t* top,
 471                            rvec              f[],
 472                            const t_block*    block,
 473                            const int         index[],
 474                            bool              bMass,
 475                            rvec              fout[])
 476 {
 477     if (bMass)
 478     {
 479         gmx_calc_com_f_block(top, f, block, index, fout);
 480     }
 481     else
 482     {
 483         gmx_calc_cog_f_block(top, f, block, index, fout);
 484     }
 485 }
 486
 487 /*!
 488  * \param[in]  top   Topology structure with masses
 489  *   (can be NULL if \p bMASS==false).
 490  * \param[in]  x     Position vectors of all atoms.
 491  * \param[in]  block Blocks for calculation.
 492  * \param[in]  bMass If true, mass weighting is used.
 493  * \param[out] xout  \p block->nr COM/COG positions.
 494  *
 495  * Calls gmx_calc_comg_block(), converting the \p t_blocka structure into
 496  * a \p t_block and an index. Other parameters are passed unmodified.
 497  *
 498  * \attention
 499  * This function assumes that a pointer to \c t_blocka can be safely typecast
 500  * into \c t_block such that the index fields can still be referenced.
 501  * With the present Gromacs defitions of these types, this is the case,
 502  * but if the layout of these structures is changed, this may lead to strange
 503  * crashes.
 504  */
 505 void gmx_calc_comg_blocka(const gmx_mtop_t* top, rvec x[], const t_blocka* block, bool bMass, rvec xout[])
 506 {
 507     /* TODO: It would probably be better to do this without the type cast */
 508     gmx_calc_comg_block(top, x, reinterpret_cast<const t_block*>(block), block->a, bMass, xout);
 509 }
 510
 511 /*!
 512  * \param[in]  top   Topology structure with masses
 513  *   (can be NULL if \p bMASS==true).
 514  * \param[in]  f     Forces on all atoms.
 515  * \param[in]  block Blocks for calculation.
 516  * \param[in]  bMass If true, force on COM is calculated.
 517  * \param[out] fout  \p block->nr forces on the COM/COG positions.
 518  *
 519  * Calls gmx_calc_comg_f_block(), converting the \p t_blocka structure into
 520  * a \p t_block and an index. Other parameters are passed unmodified.
 521  *
 522  * \attention
 523  * This function assumes that a pointer to \c t_blocka can be safely typecast
 524  * into \c t_block such that the index fields can still be referenced.
 525  * With the present Gromacs defitions of these types, this is the case,
 526  * but if the layout of these structures is changed, this may lead to strange
 527  * crashes.
 528  */
 529 void gmx_calc_comg_f_blocka(const gmx_mtop_t* top, rvec f[], const t_blocka* block, bool bMass, rvec fout[])
 530 {
 531     /* TODO: It would probably be better to do this without the type cast */
 532     gmx_calc_comg_f_block(top, f, reinterpret_cast<const t_block*>(block), block->a, bMass, fout);
 533 }