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