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,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source 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 "gmxpre.h"
  43
  44 #include "gromacs/selection/centerofmass.h"
  45
  46 #include <errno.h>
  47
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/pbcutil/pbc.h"
  50 #include "gromacs/topology/block.h"
  51 #include "gromacs/topology/topology.h"
  52
  53 int
  54 gmx_calc_cog(t_topology * /* top */, rvec x[], int nrefat, atom_id index[], rvec xout)
  55 {
  56     int                 m, ai;
  57
  58     clear_rvec(xout);
  59     for (m = 0; m < nrefat; ++m)
  60     {
  61         ai = index[m];
  62         rvec_inc(xout, x[ai]);
  63     }
  64     svmul(1.0/nrefat, xout, xout);
  65     return 0;
  66 }
  67
  68 /*!
  69  * \param[in]  top    Topology structure with masses.
  70  * \param[in]  x      Position vectors of all atoms.
  71  * \param[in]  nrefat Number of atoms in the index.
  72  * \param[in]  index  Indices of atoms.
  73  * \param[out] xout   COM position for the indexed atoms.
  74  * \returns    0 on success, EINVAL if \p top is NULL.
  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 int
  80 gmx_calc_com(t_topology *top, rvec x[], int nrefat, atom_id index[], rvec xout)
  81 {
  82     int                 m, j, ai;
  83     real                mass, mtot;
  84
  85     if (!top)
  86     {
  87         gmx_incons("no masses available while mass weighting was requested");
  88         return EINVAL;
  89     }
  90     clear_rvec(xout);
  91     mtot = 0;
  92     for (m = 0; m < nrefat; ++m)
  93     {
  94         ai   = index[m];
  95         mass = top->atoms.atom[ai].m;
  96         for (j = 0; j < DIM; ++j)
  97         {
  98             xout[j] += mass * x[ai][j];
  99         }
 100         mtot += mass;
 101     }
 102     svmul(1.0/mtot, xout, xout);
 103     return 0;
 104 }
 105
 106 /*!
 107  * \param[in]  top    Topology structure with masses.
 108  * \param[in]  f      Forces on all atoms.
 109  * \param[in]  nrefat Number of atoms in the index.
 110  * \param[in]  index  Indices of atoms.
 111  * \param[out] fout   Force on the COG position for the indexed atoms.
 112  * \returns    0 on success, EINVAL if \p top is NULL.
 113  */
 114 int
 115 gmx_calc_cog_f(t_topology *top, rvec f[], int nrefat, atom_id index[], rvec fout)
 116 {
 117     int                 m, j, ai;
 118     real                mass, mtot;
 119
 120     if (!top)
 121     {
 122         gmx_incons("no masses available while mass weighting was needed");
 123         return EINVAL;
 124     }
 125     clear_rvec(fout);
 126     mtot = 0;
 127     for (m = 0; m < nrefat; ++m)
 128     {
 129         ai   = index[m];
 130         mass = top->atoms.atom[ai].m;
 131         for (j = 0; j < DIM; ++j)
 132         {
 133             fout[j] += f[ai][j] / mass;
 134         }
 135         mtot += mass;
 136     }
 137     svmul(mtot / nrefat, fout, fout);
 138     return 0;
 139 }
 140
 141 int
 142 gmx_calc_com_f(t_topology * /* top */, rvec f[], int nrefat, atom_id index[], rvec fout)
 143 {
 144     clear_rvec(fout);
 145     for (int m = 0; m < nrefat; ++m)
 146     {
 147         const int ai = index[m];
 148         rvec_inc(fout, f[ai]);
 149     }
 150     return 0;
 151 }
 152
 153 /*!
 154  * \param[in]  top   Topology structure with masses
 155  *   (can be NULL if \p bMASS==false).
 156  * \param[in]  x     Position vectors of all atoms.
 157  * \param[in]  nrefat Number of atoms in the index.
 158  * \param[in]  index Indices of atoms.
 159  * \param[in]  bMass If true, mass weighting is used.
 160  * \param[out] xout  COM/COG position for the indexed atoms.
 161  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 162  *
 163  * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
 164  * \p bMass.
 165  * Other parameters are passed unmodified to these functions.
 166  */
 167 int
 168 gmx_calc_comg(t_topology *top, rvec x[], int nrefat, atom_id index[],
 169               bool bMass, rvec xout)
 170 {
 171     if (bMass)
 172     {
 173         return gmx_calc_com(top, x, nrefat, index, xout);
 174     }
 175     else
 176     {
 177         return gmx_calc_cog(top, x, nrefat, index, xout);
 178     }
 179 }
 180
 181 /*!
 182  * \param[in]  top   Topology structure with masses
 183  *   (can be NULL if \p bMASS==true).
 184  * \param[in]  f     Forces on all atoms.
 185  * \param[in]  nrefat Number of atoms in the index.
 186  * \param[in]  index Indices of atoms.
 187  * \param[in]  bMass If true, force on COM is calculated.
 188  * \param[out] fout  Force on the COM/COG position for the indexed atoms.
 189  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is false.
 190  *
 191  * Calls either gmx_calc_cog_f() or gmx_calc_com_f() depending on the value of
 192  * \p bMass.
 193  * Other parameters are passed unmodified to these functions.
 194  */
 195 int
 196 gmx_calc_comg_f(t_topology *top, rvec f[], int nrefat, atom_id index[],
 197                 bool bMass, rvec fout)
 198 {
 199     if (bMass)
 200     {
 201         return gmx_calc_com_f(top, f, nrefat, index, fout);
 202     }
 203     else
 204     {
 205         return gmx_calc_cog_f(top, f, nrefat, index, fout);
 206     }
 207 }
 208
 209
 210 /*!
 211  * \param[in]  top    Topology structure (unused, can be NULL).
 212  * \param[in]  x      Position vectors of all atoms.
 213  * \param[in]  pbc    Periodic boundary conditions structure.
 214  * \param[in]  nrefat Number of atoms in the index.
 215  * \param[in]  index  Indices of atoms.
 216  * \param[out] xout   COG position for the indexed atoms.
 217  * \returns    0 on success.
 218  *
 219  * Works exactly as gmx_calc_com_pbc(), but calculates the center of geometry.
 220  */
 221 int
 222 gmx_calc_cog_pbc(t_topology *top, rvec x[], t_pbc *pbc,
 223                  int nrefat, atom_id index[], rvec xout)
 224 {
 225     const real          tol = 1e-4;
 226     bool                bChanged;
 227     int                 m, j, ai, iter;
 228     rvec                dx, xtest;
 229
 230     /* First simple calculation */
 231     gmx_calc_cog(top, x, nrefat, index, xout);
 232     /* Now check if any atom is more than half the box from the COM */
 233     if (pbc)
 234     {
 235         iter = 0;
 236         do
 237         {
 238             bChanged = false;
 239             for (m = 0; m < nrefat; ++m)
 240             {
 241                 ai = index[m];
 242                 pbc_dx(pbc, x[ai], xout, dx);
 243                 rvec_add(xout, dx, xtest);
 244                 for (j = 0; j < DIM; ++j)
 245                 {
 246                     if (fabs(xtest[j] - x[ai][j]) > tol)
 247                     {
 248                         /* Here we have used the wrong image for contributing to the COM */
 249                         xout[j] += (xtest[j] - x[ai][j]) / nrefat;
 250                         x[ai][j] = xtest[j];
 251                         bChanged = true;
 252                     }
 253                 }
 254             }
 255             iter++;
 256         }
 257         while (bChanged);
 258     }
 259     return 0;
 260 }
 261
 262 /*!
 263  * \param[in]  top    Topology structure with masses.
 264  * \param[in]  x      Position vectors of all atoms.
 265  * \param[in]  pbc    Periodic boundary conditions structure.
 266  * \param[in]  nrefat Number of atoms in the index.
 267  * \param[in]  index  Indices of atoms.
 268  * \param[out] xout   COM position for the indexed atoms.
 269  * \returns    0 on success, EINVAL if \p top is NULL.
 270  *
 271  * Works as gmx_calc_com(), but takes into account periodic boundary
 272  * conditions: If any atom is more than half the box from the COM,
 273  * it is wrapped around and a new COM is calculated. This is repeated
 274  * until no atoms violate the condition.
 275  *
 276  * Modified from src/tools/gmx_sorient.c in Gromacs distribution.
 277  */
 278 int
 279 gmx_calc_com_pbc(t_topology *top, rvec x[], t_pbc *pbc,
 280                  int nrefat, atom_id index[], rvec xout)
 281 {
 282     const real          tol = 1e-4;
 283     bool                bChanged;
 284     int                 m, j, ai, iter;
 285     real                mass, mtot;
 286     rvec                dx, xtest;
 287
 288     if (!top)
 289     {
 290         gmx_incons("no masses available while mass weighting was requested");
 291         return EINVAL;
 292     }
 293     /* First simple calculation */
 294     clear_rvec(xout);
 295     mtot = 0;
 296     for (m = 0; m < nrefat; ++m)
 297     {
 298         ai   = index[m];
 299         mass = top->atoms.atom[ai].m;
 300         for (j = 0; j < DIM; ++j)
 301         {
 302             xout[j] += mass * x[ai][j];
 303         }
 304         mtot += mass;
 305     }
 306     svmul(1.0/mtot, xout, xout);
 307     /* Now check if any atom is more than half the box from the COM */
 308     if (pbc)
 309     {
 310         iter = 0;
 311         do
 312         {
 313             bChanged = false;
 314             for (m = 0; m < nrefat; ++m)
 315             {
 316                 ai   = index[m];
 317                 mass = top->atoms.atom[ai].m / mtot;
 318                 pbc_dx(pbc, x[ai], xout, dx);
 319                 rvec_add(xout, dx, xtest);
 320                 for (j = 0; j < DIM; ++j)
 321                 {
 322                     if (fabs(xtest[j] - x[ai][j]) > tol)
 323                     {
 324                         /* Here we have used the wrong image for contributing to the COM */
 325                         xout[j] += mass * (xtest[j] - x[ai][j]);
 326                         x[ai][j] = xtest[j];
 327                         bChanged = true;
 328                     }
 329                 }
 330             }
 331             iter++;
 332         }
 333         while (bChanged);
 334     }
 335     return 0;
 336 }
 337
 338 /*!
 339  * \param[in]  top   Topology structure with masses
 340  *   (can be NULL if \p bMASS==false).
 341  * \param[in]  x     Position vectors of all atoms.
 342  * \param[in]  pbc    Periodic boundary conditions structure.
 343  * \param[in]  nrefat Number of atoms in the index.
 344  * \param[in]  index Indices of atoms.
 345  * \param[in]  bMass If true, mass weighting is used.
 346  * \param[out] xout  COM/COG position for the indexed atoms.
 347  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 348  *
 349  * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
 350  * \p bMass.
 351  * Other parameters are passed unmodified to these functions.
 352  */
 353 int
 354 gmx_calc_comg_pbc(t_topology *top, rvec x[], t_pbc *pbc,
 355                   int nrefat, atom_id index[], bool bMass, rvec xout)
 356 {
 357     if (bMass)
 358     {
 359         return gmx_calc_com_pbc(top, x, pbc, nrefat, index, xout);
 360     }
 361     else
 362     {
 363         return gmx_calc_cog_pbc(top, x, pbc, nrefat, index, xout);
 364     }
 365 }
 366
 367
 368 int
 369 gmx_calc_cog_block(t_topology * /* top */, rvec x[], t_block *block, atom_id index[],
 370                    rvec xout[])
 371 {
 372     int                 b, i, ai;
 373     rvec                xb;
 374
 375     for (b = 0; b < block->nr; ++b)
 376     {
 377         clear_rvec(xb);
 378         for (i = block->index[b]; i < block->index[b+1]; ++i)
 379         {
 380             ai = index[i];
 381             rvec_inc(xb, x[ai]);
 382         }
 383         svmul(1.0/(block->index[b+1] - block->index[b]), xb, xout[b]);
 384     }
 385     return 0;
 386 }
 387
 388 /*!
 389  * \param[in]  top   Topology structure with masses.
 390  * \param[in]  x     Position vectors of all atoms.
 391  * \param[in]  block t_block structure that divides \p index into blocks.
 392  * \param[in]  index Indices of atoms.
 393  * \param[out] xout  \p block->nr COM positions.
 394  * \returns    0 on success, EINVAL if \p top is NULL.
 395  *
 396  * Works exactly as gmx_calc_cog_block() with the exception that centers of
 397  * mass are calculated, and hence a topology with masses is required.
 398  */
 399 int
 400 gmx_calc_com_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
 401                    rvec xout[])
 402 {
 403     int                 b, i, ai, d;
 404     rvec                xb;
 405     real                mass, mtot;
 406
 407     if (!top)
 408     {
 409         gmx_incons("no masses available while mass weighting was requested");
 410         return EINVAL;
 411     }
 412     for (b = 0; b < block->nr; ++b)
 413     {
 414         clear_rvec(xb);
 415         mtot = 0;
 416         for (i = block->index[b]; i < block->index[b+1]; ++i)
 417         {
 418             ai   = index[i];
 419             mass = top->atoms.atom[ai].m;
 420             for (d = 0; d < DIM; ++d)
 421             {
 422                 xb[d] += mass * x[ai][d];
 423             }
 424             mtot += mass;
 425         }
 426         svmul(1.0/mtot, xb, xout[b]);
 427     }
 428     return 0;
 429 }
 430
 431 /*!
 432  * \param[in]  top   Topology structure with masses.
 433  * \param[in]  f     Forces on all atoms.
 434  * \param[in]  block t_block structure that divides \p index into blocks.
 435  * \param[in]  index Indices of atoms.
 436  * \param[out] fout  \p block->nr Forces on COG positions.
 437  * \returns    0 on success, EINVAL if \p top is NULL.
 438  */
 439 int
 440 gmx_calc_cog_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
 441                      rvec fout[])
 442 {
 443     int                 b, i, ai, d;
 444     rvec                fb;
 445     real                mass, mtot;
 446
 447     if (!top)
 448     {
 449         gmx_incons("no masses available while mass weighting was needed");
 450         return EINVAL;
 451     }
 452     for (b = 0; b < block->nr; ++b)
 453     {
 454         clear_rvec(fb);
 455         mtot = 0;
 456         for (i = block->index[b]; i < block->index[b+1]; ++i)
 457         {
 458             ai   = index[i];
 459             mass = top->atoms.atom[ai].m;
 460             for (d = 0; d < DIM; ++d)
 461             {
 462                 fb[d] += f[ai][d] / mass;
 463             }
 464             mtot += mass;
 465         }
 466         svmul(mtot / (block->index[b+1] - block->index[b]), fb, fout[b]);
 467     }
 468     return 0;
 469 }
 470
 471 int
 472 gmx_calc_com_f_block(t_topology * /* top */, rvec f[], t_block *block, atom_id index[],
 473                      rvec fout[])
 474 {
 475     for (int b = 0; b < block->nr; ++b)
 476     {
 477         rvec fb;
 478         clear_rvec(fb);
 479         for (int i = block->index[b]; i < block->index[b+1]; ++i)
 480         {
 481             const int ai = index[i];
 482             rvec_inc(fb, f[ai]);
 483         }
 484         copy_rvec(fb, fout[b]);
 485     }
 486     return 0;
 487 }
 488
 489 /*!
 490  * \param[in]  top   Topology structure with masses
 491  *   (can be NULL if \p bMASS==false).
 492  * \param[in]  x     Position vectors of all atoms.
 493  * \param[in]  block t_block structure that divides \p index into blocks.
 494  * \param[in]  index Indices of atoms.
 495  * \param[in]  bMass If true, mass weighting is used.
 496  * \param[out] xout  \p block->nr COM/COG positions.
 497  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 498  *
 499  * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
 500  * value of \p bMass.
 501  * Other parameters are passed unmodified to these functions.
 502  */
 503 int
 504 gmx_calc_comg_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
 505                     bool bMass, rvec xout[])
 506 {
 507     if (bMass)
 508     {
 509         return gmx_calc_com_block(top, x, block, index, xout);
 510     }
 511     else
 512     {
 513         return gmx_calc_cog_block(top, x, block, index, xout);
 514     }
 515 }
 516
 517 /*!
 518  * \param[in]  top   Topology structure with masses
 519  *   (can be NULL if \p bMASS==false).
 520  * \param[in]  f     Forces on all atoms.
 521  * \param[in]  block t_block structure that divides \p index into blocks.
 522  * \param[in]  index Indices of atoms.
 523  * \param[in]  bMass If true, force on COM is calculated.
 524  * \param[out] fout  \p block->nr forces on the COM/COG positions.
 525  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 526  *
 527  * Calls either gmx_calc_com_f_block() or gmx_calc_cog_f_block() depending on
 528  * the value of \p bMass.
 529  * Other parameters are passed unmodified to these functions.
 530  */
 531 int
 532 gmx_calc_comg_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
 533                       bool bMass, rvec fout[])
 534 {
 535     if (bMass)
 536     {
 537         return gmx_calc_com_f_block(top, f, block, index, fout);
 538     }
 539     else
 540     {
 541         return gmx_calc_cog_f_block(top, f, block, index, fout);
 542     }
 543 }
 544
 545 /*!
 546  * \param[in]  top   Topology structure with masses
 547  *   (can be NULL if \p bMASS==false).
 548  * \param[in]  x     Position vectors of all atoms.
 549  * \param[in]  block Blocks for calculation.
 550  * \param[in]  bMass If true, mass weighting is used.
 551  * \param[out] xout  \p block->nr COM/COG positions.
 552  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 553  *
 554  * Calls gmx_calc_comg_block(), converting the \p t_blocka structure into
 555  * a \p t_block and an index. Other parameters are passed unmodified.
 556  *
 557  * \attention
 558  * This function assumes that a pointer to \c t_blocka can be safely typecast
 559  * into \c t_block such that the index fields can still be referenced.
 560  * With the present Gromacs defitions of these types, this is the case,
 561  * but if the layout of these structures is changed, this may lead to strange
 562  * crashes.
 563  */
 564 int
 565 gmx_calc_comg_blocka(t_topology *top, rvec x[], t_blocka *block,
 566                      bool bMass, rvec xout[])
 567 {
 568     /* TODO: It would probably be better to do this without the type cast */
 569     return gmx_calc_comg_block(top, x, (t_block *)block, block->a, bMass, xout);
 570 }
 571
 572 /*!
 573  * \param[in]  top   Topology structure with masses
 574  *   (can be NULL if \p bMASS==true).
 575  * \param[in]  f     Forces on all atoms.
 576  * \param[in]  block Blocks for calculation.
 577  * \param[in]  bMass If true, force on COM is calculated.
 578  * \param[out] fout  \p block->nr forces on the COM/COG positions.
 579  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is false.
 580  *
 581  * Calls gmx_calc_comg_f_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_f_blocka(t_topology *top, rvec f[], t_blocka *block,
 593                        bool bMass, rvec fout[])
 594 {
 595     /* TODO: It would probably be better to do this without the type cast */
 596     return gmx_calc_comg_f_block(top, f, (t_block *)block, block->a, bMass, fout);
 597 }