src/mdlib/vcm.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROwing Monsters And Cloning Shrimps
  34  */
  35 /* This file is completely threadsafe - keep it that way! */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include "macros.h"
  41 #include "vcm.h"
  42 #include "vec.h"
  43 #include "smalloc.h"
  44 #include "names.h"
  45 #include "txtdump.h"
  46 #include "network.h"
  47 #include "pbc.h"
  48
  49 t_vcm *init_vcm(FILE *fp,gmx_groups_t *groups,t_inputrec *ir)
  50 {
  51   t_vcm *vcm;
  52   int   g;
  53
  54   snew(vcm,1);
  55
  56   vcm->mode = (ir->nstcomm > 0) ? ir->comm_mode : ecmNO;
  57   vcm->ndim = ndof_com(ir);
  58
  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]);
  62
  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);
  71     }
  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];
  79
  80     /* Copy pointer to group names and print it. */
  81     if (fp) {
  82       fprintf(fp,"Center of mass motion removal mode is %s\n",
  83               ECOM(vcm->mode));
  84       fprintf(fp,"We have the following groups for center of"
  85             " mass motion removal:\n");
  86     }
  87     for(g=0; (g<vcm->nr); g++) {
  88       vcm->group_name[g] = *groups->grpname[groups->grps[egcVCM].nm_ind[g]];
  89       if (fp)
  90         fprintf(fp,"%3d:  %s\n",g,vcm->group_name[g]);
  91     }
  92   }
  93
  94   return vcm;
  95 }
  96
  97 static void update_tensor(rvec x,real m0,tensor I)
  98 {
  99   real xy,xz,yz;
 100
 101   /* Compute inertia tensor contribution due to this atom */
 102   xy         = x[XX]*x[YY]*m0;
 103   xz         = x[XX]*x[ZZ]*m0;
 104   yz         = x[YY]*x[ZZ]*m0;
 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;
 108   I[XX][YY] += xy;
 109   I[YY][XX] += xy;
 110   I[XX][ZZ] += xz;
 111   I[ZZ][XX] += xz;
 112   I[YY][ZZ] += yz;
 113   I[ZZ][YY] += yz;
 114 }
 115
 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)
 119 {
 120   int    i,g,m;
 121   real   m0,xx,xy,xz,yy,yz,zz;
 122   rvec   j0;
 123
 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]);
 130
 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]);
 137       }
 138     }
 139
 140     g = 0;
 141     for(i=start; (i<start+homenr); i++) {
 142       m0 = md->massT[i];
 143       if (md->cVCM)
 144         g = md->cVCM[i];
 145
 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];
 150
 151       if (vcm->mode == ecmANGULAR) {
 152         /* Calculate angular momentum */
 153         cprod(x[i],v[i],j0);
 154
 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];
 158         }
 159         /* Update inertia tensor */
 160         update_tensor(x[i],m0,vcm->group_i[g]);
 161       }
 162     }
 163   }
 164 }
 165
 166 void do_stopcm_grp(FILE *fp,int start,int homenr,unsigned short *group_id,
 167                    rvec x[],rvec v[],t_vcm *vcm)
 168 {
 169   int  i,g,m;
 170   real tm,tm_1;
 171   rvec dv,dx;
 172
 173   if (vcm->mode != ecmNO) {
 174     /* Subtract linear momentum */
 175     g = 0;
 176     switch (vcm->ndim) {
 177     case 1:
 178       for(i=start; (i<start+homenr); i++) {
 179         if (group_id)
 180           g = group_id[i];
 181         v[i][XX] -= vcm->group_v[g][XX];
 182       }
 183       break;
 184     case 2:
 185       for(i=start; (i<start+homenr); i++) {
 186         if (group_id)
 187           g = group_id[i];
 188         v[i][XX] -= vcm->group_v[g][XX];
 189         v[i][YY] -= vcm->group_v[g][YY];
 190       }
 191       break;
 192     case 3:
 193       for(i=start; (i<start+homenr); i++) {
 194         if (group_id)
 195           g = group_id[i];
 196         rvec_dec(v[i],vcm->group_v[g]);
 197       }
 198       break;
 199     }
 200     if (vcm->mode == ecmANGULAR) {
 201       /* Subtract angular momentum */
 202       for(i=start; (i<start+homenr); i++) {
 203         if (group_id)
 204           g = group_id[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);
 208         rvec_dec(v[i],dv);
 209       }
 210     }
 211   }
 212 }
 213
 214 static void get_minv(tensor A,tensor B)
 215 {
 216   int    m,n;
 217   double fac,rfac;
 218   tensor tmp;
 219
 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];
 229
 230   /* This is a hack to prevent very large determinants */
 231   rfac  = (tmp[XX][XX]+tmp[YY][YY]+tmp[ZZ][ZZ])/3;
 232   if (rfac == 0.0)
 233     gmx_fatal(FARGS,"Can not stop center of mass: maybe 2dimensional system");
 234   fac = 1.0/rfac;
 235   for(m=0; (m<DIM); m++)
 236     for(n=0; (n<DIM); n++)
 237       tmp[m][n] *= fac;
 238   m_inv(tmp,B);
 239   for(m=0; (m<DIM); m++)
 240     for(n=0; (n<DIM); n++)
 241       B[m][n] *= fac;
 242 }
 243
 244 void check_cm_grp(FILE *fp,t_vcm *vcm,t_inputrec *ir,real Temp_Max)
 245 {
 246   int    m,g;
 247   real   ekcm,ekrot,tm,tm_1,Temp_cm;
 248   rvec   jcm;
 249   tensor Icm,Tcm;
 250
 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];
 255       if (tm != 0) {
 256         tm_1 = 1.0/tm;
 257         svmul(tm_1,vcm->group_p[g],vcm->group_v[g]);
 258       }
 259       /* Else it's zero anyway! */
 260     }
 261     if (vcm->mode == ecmANGULAR) {
 262       for(g=0; (g<vcm->nr); g++) {
 263         tm = vcm->group_mass[g];
 264         if (tm != 0) {
 265           tm_1 = 1.0/tm;
 266
 267           /* Compute center of mass for this group */
 268           for(m=0; (m<DIM); m++)
 269             vcm->group_x[g][m] *= tm_1;
 270
 271           /* Subtract the center of mass contribution to the
 272            * angular momentum
 273            */
 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];
 277
 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).
 281            */
 282           clear_mat(Icm);
 283           update_tensor(vcm->group_x[g],tm,Icm);
 284           m_sub(vcm->group_i[g],Icm,vcm->group_i[g]);
 285
 286           /* Compute angular velocity, using matrix operation
 287            * Since J = I w
 288            * we have
 289            * w = I^-1 J
 290            */
 291           get_minv(vcm->group_i[g],Icm);
 292           mvmul(Icm,vcm->group_j[g],vcm->group_w[g]);
 293         }
 294         /* Else it's zero anyway! */
 295       }
 296     }
 297   }
 298   for(g=0; (g<vcm->nr); g++) {
 299     ekcm    = 0;
 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];
 305
 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);
 310
 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);
 329         }
 330       }
 331     }
 332   }
 333 }
 334
 335