src/gromacs/selection/centerofmass.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements functions in centerofmass.h.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_selection
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <typedefs.h>
  43 #include <pbc.h>
  44 #include <vec.h>
  45
  46 #include "gromacs/selection/centerofmass.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  * No special function is provided for calculating the force on the center of
 118  * mass, because this can be achieved with gmx_calc_cog().
 119  */
 120 int
 121 gmx_calc_cog_f(t_topology *top, rvec f[], int nrefat, atom_id index[], rvec fout)
 122 {
 123     int                 m, j, ai;
 124     real                mass, mtot;
 125
 126     if (!top)
 127     {
 128         gmx_incons("no masses available while mass weighting was needed");
 129         return EINVAL;
 130     }
 131     clear_rvec(fout);
 132     mtot = 0;
 133     for (m = 0; m < nrefat; ++m)
 134     {
 135         ai = index[m];
 136         mass = top->atoms.atom[ai].m;
 137         for (j = 0; j < DIM; ++j)
 138         {
 139             fout[j] += f[ai][j] / mass;
 140         }
 141         mtot += mass;
 142     }
 143     svmul(mtot, fout, fout);
 144     return 0;
 145 }
 146
 147 /*!
 148  * \param[in]  top   Topology structure with masses
 149  *   (can be NULL if \p bMASS==false).
 150  * \param[in]  x     Position vectors of all atoms.
 151  * \param[in]  nrefat Number of atoms in the index.
 152  * \param[in]  index Indices of atoms.
 153  * \param[in]  bMass If true, mass weighting is used.
 154  * \param[out] xout  COM/COG position for the indexed atoms.
 155  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 156  *
 157  * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
 158  * \p bMass.
 159  * Other parameters are passed unmodified to these functions.
 160  */
 161 int
 162 gmx_calc_comg(t_topology *top, rvec x[], int nrefat, atom_id index[],
 163               bool bMass, rvec xout)
 164 {
 165     if (bMass)
 166     {
 167         return gmx_calc_com(top, x, nrefat, index, xout);
 168     }
 169     else
 170     {
 171         return gmx_calc_cog(top, x, nrefat, index, xout);
 172     }
 173 }
 174
 175 /*!
 176  * \param[in]  top   Topology structure with masses
 177  *   (can be NULL if \p bMASS==true).
 178  * \param[in]  f     Forces on all atoms.
 179  * \param[in]  nrefat Number of atoms in the index.
 180  * \param[in]  index Indices of atoms.
 181  * \param[in]  bMass If true, force on COM is calculated.
 182  * \param[out] fout  Force on the COM/COG position for the indexed atoms.
 183  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is false.
 184  *
 185  * Calls either gmx_calc_cog() or gmx_calc_cog_f() depending on the value of
 186  * \p bMass.
 187  * Other parameters are passed unmodified to these functions.
 188  */
 189 int
 190 gmx_calc_comg_f(t_topology *top, rvec f[], int nrefat, atom_id index[],
 191                 bool bMass, rvec fout)
 192 {
 193     if (bMass)
 194     {
 195         return gmx_calc_cog(top, f, nrefat, index, fout);
 196     }
 197     else
 198     {
 199         return gmx_calc_cog_f(top, f, nrefat, index, fout);
 200     }
 201 }
 202
 203
 204 /*!
 205  * \param[in]  top    Topology structure (unused, can be NULL).
 206  * \param[in]  x      Position vectors of all atoms.
 207  * \param[in]  pbc    Periodic boundary conditions structure.
 208  * \param[in]  nrefat Number of atoms in the index.
 209  * \param[in]  index  Indices of atoms.
 210  * \param[out] xout   COG position for the indexed atoms.
 211  * \returns    0 on success.
 212  *
 213  * Works exactly as gmx_calc_com_pbc(), but calculates the center of geometry.
 214  */
 215 int
 216 gmx_calc_cog_pbc(t_topology *top, rvec x[], t_pbc *pbc,
 217                  int nrefat, atom_id index[], rvec xout)
 218 {
 219     const real          tol = 1e-4;
 220     bool                bChanged;
 221     int                 m, j, ai, iter;
 222     rvec                dx, xtest;
 223
 224     /* First simple calculation */
 225     gmx_calc_cog(top, x, nrefat, index, xout);
 226     /* Now check if any atom is more than half the box from the COM */
 227     if (pbc)
 228     {
 229         iter = 0;
 230         do
 231         {
 232             bChanged = false;
 233             for (m = 0; m < nrefat; ++m)
 234             {
 235                 ai = index[m];
 236                 pbc_dx(pbc, x[ai], xout, dx);
 237                 rvec_add(xout, dx, xtest);
 238                 for (j = 0; j < DIM; ++j)
 239                 {
 240                     if (fabs(xtest[j] - x[ai][j]) > tol)
 241                     {
 242                         /* Here we have used the wrong image for contributing to the COM */
 243                         xout[j] += (xtest[j] - x[ai][j]) / nrefat;
 244                         x[ai][j] = xtest[j];
 245                         bChanged = true;
 246                     }
 247                 }
 248             }
 249             iter++;
 250         }
 251         while (bChanged);
 252     }
 253     return 0;
 254 }
 255
 256 /*!
 257  * \param[in]  top    Topology structure with masses.
 258  * \param[in]  x      Position vectors of all atoms.
 259  * \param[in]  pbc    Periodic boundary conditions structure.
 260  * \param[in]  nrefat Number of atoms in the index.
 261  * \param[in]  index  Indices of atoms.
 262  * \param[out] xout   COM position for the indexed atoms.
 263  * \returns    0 on success, EINVAL if \p top is NULL.
 264  *
 265  * Works as gmx_calc_com(), but takes into account periodic boundary
 266  * conditions: If any atom is more than half the box from the COM,
 267  * it is wrapped around and a new COM is calculated. This is repeated
 268  * until no atoms violate the condition.
 269  *
 270  * Modified from src/tools/gmx_sorient.c in Gromacs distribution.
 271  */
 272 int
 273 gmx_calc_com_pbc(t_topology *top, rvec x[], t_pbc *pbc,
 274                  int nrefat, atom_id index[], rvec xout)
 275 {
 276     const real          tol = 1e-4;
 277     bool                bChanged;
 278     int                 m, j, ai, iter;
 279     real                mass, mtot;
 280     rvec                dx, xtest;
 281
 282     if (!top)
 283     {
 284         gmx_incons("no masses available while mass weighting was requested");
 285         return EINVAL;
 286     }
 287     /* First simple calculation */
 288     clear_rvec(xout);
 289     mtot = 0;
 290     for (m = 0; m < nrefat; ++m)
 291     {
 292         ai = index[m];
 293         mass = top->atoms.atom[ai].m;
 294         for (j = 0; j < DIM; ++j)
 295         {
 296             xout[j] += mass * x[ai][j];
 297         }
 298         mtot += mass;
 299     }
 300     svmul(1.0/mtot, xout, xout);
 301     /* Now check if any atom is more than half the box from the COM */
 302     if (pbc)
 303     {
 304         iter = 0;
 305         do
 306         {
 307             bChanged = false;
 308             for (m = 0; m < nrefat; ++m)
 309             {
 310                 ai = index[m];
 311                 mass = top->atoms.atom[ai].m / mtot;
 312                 pbc_dx(pbc, x[ai], xout, dx);
 313                 rvec_add(xout, dx, xtest);
 314                 for (j = 0; j < DIM; ++j)
 315                 {
 316                     if (fabs(xtest[j] - x[ai][j]) > tol)
 317                     {
 318                         /* Here we have used the wrong image for contributing to the COM */
 319                         xout[j] += mass * (xtest[j] - x[ai][j]);
 320                         x[ai][j] = xtest[j];
 321                         bChanged = true;
 322                     }
 323                 }
 324             }
 325             iter++;
 326         }
 327         while (bChanged);
 328     }
 329     return 0;
 330 }
 331
 332 /*!
 333  * \param[in]  top   Topology structure with masses
 334  *   (can be NULL if \p bMASS==false).
 335  * \param[in]  x     Position vectors of all atoms.
 336  * \param[in]  pbc    Periodic boundary conditions structure.
 337  * \param[in]  nrefat Number of atoms in the index.
 338  * \param[in]  index Indices of atoms.
 339  * \param[in]  bMass If true, mass weighting is used.
 340  * \param[out] xout  COM/COG position for the indexed atoms.
 341  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 342  *
 343  * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
 344  * \p bMass.
 345  * Other parameters are passed unmodified to these functions.
 346  */
 347 int
 348 gmx_calc_comg_pbc(t_topology *top, rvec x[], t_pbc *pbc,
 349                   int nrefat, atom_id index[], bool bMass, rvec xout)
 350 {
 351     if (bMass)
 352     {
 353         return gmx_calc_com_pbc(top, x, pbc, nrefat, index, xout);
 354     }
 355     else
 356     {
 357         return gmx_calc_cog_pbc(top, x, pbc, nrefat, index, xout);
 358     }
 359 }
 360
 361
 362 /*!
 363  * \param[in]  top   Topology structure (unused, can be NULL).
 364  * \param[in]  x     Position vectors of all atoms.
 365  * \param[in]  block t_block structure that divides \p index into blocks.
 366  * \param[in]  index Indices of atoms.
 367  * \param[out] xout  \p block->nr COG positions.
 368  * \returns    0 on success.
 369  */
 370 int
 371 gmx_calc_cog_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
 372                    rvec xout[])
 373 {
 374     int                 b, i, ai;
 375     rvec                xb;
 376
 377     for (b = 0; b < block->nr; ++b)
 378     {
 379         clear_rvec(xb);
 380         for (i = block->index[b]; i < block->index[b+1]; ++i)
 381         {
 382             ai = index[i];
 383             rvec_inc(xb, x[ai]);
 384         }
 385         svmul(1.0/(block->index[b+1] - block->index[b]), xb, xout[b]);
 386     }
 387     return 0;
 388 }
 389
 390 /*!
 391  * \param[in]  top   Topology structure with masses.
 392  * \param[in]  x     Position vectors of all atoms.
 393  * \param[in]  block t_block structure that divides \p index into blocks.
 394  * \param[in]  index Indices of atoms.
 395  * \param[out] xout  \p block->nr COM positions.
 396  * \returns    0 on success, EINVAL if \p top is NULL.
 397  *
 398  * Works exactly as gmx_calc_cog_block() with the exception that centers of
 399  * mass are calculated, and hence a topology with masses is required.
 400  */
 401 int
 402 gmx_calc_com_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
 403                    rvec xout[])
 404 {
 405     int                 b, i, ai, d;
 406     rvec                xb;
 407     real                mass, mtot;
 408
 409     if (!top)
 410     {
 411         gmx_incons("no masses available while mass weighting was requested");
 412         return EINVAL;
 413     }
 414     for (b = 0; b < block->nr; ++b)
 415     {
 416         clear_rvec(xb);
 417         mtot = 0;
 418         for (i = block->index[b]; i < block->index[b+1]; ++i)
 419         {
 420             ai = index[i];
 421             mass = top->atoms.atom[ai].m;
 422             for (d = 0; d < DIM; ++d)
 423             {
 424                 xb[d] += mass * x[ai][d];
 425             }
 426             mtot += mass;
 427         }
 428         svmul(1.0/mtot, xb, xout[b]);
 429     }
 430     return 0;
 431 }
 432
 433 /*!
 434  * \param[in]  top   Topology structure with masses.
 435  * \param[in]  f     Forces on all atoms.
 436  * \param[in]  block t_block structure that divides \p index into blocks.
 437  * \param[in]  index Indices of atoms.
 438  * \param[out] fout  \p block->nr Forces on COG positions.
 439  * \returns    0 on success, EINVAL if \p top is NULL.
 440  */
 441 int
 442 gmx_calc_cog_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
 443                      rvec fout[])
 444 {
 445     int                 b, i, ai, d;
 446     rvec                fb;
 447     real                mass, mtot;
 448
 449     if (!top)
 450     {
 451         gmx_incons("no masses available while mass weighting was needed");
 452         return EINVAL;
 453     }
 454     for (b = 0; b < block->nr; ++b)
 455     {
 456         clear_rvec(fb);
 457         mtot = 0;
 458         for (i = block->index[b]; i < block->index[b+1]; ++i)
 459         {
 460             ai = index[i];
 461             mass = top->atoms.atom[ai].m;
 462             for (d = 0; d < DIM; ++d)
 463             {
 464                 fb[d] += f[ai][d] / mass;
 465             }
 466             mtot += mass;
 467         }
 468         svmul(mtot, fb, fout[b]);
 469     }
 470     return 0;
 471 }
 472
 473 /*!
 474  * \param[in]  top   Topology structure with masses
 475  *   (can be NULL if \p bMASS==false).
 476  * \param[in]  x     Position vectors of all atoms.
 477  * \param[in]  block t_block structure that divides \p index into blocks.
 478  * \param[in]  index Indices of atoms.
 479  * \param[in]  bMass If true, mass weighting is used.
 480  * \param[out] xout  \p block->nr COM/COG positions.
 481  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 482  *
 483  * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
 484  * value of \p bMass.
 485  * Other parameters are passed unmodified to these functions.
 486  */
 487 int
 488 gmx_calc_comg_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
 489                     bool bMass, rvec xout[])
 490 {
 491     if (bMass)
 492     {
 493         return gmx_calc_com_block(top, x, block, index, xout);
 494     }
 495     else
 496     {
 497         return gmx_calc_cog_block(top, x, block, index, xout);
 498     }
 499 }
 500
 501 /*!
 502  * \param[in]  top   Topology structure with masses
 503  *   (can be NULL if \p bMASS==false).
 504  * \param[in]  f     Forces on all atoms.
 505  * \param[in]  block t_block structure that divides \p index into blocks.
 506  * \param[in]  index Indices of atoms.
 507  * \param[in]  bMass If true, force on COM is calculated.
 508  * \param[out] fout  \p block->nr forces on the COM/COG positions.
 509  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 510  *
 511  * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
 512  * value of \p bMass.
 513  * Other parameters are passed unmodified to these functions.
 514  */
 515 int
 516 gmx_calc_comg_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
 517                       bool bMass, rvec fout[])
 518 {
 519     if (bMass)
 520     {
 521         return gmx_calc_cog_block(top, f, block, index, fout);
 522     }
 523     else
 524     {
 525         return gmx_calc_cog_f_block(top, f, block, index, fout);
 526     }
 527 }
 528
 529 /*!
 530  * \param[in]  top   Topology structure with masses
 531  *   (can be NULL if \p bMASS==false).
 532  * \param[in]  x     Position vectors of all atoms.
 533  * \param[in]  block Blocks for calculation.
 534  * \param[in]  bMass If true, mass weighting is used.
 535  * \param[out] xout  \p block->nr COM/COG positions.
 536  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is true.
 537  *
 538  * Calls gmx_calc_comg_block(), converting the \p t_blocka structure into
 539  * a \p t_block and an index. Other parameters are passed unmodified.
 540  *
 541  * \attention
 542  * This function assumes that a pointer to \c t_blocka can be safely typecast
 543  * into \c t_block such that the index fields can still be referenced.
 544  * With the present Gromacs defitions of these types, this is the case,
 545  * but if the layout of these structures is changed, this may lead to strange
 546  * crashes.
 547  */
 548 int
 549 gmx_calc_comg_blocka(t_topology *top, rvec x[], t_blocka *block,
 550                      bool bMass, rvec xout[])
 551 {
 552     /* TODO: It would probably be better to do this without the type cast */
 553     return gmx_calc_comg_block(top, x, (t_block *)block, block->a, bMass, xout);
 554 }
 555
 556 /*!
 557  * \param[in]  top   Topology structure with masses
 558  *   (can be NULL if \p bMASS==true).
 559  * \param[in]  f     Forces on all atoms.
 560  * \param[in]  block Blocks for calculation.
 561  * \param[in]  bMass If true, force on COM is calculated.
 562  * \param[out] fout  \p block->nr forces on the COM/COG positions.
 563  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is false.
 564  *
 565  * Calls gmx_calc_comg_f_block(), converting the \p t_blocka structure into
 566  * a \p t_block and an index. Other parameters are passed unmodified.
 567  *
 568  * \attention
 569  * This function assumes that a pointer to \c t_blocka can be safely typecast
 570  * into \c t_block such that the index fields can still be referenced.
 571  * With the present Gromacs defitions of these types, this is the case,
 572  * but if the layout of these structures is changed, this may lead to strange
 573  * crashes.
 574  */
 575 int
 576 gmx_calc_comg_f_blocka(t_topology *top, rvec f[], t_blocka *block,
 577                        bool bMass, rvec fout[])
 578 {
 579     /* TODO: It would probably be better to do this without the type cast */
 580     return gmx_calc_comg_f_block(top, f, (t_block *)block, block->a, bMass, fout);
 581 }