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