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.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 #include <stdlib.h>
  37 #include "sysstuff.h"
  38 #include "princ.h"
  39 #include "futil.h"
  40 #include "statutil.h"
  41 #include "vec.h"
  42 #include "smalloc.h"
  43 #include "typedefs.h"
  44 #include "names.h"
  45 #include "fatal.h"
  46 #include "macros.h"
  47 #include "rdgroup.h"
  48 #include "symtab.h"
  49 #include "index.h"
  50 #include "confio.h"
  51 #include "pbc.h"
  52 #include "pull.h"
  53 #include "pull_internal.h"
  54
  55 void d_pbc_dx(matrix box,const dvec x1, const dvec x2, dvec dx)
  56 {
  57   int m,d;
  58
  59   /* Not correct for all trilinic boxes !!!
  60      Should be as pbc_dx and moved to pbc.c
  61   */
  62   dvec_sub(x1,x2,dx);
  63   for(m=DIM-1; m>=0; m--) {
  64     while (dx[m]   < -0.5*box[m][m])
  65       for(d=0; d<DIM; d++)
  66         dx[d] += box[m][d];
  67     while (dx[m]   >=  0.5*box[m][m])
  68       for(d=0; d<DIM; d++)
  69         dx[d] -= box[m][d];
  70   }
  71 }
  72
  73 void put_dvec_in_box(matrix box, dvec v)
  74 {
  75   int m,d;
  76
  77   for(m=DIM-1; m>=0; m--) {
  78     while (v[m] < 0)
  79       for(d=0; d<DIM; d++)
  80         v[d] += box[m][d];
  81     while (v[m] >= box[m][m])
  82       for(d=0; d<DIM; d++)
  83         v[d] -= box[m][d];
  84   }
  85 }
  86
  87 /* calculates center of mass of selection index from all coordinates x */
  88 void calc_com(t_pullgrp *pg, rvec x[], t_mdatoms *md, matrix box)
  89 {
  90   int  i,ii,m;
  91   real wm;
  92   dvec com;
  93
  94   clear_dvec(com);
  95   for(i=0; i<pg->ngx; i++) {
  96     ii = pg->idx[i];
  97     wm = md->massT[ii];
  98     if (pg->nweight > 0)
  99       wm *= pg->weight[i];
 100     for(m=0; m<DIM; m++)
 101       com[m] += wm*x[ii][m];
 102   }
 103   for(m=0; m<DIM; m++)
 104     com[m] *= pg->wscale/pg->tmass;
 105   put_dvec_in_box(box,com);
 106   copy_dvec(com,pg->x_unc);
 107 }
 108
 109 /* calculates com of all atoms in x[], *index has their index numbers
 110    to get the masses from atom[] */
 111 void calc_com2(t_pullgrp *pg, rvec x[], t_mdatoms *md, matrix box)
 112 {
 113   int  i,ii,m;
 114   real wm;
 115   dvec com;
 116
 117   clear_dvec(com);
 118   for(i=0; i<pg->ngx; i++) {
 119     ii = pg->idx[i];
 120     wm = md->massT[ii];
 121     if (pg->nweight > 0)
 122       wm *= pg->weight[i];
 123     for(m=0; m<DIM; m++)
 124       com[m] += wm*x[i][m];
 125   }
 126   for(m=0; m<DIM; m++)
 127     com[m] *= pg->wscale/pg->tmass;
 128   put_dvec_in_box(box,com);
 129   copy_dvec(com,pg->x_unc);
 130 }
 131
 132 void calc_running_com(t_pull *pull) {
 133   int i,j,n;
 134   dvec ave;
 135   double tm;
 136
 137   /* for group i, we have nhist[i] points of history. The maximum nr of points
 138      is pull->reflag. The array comhist[i][0..nhist[i]-1] has the positions
 139      of the center of mass over the last nhist[i] steps. x_unc[i] has the
 140      new ones. Remove the oldest one from the list, add the new one, calculate
 141      the average, put that in x_unc instead and return. We just move all coms
 142      1 down in the list and add the latest one to the top.
 143    */
 144
 145   if(pull->bCyl) {
 146     /* act on dyna groups */
 147     for(i=0;i<pull->ngrp;i++) {
 148       clear_dvec(ave);
 149       for(j=0;j<(pull->reflag-1);j++) {
 150         copy_dvec(pull->dyna[i].comhist[j+1],pull->dyna[i].comhist[j]);
 151         dvec_add(ave,pull->dyna[i].comhist[j],ave);
 152       }
 153       copy_dvec(pull->dyna[i].x_unc,pull->dyna[i].comhist[j]);
 154       dvec_add(ave,pull->dyna[i].x_unc,ave);
 155       dsvmul(1.0/pull->reflag,ave,ave);
 156
 157       /* now ave has the running com for group i, copy it to x_unc[i] */
 158       copy_dvec(ave,pull->dyna[i].x_unc);
 159
 160 #ifdef DEBUG
 161       if(pull->bVerbose)
 162         for(n=0;n<pull->reflag;n++)
 163           fprintf(stderr,"Comhist %d, grp %d: %8.3f%8.3f%8.3f\n",n,i,
 164                   pull->dyna[i].comhist[n][XX],
 165                   pull->dyna[i].comhist[n][YY],
 166                   pull->dyna[i].comhist[n][ZZ]);
 167 #endif
 168     }
 169   } else {
 170     /* act on ref group */
 171     clear_dvec(ave);
 172
 173     for(j=0;j<(pull->reflag-1);j++) {
 174       copy_dvec(pull->ref.comhist[j+1],pull->ref.comhist[j]);
 175       dvec_add(ave,pull->ref.comhist[j],ave);
 176     }
 177
 178     copy_dvec(pull->ref.x_unc,pull->ref.comhist[j]);
 179     dvec_add(ave,pull->ref.x_unc,ave);
 180     dsvmul(1.0/pull->reflag,ave,ave);
 181     /* now ave has the running com for group i, copy it to x_unc */
 182     copy_dvec(ave,pull->ref.x_unc);
 183
 184 #ifdef DEBUG
 185     if(pull->bVerbose)
 186       for(i=0;i<pull->reflag;i++)
 187         fprintf(stderr,"Comhist %d: %8.3f%8.3f%8.3f\n",i,
 188                 pull->ref.comhist[i][XX],
 189                 pull->ref.comhist[i][YY],
 190                 pull->ref.comhist[i][ZZ]);
 191 #endif
 192   }
 193 }
 194
 195 void correct_t0_pbc(t_pull *pull, rvec x[], t_mdatoms *md, matrix box) {
 196   int i,ii,j,m;
 197   rvec dx;
 198
 199   /* loop over all atoms in index for group i. Check if they moved
 200      more than half a box with respect to xp. If so add/subtract a box
 201      from x0. Copy x to xp.
 202    */
 203   for(i=0;i<pull->ref.ngx;i++) {
 204     ii = pull->ref.idx[i];
 205
 206     /* correct for jumps across the box */
 207     pbc_dx(x[ii],pull->ref.xp[i],dx);
 208     for(m=0; m<DIM; m++) {
 209       pull->ref.x0[i][m] += dx[m];
 210       pull->ref.xp[i][m]  = x[ii][m];
 211     }
 212   }
 213 }
 214
 215 /* switch function, x between r and w */
 216 real get_weight(real x, real r, real w) {
 217
 218   /* Removed static stuff. Wasn't being used anyway */
 219   /*  static bool bFirst = TRUE;
 220   static real rw, a0, a1, a2, a3; */
 221
 222   real weight;
 223
 224   /*  if (bFirst) {
 225     rw = r - w;
 226     a0 = (3*r*w*w-w*w*w)/(rw*rw*rw);
 227     a1 = -6*r*w/(rw*rw*rw);
 228     a2 = 3*(r+w)/(rw*rw*rw);
 229     a3 = -2/(rw*rw*rw);
 230     bFirst = FALSE;
 231   }
 232   */
 233
 234   if (x <= r)
 235     weight = 1;
 236   else if (x >= w)
 237     weight = 0;
 238   else
 239     weight = (w - x)/(w - r);
 240
 241   return weight;
 242 }
 243
 244 static double get_cylinder_distance(rvec x, dvec com, matrix box) {
 245   /* Triclinic compatible ??? */
 246   double dr, dx, dy, boxx, boxy;
 247
 248   boxx = box[XX][XX];
 249   boxy = box[YY][YY];
 250
 251   dx = fabs(x[XX] - com[XX]);
 252   while(dx > 0.5*boxx) dx -= boxx;
 253
 254   dy = fabs(x[YY] - com[YY]);
 255   while(dy > 0.5*boxy) dy -= boxy;
 256
 257   dr = sqrt(dx*dx+dy*dy);
 258 #ifdef CYLDEBUG
 259   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);
 260 #endif
 261   return dr;
 262 }
 263
 264 void make_refgrps(t_pull *pull,matrix box,t_mdatoms *md)
 265 {
 266   int i,ii,j,k,m;
 267
 268   double dr,mass,wmass,wwmass;
 269   dvec test;
 270   t_pullgrp *pdyna;
 271
 272   /* loop over all groups to make a reference group for each*/
 273   for(i=0;i<pull->ngrp;i++) {
 274     pdyna = &pull->dyna[i];
 275     k=0;
 276     wmass = 0;
 277     wwmass = 0;
 278     pdyna->ngx = 0;
 279
 280     /* loop over all atoms in the main ref group */
 281     for(j=0;j<pull->ref.ngx;j++) {
 282       ii = pull->ref.idx[j];
 283
 284       /* get_distance takes pbc into account */
 285       dr = get_cylinder_distance(pull->ref.x0[j],pull->grp[i].x_unc,box);
 286
 287       if(dr < pull->rc) {
 288         /* add to index, to sum of COM, to weight array */
 289         mass = md->massT[ii];
 290         pdyna->ngx++;
 291         pdyna->weight[k] = get_weight(dr,pull->r,pull->rc);
 292         pdyna->idx[k] = ii;
 293         for(m=0;m<DIM;m++)
 294           pdyna->x_unc[m] += mass*pdyna->weight[k]*pull->ref.x0[j][m];
 295         wmass += mass*pdyna->weight[k];
 296         wwmass += mass*sqr(pdyna->weight[k]);
 297         k++;
 298       }
 299     }
 300
 301     pdyna->wscale = wmass/wwmass;
 302     pdyna->tmass = pdyna->wscale*wmass;
 303
 304     /* normalize the new 'x_unc' */
 305     dsvmul(1/wmass,pdyna->x_unc,pdyna->x_unc);
 306     if(pull->bVerbose)
 307       fprintf(stderr,"Made group %d:%8.3f%8.3f%8.3f m:%8.3f\n",
 308               i,pdyna->x_unc[0],pdyna->x_unc[1],
 309               pdyna->x_unc[2],pdyna->tmass);
 310   }
 311 }
 312