src/mdlib/vcm.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 /* This file is completely threadsafe - keep it that way! */
  37
  38 #include "macros.h"
  39 #include "vcm.h"
  40 #include "vec.h"
  41 #include "smalloc.h"
  42 #include "names.h"
  43 #include "txtdump.h"
  44 #include "network.h"
  45
  46 t_vcm *init_vcm(FILE *fp,t_topology *top,t_commrec *cr,t_mdatoms *md,
  47                 int start,int homenr,int nstcomm)
  48 {
  49   t_vcm *vcm;
  50   real  *mass;
  51   int   i,g;
  52
  53   snew(vcm,1);
  54
  55   if (nstcomm < 0)
  56     vcm->mode = ecmANGULAR;
  57   else if (nstcomm > 0)
  58     vcm->mode = ecmLINEAR;
  59   else
  60     vcm->mode = ecmNO;
  61
  62   if (vcm->mode != ecmNO) {
  63     vcm->nr = top->atoms.grps[egcVCM].nr;
  64     if (vcm->mode == ecmANGULAR) {
  65       snew(vcm->group_j,vcm->nr);
  66       snew(vcm->group_x,vcm->nr);
  67       snew(vcm->group_i,vcm->nr);
  68       snew(vcm->group_w,vcm->nr);
  69     }
  70     snew(vcm->group_p,vcm->nr);
  71     snew(vcm->group_v,vcm->nr);
  72     snew(vcm->group_mass,vcm->nr);
  73     snew(vcm->group_name,vcm->nr);
  74     vcm->group_id = md->cVCM;
  75     for(i=start; (i<start+homenr); i++) {
  76       g = vcm->group_id[i];
  77       vcm->group_mass[g] += md->massT[i];
  78     }
  79     snew(mass,vcm->nr);
  80     for(g=0; (g<vcm->nr); g++)
  81       mass[g] = vcm->group_mass[g];
  82     if(PAR(cr))
  83       gmx_sum(vcm->nr,mass,cr);
  84
  85     /* Copy pointer to group names and print it. */
  86     fprintf(fp,"Center of mass motion removal mode is %s\n",ECOM(vcm->mode));
  87     fprintf(fp,"We have the following groups for center of"
  88             " mass motion removal:\n");
  89     for(g=0; (g<vcm->nr); g++) {
  90       vcm->group_name[g] =
  91         *top->atoms.grpname[top->atoms.grps[egcVCM].nm_ind[g]];
  92       fprintf(fp,"%3d:  %s, initial mass: %g\n",
  93               g,vcm->group_name[g],mass[g]);
  94     }
  95     sfree(mass);
  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,real mass[],
 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     for(g=0; (g<vcm->nr); g++) {
 129       /* Reset linear momentum */
 130       vcm->group_mass[g] = 0;
 131       clear_rvec(vcm->group_p[g]);
 132
 133       if (vcm->mode == ecmANGULAR) {
 134         /* Reset anular momentum */
 135         clear_rvec(vcm->group_j[g]);
 136         clear_rvec(vcm->group_x[g]);
 137         clear_rvec(vcm->group_w[g]);
 138         clear_mat(vcm->group_i[g]);
 139       }
 140     }
 141
 142     for(i=start; (i<start+homenr); i++) {
 143       m0 = mass[i];
 144       g  = vcm->group_id[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         oprod(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,rvec x[],rvec v[],
 167                    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     for(i=start; (i<start+homenr); i++) {
 176       g = vcm->group_id[i];
 177       rvec_dec(v[i],vcm->group_v[g]);
 178     }
 179     if (vcm->mode == ecmANGULAR) {
 180       /* Subtract angular momentum */
 181       for(i=start; (i<start+homenr); i++) {
 182         g = vcm->group_id[i];
 183         /* Compute the correction to the velocity for each atom */
 184         rvec_sub(x[i],vcm->group_x[g],dx);
 185         oprod(vcm->group_w[g],dx,dv);
 186         rvec_dec(v[i],dv);
 187       }
 188     }
 189   }
 190 }
 191
 192 static void get_minv(tensor A,tensor B)
 193 {
 194   int    m,n;
 195   double fac,rfac;
 196   tensor tmp;
 197
 198   tmp[XX][XX] =  A[YY][YY] + A[ZZ][ZZ];
 199   tmp[YY][XX] = -A[XX][YY];
 200   tmp[ZZ][XX] = -A[XX][ZZ];
 201   tmp[XX][YY] = -A[XX][YY];
 202   tmp[YY][YY] =  A[XX][XX] + A[ZZ][ZZ];
 203   tmp[ZZ][YY] = -A[YY][ZZ];
 204   tmp[XX][ZZ] = -A[XX][ZZ];
 205   tmp[YY][ZZ] = -A[YY][ZZ];
 206   tmp[ZZ][ZZ] =  A[XX][XX] + A[YY][YY];
 207
 208   /* This is a hack to prevent very large determinants */
 209   rfac  = (tmp[XX][XX]+tmp[YY][YY]+tmp[ZZ][ZZ])/3;
 210   if (rfac == 0.0)
 211     fatal_error(0,"Can not stop center of mass: maybe 2dimensional system");
 212   fac = 1.0/rfac;
 213   for(m=0; (m<DIM); m++)
 214     for(n=0; (n<DIM); n++)
 215       tmp[m][n] *= fac;
 216   m_inv(tmp,B);
 217   for(m=0; (m<DIM); m++)
 218     for(n=0; (n<DIM); n++)
 219       B[m][n] *= fac;
 220 }
 221
 222 void check_cm_grp(FILE *fp,t_vcm *vcm)
 223 {
 224   int    m,g;
 225   real   ekcm,ekrot,tm,tm_1;
 226   rvec   jcm;
 227   tensor Icm,Tcm;
 228
 229   /* First analyse the total results */
 230   if (vcm->mode != ecmNO) {
 231     for(g=0; (g<vcm->nr); g++) {
 232       tm = vcm->group_mass[g];
 233       if (tm != 0) {
 234         tm_1 = 1.0/tm;
 235         svmul(tm_1,vcm->group_p[g],vcm->group_v[g]);
 236       }
 237       /* Else it's zero anyway! */
 238     }
 239     if (vcm->mode == ecmANGULAR) {
 240       for(g=0; (g<vcm->nr); g++) {
 241         tm = vcm->group_mass[g];
 242         if (tm != 0) {
 243           tm_1 = 1.0/tm;
 244
 245           /* Compute center of mass for this group */
 246           for(m=0; (m<DIM); m++)
 247             vcm->group_x[g][m] *= tm_1;
 248
 249           /* Subtract the center of mass contribution to the
 250            * angular momentum
 251            */
 252           oprod(vcm->group_x[g],vcm->group_v[g],jcm);
 253           for(m=0; (m<DIM); m++)
 254             vcm->group_j[g][m] -= tm*jcm[m];
 255
 256           /* Subtract the center of mass contribution from the inertia
 257            * tensor (this is not as trivial as it seems, but due to
 258            * some cancellation we can still do it, even in parallel).
 259            */
 260           clear_mat(Icm);
 261           update_tensor(vcm->group_x[g],tm,Icm);
 262           m_sub(vcm->group_i[g],Icm,vcm->group_i[g]);
 263
 264           /* Compute angular velocity, using matrix operation
 265            * Since J = I w
 266            * we have
 267            * w = I^-1 J
 268            */
 269           get_minv(vcm->group_i[g],Icm);
 270           mvmul(Icm,vcm->group_j[g],vcm->group_w[g]);
 271         }
 272         /* Else it's zero anyway! */
 273       }
 274     }
 275   }
 276   for(g=0; (g<vcm->nr); g++) {
 277     ekcm    = 0;
 278     if (vcm->group_mass[g] != 0) {
 279       for(m=0; (m<DIM); m++)
 280         ekcm += sqr(vcm->group_v[g][m]);
 281       ekcm *= 0.5*vcm->group_mass[g];
 282
 283       if ((ekcm > 1) || debug)
 284         fprintf(fp,"Large VCM(group %s): %12.5f, %12.5f, %12.5f, ekin-cm: %12.5e\n",
 285                 vcm->group_name[g],vcm->group_v[g][XX],
 286                 vcm->group_v[g][YY],vcm->group_v[g][ZZ],ekcm);
 287
 288       if (vcm->mode == ecmANGULAR) {
 289         ekrot = 0.5*iprod(vcm->group_j[g],vcm->group_w[g]);
 290         if ((ekrot > 1) || debug) {
 291           tm    = vcm->group_mass[g];
 292           fprintf(fp,"Group %s with mass %12.5e, Ekrot %12.5e Det(I) = %12.5e\n",
 293                   vcm->group_name[g],tm,ekrot,det(vcm->group_i[g]));
 294           fprintf(fp,"  COM: %12.5f  %12.5f  %12.5f\n",
 295                   vcm->group_x[g][XX],vcm->group_x[g][YY],vcm->group_x[g][ZZ]);
 296           fprintf(fp,"  P:   %12.5f  %12.5f  %12.5f\n",
 297                   vcm->group_p[g][XX],vcm->group_p[g][YY],vcm->group_p[g][ZZ]);
 298           fprintf(fp,"  V:   %12.5f  %12.5f  %12.5f\n",
 299                   vcm->group_v[g][XX],vcm->group_v[g][YY],vcm->group_v[g][ZZ]);
 300           fprintf(fp,"  J:   %12.5f  %12.5f  %12.5f\n",
 301                   vcm->group_j[g][XX],vcm->group_j[g][YY],vcm->group_j[g][ZZ]);
 302           fprintf(fp,"  w:   %12.5f  %12.5f  %12.5f\n",
 303                 vcm->group_w[g][XX],vcm->group_w[g][YY],vcm->group_w[g][ZZ]);
 304           pr_rvecs(fp,0,"Inertia tensor",vcm->group_i[g],DIM);
 305         }
 306       }
 307     }
 308   }
 309 }
 310
 311