2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
37 /* This file is completely threadsafe - keep it that way! */
42 #include "gromacs/gmxlib/network.h"
43 #include "gromacs/math/functions.h"
44 #include "gromacs/math/invertmatrix.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/math/vecdump.h"
47 #include "gromacs/mdlib/gmx_omp_nthreads.h"
48 #include "gromacs/mdtypes/inputrec.h"
49 #include "gromacs/mdtypes/md_enums.h"
50 #include "gromacs/pbcutil/pbc.h"
51 #include "gromacs/topology/topology.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/gmxassert.h"
54 #include "gromacs/utility/gmxomp.h"
55 #include "gromacs/utility/smalloc.h"
57 t_vcm
*init_vcm(FILE *fp
, gmx_groups_t
*groups
, const t_inputrec
*ir
)
64 vcm
->mode
= (ir
->nstcomm
> 0) ? ir
->comm_mode
: ecmNO
;
65 vcm
->ndim
= ndof_com(ir
);
66 vcm
->timeStep
= ir
->nstcomm
*ir
->delta_t
;
68 if (vcm
->mode
== ecmANGULAR
&& vcm
->ndim
< 3)
70 gmx_fatal(FARGS
, "Can not have angular comm removal with pbc=%s",
71 epbc_names
[ir
->ePBC
]);
74 if (vcm
->mode
!= ecmNO
)
76 vcm
->nr
= groups
->grps
[egcVCM
].nr
;
77 /* Allocate one extra for a possible rest group */
78 vcm
->size
= vcm
->nr
+ 1;
79 /* We need vcm->nr+1 elements per thread, but to avoid cache
80 * invalidation we add 2 elements to get a 152 byte separation.
82 vcm
->stride
= vcm
->nr
+ 3;
83 if (vcm
->mode
== ecmANGULAR
)
85 snew(vcm
->group_j
, vcm
->size
);
86 snew(vcm
->group_x
, vcm
->size
);
87 snew(vcm
->group_i
, vcm
->size
);
88 snew(vcm
->group_w
, vcm
->size
);
90 snew(vcm
->group_p
, vcm
->size
);
91 snew(vcm
->group_v
, vcm
->size
);
92 snew(vcm
->group_mass
, vcm
->size
);
93 snew(vcm
->group_name
, vcm
->size
);
94 snew(vcm
->group_ndf
, vcm
->size
);
95 for (g
= 0; (g
< vcm
->nr
); g
++)
97 vcm
->group_ndf
[g
] = ir
->opts
.nrdf
[g
];
100 /* Copy pointer to group names and print it. */
103 fprintf(fp
, "Center of mass motion removal mode is %s\n",
105 fprintf(fp
, "We have the following groups for center of"
106 " mass motion removal:\n");
108 for (g
= 0; (g
< vcm
->nr
); g
++)
110 vcm
->group_name
[g
] = *groups
->grpname
[groups
->grps
[egcVCM
].nm_ind
[g
]];
113 fprintf(fp
, "%3d: %s\n", g
, vcm
->group_name
[g
]);
117 snew(vcm
->thread_vcm
, gmx_omp_nthreads_get(emntDefault
) * vcm
->stride
);
123 static void update_tensor(rvec x
, real m0
, tensor I
)
127 /* Compute inertia tensor contribution due to this atom */
131 I
[XX
][XX
] += x
[XX
]*x
[XX
]*m0
;
132 I
[YY
][YY
] += x
[YY
]*x
[YY
]*m0
;
133 I
[ZZ
][ZZ
] += x
[ZZ
]*x
[ZZ
]*m0
;
142 /* Center of mass code for groups */
143 void calc_vcm_grp(int start
, int homenr
, t_mdatoms
*md
,
144 rvec x
[], rvec v
[], t_vcm
*vcm
)
146 int nthreads
= gmx_omp_nthreads_get(emntDefault
);
147 if (vcm
->mode
!= ecmNO
)
149 #pragma omp parallel num_threads(nthreads)
151 int t
= gmx_omp_get_thread_num();
152 for (int g
= 0; g
< vcm
->size
; g
++)
154 /* Reset linear momentum */
155 t_vcm_thread
*vcm_t
= &vcm
->thread_vcm
[t
*vcm
->stride
+ g
];
157 clear_rvec(vcm_t
->p
);
158 if (vcm
->mode
== ecmANGULAR
)
160 /* Reset angular momentum */
161 clear_rvec(vcm_t
->j
);
162 clear_rvec(vcm_t
->x
);
167 #pragma omp for schedule(static)
168 for (int i
= start
; i
< start
+homenr
; i
++)
171 real m0
= md
->massT
[i
];
176 t_vcm_thread
*vcm_t
= &vcm
->thread_vcm
[t
*vcm
->stride
+ g
];
177 /* Calculate linear momentum */
180 for (m
= 0; (m
< DIM
); m
++)
182 vcm_t
->p
[m
] += m0
*v
[i
][m
];
185 if (vcm
->mode
== ecmANGULAR
)
187 /* Calculate angular momentum */
189 cprod(x
[i
], v
[i
], j0
);
191 for (m
= 0; (m
< DIM
); m
++)
193 vcm_t
->j
[m
] += m0
*j0
[m
];
194 vcm_t
->x
[m
] += m0
*x
[i
][m
];
196 /* Update inertia tensor */
197 update_tensor(x
[i
], m0
, vcm_t
->i
);
201 for (int g
= 0; g
< vcm
->size
; g
++)
203 /* Reset linear momentum */
204 vcm
->group_mass
[g
] = 0;
205 clear_rvec(vcm
->group_p
[g
]);
206 if (vcm
->mode
== ecmANGULAR
)
208 /* Reset angular momentum */
209 clear_rvec(vcm
->group_j
[g
]);
210 clear_rvec(vcm
->group_x
[g
]);
211 clear_rvec(vcm
->group_w
[g
]);
212 clear_mat(vcm
->group_i
[g
]);
215 for (int t
= 0; t
< nthreads
; t
++)
217 t_vcm_thread
*vcm_t
= &vcm
->thread_vcm
[t
*vcm
->stride
+ g
];
218 vcm
->group_mass
[g
] += vcm_t
->mass
;
219 rvec_inc(vcm
->group_p
[g
], vcm_t
->p
);
220 if (vcm
->mode
== ecmANGULAR
)
222 rvec_inc(vcm
->group_j
[g
], vcm_t
->j
);
223 rvec_inc(vcm
->group_x
[g
], vcm_t
->x
);
224 m_add(vcm_t
->i
, vcm
->group_i
[g
], vcm
->group_i
[g
]);
232 /*! \brief Remove the COM motion velocity from the velocities
234 * \note This routine should be called from within an OpenMP parallel region.
236 * \tparam numDimensions Correct dimensions 0 to \p numDimensions-1
237 * \param[in] homenr The number of atoms to correct
238 * \param[in] group_id List of VCM group ids, when nullptr is passed all atoms are assumed to be in group 0
239 * \param[in,out] v The velocities to correct
240 * \param[in] vcm VCM data
242 template<int numDimensions
>
244 doStopComMotionLinear(int homenr
,
245 const unsigned short *group_id
,
249 if (group_id
== nullptr)
251 #pragma omp for schedule(static)
252 for (int i
= 0; i
< homenr
; i
++)
254 for (int d
= 0; d
< numDimensions
; d
++)
256 v
[i
][d
] -= vcm
.group_v
[0][d
];
262 #pragma omp for schedule(static)
263 for (int i
= 0; i
< homenr
; i
++)
265 const int g
= group_id
[i
];
266 for (int d
= 0; d
< numDimensions
; d
++)
268 v
[i
][d
] -= vcm
.group_v
[g
][d
];
274 /*! \brief Remove the COM motion velocity from the velocities, correct the coordinates assuming constant acceleration
276 * \note This routine should be called from within an OpenMP parallel region.
278 * \tparam numDimensions Correct dimensions 0 to \p numDimensions-1
279 * \param[in] homenr The number of atoms to correct
280 * \param[in] group_id List of VCM group ids, when nullptr is passed all atoms are assumed to be in group 0
281 * \param[in,out] x The coordinates to correct
282 * \param[in,out] v The velocities to correct
283 * \param[in] vcm VCM data
285 template<int numDimensions
>
287 doStopComMotionAccelerationCorrection(int homenr
,
288 const unsigned short *group_id
,
289 rvec
* gmx_restrict x
,
290 rvec
* gmx_restrict v
,
293 const real xCorrectionFactor
= 0.5*vcm
.timeStep
;
295 if (group_id
== nullptr)
297 #pragma omp for schedule(static)
298 for (int i
= 0; i
< homenr
; i
++)
300 for (int d
= 0; d
< numDimensions
; d
++)
302 x
[i
][d
] -= vcm
.group_v
[0][d
]*xCorrectionFactor
;
303 v
[i
][d
] -= vcm
.group_v
[0][d
];
309 #pragma omp for schedule(static)
310 for (int i
= 0; i
< homenr
; i
++)
312 const int g
= group_id
[i
];
313 for (int d
= 0; d
< numDimensions
; d
++)
315 x
[i
][d
] -= vcm
.group_v
[g
][d
]*xCorrectionFactor
;
316 v
[i
][d
] -= vcm
.group_v
[g
][d
];
322 void do_stopcm_grp(int homenr
, const unsigned short *group_id
,
323 rvec x
[], rvec v
[], const t_vcm
&vcm
)
325 if (vcm
.mode
!= ecmNO
)
327 // cppcheck-suppress unreadVariable
328 int gmx_unused nth
= gmx_omp_nthreads_get(emntDefault
);
329 #pragma omp parallel num_threads(nth)
331 if (vcm
.mode
== ecmLINEAR
||
332 vcm
.mode
== ecmANGULAR
||
333 (vcm
.mode
== ecmLINEAR_ACCELERATION_CORRECTION
&& x
== nullptr))
335 /* Subtract linear momentum for v */
339 doStopComMotionLinear
<1>(homenr
, group_id
, v
, vcm
);
342 doStopComMotionLinear
<2>(homenr
, group_id
, v
, vcm
);
345 doStopComMotionLinear
<3>(homenr
, group_id
, v
, vcm
);
351 GMX_ASSERT(vcm
.mode
== ecmLINEAR_ACCELERATION_CORRECTION
, "When the mode is not linear or angular, it should be acceleration correction");
352 /* Subtract linear momentum for v and x*/
356 doStopComMotionAccelerationCorrection
<1>(homenr
, group_id
, x
, v
, vcm
);
359 doStopComMotionAccelerationCorrection
<2>(homenr
, group_id
, x
, v
, vcm
);
362 doStopComMotionAccelerationCorrection
<3>(homenr
, group_id
, x
, v
, vcm
);
367 if (vcm
.mode
== ecmANGULAR
)
369 /* Subtract angular momentum */
370 GMX_ASSERT(x
, "Need x to compute angular momentum correction");
373 #pragma omp for schedule(static)
374 for (int i
= 0; i
< homenr
; i
++)
380 /* Compute the correction to the velocity for each atom */
382 rvec_sub(x
[i
], vcm
.group_x
[g
], dx
);
383 cprod(vcm
.group_w
[g
], dx
, dv
);
391 static void get_minv(tensor A
, tensor B
)
397 tmp
[XX
][XX
] = A
[YY
][YY
] + A
[ZZ
][ZZ
];
398 tmp
[YY
][XX
] = -A
[XX
][YY
];
399 tmp
[ZZ
][XX
] = -A
[XX
][ZZ
];
400 tmp
[XX
][YY
] = -A
[XX
][YY
];
401 tmp
[YY
][YY
] = A
[XX
][XX
] + A
[ZZ
][ZZ
];
402 tmp
[ZZ
][YY
] = -A
[YY
][ZZ
];
403 tmp
[XX
][ZZ
] = -A
[XX
][ZZ
];
404 tmp
[YY
][ZZ
] = -A
[YY
][ZZ
];
405 tmp
[ZZ
][ZZ
] = A
[XX
][XX
] + A
[YY
][YY
];
407 /* This is a hack to prevent very large determinants */
408 rfac
= (tmp
[XX
][XX
]+tmp
[YY
][YY
]+tmp
[ZZ
][ZZ
])/3;
411 gmx_fatal(FARGS
, "Can not stop center of mass: maybe 2dimensional system");
414 for (m
= 0; (m
< DIM
); m
++)
416 for (n
= 0; (n
< DIM
); n
++)
421 gmx::invertMatrix(tmp
, B
);
422 for (m
= 0; (m
< DIM
); m
++)
424 for (n
= 0; (n
< DIM
); n
++)
431 void check_cm_grp(FILE *fp
, t_vcm
*vcm
, t_inputrec
*ir
, real Temp_Max
)
434 real ekcm
, ekrot
, tm
, tm_1
, Temp_cm
;
438 /* First analyse the total results */
439 if (vcm
->mode
!= ecmNO
)
441 for (g
= 0; (g
< vcm
->nr
); g
++)
443 tm
= vcm
->group_mass
[g
];
447 svmul(tm_1
, vcm
->group_p
[g
], vcm
->group_v
[g
]);
449 /* Else it's zero anyway! */
451 if (vcm
->mode
== ecmANGULAR
)
453 for (g
= 0; (g
< vcm
->nr
); g
++)
455 tm
= vcm
->group_mass
[g
];
460 /* Compute center of mass for this group */
461 for (m
= 0; (m
< DIM
); m
++)
463 vcm
->group_x
[g
][m
] *= tm_1
;
466 /* Subtract the center of mass contribution to the
469 cprod(vcm
->group_x
[g
], vcm
->group_v
[g
], jcm
);
470 for (m
= 0; (m
< DIM
); m
++)
472 vcm
->group_j
[g
][m
] -= tm
*jcm
[m
];
475 /* Subtract the center of mass contribution from the inertia
476 * tensor (this is not as trivial as it seems, but due to
477 * some cancellation we can still do it, even in parallel).
480 update_tensor(vcm
->group_x
[g
], tm
, Icm
);
481 m_sub(vcm
->group_i
[g
], Icm
, vcm
->group_i
[g
]);
483 /* Compute angular velocity, using matrix operation
488 get_minv(vcm
->group_i
[g
], Icm
);
489 mvmul(Icm
, vcm
->group_j
[g
], vcm
->group_w
[g
]);
491 /* Else it's zero anyway! */
495 for (g
= 0; (g
< vcm
->nr
); g
++)
498 if (vcm
->group_mass
[g
] != 0 && vcm
->group_ndf
[g
] > 0)
500 for (m
= 0; m
< vcm
->ndim
; m
++)
502 ekcm
+= gmx::square(vcm
->group_v
[g
][m
]);
504 ekcm
*= 0.5*vcm
->group_mass
[g
];
505 Temp_cm
= 2*ekcm
/vcm
->group_ndf
[g
];
507 if ((Temp_cm
> Temp_Max
) && fp
)
509 fprintf(fp
, "Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
510 vcm
->group_name
[g
], vcm
->group_v
[g
][XX
],
511 vcm
->group_v
[g
][YY
], vcm
->group_v
[g
][ZZ
], Temp_cm
);
514 if (vcm
->mode
== ecmANGULAR
)
516 ekrot
= 0.5*iprod(vcm
->group_j
[g
], vcm
->group_w
[g
]);
517 if ((ekrot
> 1) && fp
&& !EI_RANDOM(ir
->eI
))
519 /* if we have an integrator that may not conserve momenta, skip */
520 tm
= vcm
->group_mass
[g
];
521 fprintf(fp
, "Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
522 vcm
->group_name
[g
], tm
, ekrot
, det(vcm
->group_i
[g
]));
523 fprintf(fp
, " COM: %12.5f %12.5f %12.5f\n",
524 vcm
->group_x
[g
][XX
], vcm
->group_x
[g
][YY
], vcm
->group_x
[g
][ZZ
]);
525 fprintf(fp
, " P: %12.5f %12.5f %12.5f\n",
526 vcm
->group_p
[g
][XX
], vcm
->group_p
[g
][YY
], vcm
->group_p
[g
][ZZ
]);
527 fprintf(fp
, " V: %12.5f %12.5f %12.5f\n",
528 vcm
->group_v
[g
][XX
], vcm
->group_v
[g
][YY
], vcm
->group_v
[g
][ZZ
]);
529 fprintf(fp
, " J: %12.5f %12.5f %12.5f\n",
530 vcm
->group_j
[g
][XX
], vcm
->group_j
[g
][YY
], vcm
->group_j
[g
][ZZ
]);
531 fprintf(fp
, " w: %12.5f %12.5f %12.5f\n",
532 vcm
->group_w
[g
][XX
], vcm
->group_w
[g
][YY
], vcm
->group_w
[g
][ZZ
]);
533 pr_rvecs(fp
, 0, "Inertia tensor", vcm
->group_i
[g
], DIM
);