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)
60 gmx_fatal(FARGS,"Can not have angular comm removal with pbc=%s",
61 epbc_names[ir->ePBC]);
63 if (vcm->mode != ecmNO) {
64 vcm->nr = groups->grps[egcVCM].nr;
65 /* Allocate one extra for a possible rest group */
66 if (vcm->mode == ecmANGULAR) {
67 snew(vcm->group_j,vcm->nr+1);
68 snew(vcm->group_x,vcm->nr+1);
69 snew(vcm->group_i,vcm->nr+1);
70 snew(vcm->group_w,vcm->nr+1);
72 snew(vcm->group_p,vcm->nr+1);
73 snew(vcm->group_v,vcm->nr+1);
74 snew(vcm->group_mass,vcm->nr+1);
75 snew(vcm->group_name,vcm->nr);
76 snew(vcm->group_ndf,vcm->nr);
77 for(g=0; (g<vcm->nr); g++)
78 vcm->group_ndf[g] = ir->opts.nrdf[g];
80 /* Copy pointer to group names and print it. */
82 fprintf(fp,"Center of mass motion removal mode is %s\n",
84 fprintf(fp,"We have the following groups for center of"
85 " mass motion removal:\n");
87 for(g=0; (g<vcm->nr); g++) {
88 vcm->group_name[g] = *groups->grpname[groups->grps[egcVCM].nm_ind[g]];
90 fprintf(fp,"%3d: %s\n",g,vcm->group_name[g]);
97 static void update_tensor(rvec x,real m0,tensor I)
101 /* Compute inertia tensor contribution due to this atom */
105 I[XX][XX] += x[XX]*x[XX]*m0;
106 I[YY][YY] += x[YY]*x[YY]*m0;
107 I[ZZ][ZZ] += x[ZZ]*x[ZZ]*m0;
116 /* Center of mass code for groups */
117 void calc_vcm_grp(FILE *fp,int start,int homenr,t_mdatoms *md,
118 rvec x[],rvec v[],t_vcm *vcm)
121 real m0,xx,xy,xz,yy,yz,zz;
124 if (vcm->mode != ecmNO) {
125 /* Also clear a possible rest group */
126 for(g=0; (g<vcm->nr+1); g++) {
127 /* Reset linear momentum */
128 vcm->group_mass[g] = 0;
129 clear_rvec(vcm->group_p[g]);
131 if (vcm->mode == ecmANGULAR) {
132 /* Reset anular momentum */
133 clear_rvec(vcm->group_j[g]);
134 clear_rvec(vcm->group_x[g]);
135 clear_rvec(vcm->group_w[g]);
136 clear_mat(vcm->group_i[g]);
141 for(i=start; (i<start+homenr); i++) {
146 /* Calculate linear momentum */
147 vcm->group_mass[g] += m0;
148 for(m=0; (m<DIM);m++)
149 vcm->group_p[g][m] += m0*v[i][m];
151 if (vcm->mode == ecmANGULAR) {
152 /* Calculate angular momentum */
155 for(m=0; (m<DIM); m++) {
156 vcm->group_j[g][m] += m0*j0[m];
157 vcm->group_x[g][m] += m0*x[i][m];
159 /* Update inertia tensor */
160 update_tensor(x[i],m0,vcm->group_i[g]);
166 void do_stopcm_grp(FILE *fp,int start,int homenr,unsigned short *group_id,
167 rvec x[],rvec v[],t_vcm *vcm)
173 if (vcm->mode != ecmNO) {
174 /* Subtract linear momentum */
178 for(i=start; (i<start+homenr); i++) {
181 v[i][XX] -= vcm->group_v[g][XX];
185 for(i=start; (i<start+homenr); i++) {
188 v[i][XX] -= vcm->group_v[g][XX];
189 v[i][YY] -= vcm->group_v[g][YY];
193 for(i=start; (i<start+homenr); i++) {
196 rvec_dec(v[i],vcm->group_v[g]);
200 if (vcm->mode == ecmANGULAR) {
201 /* Subtract angular momentum */
202 for(i=start; (i<start+homenr); i++) {
205 /* Compute the correction to the velocity for each atom */
206 rvec_sub(x[i],vcm->group_x[g],dx);
207 cprod(vcm->group_w[g],dx,dv);
214 static void get_minv(tensor A,tensor B)
220 tmp[XX][XX] = A[YY][YY] + A[ZZ][ZZ];
221 tmp[YY][XX] = -A[XX][YY];
222 tmp[ZZ][XX] = -A[XX][ZZ];
223 tmp[XX][YY] = -A[XX][YY];
224 tmp[YY][YY] = A[XX][XX] + A[ZZ][ZZ];
225 tmp[ZZ][YY] = -A[YY][ZZ];
226 tmp[XX][ZZ] = -A[XX][ZZ];
227 tmp[YY][ZZ] = -A[YY][ZZ];
228 tmp[ZZ][ZZ] = A[XX][XX] + A[YY][YY];
230 /* This is a hack to prevent very large determinants */
231 rfac = (tmp[XX][XX]+tmp[YY][YY]+tmp[ZZ][ZZ])/3;
233 gmx_fatal(FARGS,"Can not stop center of mass: maybe 2dimensional system");
235 for(m=0; (m<DIM); m++)
236 for(n=0; (n<DIM); n++)
239 for(m=0; (m<DIM); m++)
240 for(n=0; (n<DIM); n++)
244 void check_cm_grp(FILE *fp,t_vcm *vcm,t_inputrec *ir,real Temp_Max)
247 real ekcm,ekrot,tm,tm_1,Temp_cm;
251 /* First analyse the total results */
252 if (vcm->mode != ecmNO) {
253 for(g=0; (g<vcm->nr); g++) {
254 tm = vcm->group_mass[g];
257 svmul(tm_1,vcm->group_p[g],vcm->group_v[g]);
259 /* Else it's zero anyway! */
261 if (vcm->mode == ecmANGULAR) {
262 for(g=0; (g<vcm->nr); g++) {
263 tm = vcm->group_mass[g];
267 /* Compute center of mass for this group */
268 for(m=0; (m<DIM); m++)
269 vcm->group_x[g][m] *= tm_1;
271 /* Subtract the center of mass contribution to the
274 cprod(vcm->group_x[g],vcm->group_v[g],jcm);
275 for(m=0; (m<DIM); m++)
276 vcm->group_j[g][m] -= tm*jcm[m];
278 /* Subtract the center of mass contribution from the inertia
279 * tensor (this is not as trivial as it seems, but due to
280 * some cancellation we can still do it, even in parallel).
283 update_tensor(vcm->group_x[g],tm,Icm);
284 m_sub(vcm->group_i[g],Icm,vcm->group_i[g]);
286 /* Compute angular velocity, using matrix operation
291 get_minv(vcm->group_i[g],Icm);
292 mvmul(Icm,vcm->group_j[g],vcm->group_w[g]);
294 /* Else it's zero anyway! */
298 for(g=0; (g<vcm->nr); g++) {
300 if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0) {
301 for(m=0; m<vcm->ndim; m++)
302 ekcm += sqr(vcm->group_v[g][m]);
303 ekcm *= 0.5*vcm->group_mass[g];
304 Temp_cm = 2*ekcm/vcm->group_ndf[g];
306 if ((Temp_cm > Temp_Max) && fp)
307 fprintf(fp,"Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
308 vcm->group_name[g],vcm->group_v[g][XX],
309 vcm->group_v[g][YY],vcm->group_v[g][ZZ],Temp_cm);
311 if (vcm->mode == ecmANGULAR) {
312 ekrot = 0.5*iprod(vcm->group_j[g],vcm->group_w[g]);
313 if ((ekrot > 1) && fp && !EI_RANDOM(ir->eI)) {
314 /* if we have an integrator that may not conserve momenta, skip */
315 tm = vcm->group_mass[g];
316 fprintf(fp,"Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
317 vcm->group_name[g],tm,ekrot,det(vcm->group_i[g]));
318 fprintf(fp," COM: %12.5f %12.5f %12.5f\n",
319 vcm->group_x[g][XX],vcm->group_x[g][YY],vcm->group_x[g][ZZ]);
320 fprintf(fp," P: %12.5f %12.5f %12.5f\n",
321 vcm->group_p[g][XX],vcm->group_p[g][YY],vcm->group_p[g][ZZ]);
322 fprintf(fp," V: %12.5f %12.5f %12.5f\n",
323 vcm->group_v[g][XX],vcm->group_v[g][YY],vcm->group_v[g][ZZ]);
324 fprintf(fp," J: %12.5f %12.5f %12.5f\n",
325 vcm->group_j[g][XX],vcm->group_j[g][YY],vcm->group_j[g][ZZ]);
326 fprintf(fp," w: %12.5f %12.5f %12.5f\n",
327 vcm->group_w[g][XX],vcm->group_w[g][YY],vcm->group_w[g][ZZ]);
328 pr_rvecs(fp,0,"Inertia tensor",vcm->group_i[g],DIM);