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