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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_vcm_c = "$Id$";
  30
  31 #include "macros.h"
  32 #include "vcm.h"
  33 #include "vec.h"
  34 #include "smalloc.h"
  35 #include "do_fit.h"
  36
  37 void calc_vcm(FILE *log,int homenr,int start,real mass[],rvec v[],rvec vcm)
  38 {
  39   int    i;
  40   real   m0;
  41   real   x,y,z;
  42
  43   /* Calculate */
  44   x=y=z=0;
  45   for(i=start; (i<start+homenr); i++) {
  46     m0=mass[i];
  47     x+=m0*v[i][XX];
  48     y+=m0*v[i][YY];
  49     z+=m0*v[i][ZZ];
  50   }
  51   vcm[XX]=x;
  52   vcm[YY]=y;
  53   vcm[ZZ]=z;
  54 }
  55
  56 void do_stopcm(FILE *log,int homenr,int start,rvec v[],rvec mvcm,
  57                real tm,real invmass[])
  58 {
  59   int  i;
  60   rvec vcm;
  61
  62   vcm[XX]=mvcm[XX]/tm;
  63   vcm[YY]=mvcm[YY]/tm;
  64   vcm[ZZ]=mvcm[ZZ]/tm;
  65
  66   for(i=start; (i<start+homenr); i++)
  67     if (invmass[i] != 0)
  68       rvec_dec(v[i],vcm);
  69 }
  70
  71 void check_cm(FILE *log,rvec mvcm,real tm)
  72 {
  73   int    m;
  74   rvec   vcm;
  75   real   ekcm=0,max_vcm=0;
  76
  77   for(m=0; (m<DIM); m++) {
  78     vcm[m]=mvcm[m]/tm;
  79     max_vcm=max(max_vcm,fabs(vcm[m]));
  80     ekcm+=vcm[m]*vcm[m];
  81   }
  82   if (max_vcm > 0.1) {
  83     ekcm*=0.5*tm;
  84     fprintf(log,"Large VCM: (%12.5f,  %12.5f,  %12.5f), ekin-cm: %12.5f\n",
  85             vcm[XX],vcm[YY],vcm[ZZ],ekcm);
  86   }
  87 }
  88
  89 void do_stoprot(FILE *log, int natoms, rvec box, rvec x[], real mass[])
  90 {
  91   static atom_id *index=NULL;
  92   static rvec *old_x=NULL;
  93   rvec   half_box;
  94   int    i;
  95
  96   for (i=0; (i<DIM); i++)
  97     half_box[i]=box[i]*0.5;
  98   if (!index) {
  99     snew(index,natoms);
 100     for (i=0; (i<natoms); i++)
 101       index[i]=i;
 102   }
 103
 104   /* first center on origin (do_fit does not do this!) */
 105   reset_x(natoms,index,natoms,index,x,mass);
 106
 107   if (old_x) {
 108     /* do least squares fit to previous structure */
 109     do_fit(natoms,mass,old_x,x);
 110   } else {
 111     /* no previous structure, make room to copy this one */
 112     snew(old_x, natoms);
 113   }
 114   /* save structure for next fit and move structure back to center of box */
 115   for (i=0; (i<natoms); i++) {
 116     copy_rvec(x[i],old_x[i]);
 117     rvec_inc(x[i],half_box);
 118   }
 119 }
 120
 121
 122