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