src/gromacs/selection/centerofmass.cpp

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