3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
35 /* This file is completely threadsafe - keep it that way! */
49 t_vcm *init_vcm(FILE *fp, gmx_groups_t *groups, t_inputrec *ir)
56 vcm->mode = (ir->nstcomm > 0) ? ir->comm_mode : ecmNO;
57 vcm->ndim = ndof_com(ir);
59 if (vcm->mode == ecmANGULAR && vcm->ndim < 3)
61 gmx_fatal(FARGS, "Can not have angular comm removal with pbc=%s",
62 epbc_names[ir->ePBC]);
65 if (vcm->mode != ecmNO)
67 vcm->nr = groups->grps[egcVCM].nr;
68 /* Allocate one extra for a possible rest group */
69 if (vcm->mode == ecmANGULAR)
71 snew(vcm->group_j, vcm->nr+1);
72 snew(vcm->group_x, vcm->nr+1);
73 snew(vcm->group_i, vcm->nr+1);
74 snew(vcm->group_w, vcm->nr+1);
76 snew(vcm->group_p, vcm->nr+1);
77 snew(vcm->group_v, vcm->nr+1);
78 snew(vcm->group_mass, vcm->nr+1);
79 snew(vcm->group_name, vcm->nr);
80 snew(vcm->group_ndf, vcm->nr);
81 for (g = 0; (g < vcm->nr); g++)
83 vcm->group_ndf[g] = ir->opts.nrdf[g];
86 /* Copy pointer to group names and print it. */
89 fprintf(fp, "Center of mass motion removal mode is %s\n",
91 fprintf(fp, "We have the following groups for center of"
92 " mass motion removal:\n");
94 for (g = 0; (g < vcm->nr); g++)
96 vcm->group_name[g] = *groups->grpname[groups->grps[egcVCM].nm_ind[g]];
99 fprintf(fp, "%3d: %s\n", g, vcm->group_name[g]);
107 static void update_tensor(rvec x, real m0, tensor I)
111 /* Compute inertia tensor contribution due to this atom */
115 I[XX][XX] += x[XX]*x[XX]*m0;
116 I[YY][YY] += x[YY]*x[YY]*m0;
117 I[ZZ][ZZ] += x[ZZ]*x[ZZ]*m0;
126 /* Center of mass code for groups */
127 void calc_vcm_grp(FILE *fp, int start, int homenr, t_mdatoms *md,
128 rvec x[], rvec v[], t_vcm *vcm)
131 real m0, xx, xy, xz, yy, yz, zz;
134 if (vcm->mode != ecmNO)
136 /* Also clear a possible rest group */
137 for (g = 0; (g < vcm->nr+1); g++)
139 /* Reset linear momentum */
140 vcm->group_mass[g] = 0;
141 clear_rvec(vcm->group_p[g]);
143 if (vcm->mode == ecmANGULAR)
145 /* Reset anular momentum */
146 clear_rvec(vcm->group_j[g]);
147 clear_rvec(vcm->group_x[g]);
148 clear_rvec(vcm->group_w[g]);
149 clear_mat(vcm->group_i[g]);
154 for (i = start; (i < start+homenr); i++)
162 /* Calculate linear momentum */
163 vcm->group_mass[g] += m0;
164 for (m = 0; (m < DIM); m++)
166 vcm->group_p[g][m] += m0*v[i][m];
169 if (vcm->mode == ecmANGULAR)
171 /* Calculate angular momentum */
172 cprod(x[i], v[i], j0);
174 for (m = 0; (m < DIM); m++)
176 vcm->group_j[g][m] += m0*j0[m];
177 vcm->group_x[g][m] += m0*x[i][m];
179 /* Update inertia tensor */
180 update_tensor(x[i], m0, vcm->group_i[g]);
186 void do_stopcm_grp(FILE *fp, int start, int homenr, unsigned short *group_id,
187 rvec x[], rvec v[], t_vcm *vcm)
193 if (vcm->mode != ecmNO)
195 /* Subtract linear momentum */
200 for (i = start; (i < start+homenr); i++)
206 v[i][XX] -= vcm->group_v[g][XX];
210 for (i = start; (i < start+homenr); i++)
216 v[i][XX] -= vcm->group_v[g][XX];
217 v[i][YY] -= vcm->group_v[g][YY];
221 for (i = start; (i < start+homenr); i++)
227 rvec_dec(v[i], vcm->group_v[g]);
231 if (vcm->mode == ecmANGULAR)
233 /* Subtract angular momentum */
234 for (i = start; (i < start+homenr); i++)
240 /* Compute the correction to the velocity for each atom */
241 rvec_sub(x[i], vcm->group_x[g], dx);
242 cprod(vcm->group_w[g], dx, dv);
249 static void get_minv(tensor A, tensor B)
255 tmp[XX][XX] = A[YY][YY] + A[ZZ][ZZ];
256 tmp[YY][XX] = -A[XX][YY];
257 tmp[ZZ][XX] = -A[XX][ZZ];
258 tmp[XX][YY] = -A[XX][YY];
259 tmp[YY][YY] = A[XX][XX] + A[ZZ][ZZ];
260 tmp[ZZ][YY] = -A[YY][ZZ];
261 tmp[XX][ZZ] = -A[XX][ZZ];
262 tmp[YY][ZZ] = -A[YY][ZZ];
263 tmp[ZZ][ZZ] = A[XX][XX] + A[YY][YY];
265 /* This is a hack to prevent very large determinants */
266 rfac = (tmp[XX][XX]+tmp[YY][YY]+tmp[ZZ][ZZ])/3;
269 gmx_fatal(FARGS, "Can not stop center of mass: maybe 2dimensional system");
272 for (m = 0; (m < DIM); m++)
274 for (n = 0; (n < DIM); n++)
280 for (m = 0; (m < DIM); m++)
282 for (n = 0; (n < DIM); n++)
289 void check_cm_grp(FILE *fp, t_vcm *vcm, t_inputrec *ir, real Temp_Max)
292 real ekcm, ekrot, tm, tm_1, Temp_cm;
296 /* First analyse the total results */
297 if (vcm->mode != ecmNO)
299 for (g = 0; (g < vcm->nr); g++)
301 tm = vcm->group_mass[g];
305 svmul(tm_1, vcm->group_p[g], vcm->group_v[g]);
307 /* Else it's zero anyway! */
309 if (vcm->mode == ecmANGULAR)
311 for (g = 0; (g < vcm->nr); g++)
313 tm = vcm->group_mass[g];
318 /* Compute center of mass for this group */
319 for (m = 0; (m < DIM); m++)
321 vcm->group_x[g][m] *= tm_1;
324 /* Subtract the center of mass contribution to the
327 cprod(vcm->group_x[g], vcm->group_v[g], jcm);
328 for (m = 0; (m < DIM); m++)
330 vcm->group_j[g][m] -= tm*jcm[m];
333 /* Subtract the center of mass contribution from the inertia
334 * tensor (this is not as trivial as it seems, but due to
335 * some cancellation we can still do it, even in parallel).
338 update_tensor(vcm->group_x[g], tm, Icm);
339 m_sub(vcm->group_i[g], Icm, vcm->group_i[g]);
341 /* Compute angular velocity, using matrix operation
346 get_minv(vcm->group_i[g], Icm);
347 mvmul(Icm, vcm->group_j[g], vcm->group_w[g]);
349 /* Else it's zero anyway! */
353 for (g = 0; (g < vcm->nr); g++)
356 if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0)
358 for (m = 0; m < vcm->ndim; m++)
360 ekcm += sqr(vcm->group_v[g][m]);
362 ekcm *= 0.5*vcm->group_mass[g];
363 Temp_cm = 2*ekcm/vcm->group_ndf[g];
365 if ((Temp_cm > Temp_Max) && fp)
367 fprintf(fp, "Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
368 vcm->group_name[g], vcm->group_v[g][XX],
369 vcm->group_v[g][YY], vcm->group_v[g][ZZ], Temp_cm);
372 if (vcm->mode == ecmANGULAR)
374 ekrot = 0.5*iprod(vcm->group_j[g], vcm->group_w[g]);
375 if ((ekrot > 1) && fp && !EI_RANDOM(ir->eI))
377 /* if we have an integrator that may not conserve momenta, skip */
378 tm = vcm->group_mass[g];
379 fprintf(fp, "Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
380 vcm->group_name[g], tm, ekrot, det(vcm->group_i[g]));
381 fprintf(fp, " COM: %12.5f %12.5f %12.5f\n",
382 vcm->group_x[g][XX], vcm->group_x[g][YY], vcm->group_x[g][ZZ]);
383 fprintf(fp, " P: %12.5f %12.5f %12.5f\n",
384 vcm->group_p[g][XX], vcm->group_p[g][YY], vcm->group_p[g][ZZ]);
385 fprintf(fp, " V: %12.5f %12.5f %12.5f\n",
386 vcm->group_v[g][XX], vcm->group_v[g][YY], vcm->group_v[g][ZZ]);
387 fprintf(fp, " J: %12.5f %12.5f %12.5f\n",
388 vcm->group_j[g][XX], vcm->group_j[g][YY], vcm->group_j[g][ZZ]);
389 fprintf(fp, " w: %12.5f %12.5f %12.5f\n",
390 vcm->group_w[g][XX], vcm->group_w[g][YY], vcm->group_w[g][ZZ]);
391 pr_rvecs(fp, 0, "Inertia tensor", vcm->group_i[g], DIM);