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