src/mdlib/vcm.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include "macros.h"
  44 #include "vcm.h"
  45 #include "vec.h"
  46 #include "smalloc.h"
  47 #include "names.h"
  48 #include "txtdump.h"
  49 #include "network.h"
  50 #include "pbc.h"
  51
  52 t_vcm *init_vcm(FILE *fp,gmx_groups_t *groups,t_inputrec *ir)
  53 {
  54   t_vcm *vcm;
  55   int   g;
  56
  57   snew(vcm,1);
  58
  59   vcm->mode = (ir->nstcomm > 0) ? ir->comm_mode : ecmNO;
  60   vcm->ndim = ndof_com(ir);
  61
  62   if (vcm->mode == ecmANGULAR && vcm->ndim < 3)
  63     gmx_fatal(FARGS,"Can not have angular comm removal with pbc=%s",
  64               epbc_names[ir->ePBC]);
  65
  66   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) {
  70       snew(vcm->group_j,vcm->nr+1);
  71       snew(vcm->group_x,vcm->nr+1);
  72       snew(vcm->group_i,vcm->nr+1);
  73       snew(vcm->group_w,vcm->nr+1);
  74     }
  75     snew(vcm->group_p,vcm->nr+1);
  76     snew(vcm->group_v,vcm->nr+1);
  77     snew(vcm->group_mass,vcm->nr+1);
  78     snew(vcm->group_name,vcm->nr);
  79     snew(vcm->group_ndf,vcm->nr);
  80     for(g=0; (g<vcm->nr); g++)
  81       vcm->group_ndf[g] = ir->opts.nrdf[g];
  82
  83     /* Copy pointer to group names and print it. */
  84     if (fp) {
  85       fprintf(fp,"Center of mass motion removal mode is %s\n",
  86               ECOM(vcm->mode));
  87       fprintf(fp,"We have the following groups for center of"
  88             " mass motion removal:\n");
  89     }
  90     for(g=0; (g<vcm->nr); g++) {
  91       vcm->group_name[g] = *groups->grpname[groups->grps[egcVCM].nm_ind[g]];
  92       if (fp)
  93         fprintf(fp,"%3d:  %s\n",g,vcm->group_name[g]);
  94     }
  95   }
  96
  97   return vcm;
  98 }
  99
 100 static void update_tensor(rvec x,real m0,tensor I)
 101 {
 102   real xy,xz,yz;
 103
 104   /* Compute inertia tensor contribution due to this atom */
 105   xy         = x[XX]*x[YY]*m0;
 106   xz         = x[XX]*x[ZZ]*m0;
 107   yz         = x[YY]*x[ZZ]*m0;
 108   I[XX][XX] += x[XX]*x[XX]*m0;
 109   I[YY][YY] += x[YY]*x[YY]*m0;
 110   I[ZZ][ZZ] += x[ZZ]*x[ZZ]*m0;
 111   I[XX][YY] += xy;
 112   I[YY][XX] += xy;
 113   I[XX][ZZ] += xz;
 114   I[ZZ][XX] += xz;
 115   I[YY][ZZ] += yz;
 116   I[ZZ][YY] += yz;
 117 }
 118
 119 /* Center of mass code for groups */
 120 void calc_vcm_grp(FILE *fp,int start,int homenr,t_mdatoms *md,
 121                   rvec x[],rvec v[],t_vcm *vcm)
 122 {
 123   int    i,g,m;
 124   real   m0,xx,xy,xz,yy,yz,zz;
 125   rvec   j0;
 126
 127   if (vcm->mode != ecmNO) {
 128     /* Also clear a possible rest group */
 129     for(g=0; (g<vcm->nr+1); g++) {
 130       /* Reset linear momentum */
 131       vcm->group_mass[g] = 0;
 132       clear_rvec(vcm->group_p[g]);
 133
 134       if (vcm->mode == ecmANGULAR) {
 135         /* Reset anular momentum */
 136         clear_rvec(vcm->group_j[g]);
 137         clear_rvec(vcm->group_x[g]);
 138         clear_rvec(vcm->group_w[g]);
 139         clear_mat(vcm->group_i[g]);
 140       }
 141     }
 142
 143     g = 0;
 144     for(i=start; (i<start+homenr); i++) {
 145       m0 = md->massT[i];
 146       if (md->cVCM)
 147         g = md->cVCM[i];
 148
 149       /* Calculate linear momentum */
 150       vcm->group_mass[g]  += m0;
 151       for(m=0; (m<DIM);m++)
 152         vcm->group_p[g][m] += m0*v[i][m];
 153
 154       if (vcm->mode == ecmANGULAR) {
 155         /* Calculate angular momentum */
 156         cprod(x[i],v[i],j0);
 157
 158         for(m=0; (m<DIM); m++) {
 159           vcm->group_j[g][m] += m0*j0[m];
 160           vcm->group_x[g][m] += m0*x[i][m];
 161         }
 162         /* Update inertia tensor */
 163         update_tensor(x[i],m0,vcm->group_i[g]);
 164       }
 165     }
 166   }
 167 }
 168
 169 void do_stopcm_grp(FILE *fp,int start,int homenr,unsigned short *group_id,
 170                    rvec x[],rvec v[],t_vcm *vcm)
 171 {
 172   int  i,g,m;
 173   real tm,tm_1;
 174   rvec dv,dx;
 175
 176   if (vcm->mode != ecmNO) {
 177     /* Subtract linear momentum */
 178     g = 0;
 179     switch (vcm->ndim) {
 180     case 1:
 181       for(i=start; (i<start+homenr); i++) {
 182         if (group_id)
 183           g = group_id[i];
 184         v[i][XX] -= vcm->group_v[g][XX];
 185       }
 186       break;
 187     case 2:
 188       for(i=start; (i<start+homenr); i++) {
 189         if (group_id)
 190           g = group_id[i];
 191         v[i][XX] -= vcm->group_v[g][XX];
 192         v[i][YY] -= vcm->group_v[g][YY];
 193       }
 194       break;
 195     case 3:
 196       for(i=start; (i<start+homenr); i++) {
 197         if (group_id)
 198           g = group_id[i];
 199         rvec_dec(v[i],vcm->group_v[g]);
 200       }
 201       break;
 202     }
 203     if (vcm->mode == ecmANGULAR) {
 204       /* Subtract angular momentum */
 205       for(i=start; (i<start+homenr); i++) {
 206         if (group_id)
 207           g = group_id[i];
 208         /* Compute the correction to the velocity for each atom */
 209         rvec_sub(x[i],vcm->group_x[g],dx);
 210         cprod(vcm->group_w[g],dx,dv);
 211         rvec_dec(v[i],dv);
 212       }
 213     }
 214   }
 215 }
 216
 217 static void get_minv(tensor A,tensor B)
 218 {
 219   int    m,n;
 220   double fac,rfac;
 221   tensor tmp;
 222
 223   tmp[XX][XX] =  A[YY][YY] + A[ZZ][ZZ];
 224   tmp[YY][XX] = -A[XX][YY];
 225   tmp[ZZ][XX] = -A[XX][ZZ];
 226   tmp[XX][YY] = -A[XX][YY];
 227   tmp[YY][YY] =  A[XX][XX] + A[ZZ][ZZ];
 228   tmp[ZZ][YY] = -A[YY][ZZ];
 229   tmp[XX][ZZ] = -A[XX][ZZ];
 230   tmp[YY][ZZ] = -A[YY][ZZ];
 231   tmp[ZZ][ZZ] =  A[XX][XX] + A[YY][YY];
 232
 233   /* This is a hack to prevent very large determinants */
 234   rfac  = (tmp[XX][XX]+tmp[YY][YY]+tmp[ZZ][ZZ])/3;
 235   if (rfac == 0.0)
 236     gmx_fatal(FARGS,"Can not stop center of mass: maybe 2dimensional system");
 237   fac = 1.0/rfac;
 238   for(m=0; (m<DIM); m++)
 239     for(n=0; (n<DIM); n++)
 240       tmp[m][n] *= fac;
 241   m_inv(tmp,B);
 242   for(m=0; (m<DIM); m++)
 243     for(n=0; (n<DIM); n++)
 244       B[m][n] *= fac;
 245 }
 246
 247 void check_cm_grp(FILE *fp,t_vcm *vcm,t_inputrec *ir,real Temp_Max)
 248 {
 249   int    m,g;
 250   real   ekcm,ekrot,tm,tm_1,Temp_cm;
 251   rvec   jcm;
 252   tensor Icm,Tcm;
 253
 254   /* First analyse the total results */
 255   if (vcm->mode != ecmNO) {
 256     for(g=0; (g<vcm->nr); g++) {
 257       tm = vcm->group_mass[g];
 258       if (tm != 0) {
 259         tm_1 = 1.0/tm;
 260         svmul(tm_1,vcm->group_p[g],vcm->group_v[g]);
 261       }
 262       /* Else it's zero anyway! */
 263     }
 264     if (vcm->mode == ecmANGULAR) {
 265       for(g=0; (g<vcm->nr); g++) {
 266         tm = vcm->group_mass[g];
 267         if (tm != 0) {
 268           tm_1 = 1.0/tm;
 269
 270           /* Compute center of mass for this group */
 271           for(m=0; (m<DIM); m++)
 272             vcm->group_x[g][m] *= tm_1;
 273
 274           /* Subtract the center of mass contribution to the
 275            * angular momentum
 276            */
 277           cprod(vcm->group_x[g],vcm->group_v[g],jcm);
 278           for(m=0; (m<DIM); m++)
 279             vcm->group_j[g][m] -= tm*jcm[m];
 280
 281           /* Subtract the center of mass contribution from the inertia
 282            * tensor (this is not as trivial as it seems, but due to
 283            * some cancellation we can still do it, even in parallel).
 284            */
 285           clear_mat(Icm);
 286           update_tensor(vcm->group_x[g],tm,Icm);
 287           m_sub(vcm->group_i[g],Icm,vcm->group_i[g]);
 288
 289           /* Compute angular velocity, using matrix operation
 290            * Since J = I w
 291            * we have
 292            * w = I^-1 J
 293            */
 294           get_minv(vcm->group_i[g],Icm);
 295           mvmul(Icm,vcm->group_j[g],vcm->group_w[g]);
 296         }
 297         /* Else it's zero anyway! */
 298       }
 299     }
 300   }
 301   for(g=0; (g<vcm->nr); g++) {
 302     ekcm    = 0;
 303     if (vcm->group_mass[g] != 0 && vcm->group_ndf[g] > 0) {
 304       for(m=0; m<vcm->ndim; m++)
 305         ekcm += sqr(vcm->group_v[g][m]);
 306       ekcm *= 0.5*vcm->group_mass[g];
 307       Temp_cm = 2*ekcm/vcm->group_ndf[g];
 308
 309       if ((Temp_cm > Temp_Max) && fp)
 310         fprintf(fp,"Large VCM(group %s): %12.5f, %12.5f, %12.5f, Temp-cm: %12.5e\n",
 311                 vcm->group_name[g],vcm->group_v[g][XX],
 312                 vcm->group_v[g][YY],vcm->group_v[g][ZZ],Temp_cm);
 313
 314       if (vcm->mode == ecmANGULAR) {
 315         ekrot = 0.5*iprod(vcm->group_j[g],vcm->group_w[g]);
 316         if ((ekrot > 1) && fp && !EI_RANDOM(ir->eI)) {
 317           /* if we have an integrator that may not conserve momenta, skip */
 318           tm    = vcm->group_mass[g];
 319           fprintf(fp,"Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
 320                   vcm->group_name[g],tm,ekrot,det(vcm->group_i[g]));
 321           fprintf(fp,"  COM: %12.5f  %12.5f  %12.5f\n",
 322                   vcm->group_x[g][XX],vcm->group_x[g][YY],vcm->group_x[g][ZZ]);
 323           fprintf(fp,"  P:   %12.5f  %12.5f  %12.5f\n",
 324                   vcm->group_p[g][XX],vcm->group_p[g][YY],vcm->group_p[g][ZZ]);
 325           fprintf(fp,"  V:   %12.5f  %12.5f  %12.5f\n",
 326                   vcm->group_v[g][XX],vcm->group_v[g][YY],vcm->group_v[g][ZZ]);
 327           fprintf(fp,"  J:   %12.5f  %12.5f  %12.5f\n",
 328                   vcm->group_j[g][XX],vcm->group_j[g][YY],vcm->group_j[g][ZZ]);
 329           fprintf(fp,"  w:   %12.5f  %12.5f  %12.5f\n",
 330                 vcm->group_w[g][XX],vcm->group_w[g][YY],vcm->group_w[g][ZZ]);
 331           pr_rvecs(fp,0,"Inertia tensor",vcm->group_i[g],DIM);
 332         }
 333       }
 334     }
 335   }
 336 }
 337
 338