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, 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! */
40 #include "gromacs/legacyheaders/vcm.h"
42 #include "gromacs/legacyheaders/macros.h"
43 #include "gromacs/legacyheaders/names.h"
44 #include "gromacs/legacyheaders/network.h"
45 #include "gromacs/legacyheaders/txtdump.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/pbcutil/pbc.h"
48 #include "gromacs/utility/smalloc.h"
50 t_vcm *init_vcm(FILE *fp, gmx_groups_t *groups, t_inputrec *ir)
57 vcm->mode = (ir->nstcomm > 0) ? ir->comm_mode : ecmNO;
58 vcm->ndim = ndof_com(ir);
60 if (vcm->mode == ecmANGULAR && vcm->ndim < 3)
62 gmx_fatal(FARGS, "Can not have angular comm removal with pbc=%s",
63 epbc_names[ir->ePBC]);
66 if (vcm->mode != ecmNO)
68 vcm->nr = groups->grps[egcVCM].nr;
69 /* Allocate one extra for a possible rest group */
70 if (vcm->mode == ecmANGULAR)
72 snew(vcm->group_j, vcm->nr+1);
73 snew(vcm->group_x, vcm->nr+1);
74 snew(vcm->group_i, vcm->nr+1);
75 snew(vcm->group_w, vcm->nr+1);
77 snew(vcm->group_p, vcm->nr+1);
78 snew(vcm->group_v, vcm->nr+1);
79 snew(vcm->group_mass, vcm->nr+1);
80 snew(vcm->group_name, vcm->nr);
81 snew(vcm->group_ndf, vcm->nr);
82 for (g = 0; (g < vcm->nr); g++)
84 vcm->group_ndf[g] = ir->opts.nrdf[g];
87 /* Copy pointer to group names and print it. */
90 fprintf(fp, "Center of mass motion removal mode is %s\n",
92 fprintf(fp, "We have the following groups for center of"
93 " mass motion removal:\n");
95 for (g = 0; (g < vcm->nr); g++)
97 vcm->group_name[g] = *groups->grpname[groups->grps[egcVCM].nm_ind[g]];
100 fprintf(fp, "%3d: %s\n", g, vcm->group_name[g]);
108 static void update_tensor(rvec x, real m0, tensor I)
112 /* Compute inertia tensor contribution due to this atom */
116 I[XX][XX] += x[XX]*x[XX]*m0;
117 I[YY][YY] += x[YY]*x[YY]*m0;
118 I[ZZ][ZZ] += x[ZZ]*x[ZZ]*m0;
127 /* Center of mass code for groups */
128 void calc_vcm_grp(int start, int homenr, t_mdatoms *md,
129 rvec x[], rvec v[], t_vcm *vcm)
132 real m0, xx, xy, xz, yy, yz, zz;
135 if (vcm->mode != ecmNO)
137 /* Also clear a possible rest group */
138 for (g = 0; (g < vcm->nr+1); g++)
140 /* Reset linear momentum */
141 vcm->group_mass[g] = 0;
142 clear_rvec(vcm->group_p[g]);
144 if (vcm->mode == ecmANGULAR)
146 /* Reset anular momentum */
147 clear_rvec(vcm->group_j[g]);
148 clear_rvec(vcm->group_x[g]);
149 clear_rvec(vcm->group_w[g]);
150 clear_mat(vcm->group_i[g]);
155 for (i = start; (i < start+homenr); i++)
163 /* Calculate linear momentum */
164 vcm->group_mass[g] += m0;
165 for (m = 0; (m < DIM); m++)
167 vcm->group_p[g][m] += m0*v[i][m];
170 if (vcm->mode == ecmANGULAR)
172 /* Calculate angular momentum */
173 cprod(x[i], v[i], j0);
175 for (m = 0; (m < DIM); m++)
177 vcm->group_j[g][m] += m0*j0[m];
178 vcm->group_x[g][m] += m0*x[i][m];
180 /* Update inertia tensor */
181 update_tensor(x[i], m0, vcm->group_i[g]);
187 void do_stopcm_grp(int start, int homenr, unsigned short *group_id,
188 rvec x[], rvec v[], t_vcm *vcm)
194 if (vcm->mode != ecmNO)
196 /* Subtract linear momentum */
201 for (i = start; (i < start+homenr); i++)
207 v[i][XX] -= vcm->group_v[g][XX];
211 for (i = start; (i < start+homenr); i++)
217 v[i][XX] -= vcm->group_v[g][XX];
218 v[i][YY] -= vcm->group_v[g][YY];
222 for (i = start; (i < start+homenr); i++)
228 rvec_dec(v[i], vcm->group_v[g]);
232 if (vcm->mode == ecmANGULAR)
234 /* Subtract angular momentum */
235 for (i = start; (i < start+homenr); i++)
241 /* Compute the correction to the velocity for each atom */
242 rvec_sub(x[i], vcm->group_x[g], dx);
243 cprod(vcm->group_w[g], dx, dv);
250 static void get_minv(tensor A, tensor B)
256 tmp[XX][XX] = A[YY][YY] + A[ZZ][ZZ];
257 tmp[YY][XX] = -A[XX][YY];
258 tmp[ZZ][XX] = -A[XX][ZZ];
259 tmp[XX][YY] = -A[XX][YY];
260 tmp[YY][YY] = A[XX][XX] + A[ZZ][ZZ];
261 tmp[ZZ][YY] = -A[YY][ZZ];
262 tmp[XX][ZZ] = -A[XX][ZZ];
263 tmp[YY][ZZ] = -A[YY][ZZ];
264 tmp[ZZ][ZZ] = A[XX][XX] + A[YY][YY];
266 /* This is a hack to prevent very large determinants */
267 rfac = (tmp[XX][XX]+tmp[YY][YY]+tmp[ZZ][ZZ])/3;
270 gmx_fatal(FARGS, "Can not stop center of mass: maybe 2dimensional system");
273 for (m = 0; (m < DIM); m++)
275 for (n = 0; (n < DIM); n++)
281 for (m = 0; (m < DIM); m++)
283 for (n = 0; (n < DIM); n++)
290 void check_cm_grp(FILE *fp, t_vcm *vcm, t_inputrec *ir, real Temp_Max)
293 real ekcm, ekrot, tm, tm_1, Temp_cm;
297 /* First analyse the total results */
298 if (vcm->mode != ecmNO)
300 for (g = 0; (g < vcm->nr); g++)
302 tm = vcm->group_mass[g];
306 svmul(tm_1, vcm->group_p[g], vcm->group_v[g]);
308 /* Else it's zero anyway! */
310 if (vcm->mode == ecmANGULAR)
312 for (g = 0; (g < vcm->nr); g++)
314 tm = vcm->group_mass[g];
319 /* Compute center of mass for this group */
320 for (m = 0; (m < DIM); m++)
322 vcm->group_x[g][m] *= tm_1;
325 /* Subtract the center of mass contribution to the
328 cprod(vcm->group_x[g], vcm->group_v[g], jcm);
329 for (m = 0; (m < DIM); m++)
331 vcm->group_j[g][m] -= tm*jcm[m];
334 /* Subtract the center of mass contribution from the inertia
335 * tensor (this is not as trivial as it seems, but due to
336 * some cancellation we can still do it, even in parallel).
339 update_tensor(vcm->group_x[g], tm, Icm);
340 m_sub(vcm->group_i[g], Icm, vcm->group_i[g]);
342 /* Compute angular velocity, using matrix operation
347 get_minv(vcm->group_i[g], Icm);
348 mvmul(Icm, vcm->group_j[g], vcm->group_w[g]);
350 /* Else it's zero anyway! */
354 for (g = 0; (g < vcm->nr); g++)
357 if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0)
359 for (m = 0; m < vcm->ndim; m++)
361 ekcm += sqr(vcm->group_v[g][m]);
363 ekcm *= 0.5*vcm->group_mass[g];
364 Temp_cm = 2*ekcm/vcm->group_ndf[g];
366 if ((Temp_cm > Temp_Max) && fp)
368 fprintf(fp, "Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
369 vcm->group_name[g], vcm->group_v[g][XX],
370 vcm->group_v[g][YY], vcm->group_v[g][ZZ], Temp_cm);
373 if (vcm->mode == ecmANGULAR)
375 ekrot = 0.5*iprod(vcm->group_j[g], vcm->group_w[g]);
376 if ((ekrot > 1) && fp && !EI_RANDOM(ir->eI))
378 /* if we have an integrator that may not conserve momenta, skip */
379 tm = vcm->group_mass[g];
380 fprintf(fp, "Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
381 vcm->group_name[g], tm, ekrot, det(vcm->group_i[g]));
382 fprintf(fp, " COM: %12.5f %12.5f %12.5f\n",
383 vcm->group_x[g][XX], vcm->group_x[g][YY], vcm->group_x[g][ZZ]);
384 fprintf(fp, " P: %12.5f %12.5f %12.5f\n",
385 vcm->group_p[g][XX], vcm->group_p[g][YY], vcm->group_p[g][ZZ]);
386 fprintf(fp, " V: %12.5f %12.5f %12.5f\n",
387 vcm->group_v[g][XX], vcm->group_v[g][YY], vcm->group_v[g][ZZ]);
388 fprintf(fp, " J: %12.5f %12.5f %12.5f\n",
389 vcm->group_j[g][XX], vcm->group_j[g][YY], vcm->group_j[g][ZZ]);
390 fprintf(fp, " w: %12.5f %12.5f %12.5f\n",
391 vcm->group_w[g][XX], vcm->group_w[g][YY], vcm->group_w[g][ZZ]);
392 pr_rvecs(fp, 0, "Inertia tensor", vcm->group_i[g], DIM);