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
38 #include "gromacs/legacyheaders/typedefs.h"
39 #include "gromacs/legacyheaders/pbc.h"
40 #include "gromacs/legacyheaders/vec.h"
42 #include "gromacs/selection/centerofmass.h"
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.
53 gmx_calc_cog(t_topology *top, rvec x[], int nrefat, atom_id index[], rvec xout)
58 for (m = 0; m < nrefat; ++m)
61 rvec_inc(xout, x[ai]);
63 svmul(1.0/nrefat, xout, xout);
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.
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.
79 gmx_calc_com(t_topology *top, rvec x[], int nrefat, atom_id index[], rvec xout)
86 gmx_incons("no masses available while mass weighting was requested");
91 for (m = 0; m < nrefat; ++m)
94 mass = top->atoms.atom[ai].m;
95 for (j = 0; j < DIM; ++j)
97 xout[j] += mass * x[ai][j];
101 svmul(1.0/mtot, xout, xout);
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.
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().
117 gmx_calc_cog_f(t_topology *top, rvec f[], int nrefat, atom_id index[], rvec fout)
124 gmx_incons("no masses available while mass weighting was needed");
129 for (m = 0; m < nrefat; ++m)
132 mass = top->atoms.atom[ai].m;
133 for (j = 0; j < DIM; ++j)
135 fout[j] += f[ai][j] / mass;
139 svmul(mtot, fout, fout);
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.
153 * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
155 * Other parameters are passed unmodified to these functions.
158 gmx_calc_comg(t_topology *top, rvec x[], int nrefat, atom_id index[],
159 bool bMass, rvec xout)
163 return gmx_calc_com(top, x, nrefat, index, xout);
167 return gmx_calc_cog(top, x, nrefat, index, xout);
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.
181 * Calls either gmx_calc_cog() or gmx_calc_cog_f() depending on the value of
183 * Other parameters are passed unmodified to these functions.
186 gmx_calc_comg_f(t_topology *top, rvec f[], int nrefat, atom_id index[],
187 bool bMass, rvec fout)
191 return gmx_calc_cog(top, f, nrefat, index, fout);
195 return gmx_calc_cog_f(top, f, nrefat, index, fout);
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.
209 * Works exactly as gmx_calc_com_pbc(), but calculates the center of geometry.
212 gmx_calc_cog_pbc(t_topology *top, rvec x[], t_pbc *pbc,
213 int nrefat, atom_id index[], rvec xout)
215 const real tol = 1e-4;
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 */
229 for (m = 0; m < nrefat; ++m)
232 pbc_dx(pbc, x[ai], xout, dx);
233 rvec_add(xout, dx, xtest);
234 for (j = 0; j < DIM; ++j)
236 if (fabs(xtest[j] - x[ai][j]) > tol)
238 /* Here we have used the wrong image for contributing to the COM */
239 xout[j] += (xtest[j] - x[ai][j]) / nrefat;
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.
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.
266 * Modified from src/tools/gmx_sorient.c in Gromacs distribution.
269 gmx_calc_com_pbc(t_topology *top, rvec x[], t_pbc *pbc,
270 int nrefat, atom_id index[], rvec xout)
272 const real tol = 1e-4;
280 gmx_incons("no masses available while mass weighting was requested");
283 /* First simple calculation */
286 for (m = 0; m < nrefat; ++m)
289 mass = top->atoms.atom[ai].m;
290 for (j = 0; j < DIM; ++j)
292 xout[j] += mass * x[ai][j];
296 svmul(1.0/mtot, xout, xout);
297 /* Now check if any atom is more than half the box from the COM */
304 for (m = 0; m < nrefat; ++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)
312 if (fabs(xtest[j] - x[ai][j]) > tol)
314 /* Here we have used the wrong image for contributing to the COM */
315 xout[j] += mass * (xtest[j] - x[ai][j]);
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.
339 * Calls either gmx_calc_com() or gmx_calc_cog() depending on the value of
341 * Other parameters are passed unmodified to these functions.
344 gmx_calc_comg_pbc(t_topology *top, rvec x[], t_pbc *pbc,
345 int nrefat, atom_id index[], bool bMass, rvec xout)
349 return gmx_calc_com_pbc(top, x, pbc, nrefat, index, xout);
353 return gmx_calc_cog_pbc(top, x, pbc, nrefat, index, xout);
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.
367 gmx_calc_cog_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
373 for (b = 0; b < block->nr; ++b)
376 for (i = block->index[b]; i < block->index[b+1]; ++i)
381 svmul(1.0/(block->index[b+1] - block->index[b]), xb, xout[b]);
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.
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.
398 gmx_calc_com_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
407 gmx_incons("no masses available while mass weighting was requested");
410 for (b = 0; b < block->nr; ++b)
414 for (i = block->index[b]; i < block->index[b+1]; ++i)
417 mass = top->atoms.atom[ai].m;
418 for (d = 0; d < DIM; ++d)
420 xb[d] += mass * x[ai][d];
424 svmul(1.0/mtot, xb, xout[b]);
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.
438 gmx_calc_cog_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
447 gmx_incons("no masses available while mass weighting was needed");
450 for (b = 0; b < block->nr; ++b)
454 for (i = block->index[b]; i < block->index[b+1]; ++i)
457 mass = top->atoms.atom[ai].m;
458 for (d = 0; d < DIM; ++d)
460 fb[d] += f[ai][d] / mass;
464 svmul(mtot, fb, fout[b]);
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.
479 * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
481 * Other parameters are passed unmodified to these functions.
484 gmx_calc_comg_block(t_topology *top, rvec x[], t_block *block, atom_id index[],
485 bool bMass, rvec xout[])
489 return gmx_calc_com_block(top, x, block, index, xout);
493 return gmx_calc_cog_block(top, x, block, index, xout);
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.
507 * Calls either gmx_calc_com_block() or gmx_calc_cog_block() depending on the
509 * Other parameters are passed unmodified to these functions.
512 gmx_calc_comg_f_block(t_topology *top, rvec f[], t_block *block, atom_id index[],
513 bool bMass, rvec fout[])
517 return gmx_calc_cog_block(top, f, block, index, fout);
521 return gmx_calc_cog_f_block(top, f, block, index, fout);
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.
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.
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
545 gmx_calc_comg_blocka(t_topology *top, rvec x[], t_blocka *block,
546 bool bMass, rvec xout[])
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);
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.
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.
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
572 gmx_calc_comg_f_blocka(t_topology *top, rvec f[], t_blocka *block,
573 bool bMass, rvec fout[])
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);