3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
29 * For more info, check our website at http://www.gromacs.org
33 * Implements functions in centerofmass.h.
35 * \author Teemu Murtola <teemu.murtola@cbr.su.se>
36 * \ingroup module_selection
46 #include "gromacs/selection/centerofmass.h"
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.
57 gmx_calc_cog(t_topology *top, rvec x[], int nrefat, atom_id index[], rvec xout)
62 for (m = 0; m < nrefat; ++m)
65 rvec_inc(xout, x[ai]);
67 svmul(1.0/nrefat, xout, xout);
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.
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.
83 gmx_calc_com(t_topology *top, rvec x[], int nrefat, atom_id index[], rvec xout)
90 gmx_incons("no masses available while mass weighting was requested");
95 for (m = 0; m < nrefat; ++m)
98 mass = top->atoms.atom[ai].m;
99 for (j = 0; j < DIM; ++j)
101 xout[j] += mass * x[ai][j];
105 svmul(1.0/mtot, xout, xout);
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.
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().
121 gmx_calc_cog_f(t_topology *top, rvec f[], int nrefat, atom_id index[], rvec fout)
128 gmx_incons("no masses available while mass weighting was needed");
133 for (m = 0; m < nrefat; ++m)
136 mass = top->atoms.atom[ai].m;
137 for (j = 0; j < DIM; ++j)
139 fout[j] += f[ai][j] / mass;
143 svmul(mtot, fout, fout);
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.
157 * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
159 * Other parameters are passed unmodified to these functions.
162 gmx_calc_comg(t_topology *top, rvec x[], int nrefat, atom_id index[],
163 bool bMass, rvec xout)
167 return gmx_calc_com(top, x, nrefat, index, xout);
171 return gmx_calc_cog(top, x, nrefat, index, xout);
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.
185 * Calls either gmx_calc_cog() or gmx_calc_cog_f() depending on the value of
187 * Other parameters are passed unmodified to these functions.
190 gmx_calc_comg_f(t_topology *top, rvec f[], int nrefat, atom_id index[],
191 bool bMass, rvec fout)
195 return gmx_calc_cog(top, f, nrefat, index, fout);
199 return gmx_calc_cog_f(top, f, nrefat, index, fout);
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.
213 * Works exactly as gmx_calc_com_pbc(), but calculates the center of geometry.
216 gmx_calc_cog_pbc(t_topology *top, rvec x[], t_pbc *pbc,
217 int nrefat, atom_id index[], rvec xout)
219 const real tol = 1e-4;
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 */
233 for (m = 0; m < nrefat; ++m)
236 pbc_dx(pbc, x[ai], xout, dx);
237 rvec_add(xout, dx, xtest);
238 for (j = 0; j < DIM; ++j)
240 if (fabs(xtest[j] - x[ai][j]) > tol)
242 /* Here we have used the wrong image for contributing to the COM */
243 xout[j] += (xtest[j] - x[ai][j]) / nrefat;
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.
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.
270 * Modified from src/tools/gmx_sorient.c in Gromacs distribution.
273 gmx_calc_com_pbc(t_topology *top, rvec x[], t_pbc *pbc,
274 int nrefat, atom_id index[], rvec xout)
276 const real tol = 1e-4;
284 gmx_incons("no masses available while mass weighting was requested");
287 /* First simple calculation */
290 for (m = 0; m < nrefat; ++m)
293 mass = top->atoms.atom[ai].m;
294 for (j = 0; j < DIM; ++j)
296 xout[j] += mass * x[ai][j];
300 svmul(1.0/mtot, xout, xout);
301 /* Now check if any atom is more than half the box from the COM */
308 for (m = 0; m < nrefat; ++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)
316 if (fabs(xtest[j] - x[ai][j]) > tol)
318 /* Here we have used the wrong image for contributing to the COM */
319 xout[j] += mass * (xtest[j] - x[ai][j]);
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.
343 * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
345 * Other parameters are passed unmodified to these functions.
348 gmx_calc_comg_pbc(t_topology *top, rvec x[], t_pbc *pbc,
349 int nrefat, atom_id index[], bool bMass, rvec xout)
353 return gmx_calc_com_pbc(top, x, pbc, nrefat, index, xout);
357 return gmx_calc_cog_pbc(top, x, pbc, nrefat, index, xout);
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.
371 gmx_calc_cog_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
377 for (b = 0; b < block->nr; ++b)
380 for (i = block->index[b]; i < block->index[b+1]; ++i)
385 svmul(1.0/(block->index[b+1] - block->index[b]), xb, xout[b]);
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.
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.
402 gmx_calc_com_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
411 gmx_incons("no masses available while mass weighting was requested");
414 for (b = 0; b < block->nr; ++b)
418 for (i = block->index[b]; i < block->index[b+1]; ++i)
421 mass = top->atoms.atom[ai].m;
422 for (d = 0; d < DIM; ++d)
424 xb[d] += mass * x[ai][d];
428 svmul(1.0/mtot, xb, xout[b]);
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.
442 gmx_calc_cog_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
451 gmx_incons("no masses available while mass weighting was needed");
454 for (b = 0; b < block->nr; ++b)
458 for (i = block->index[b]; i < block->index[b+1]; ++i)
461 mass = top->atoms.atom[ai].m;
462 for (d = 0; d < DIM; ++d)
464 fb[d] += f[ai][d] / mass;
468 svmul(mtot, fb, fout[b]);
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.
483 * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
485 * Other parameters are passed unmodified to these functions.
488 gmx_calc_comg_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
489 bool bMass, rvec xout[])
493 return gmx_calc_com_block(top, x, block, index, xout);
497 return gmx_calc_cog_block(top, x, block, index, xout);
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.
511 * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
513 * Other parameters are passed unmodified to these functions.
516 gmx_calc_comg_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
517 bool bMass, rvec fout[])
521 return gmx_calc_cog_block(top, f, block, index, fout);
525 return gmx_calc_cog_f_block(top, f, block, index, fout);
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.
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.
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
549 gmx_calc_comg_blocka(t_topology *top, rvec x[], t_blocka *block,
550 bool bMass, rvec xout[])
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);
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.
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.
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
576 gmx_calc_comg_f_blocka(t_topology *top, rvec f[], t_blocka *block,
577 bool bMass, rvec fout[])
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);