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