src/mdlib/pullinit.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_pullinit_c = "$Id$";
  30
  31 #include <string.h>
  32 #include "princ.h"
  33 #include <stdlib.h>
  34 #include "readinp.h"
  35 #include "sysstuff.h"
  36 #include "futil.h"
  37 #include "statutil.h"
  38 #include "vec.h"
  39 #include "smalloc.h"
  40 #include "typedefs.h"
  41 #include "names.h"
  42 #include "fatal.h"
  43 #include "macros.h"
  44 #include "rdgroup.h"
  45 #include "symtab.h"
  46 #include "index.h"
  47 #include "confio.h"
  48 #include "pull.h"
  49 #include "string.h"
  50
  51 static void get_pullmemory(t_pullgrps *pull, int ngrps)
  52 {
  53   snew(pull->ngx,ngrps);
  54   snew(pull->x_con,ngrps);
  55   snew(pull->xprev,ngrps);
  56   snew(pull->tmass,ngrps);
  57   snew(pull->idx,ngrps);
  58   snew(pull->f,ngrps);
  59   snew(pull->spring,ngrps);
  60   snew(pull->dir,ngrps);
  61   snew(pull->x_unc,ngrps);
  62   snew(pull->x_ref,ngrps);
  63   snew(pull->d_ref,ngrps);
  64   snew(pull->x0,ngrps);
  65   snew(pull->xp,ngrps);
  66   snew(pull->comhist,ngrps);
  67 }
  68
  69 static void print_info(FILE *log,t_pull *pull)
  70 {
  71
  72   fprintf(log,"\n**************************************************\n"
  73           "                         PULL INFO                      \n"
  74           "**************************************************\n");
  75
  76   switch (pull->runtype) {
  77   case eStart:
  78     fprintf(log,"RUN TYPE: Generating Starting structures\n");
  79     break;
  80   case eAfm:
  81     fprintf(log,"RUN TYPE: Afm\n");
  82     break;
  83   case eConstraint:
  84     fprintf(log,"RUN TYPE: Constraint\n");
  85     break;
  86   case eUmbrella:
  87     fprintf(log,"RUN TYPE: Umbrella sampling\n");
  88     break;
  89   case eTest:
  90     fprintf(log,"RUN TYPE: Test run\n");
  91     break;
  92   default:
  93     fprintf(log,"RUN TYPE: WARNING! pullinit does not know this runtype\n");
  94   }
  95
  96   if (pull->runtype == eConstraint || pull->runtype == eStart)
  97     switch (pull->reftype) {
  98     case eCom:
  99       fprintf(log,"REFERENCE TYPE: center of mass of reference group\n");
 100       break;
 101     case eComT0:
 102       fprintf(log,
 103               "REFERENCE TYPE: center of mass of reference group at t=0\n");
 104       break;
 105     case eDyn:
 106       fprintf(log,
 107               "REFERENCE TYPE: center of mass of dynamically made groups\n");
 108       fprintf(log,"Using dynamic reference groups: r=%8.3f, rc=%8.3f\n",
 109               pull->r,pull->rc);
 110       break;
 111     case eDynT0:
 112       fprintf(log,
 113               "REFERENCE TYPE: center of mass of dynamically made groups,\n"
 114               "based on the positions of its atoms at t=0\n");
 115       fprintf(log,"Using dynamic reference groups: r=%8.3f, rc=%8.3f\n",
 116               pull->r,pull->rc);
 117       break;
 118     default:
 119       fprintf(log,"REFERENCE TYPE: no clue! What did you do wrong?\n");
 120     }
 121 }
 122
 123 static void get_named_indexgroup(FILE *log,atom_id **target,int *isize,
 124                                  char *name,atom_id **index,int ngx[],
 125                                  char **grpnames,int ngrps)
 126 {
 127   int i,j;
 128   bool bFound = FALSE;
 129   atom_id *tmp = NULL;
 130
 131   fprintf(log,"Looking for group %s: ",name);
 132   for (i=0;i<ngrps;i++) {
 133     if (!strcmp(name,grpnames[i])) {
 134       /* found the group we're looking for */
 135       snew(tmp,ngx[i]);
 136       for (j=0;j<ngx[i];j++)
 137         tmp[j]=index[i][j];
 138       *isize=ngx[i];
 139       bFound = TRUE;
 140       fprintf(log,"found group %s: %d elements. First: %d\n",name,ngx[i],
 141               tmp[0]);
 142     }
 143   }
 144
 145   *target=tmp;
 146   if (!bFound)
 147     fatal_error(0,"Can't find group %s in the index file",name);
 148 }
 149
 150 static void read_whole_index(char *indexfile,char ***grpnames,
 151                              atom_id ***index, int **ngx,int *totalgrps)
 152 {
 153   t_block *grps;
 154   char    **gnames;
 155   int i,j;
 156
 157   if (!indexfile)
 158     fatal_error(0,"No index file specified");
 159
 160   grps = init_index(indexfile,&gnames);
 161   if (grps->nr==0)
 162     fatal_error(0,"No groups in indexfile\n");
 163
 164   *totalgrps = grps->nr;
 165   snew(*index,grps->nr);
 166   snew(*grpnames,grps->nr);
 167   snew((*ngx),grps->nr);
 168   /* memory leak, can't free ngx anymore. 4bytes/group */
 169
 170   for(i=0; (i<*totalgrps); i++) {
 171     (*grpnames)[i]=strdup(gnames[i]);
 172     (*ngx)[i]=grps->index[i+1]-grps->index[i];
 173     snew((*index)[i],(*ngx)[i]);
 174     for(j=0; (j<(*ngx)[i]); j++)
 175       (*index)[i][j]=grps->a[grps->index[i]+j];
 176   }
 177
 178   for(i=0; (i<grps->nr); i++)
 179     sfree(gnames[i]);
 180   sfree(gnames);
 181   done_block(grps);
 182   sfree(grps);
 183 }
 184
 185 static void print_whole_index(char **grpnames, atom_id **index, int *ngx,
 186                               int ngrps)
 187 {
 188   int i,j;
 189   FILE *tmp;
 190
 191   tmp = ffopen("indexcheck","w");
 192   for (i=0;i<ngrps;i++) {
 193     fprintf(tmp,"\nGrp %d: %s. %d elements\n",i,grpnames[i],ngx[i]);
 194     for (j=0;j<ngx[i];j++)
 195       fprintf(tmp," %d ",index[i][j]);
 196   }
 197   fflush(tmp);
 198 }
 199
 200 void init_pull(FILE *log,int nfile,t_filenm fnm[],t_pull *pull,rvec *x,
 201                t_mdatoms *md,rvec boxsize)
 202 {
 203   int i,j,m,ii;
 204   int ngrps;
 205   real tm;
 206   rvec tmp;
 207   matrix box;
 208   char **grpnames;
 209   atom_id **index;
 210   int *ngx;
 211   int totalgrps;    /* total number of groups in the index file */
 212
 213   clear_mat(box);
 214   for (m=0;m<DIM;m++)
 215     box[m][m]=boxsize[m];
 216
 217   /* do we have to do any pulling at all? If not return */
 218   pull->bPull = opt2bSet("-pi",nfile,fnm);
 219   if (!pull->bPull) return;
 220
 221   pull->out = ffopen(opt2fn("-pd",nfile,fnm),"w");
 222   read_pullparams(pull, opt2fn("-pi",nfile,fnm), opt2fn("-po",nfile,fnm));
 223   ngrps = pull->pull.n;
 224
 225   if (pull->reftype == eDyn || pull->reftype == eDynT0)
 226     pull->bCyl = TRUE;
 227   else
 228     pull->bCyl = FALSE;
 229
 230   if (pull->bCyl && (pull->rc < 0.01 || pull->r < 0.01))
 231     fatal_error(0,"rc or r is too small or zero.");
 232
 233   print_info(log,pull);
 234
 235   get_pullmemory(&pull->pull,ngrps);
 236   get_pullmemory(&pull->dyna,ngrps);
 237   get_pullmemory(&pull->ref,1);
 238
 239   /* read the whole index file */
 240   read_whole_index(opt2fn("-pn",nfile,fnm),&grpnames,&index,&ngx,&totalgrps);
 241
 242   if (pull->bVerbose) {
 243     fprintf(stderr,"read_whole_index: %d groups total\n",totalgrps);
 244     for(i=0;i<totalgrps;i++)
 245       fprintf(stderr,"group %i (%s) %d elements\n",
 246               i,grpnames[i],ngx[i]);
 247     /*    print_whole_index(grpnames,index,ngx,totalgrps); */
 248   }
 249
 250   /* grab the groups that are specified in the param file */
 251   for (i=0;i<pull->pull.n;i++)
 252     get_named_indexgroup(log,&pull->pull.idx[i],&pull->pull.ngx[i],
 253                          pull->pull.grps[i],index,ngx,grpnames,totalgrps) ;
 254   get_named_indexgroup(log,&pull->ref.idx[0],&pull->ref.ngx[0],
 255                        pull->ref.grps[0],index,ngx,grpnames,totalgrps);
 256
 257   /* get more memory! Don't we love C? */
 258   snew(pull->ref.x0[0],pull->ref.ngx[0]);
 259   snew(pull->ref.xp[0],pull->ref.ngx[0]);
 260   snew(pull->ref.comhist[0],pull->reflag);
 261   for (i=0;i<pull->pull.n;i++)
 262     snew(pull->dyna.comhist[i],pull->reflag);
 263
 264   for (i=0;i<ngrps;i++) {
 265     tm = calc_com(x,pull->pull.ngx[i],pull->pull.idx[i],
 266                   md,tmp,box);
 267     copy_rvec(tmp,pull->pull.x_con[i]);
 268     copy_rvec(tmp,pull->pull.x_unc[i]);
 269     copy_rvec(tmp,pull->pull.x_ref[i]);
 270     copy_rvec(tmp,pull->pull.spring[i]);
 271     fprintf(log,"Initializing pull groups. Mass of group %d: %8.3f\n"
 272             "Initial coordinates center of mass: %8.3f %8.3f %8.3f\n",
 273             i,tm,tmp[XX],tmp[YY],tmp[ZZ]);
 274     pull->pull.tmass[i] = tm;
 275   }
 276
 277   /* initialize the reference group, in all cases */
 278   tm = calc_com(x,pull->ref.ngx[0],pull->ref.idx[0],md,
 279                 tmp,box);
 280
 281   copy_rvec(tmp,pull->ref.x_unc[0]);
 282   copy_rvec(tmp,pull->ref.x_con[0]);
 283   copy_rvec(tmp,pull->ref.x_ref[0]);
 284
 285   for (j=0;j<pull->reflag;j++)
 286     copy_rvec(pull->ref.x_unc[0],pull->ref.comhist[0][j]);
 287
 288   fprintf(log,"Initializing reference group. Mass: %8.3f\n"
 289           "Initial coordinates center of mass: %8.3f %8.3f %8.3f\n",
 290           tm,tmp[XX],tmp[YY],tmp[ZZ]);
 291   /* keep the initial coordinates for center of mass at t0 */
 292   for (j=0;j<pull->ref.ngx[0];j++) {
 293     copy_rvec(x[pull->ref.idx[0][j]],pull->ref.x0[0][j]);
 294     copy_rvec(x[pull->ref.idx[0][j]],pull->ref.xp[0][j]);
 295   }
 296   pull->ref.tmass[0] = tm;
 297
 298   /* if we use dynamic reference groups, do some initialising for them */
 299   if (pull->bCyl) {
 300     make_refgrps(pull,x,md);
 301     for (i=0;i<ngrps;i++) {
 302       copy_rvec(pull->dyna.x_unc[i],pull->dyna.x_con[i]);
 303       copy_rvec(pull->dyna.x_unc[i],pull->dyna.x_ref[i]);
 304
 305       /* initialize comhist values for running averages */
 306       for (j=0;j<pull->reflag;j++)
 307         copy_rvec(pull->dyna.x_unc[i],pull->dyna.comhist[i][j]);
 308
 309       fprintf(log,"Initializing dynamic groups. %d atoms. Weighted mass"
 310               "for %d:%8.3f\n"
 311               "Initial coordinates center of mass: %8.3f %8.3f %8.3f\n",
 312               pull->dyna.ngx[i],i,pull->dyna.tmass[i],pull->dyna.x_ref[i][XX],
 313               pull->dyna.x_ref[i][YY],pull->dyna.x_ref[i][ZZ]);
 314     }
 315   }
 316
 317   /* set the reference distances and directions, taking into account pbc */
 318   for (i=0;i<ngrps;i++) {
 319     if (pull->bCyl) {
 320       rvec_sub(pull->pull.x_ref[i],pull->dyna.x_ref[i],tmp);
 321       for (m=0;m<DIM;m++) {
 322         if (tmp[m] < -0.5*boxsize[m]) tmp[m] += boxsize[m];
 323         if (tmp[m] >  0.5*boxsize[m]) tmp[m] -= boxsize[m];
 324       }
 325     } else {
 326       rvec_sub(pull->pull.x_ref[i],pull->ref.x_ref[0],tmp);
 327       for (m=0;m<DIM;m++) {
 328         if (tmp[m] < -0.5*boxsize[m]) tmp[m] += boxsize[m];
 329         if (tmp[m] >  0.5*boxsize[m]) tmp[m] -= boxsize[m];
 330       }
 331     }
 332
 333     /* reference distance for constraint run */
 334     pull->pull.d_ref[i] = norm(tmp);
 335
 336     /* select elements of direction vector to use for Afm and Start runs */
 337     for (m=0;m<DIM;m++)
 338       tmp[m] = tmp[m] * pull->dims[m];
 339     svmul(1/norm(tmp),tmp,pull->pull.dir[i]);
 340     if (pull->bReverse)
 341       svmul(-1.0,pull->pull.dir[i],pull->pull.dir[i]);
 342
 343     if (pull->runtype == eAfm)
 344       fprintf(log,"\nPull rate: %e nm/step. Force constant: %e kJ/(mol nm)",
 345               pull->rate,pull->k);
 346     if (pull->runtype == eAfm || pull->runtype == eStart)
 347       fprintf(log,"\nPull direction: %8.3f %8.3f %8.3f bReverse = %d\n",
 348               pull->pull.dir[i][XX],pull->pull.dir[i][YY],
 349               pull->pull.dir[i][ZZ],pull->bReverse);
 350
 351   }
 352   fprintf(log,"**************************************************\n"
 353           "                      END   PULL INFO                    \n"
 354           "**************************************************\n\n");
 355 }
 356
 357
 358
 359
 360