src/gmxlib/trajana/centerofmass.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2009, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /*! \internal \file
  39  * \brief Implementation of functions in centerofmass.h.
  40  */
  41 #ifdef HAVE_CONFIG_H
  42 #include <config.h>
  43 #endif
  44
  45 #include <typedefs.h>
  46 #include <pbc.h>
  47 #include <vec.h>
  48
  49 #include <centerofmass.h>
  50
  51 /*!
  52  * \param[in]  top    Topology structure (unused, can be NULL).
  53  * \param[in]  x      Position vectors of all atoms.
  54  * \param[in]  nrefat Number of atoms in the index.
  55  * \param[in]  index  Indices of atoms.
  56  * \param[out] xout   COG position for the indexed atoms.
  57  * \returns    0 on success.
  58  */
  59 int
  60 gmx_calc_cog(t_topology *top, rvec x[], int nrefat, atom_id index[], rvec xout)
  61 {
  62     int                 m, j, ai;
  63
  64     clear_rvec(xout);
  65     for (m = 0; m < nrefat; ++m)
  66     {
  67         ai = index[m];
  68         rvec_inc(xout, x[ai]);
  69     }
  70     svmul(1.0/nrefat, xout, xout);
  71     return 0;
  72 }
  73
  74 /*!
  75  * \param[in]  top    Topology structure with masses.
  76  * \param[in]  x      Position vectors of all atoms.
  77  * \param[in]  nrefat Number of atoms in the index.
  78  * \param[in]  index  Indices of atoms.
  79  * \param[out] xout   COM position for the indexed atoms.
  80  * \returns    0 on success, EINVAL if \p top is NULL.
  81  *
  82  * Works exactly as gmx_calc_cog() with the exception that a center of
  83  * mass are calculated, and hence a topology with masses is required.
  84  */
  85 int
  86 gmx_calc_com(t_topology *top, rvec x[], int nrefat, atom_id index[], rvec xout)
  87 {
  88     int                 m, j, ai;
  89     real                mass, mtot;
  90
  91     if (!top)
  92     {
  93         gmx_incons("no masses available while mass weighting was requested");
  94         return EINVAL;
  95     }
  96     clear_rvec(xout);
  97     mtot = 0;
  98     for (m = 0; m < nrefat; ++m)
  99     {
 100         ai = index[m];
 101         mass = top->atoms.atom[ai].m;
 102         for (j = 0; j < DIM; ++j)
 103         {
 104             xout[j] += mass * x[ai][j];
 105         }
 106         mtot += mass;
 107     }
 108     svmul(1.0/mtot, xout, xout);
 109     return 0;
 110 }
 111
 112 /*!
 113  * \param[in]  top    Topology structure with masses.
 114  * \param[in]  f      Forces on all atoms.
 115  * \param[in]  nrefat Number of atoms in the index.
 116  * \param[in]  index  Indices of atoms.
 117  * \param[out] fout   Force on the COG position for the indexed atoms.
 118  * \returns    0 on success, EINVAL if \p top is NULL.
 119  *
 120  * No special function is provided for calculating the force on the center of
 121  * mass, because this can be achieved with gmx_calc_cog().
 122  */
 123 int
 124 gmx_calc_cog_f(t_topology *top, rvec f[], int nrefat, atom_id index[], rvec fout)
 125 {
 126     int                 m, j, ai;
 127     real                mass, mtot;
 128
 129     if (!top)
 130     {
 131         gmx_incons("no masses available while mass weighting was needed");
 132         return EINVAL;
 133     }
 134     clear_rvec(fout);
 135     mtot = 0;
 136     for (m = 0; m < nrefat; ++m)
 137     {
 138         ai = index[m];
 139         mass = top->atoms.atom[ai].m;
 140         for (j = 0; j < DIM; ++j)
 141         {
 142             fout[j] += f[ai][j] / mass;
 143         }
 144         mtot += mass;
 145     }
 146     svmul(mtot, fout, fout);
 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               gmx_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]  x     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] xout  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() or gmx_calc_cog_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                 gmx_bool bMass, rvec fout)
 195 {
 196     if (bMass)
 197     {
 198         return gmx_calc_cog(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     gmx_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     gmx_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[], gmx_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 /*!
 366  * \param[in]  top   Topology structure (unused, can be NULL).
 367  * \param[in]  x     Position vectors of all atoms.
 368  * \param[in]  block t_block structure that divides \p index into blocks.
 369  * \param[in]  index Indices of atoms.
 370  * \param[out] xout  \p block->nr COG positions.
 371  * \returns    0 on success.
 372  */
 373 int
 374 gmx_calc_cog_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
 375                    rvec xout[])
 376 {
 377     int                 b, i, ai;
 378     rvec                xb;
 379
 380     for (b = 0; b < block->nr; ++b)
 381     {
 382         clear_rvec(xb);
 383         for (i = block->index[b]; i < block->index[b+1]; ++i)
 384         {
 385             ai = index[i];
 386             rvec_inc(xb, x[ai]);
 387         }
 388         svmul(1.0/(block->index[b+1] - block->index[b]), xb, xout[b]);
 389     }
 390     return 0;
 391 }
 392
 393 /*!
 394  * \param[in]  top   Topology structure with masses.
 395  * \param[in]  x     Position vectors of all atoms.
 396  * \param[in]  block t_block structure that divides \p index into blocks.
 397  * \param[in]  index Indices of atoms.
 398  * \param[out] xout  \p block->nr COM positions.
 399  * \returns    0 on success, EINVAL if \p top is NULL.
 400  *
 401  * Works exactly as gmx_calc_cog_block() with the exception that centers of
 402  * mass are calculated, and hence a topology with masses is required.
 403  */
 404 int
 405 gmx_calc_com_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
 406                    rvec xout[])
 407 {
 408     int                 b, i, ai, d;
 409     rvec                xb;
 410     real                mass, mtot;
 411
 412     if (!top)
 413     {
 414         gmx_incons("no masses available while mass weighting was requested");
 415         return EINVAL;
 416     }
 417     for (b = 0; b < block->nr; ++b)
 418     {
 419         clear_rvec(xb);
 420         mtot = 0;
 421         for (i = block->index[b]; i < block->index[b+1]; ++i)
 422         {
 423             ai = index[i];
 424             mass = top->atoms.atom[ai].m;
 425             for (d = 0; d < DIM; ++d)
 426             {
 427                 xb[d] += mass * x[ai][d];
 428             }
 429             mtot += mass;
 430         }
 431         svmul(1.0/mtot, xb, xout[b]);
 432     }
 433     return 0;
 434 }
 435
 436 /*!
 437  * \param[in]  top   Topology structure with masses.
 438  * \param[in]  f     Forces on all atoms.
 439  * \param[in]  block t_block structure that divides \p index into blocks.
 440  * \param[in]  index Indices of atoms.
 441  * \param[out] xout  \p block->nr Forces on COG positions.
 442  * \returns    0 on success, EINVAL if \p top is NULL.
 443  */
 444 int
 445 gmx_calc_cog_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
 446                      rvec fout[])
 447 {
 448     int                 b, i, ai, d;
 449     rvec                fb;
 450     real                mass, mtot;
 451
 452     if (!top)
 453     {
 454         gmx_incons("no masses available while mass weighting was needed");
 455         return EINVAL;
 456     }
 457     for (b = 0; b < block->nr; ++b)
 458     {
 459         clear_rvec(fb);
 460         mtot = 0;
 461         for (i = block->index[b]; i < block->index[b+1]; ++i)
 462         {
 463             ai = index[i];
 464             mass = top->atoms.atom[ai].m;
 465             for (d = 0; d < DIM; ++d)
 466             {
 467                 fb[d] += f[ai][d] / mass;
 468             }
 469             mtot += mass;
 470         }
 471         svmul(mtot, fb, fout[b]);
 472     }
 473     return 0;
 474 }
 475
 476 /*!
 477  * \param[in]  top   Topology structure with masses
 478  *   (can be NULL if \p bMASS==FALSE).
 479  * \param[in]  x     Position vectors of all atoms.
 480  * \param[in]  block t_block structure that divides \p index into blocks.
 481  * \param[in]  index Indices of atoms.
 482  * \param[in]  bMass If TRUE, mass weighting is used.
 483  * \param[out] xout  \p block->nr COM/COG positions.
 484  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is TRUE.
 485  *
 486  * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
 487  * value of \p bMass.
 488  * Other parameters are passed unmodified to these functions.
 489  */
 490 int
 491 gmx_calc_comg_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
 492                     gmx_bool bMass, rvec xout[])
 493 {
 494     if (bMass)
 495     {
 496         return gmx_calc_com_block(top, x, block, index, xout);
 497     }
 498     else
 499     {
 500         return gmx_calc_cog_block(top, x, block, index, xout);
 501     }
 502 }
 503
 504 /*!
 505  * \param[in]  top   Topology structure with masses
 506  *   (can be NULL if \p bMASS==FALSE).
 507  * \param[in]  f     Forces on all atoms.
 508  * \param[in]  block t_block structure that divides \p index into blocks.
 509  * \param[in]  index Indices of atoms.
 510  * \param[in]  bMass If TRUE, force on COM is calculated.
 511  * \param[out] xout  \p block->nr forces on the COM/COG positions.
 512  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is TRUE.
 513  *
 514  * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
 515  * value of \p bMass.
 516  * Other parameters are passed unmodified to these functions.
 517  */
 518 int
 519 gmx_calc_comg_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
 520                       gmx_bool bMass, rvec fout[])
 521 {
 522     if (bMass)
 523     {
 524         return gmx_calc_cog_block(top, f, block, index, fout);
 525     }
 526     else
 527     {
 528         return gmx_calc_cog_f_block(top, f, block, index, fout);
 529     }
 530 }
 531
 532 /*!
 533  * \param[in]  top   Topology structure with masses
 534  *   (can be NULL if \p bMASS==FALSE).
 535  * \param[in]  x     Position vectors of all atoms.
 536  * \param[in]  block Blocks for calculation.
 537  * \param[in]  bMass If TRUE, mass weighting is used.
 538  * \param[out] xout  \p block->nr COM/COG positions.
 539  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is TRUE.
 540  *
 541  * Calls gmx_calc_comg_block(), converting the \p t_blocka structure into
 542  * a \p t_block and an index. Other parameters are passed unmodified.
 543  *
 544  * \attention
 545  * This function assumes that a pointer to \c t_blocka can be safely typecast
 546  * into \c t_block such that the index fields can still be referenced.
 547  * With the present Gromacs defitions of these types, this is the case,
 548  * but if the layout of these structures is changed, this may lead to strange
 549  * crashes.
 550  */
 551 int
 552 gmx_calc_comg_blocka(t_topology *top, rvec x[], t_blocka *block,
 553                      gmx_bool bMass, rvec xout[])
 554 {
 555     /* TODO: It would probably be better to do this without the type cast */
 556     return gmx_calc_comg_block(top, x, (t_block *)block, block->a, bMass, xout);
 557 }
 558
 559 /*!
 560  * \param[in]  top   Topology structure with masses
 561  *   (can be NULL if \p bMASS==TRUE).
 562  * \param[in]  f     Forces on all atoms.
 563  * \param[in]  block Blocks for calculation.
 564  * \param[in]  bMass If TRUE, force on COM is calculated.
 565  * \param[out] fout  \p block->nr forces on the COM/COG positions.
 566  * \returns    0 on success, EINVAL if \p top is NULL and \p bMass is FALSE.
 567  *
 568  * Calls gmx_calc_comg_f_block(), converting the \p t_blocka structure into
 569  * a \p t_block and an index. Other parameters are passed unmodified.
 570  *
 571  * \attention
 572  * This function assumes that a pointer to \c t_blocka can be safely typecast
 573  * into \c t_block such that the index fields can still be referenced.
 574  * With the present Gromacs defitions of these types, this is the case,
 575  * but if the layout of these structures is changed, this may lead to strange
 576  * crashes.
 577  */
 578 int
 579 gmx_calc_comg_f_blocka(t_topology *top, rvec f[], t_blocka *block,
 580                        gmx_bool bMass, rvec fout[])
 581 {
 582     /* TODO: It would probably be better to do this without the type cast */
 583     return gmx_calc_comg_f_block(top, f, (t_block *)block, block->a, bMass, fout);
 584 }