src/mdlib/pullutil.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.0
  11  *
  12  * Copyright (c) 1991-2001
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  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  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 static char *SRCID_pullutil_c = "$Id$";
  37 #include <stdlib.h>
  38 #include "sysstuff.h"
  39 #include "princ.h"
  40 #include "futil.h"
  41 #include "statutil.h"
  42 #include "vec.h"
  43 #include "smalloc.h"
  44 #include "typedefs.h"
  45 #include "names.h"
  46 #include "fatal.h"
  47 #include "macros.h"
  48 #include "rdgroup.h"
  49 #include "symtab.h"
  50 #include "index.h"
  51 #include "confio.h"
  52 #include "pull.h"
  53 #include "pull_internal.h"
  54
  55 /* calculates center of mass of selection index from all coordines x */
  56 real calc_com(rvec x[],int gnx,atom_id *index,t_mdatoms *md,
  57                      rvec com,matrix box)
  58 {
  59   int  i,ii,m;
  60   real m0,tm;
  61
  62   clear_rvec(com);
  63   tm=0;
  64   for(i=0; (i<gnx); i++) {
  65     ii=index[i];
  66     m0=md->massT[ii];
  67     tm+=m0;
  68     for(m=0; (m<DIM); m++)
  69       com[m]+=m0*x[ii][m];
  70   }
  71   svmul(1/tm,com,com);
  72   for(m=DIM-1; m>=0; m--) {
  73     if (com[m] < 0        ) rvec_inc(com,box[m]);
  74     if (com[m] > box[m][m]) rvec_dec(com,box[m]);
  75   }
  76
  77   return tm;
  78 }
  79
  80 /* calculates com of all atoms in x[], *index has their index numbers
  81    to get the masses from atom[] */
  82 real calc_com2(rvec x[],int gnx,atom_id *index,t_mdatoms *md,rvec com,
  83                matrix box)
  84 {
  85   int  i,ii,m;
  86   real m0,tm;
  87
  88   clear_rvec(com);
  89   tm=0;
  90   for(i=0; (i<gnx); i++) {
  91     ii=index[i];
  92     m0=md->massT[ii];
  93     tm+=m0;
  94     for(m=0; (m<DIM); m++)
  95       com[m]+=m0*x[i][m];
  96   }
  97   svmul(1/tm,com,com);
  98   for(m=DIM-1; m>=0; m--) {
  99     /* next two lines used to be commented out */
 100     if (com[m] < 0        ) rvec_inc(com,box[m]);
 101     if (com[m] > box[m][m]) rvec_dec(com,box[m]);
 102   }
 103   return tm;
 104 }
 105
 106 void calc_running_com(t_pull *pull) {
 107   int i,j,n;
 108   rvec ave;
 109   real tm;
 110
 111   /* for group i, we have nhist[i] points of history. The maximum nr of points
 112      is pull->reflag. The array comhist[i][0..nhist[i]-1] has the positions
 113      of the center of mass over the last nhist[i] steps. x_unc[i] has the
 114      new ones. Remove the oldest one from the list, add the new one, calculate
 115      the average, put that in x_unc instead and return. We just move all coms
 116      1 down in the list and add the latest one to the top.
 117    */
 118
 119   if (pull->bCyl) {
 120     /* act on dyna groups */
 121     for (i=0;i<pull->pull.n;i++) {
 122       clear_rvec(ave);
 123       for (j=0;j<(pull->reflag-1);j++) {
 124         copy_rvec(pull->dyna.comhist[i][j+1],pull->dyna.comhist[i][j]);
 125         rvec_add(ave,pull->dyna.comhist[i][j],ave);
 126       }
 127       copy_rvec(pull->dyna.x_unc[i],pull->dyna.comhist[i][j]);
 128       rvec_add(ave,pull->dyna.x_unc[i],ave);
 129       svmul(1.0/pull->reflag,ave,ave);
 130
 131       /* now ave has the running com for group i, copy it to x_unc[i] */
 132       copy_rvec(ave,pull->dyna.x_unc[i]);
 133
 134 #ifdef DEBUG
 135       if (pull->bVerbose)
 136         for (n=0;n<pull->reflag;n++)
 137           fprintf(stderr,"Comhist %d, grp %d: %8.3f%8.3f%8.3f\n",n,i,
 138                   pull->dyna.comhist[i][n][XX],
 139                   pull->dyna.comhist[i][n][YY],
 140                   pull->dyna.comhist[i][n][ZZ]);
 141 #endif
 142     }
 143   } else {
 144     /* act on ref group */
 145     clear_rvec(ave);
 146
 147     for (j=0;j<(pull->reflag-1);j++) {
 148       copy_rvec(pull->ref.comhist[0][j+1],pull->ref.comhist[0][j]);
 149       rvec_add(ave,pull->ref.comhist[0][j],ave);
 150     }
 151
 152     copy_rvec(pull->ref.x_unc[0],pull->ref.comhist[0][j]);
 153     rvec_add(ave,pull->ref.x_unc[0],ave);
 154     svmul(1.0/pull->reflag,ave,ave);
 155     /* now ave has the running com for group i, copy it to x_unc[0] */
 156     copy_rvec(ave,pull->ref.x_unc[0]);
 157
 158 #ifdef DEBUG
 159     if (pull->bVerbose)
 160       for (i=0;i<pull->reflag;i++)
 161         fprintf(stderr,"Comhist %d: %8.3f%8.3f%8.3f\n",i,
 162                 pull->ref.comhist[0][i][XX],
 163                 pull->ref.comhist[0][i][YY],
 164                 pull->ref.comhist[0][i][ZZ]);
 165 #endif
 166   }
 167 }
 168
 169 void correct_t0_pbc(t_pull *pull, rvec x[], t_mdatoms *md, matrix box) {
 170   int i,ii,j,m;
 171   real tm;
 172   rvec com,dx;
 173
 174   /* loop over all atoms in index for group i. Check if they moved
 175      more than half a box with respect to xp. If so add/subtract a box
 176      from x0. Copy x to xp.
 177    */
 178   for (i=0;i<pull->ref.ngx[0];i++) {
 179     ii = pull->ref.idx[0][i];
 180
 181     /* correct for jumps across the box */
 182     rvec_sub(x[ii],pull->ref.xp[0][i],dx);
 183     for (m=DIM-1; m>=0; m--) {
 184       if (dx[m] < -0.5*box[m][m]) {
 185         rvec_inc(dx,box[m]);
 186         if (pull->bVerbose && pull->dims[m])
 187           fprintf(stderr,"Jumped +box: nr %d dir: %d old:%8.3f\n",ii,m,
 188                   pull->ref.x0[0][i][m]);
 189       }
 190
 191       if (dx[m] >  0.5*box[m][m]) {
 192         rvec_dec(dx,box[m]);
 193         if (pull->bVerbose && pull->dims[m])
 194           fprintf(stderr,"Jumped -box: nr %d dir: %d old:%8.3f\n",ii,m,
 195                   pull->ref.x0[0][i][m]);
 196       }
 197
 198       pull->ref.x0[0][i][m] += dx[m];
 199       pull->ref.xp[0][i][m]  = x[ii][m];
 200     }
 201   }
 202   tm = calc_com2(pull->ref.x0[0],pull->ref.ngx[0],pull->ref.idx[0],
 203                  md,com,box);
 204   if (pull->bVerbose)
 205     fprintf(stderr,"correct_t0: Group %s: mass:%8.3f com:%8.3f%8.3f%8.3f\n",
 206             pull->ref.grps[0],tm,com[0],com[1],com[2]);
 207 }
 208
 209 void string2rvec(char buf[], rvec nums) {
 210   float a,b,c;
 211   if (sscanf(buf,"%f%f%f",&a,&b,&c) != 3)
 212     fatal_error(0,"Expected three numbers at input line %s",buf);
 213   nums[0]=a; nums[1]=b; nums[2]=c;
 214 }
 215
 216 /* switch function, x between r and w */
 217 real get_weight(real x, real r, real w) {
 218   static bool bFirst = TRUE;
 219   static real rw, a0, a1, a2, a3;
 220   real weight;
 221
 222   /*  if (bFirst) {
 223     rw = r - w;
 224     a0 = (3*r*w*w-w*w*w)/(rw*rw*rw);
 225     a1 = -6*r*w/(rw*rw*rw);
 226     a2 = 3*(r+w)/(rw*rw*rw);
 227     a3 = -2/(rw*rw*rw);
 228     bFirst = FALSE;
 229   }
 230   */
 231
 232   if (x < r)
 233     weight = 1;
 234   else if (x > w)
 235     weight = 0;
 236   else
 237     weight = -w/(r-w) + x/(r-w);
 238
 239   return weight;
 240 }
 241
 242 static real get_cylinder_distance(rvec x, rvec com, matrix box) {
 243   /* Triclinic compatible ??? */
 244   real dr, dx, dy, boxx, boxy;
 245
 246   boxx = box[XX][XX];
 247   boxy = box[YY][YY];
 248
 249   dx = fabs(x[XX] - com[XX]);
 250   while (dx > 0.5*boxx) dx -= boxx;
 251
 252   dy = fabs(x[YY] - com[YY]);
 253   while (dy > 0.5*boxy) dy -= boxy;
 254
 255   dr = sqrt(dx*dx+dy*dy);
 256 #ifdef CYLDEBUG
 257   fprintf(stderr,"x:%8.3f%8.3f%8.3f com:%8.3f%8.3f%8.3f dx,dy,dr:%8.3f%8.3f%8.3f\n",x[0],x[1],x[2],com[0],com[1],com[2],dx,dy,dr);
 258 #endif
 259   return dr;
 260 }
 261
 262 void make_refgrps(t_pull *pull,matrix box,t_mdatoms *md)
 263 {
 264   int ngrps,i,ii,j,k,m;
 265   static bool bFirst = TRUE;
 266   real dr,mass;
 267   real truemass;
 268   rvec test;
 269
 270   ngrps = pull->pull.n;
 271   if (bFirst) {
 272     snew(pull->dyna.ngx,ngrps);
 273     snew(pull->dyna.idx,ngrps);
 274     snew(pull->dyna.weights,ngrps);
 275     for (i=0;i<ngrps;i++) {
 276       snew(pull->dyna.idx[i],pull->ref.ngx[0]);    /* more than nessary */
 277       snew(pull->dyna.weights[i],pull->ref.ngx[0]);
 278     }
 279     bFirst = FALSE;
 280   }
 281
 282   /* loop over all groups to make a reference group for each*/
 283   for (i=0;i<ngrps;i++) {
 284     k=0;
 285     truemass=0;
 286     pull->dyna.tmass[i] = 0;
 287     pull->dyna.ngx[i] = 0;
 288
 289     /* loop over all atoms in the main ref group */
 290     for (j=0;j<pull->ref.ngx[0];j++) {
 291       ii = pull->ref.idx[0][j];
 292
 293       /* get_distance takes pbc into account */
 294       dr = get_cylinder_distance(pull->ref.x0[0][j],pull->pull.x_unc[i],box);
 295
 296       if (dr < pull->rc) {
 297         /* add to index, to sum of COM, to weight array */
 298         mass = md->massT[ii];
 299         truemass += mass;
 300         pull->dyna.ngx[i]++;
 301         pull->dyna.weights[i][k] = get_weight(dr,pull->r,pull->rc);
 302         pull->dyna.idx[i][k] = ii;
 303         for (m=0;m<DIM;m++)
 304           pull->dyna.x_unc[i][m] += mass*pull->dyna.weights[i][k]*
 305             pull->ref.x0[0][j][m];
 306         pull->dyna.tmass[i] += mass*pull->dyna.weights[i][k];
 307         k++;
 308       }
 309     }
 310
 311     /* normalize the new 'x_unc' */
 312     svmul(1/pull->dyna.tmass[i],pull->dyna.x_unc[i],pull->dyna.x_unc[i]);
 313     if (pull->bVerbose)
 314       fprintf(stderr,"Made group %d:%8.3f%8.3f%8.3f wm:%8.3f m:%8.3f\n",
 315               i,pull->dyna.x_unc[i][0],pull->dyna.x_unc[i][1],
 316               pull->dyna.x_unc[i][2],pull->dyna.tmass[i],truemass);
 317   }
 318 }
 319