src/tools/gmx_make_edi.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  *
  11  * The make_edi program was generously contributed by Oliver Lange, based
  12  * on the code from g_anaeig. You can reach him as olange@gwdg.de.
  13  *
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  *
  31  * And Hey:
  32  * Gromacs Runs One Microsecond At Cannonball Speeds
  33  */
  34 #ifdef HAVE_CONFIG_H
  35 #include <config.h>
  36 #endif
  37
  38 #include <math.h>
  39 #include <stdlib.h>
  40 #include <ctype.h>
  41 #include <string.h>
  42 #include "readinp.h"
  43 #include "statutil.h"
  44 #include "sysstuff.h"
  45 #include "typedefs.h"
  46 #include "smalloc.h"
  47 #include "macros.h"
  48 #include "gmx_fatal.h"
  49 #include "vec.h"
  50 #include "pbc.h"
  51 #include "copyrite.h"
  52 #include "futil.h"
  53 #include "statutil.h"
  54 #include "pdbio.h"
  55 #include "confio.h"
  56 #include "tpxio.h"
  57 #include "matio.h"
  58 #include "mshift.h"
  59 #include "xvgr.h"
  60 #include "do_fit.h"
  61 #include "rmpbc.h"
  62 #include "txtdump.h"
  63 #include "eigio.h"
  64 #include "index.h"
  65 #include "string2.h"
  66
  67 typedef struct
  68 {
  69     real        deltaF0;
  70     gmx_bool    bHarmonic;
  71     gmx_bool    bConstForce;   /* Do constant force flooding instead of
  72                                   evaluating a flooding potential             */
  73     real        tau;
  74     real        deltaF;
  75     real        kT;
  76     real        constEfl;
  77     real        alpha2;
  78 } t_edflood;
  79
  80
  81 /* This type is for the average, reference, target, and origin structure   */
  82 typedef struct edix
  83 {
  84     int         nr;             /* number of atoms this structure contains */
  85     int         *anrs;          /* atom index numbers                      */
  86     rvec        *x;             /* positions                               */
  87     real        *sqrtm;         /* sqrt of the masses used for mass-
  88                                  * weighting of analysis                   */
  89 } t_edix;
  90
  91
  92 typedef struct edipar
  93 {
  94     int         nini;           /* total Nr of atoms                    */
  95     gmx_bool    fitmas;         /* true if trans fit with cm            */
  96     gmx_bool    pcamas;         /* true if mass-weighted PCA            */
  97     int         presteps;       /* number of steps to run without any
  98                                  *    perturbations ... just monitoring */
  99     int         outfrq;         /* freq (in steps) of writing to edo    */
 100     int         maxedsteps;     /* max nr of steps per cycle            */
 101     struct edix sref;           /* reference positions, to these fitting
 102                                  * will be done                         */
 103     struct edix sav;            /* average positions                    */
 104     struct edix star;           /* target positions                     */
 105     struct edix sori;           /* origin positions                     */
 106     real        slope;          /* minimal slope in acceptance radexp   */
 107     int         ned;            /* Nr of atoms in essdyn buffer         */
 108     t_edflood   flood;          /* parameters especially for flooding   */
 109 } t_edipar;
 110
 111
 112
 113 void make_t_edx(struct edix *edx, int natoms, rvec *pos, atom_id index[]) {
 114   edx->nr=natoms;
 115   edx->anrs=index;
 116   edx->x=pos;
 117 }
 118
 119 void write_t_edx(FILE *fp, struct edix edx, const char *comment) {
 120  /*here we copy only the pointers into the t_edx struct
 121   no data is copied and edx.box is ignored  */
 122  int i;
 123   fprintf(fp,"#%s \n %d \n",comment,edx.nr);
 124   for (i=0;i<edx.nr;i++) {
 125     fprintf(fp,"%d  %f  %f  %f\n",(edx.anrs)[i]+1,(edx.x)[i][XX],(edx.x)[i][YY],(edx.x)[i][ZZ]);
 126   }
 127 }
 128
 129 int sscan_list(int *list[], const char *str, const char *listname) {
 130  /*this routine scans a string of the form 1,3-6,9 and returns the
 131     selected numbers (in this case 1 3 4 5 6 9) in NULL-terminated array of integers.
 132     memory for this list will be allocated  in this routine -- sscan_list expects *list to
 133     be a NULL-Pointer
 134
 135     listname is a string used in the errormessage*/
 136
 137
 138    int i,istep;
 139    char c;
 140    char *pos,*startpos,*step;
 141    int n=strlen(str);
 142
 143    /*enums to define the different lexical stati */
 144    enum { sBefore, sNumber, sMinus, sRange, sZero, sSmaller, sError, sSteppedRange };
 145
 146    int status=sBefore; /*status of the deterministic automat to scan str   */
 147    int number=0;
 148    int end_number=0;
 149
 150    char *start=NULL; /*holds the string of the number behind a ','*/
 151    char *end=NULL; /*holds the string of the number behind a '-' */
 152
 153    int nvecs=0; /* counts the number of vectors in the list*/
 154
 155    step=NULL;
 156    snew(pos,n+4);
 157    startpos=pos;
 158    strcpy(pos,str);
 159    pos[n]=',';
 160    pos[n+1]='1';
 161    pos[n+2]='\0';
 162
 163    *list = NULL;
 164
 165    while ((c=*pos)!=0) {
 166      switch(status) {
 167         /* expect a number */
 168         case sBefore: if (isdigit(c)) {
 169               start=pos;
 170               status=sNumber;
 171               break;
 172            } else status=sError; break;
 173
 174         /* have read a number, expect ',' or '-' */
 175         case sNumber: if (c==',') {
 176              /*store number*/
 177              srenew(*list,nvecs+1);
 178              (*list)[nvecs++]=number=strtol(start,NULL,10);
 179              status=sBefore;
 180              if (number==0)
 181                  status=sZero;
 182               break;
 183              } else if (c=='-') { status=sMinus; break; }
 184              else if (isdigit(c)) break;
 185              else status=sError; break;
 186
 187         /* have read a '-' -> expect a number */
 188      case sMinus:
 189        if (isdigit(c)) {
 190          end=pos;
 191          status=sRange; break;
 192        } else status=sError; break;
 193
 194      case sSteppedRange:
 195        if (isdigit(c)) {
 196          if (step) {
 197            status=sError; break;
 198          } else
 199            step=pos;
 200          status=sRange;
 201          break;
 202        } else status=sError; break;
 203
 204         /* have read the number after a minus, expect ',' or ':' */
 205         case sRange:
 206             if (c==',') {
 207                /*store numbers*/
 208                end_number=strtol(end,NULL,10);
 209                number=strtol(start,NULL,10);
 210                status=sBefore;
 211                if (number==0) {
 212                   status=sZero; break;
 213                }
 214                if (end_number<=number) {
 215                   status=sSmaller; break;
 216                }
 217                srenew(*list,nvecs+end_number-number+1);
 218                if (step) {
 219                  istep=strtol(step,NULL,10);
 220                  step=NULL;
 221                } else istep=1;
 222                for (i=number;i<=end_number;i+=istep)
 223                     (*list)[nvecs++]=i;
 224                break;
 225             } else if (c==':') {
 226               status = sSteppedRange;
 227               break;
 228             } else if (isdigit(c)) break; else status=sError;  break;
 229
 230        /* format error occured */
 231        case sError:
 232            gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d with char %c",listname,pos-startpos,*(pos-1));
 233            break;
 234        /* logical error occured */
 235        case sZero:
 236                    gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d: eigenvector 0 is not valid",listname,pos-startpos);
 237                    break;
 238        case sSmaller:
 239                    gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d: second index %d is not bigger than %d",listname,pos-startpos,end_number,number);
 240                    break;
 241      }
 242    ++pos; /* read next character */
 243    } /*scanner has finished */
 244
 245    /* append zero to list of eigenvectors */
 246    srenew(*list,nvecs+1);
 247    (*list)[nvecs]=0;
 248    sfree(startpos);
 249    return nvecs;
 250 } /*sscan_list*/
 251
 252 void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs,int nvec, const char *grouptitle,real steps[]) {
 253 /* eig_list is a zero-terminated list of indices into the eigvecs array.
 254    eigvecs are coordinates of eigenvectors
 255    grouptitle to write in the comment line
 256    steps  -- array with stepsizes for evLINFIX, evLINACC and evRADACC
 257 */
 258
 259   int n=0,i; rvec x;
 260   real sum;
 261   while (eig_list[n++]);  /*count selected eigenvecs*/
 262
 263   fprintf(fp,"# NUMBER OF EIGENVECTORS + %s\n %d\n",grouptitle,n-1);
 264
 265   /* write list of eigenvector indicess */
 266   for(n=0;eig_list[n];n++) {
 267     if (steps)
 268       fprintf(fp,"%8d   %g\n",eig_list[n],steps[n]);
 269     else
 270       fprintf(fp,"%8d   %g\n",eig_list[n],1.0);
 271   }
 272   n=0;
 273
 274   /* dump coordinates of the selected eigenvectors */
 275   while (eig_list[n]) {
 276     sum=0;
 277     for (i=0; i<natoms; i++) {
 278       if (eig_list[n]>nvec)
 279         gmx_fatal(FARGS,"Selected eigenvector %d is higher than maximum number %d of available eigenvectors",eig_list[n],nvec);
 280       copy_rvec(eigvecs[eig_list[n]-1][i],x);
 281       sum+=norm2(x);
 282       fprintf(fp,"%8.5f %8.5f %8.5f\n",x[XX],x[YY],x[ZZ]);
 283     }
 284     n++;
 285   }
 286 }
 287
 288
 289 /*enum referring to the different lists of eigenvectors*/
 290 enum { evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON , evMON,  evNr };
 291 #define oldMAGIC 666
 292 #define MAGIC 670
 293
 294
 295 void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
 296                            int nvec, int *eig_listen[], real* evStepList[])
 297 {
 298 /* write edi-file */
 299
 300     /*Header*/
 301     fprintf(fp,"#MAGIC\n %d \n#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
 302         MAGIC,edpars->nini,edpars->fitmas,edpars->pcamas);
 303     fprintf(fp,"#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
 304         edpars->outfrq,edpars->maxedsteps,edpars->slope);
 305     fprintf(fp,"#PRESTEPS\n %d\n#DELTA_F0\n %f\n#INIT_DELTA_F\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n#KT\n %f\n#HARMONIC\n %d\n#CONST_FORCE_FLOODING\n %d\n",
 306         edpars->presteps,edpars->flood.deltaF0,edpars->flood.deltaF,edpars->flood.tau,edpars->flood.constEfl,
 307         edpars->flood.alpha2,edpars->flood.kT,edpars->flood.bHarmonic,edpars->flood.bConstForce);
 308
 309     /* Average and reference positions */
 310     write_t_edx(fp,edpars->sref,"NREF, XREF");
 311     write_t_edx(fp,edpars->sav,"NAV, XAV");
 312
 313     /*Eigenvectors */
 314
 315     write_eigvec(fp, edpars->ned, eig_listen[evMON],eigvecs,nvec,"COMPONENTS GROUP 1",NULL);
 316     write_eigvec(fp, edpars->ned, eig_listen[evLINFIX],eigvecs,nvec,"COMPONENTS GROUP 2",evStepList[evLINFIX]);
 317     write_eigvec(fp, edpars->ned, eig_listen[evLINACC],eigvecs,nvec,"COMPONENTS GROUP 3",evStepList[evLINACC]);
 318     write_eigvec(fp, edpars->ned, eig_listen[evRADFIX],eigvecs,nvec,"COMPONENTS GROUP 4",evStepList[evRADFIX]);
 319     write_eigvec(fp, edpars->ned, eig_listen[evRADACC],eigvecs,nvec,"COMPONENTS GROUP 5",NULL);
 320     write_eigvec(fp, edpars->ned, eig_listen[evRADCON],eigvecs,nvec,"COMPONENTS GROUP 6",NULL);
 321     write_eigvec(fp, edpars->ned, eig_listen[evFLOOD],eigvecs,nvec,"COMPONENTS GROUP 7",evStepList[evFLOOD]);
 322
 323
 324     /*Target and Origin positions */
 325     write_t_edx(fp,edpars->star,"NTARGET, XTARGET");
 326     write_t_edx(fp,edpars->sori,"NORIGIN, XORIGIN");
 327 }
 328
 329 int read_conffile(const char *confin,char *title,rvec *x[])
 330 {
 331 /* read coordinates out of STX file  */
 332   int natoms;
 333   t_atoms  confat;
 334   matrix box;
 335   printf("read coordnumber from file %s\n",confin);
 336   get_stx_coordnum(confin,&natoms);
 337   printf("number of coordinates in file %d\n",natoms);
 338 /*  if (natoms != ncoords)
 339      gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n"
 340            "             does not match topology (= %d)",
 341            confin,natoms,ncoords);
 342   else {*/
 343     /* make space for coordinates and velocities */
 344     init_t_atoms(&confat,natoms,FALSE);
 345     snew(*x,natoms);
 346     read_stx_conf(confin,title,&confat,*x,NULL,NULL,box);
 347     return natoms;
 348 }
 349
 350
 351 void read_eigenvalues(int vecs[],const char *eigfile, real values[],
 352                       gmx_bool bHesse, real kT)
 353 {
 354   int  neig,nrow,i;
 355   double **eigval;
 356
 357   neig = read_xvg(eigfile,&eigval,&nrow);
 358
 359   fprintf(stderr,"Read %d eigenvalues\n",neig);
 360   for(i=bHesse ? 6 : 0 ; i<neig; i++) {
 361     if (eigval[1][i] < -0.001 && bHesse)
 362       fprintf(stderr,
 363               "WARNING: The Hessian Matrix has negative eigenvalue %f, we set it to zero (no flooding in this direction)\n\n",eigval[1][i]);
 364
 365     if (eigval[1][i] < 0)
 366       eigval[1][i] = 0;
 367   }
 368   if (bHesse)
 369     for (i=0; vecs[i]; i++) {
 370       if (vecs[i]<7)
 371         gmx_fatal(FARGS,"ERROR: You have chosen one of the first 6 eigenvectors of the HESSE Matrix. That does not make sense, since they correspond to the 6 rotational and translational degrees of freedom.\n\n");
 372       values[i]=eigval[1][vecs[i]-1]/kT;
 373     }
 374   else
 375     for (i=0; vecs[i]; i++) {
 376       if (vecs[i]>(neig-6))
 377         gmx_fatal(FARGS,"ERROR: You have chosen one of the last 6 eigenvectors of the COVARIANCE Matrix. That does not make sense, since they correspond to the 6 rotational and translational degrees of freedom.\n\n");
 378       values[i]=1/eigval[1][vecs[i]-1];
 379     }
 380   /* free memory */
 381   for (i=0; i<nrow; i++)
 382     sfree(eigval[i]);
 383   sfree(eigval);
 384 }
 385
 386
 387 static real *scan_vecparams(const char *str,const char * par, int nvecs)
 388 {
 389   char   f0[256],f1[256];             /*format strings adapted every pass of the loop*/
 390   double d,tcap=0;
 391   int    i;
 392   real   *vec_params;
 393
 394   snew(vec_params,nvecs);
 395   if (str) {
 396     f0[0] = '\0';
 397     for(i=0; (i<nvecs); i++) {
 398       strcpy(f1,f0);  /*f0 is the format string for the "to-be-ignored" numbers*/
 399       strcat(f1,"%lf"); /*and f1 to read the actual number in this pass of the loop*/
 400       if (sscanf(str,f1,&d) != 1)
 401         gmx_fatal(FARGS,"Not enough elements for %s parameter (I need %d)",par,nvecs);
 402       vec_params[i] = d;
 403       tcap += d;
 404       strcat(f0,"%*s");
 405     }
 406   }
 407   return vec_params;
 408 }
 409
 410
 411 void init_edx(struct edix *edx) {
 412   edx->nr=0;
 413   snew(edx->x,1);
 414   snew(edx->anrs,1);
 415 }
 416
 417 void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
 418                 atom_id igro[],rvec *x,const char* structure)
 419 {
 420 /* filter2edx copies coordinates from x to edx which are given in index
 421 */
 422
 423    int pos,i;
 424    int ix=edx->nr;
 425    edx->nr+=nindex;
 426    srenew(edx->x,edx->nr);
 427    srenew(edx->anrs,edx->nr);
 428    for (i=0;i<nindex;i++,ix++) {
 429          for (pos=0; pos<ngro-1 && igro[pos]!=index[i] ; ++pos) {};  /*search element in igro*/
 430          if (igro[pos]!=index[i])
 431               gmx_fatal(FARGS,"Couldn't find atom with index %d in structure %s",index[i],structure);
 432          edx->anrs[ix]=index[i];
 433          copy_rvec(x[pos],edx->x[ix]);
 434    }
 435 }
 436
 437 void get_structure(t_atoms *atoms,const char *IndexFile,
 438                    const char *StructureFile,struct edix *edx,int nfit,
 439                    atom_id ifit[],int nav, atom_id index[])
 440 {
 441   atom_id *igro;  /*index corresponding to target or origin structure*/
 442   int ngro;
 443   int ntar;
 444   rvec *xtar;
 445   char  title[STRLEN];
 446   char* grpname;
 447
 448
 449   ntar=read_conffile(StructureFile,title,&xtar);
 450   printf("Select an index group of %d elements that corresponds to the atoms in the structure file %s\n",
 451             ntar,StructureFile);
 452   get_index(atoms,IndexFile,1,&ngro,&igro,&grpname);
 453   if (ngro!=ntar)
 454      gmx_fatal(FARGS,"You selected an index group with %d elements instead of %d",ngro,ntar);
 455   init_edx(edx);
 456   filter2edx(edx,nfit,ifit,ngro,igro,xtar,StructureFile);
 457
 458   /* If average and reference/fitting structure differ, append the average structure as well */
 459   if (ifit!=index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/
 460      filter2edx(edx,nav,index,ngro,igro,xtar,StructureFile);
 461 }
 462
 463 int gmx_make_edi(int argc,char *argv[])
 464 {
 465
 466   static const char *desc[] = {
 467       "[TT]make_edi[tt] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]",
 468       "based on eigenvectors of a covariance matrix ([TT]g_covar[tt]) or from a",
 469       "normal modes anaysis ([TT]g_nmeig[tt]).",
 470       "ED sampling can be used to manipulate the position along collective coordinates",
 471       "(eigenvectors) of (biological) macromolecules during a simulation. Particularly,",
 472       "it may be used to enhance the sampling efficiency of MD simulations by stimulating",
 473       "the system to explore new regions along these collective coordinates. A number",
 474       "of different algorithms are implemented to drive the system along the eigenvectors",
 475       "([TT]-linfix[tt], [TT]-linacc[tt], [TT]-radfix[tt], [TT]-radacc[tt], [TT]-radcon[tt]),",
 476       "to keep the position along a certain (set of) coordinate(s) fixed ([TT]-linfix[tt]),",
 477       "or to only monitor the projections of the positions onto",
 478       "these coordinates ([TT]-mon[tt]).[PAR]",
 479       "References:[BR]",
 480       "A. Amadei, A.B.M. Linssen, B.L. de Groot, D.M.F. van Aalten and ",
 481       "H.J.C. Berendsen; An efficient method for sampling the essential subspace ",
 482       "of proteins., J. Biomol. Struct. Dyn. 13:615-626 (1996)[BR]",
 483       "B.L. de Groot, A. Amadei, D.M.F. van Aalten and H.J.C. Berendsen; ",
 484       "Towards an exhaustive sampling of the configurational spaces of the ",
 485       "two forms of the peptide hormone guanylin,",
 486       "J. Biomol. Struct. Dyn. 13 : 741-751 (1996)[BR]",
 487       "B.L. de Groot, A.Amadei, R.M. Scheek, N.A.J. van Nuland and H.J.C. Berendsen; ",
 488       "An extended sampling of the configurational space of HPr from E. coli",
 489       "Proteins: Struct. Funct. Gen. 26: 314-322 (1996)",
 490       "[PAR]You will be prompted for one or more index groups that correspond to the eigenvectors,",
 491       "reference structure, target positions, etc.[PAR]",
 492
 493       "[TT]-mon[tt]: monitor projections of the coordinates onto selected eigenvectors.[PAR]",
 494       "[TT]-linfix[tt]: perform fixed-step linear expansion along selected eigenvectors.[PAR]",
 495       "[TT]-linacc[tt]: perform acceptance linear expansion along selected eigenvectors.",
 496       "(steps in the desired directions will be accepted, others will be rejected).[PAR]",
 497       "[TT]-radfix[tt]: perform fixed-step radius expansion along selected eigenvectors.[PAR]",
 498       "[TT]-radacc[tt]: perform acceptance radius expansion along selected eigenvectors.",
 499       "(steps in the desired direction will be accepted, others will be rejected).",
 500       "[BB]Note:[bb] by default the starting MD structure will be taken as origin of the first",
 501       "expansion cycle for radius expansion. If [TT]-ori[tt] is specified, you will be able",
 502       "to read in a structure file that defines an external origin.[PAR]",
 503       "[TT]-radcon[tt]: perform acceptance radius contraction along selected eigenvectors",
 504       "towards a target structure specified with [TT]-tar[tt].[PAR]",
 505       "NOTE: each eigenvector can be selected only once. [PAR]",
 506       "[TT]-outfrq[tt]: frequency (in steps) of writing out projections etc. to [TT].edo[tt] file[PAR]",
 507       "[TT]-slope[tt]: minimal slope in acceptance radius expansion. A new expansion",
 508       "cycle will be started if the spontaneous increase of the radius (in nm/step)",
 509       "is less than the value specified.[PAR]",
 510       "[TT]-maxedsteps[tt]: maximum number of steps per cycle in radius expansion",
 511       "before a new cycle is started.[PAR]",
 512       "Note on the parallel implementation: since ED sampling is a 'global' thing",
 513       "(collective coordinates etc.), at least on the 'protein' side, ED sampling",
 514       "is not very parallel-friendly from an implentation point of view. Because",
 515       "parallel ED requires some extra communication, expect the performance to be",
 516       "lower as in a free MD simulation, especially on a large number of nodes. [PAR]",
 517       "All output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a .edo file. In the output",
 518       "file, per OUTFRQ step the following information is present: [PAR]",
 519       "[TT]*[tt] the step number[BR]",
 520       "[TT]*[tt] the number of the ED dataset. ([BB]Note[bb] that you can impose multiple ED constraints in",
 521       "a single simulation (on different molecules) if several [TT].edi[tt] files were concatenated",
 522       "first. The constraints are applied in the order they appear in the [TT].edi[tt] file.) [BR]",
 523       "[TT]*[tt] RMSD (for atoms involved in fitting prior to calculating the ED constraints)[BR]",
 524       "* projections of the positions onto selected eigenvectors[BR]",
 525       "[PAR][PAR]",
 526       "FLOODING:[PAR]",
 527       "with [TT]-flood[tt], you can specify which eigenvectors are used to compute a flooding potential,",
 528       "which will lead to extra forces expelling the structure out of the region described",
 529       "by the covariance matrix. If you switch -restrain the potential is inverted and the structure",
 530       "is kept in that region.",
 531       "[PAR]",
 532       "The origin is normally the average structure stored in the [TT]eigvec.trr[tt] file.",
 533       "It can be changed with [TT]-ori[tt] to an arbitrary position in configurational space.",
 534       "With [TT]-tau[tt], [TT]-deltaF0[tt], and [TT]-Eflnull[tt] you control the flooding behaviour.",
 535       "Efl is the flooding strength, it is updated according to the rule of adaptive flooding.",
 536       "Tau is the time constant of adaptive flooding, high [GRK]tau[grk] means slow adaption (i.e. growth). ",
 537       "DeltaF0 is the flooding strength you want to reach after tau ps of simulation.",
 538       "To use constant Efl set [TT]-tau[tt] to zero.",
 539       "[PAR]",
 540       "[TT]-alpha[tt] is a fudge parameter to control the width of the flooding potential. A value of 2 has been found",
 541       "to give good results for most standard cases in flooding of proteins.",
 542       "[GRK]alpha[grk] basically accounts for incomplete sampling, if you sampled further the width of the ensemble would",
 543       "increase, this is mimicked by [GRK]alpha[grk] > 1.",
 544       "For restraining, [GRK]alpha[grk] < 1 can give you smaller width in the restraining potential.",
 545       "[PAR]",
 546       "RESTART and FLOODING:",
 547       "If you want to restart a crashed flooding simulation please find the values deltaF and Efl in",
 548       "the output file and manually put them into the [TT].edi[tt] file under DELTA_F0 and EFL_NULL."
 549   };
 550
 551     /* Save all the params in this struct and then save it in an edi file.
 552     * ignoring fields nmass,massnrs,mass,tmass,nfit,fitnrs,edo
 553     */
 554     static t_edipar edi_params;
 555
 556     enum  { evStepNr = evRADFIX + 1};
 557     static const char* evSelections[evNr]= {NULL,NULL,NULL,NULL,NULL,NULL};
 558     static const char* evOptions[evNr] = {"-linfix","-linacc","-flood","-radfix","-radacc","-radcon","-mon"};
 559     static const char* evParams[evStepNr] ={NULL,NULL};
 560     static const char* evStepOptions[evStepNr] = {"-linstep","-accdir","-not_used","-radstep"};
 561     static const char* ConstForceStr;
 562     static real* evStepList[evStepNr];
 563     static real  radfix=0.0;
 564     static real deltaF0=150;
 565     static real deltaF=0;
 566     static real tau=.1;
 567     static real constEfl=0.0;
 568     static real alpha=1;
 569     static int eqSteps=0;
 570     static int* listen[evNr];
 571     static real T=300.0;
 572     const real kB = 2.5 / 300.0; /* k_boltzmann in MD units */
 573     static gmx_bool bRestrain = FALSE;
 574     static gmx_bool bHesse=FALSE;
 575     static gmx_bool bHarmonic=FALSE;
 576     t_pargs pa[] = {
 577     { "-mon", FALSE, etSTR, {&evSelections[evMON]},
 578         "Indices of eigenvectors for projections of x (e.g. 1,2-5,9) or 1-100:10 means 1 11 21 31 ... 91" },
 579     { "-linfix", FALSE, etSTR, {&evSelections[0]},
 580         "Indices of eigenvectors for fixed increment linear sampling" },
 581     { "-linacc", FALSE, etSTR, {&evSelections[1]},
 582         "Indices of eigenvectors for acceptance linear sampling" },
 583     { "-radfix", FALSE, etSTR, {&evSelections[3]},
 584         "Indices of eigenvectors for fixed increment radius expansion" },
 585     { "-radacc", FALSE, etSTR, {&evSelections[4]},
 586         "Indices of eigenvectors for acceptance radius expansion" },
 587     { "-radcon", FALSE, etSTR, {&evSelections[5]},
 588         "Indices of eigenvectors for acceptance radius contraction" },
 589     { "-flood",  FALSE, etSTR, {&evSelections[2]},
 590         "Indices of eigenvectors for flooding"},
 591     { "-outfrq", FALSE, etINT, {&edi_params.outfrq},
 592         "Freqency (in steps) of writing output in [TT].edo[tt] file" },
 593     { "-slope", FALSE, etREAL, { &edi_params.slope},
 594         "Minimal slope in acceptance radius expansion"},
 595     { "-linstep", FALSE, etSTR, {&evParams[0]},
 596         "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"},
 597     { "-accdir", FALSE, etSTR, {&evParams[1]},
 598         "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"},
 599     { "-radstep", FALSE, etREAL, {&radfix},
 600         "Stepsize (nm/step) for fixed increment radius expansion"},
 601     { "-maxedsteps", FALSE, etINT, {&edi_params.maxedsteps},
 602         "Maximum number of steps per cycle" },
 603     { "-eqsteps", FALSE, etINT, {&eqSteps},
 604         "Number of steps to run without any perturbations "},
 605     { "-deltaF0", FALSE,etREAL, {&deltaF0},
 606         "Target destabilization energy for flooding"},
 607     { "-deltaF", FALSE, etREAL, {&deltaF},
 608         "Start deltaF with this parameter - default 0, nonzero values only needed for restart"},
 609     { "-tau", FALSE, etREAL, {&tau},
 610         "Coupling constant for adaption of flooding strength according to deltaF0, 0 = infinity i.e. constant flooding strength"},
 611     { "-Eflnull", FALSE, etREAL, {&constEfl},
 612         "The starting value of the flooding strength. The flooding strength is updated "
 613         "according to the adaptive flooding scheme. For a constant flooding strength use [TT]-tau[tt] 0. "},
 614     { "-T", FALSE, etREAL, {&T},
 615         "T is temperature, the value is needed if you want to do flooding "},
 616     { "-alpha",FALSE,etREAL,{&alpha},
 617         "Scale width of gaussian flooding potential with alpha^2 "},
 618     { "-restrain",FALSE, etBOOL, {&bRestrain},
 619         "Use the flooding potential with inverted sign -> effects as quasiharmonic restraining potential"},
 620     { "-hessian",FALSE, etBOOL, {&bHesse},
 621         "The eigenvectors and eigenvalues are from a Hessian matrix"},
 622     { "-harmonic",FALSE, etBOOL, {&bHarmonic},
 623         "The eigenvalues are interpreted as spring constant"},
 624     { "-constF", FALSE, etSTR, {&ConstForceStr},
 625         "Constant force flooding: manually set the forces for the eigenvectors selected with -flood "
 626         "(put in quotes! \"1.0 2.3 5.1 -3.1\"). No other flooding parameters are needed when specifying the forces directly."}
 627     };
 628 #define NPA asize(pa)
 629
 630     rvec       *xref1;
 631     int        nvec1,*eignr1=NULL;
 632     rvec       *xav1,**eigvec1=NULL;
 633     t_atoms    *atoms=NULL;
 634     int        nav;  /* Number of atoms in the average structure */
 635     char       *grpname;
 636     const char *indexfile;
 637     int        i;
 638     atom_id    *index,*ifit;
 639     int        nfit; /* Number of atoms in the reference/fit structure */
 640     int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */
 641     int nvecs;
 642     real *eigval1=NULL; /* in V3.3 this is parameter of read_eigenvectors */
 643
 644     const char *EdiFile;
 645     const char *TargetFile;
 646     const char *OriginFile;
 647     const char *EigvecFile;
 648
 649     output_env_t oenv;
 650
 651     /*to read topology file*/
 652     t_topology top;
 653     int        ePBC;
 654     char       title[STRLEN];
 655     matrix     topbox;
 656     rvec       *xtop;
 657     gmx_bool bTop, bFit1;
 658
 659     t_filenm fnm[] = {
 660     { efTRN, "-f",    "eigenvec",    ffREAD  },
 661     { efXVG, "-eig",  "eigenval",    ffOPTRD },
 662     { efTPS, NULL,    NULL,          ffREAD },
 663     { efNDX, NULL,    NULL,  ffOPTRD },
 664     { efSTX, "-tar", "target", ffOPTRD},
 665     { efSTX, "-ori", "origin", ffOPTRD},
 666     { efEDI, "-o", "sam", ffWRITE }
 667     };
 668 #define NFILE asize(fnm)
 669     edi_params.outfrq=100; edi_params.slope=0.0; edi_params.maxedsteps=0;
 670     CopyRight(stderr,argv[0]);
 671     parse_common_args(&argc,argv, 0 ,
 672                       NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);
 673
 674     indexfile=ftp2fn_null(efNDX,NFILE,fnm);
 675     EdiFile=ftp2fn(efEDI,NFILE,fnm);
 676     TargetFile      = opt2fn_null("-tar",NFILE,fnm);
 677     OriginFile      = opt2fn_null("-ori",NFILE,fnm);
 678
 679
 680     for (ev_class=0; ev_class<evNr; ++ev_class) {
 681         if (opt2parg_bSet(evOptions[ev_class],NPA,pa)) {
 682             /*get list of eigenvectors*/
 683             nvecs=sscan_list(&(listen[ev_class]),opt2parg_str(evOptions[ev_class],NPA,pa),evOptions[ev_class]);
 684             if (ev_class<evStepNr-2) {
 685                 /*if apropriate get list of stepsizes for these eigenvectors*/
 686                 if (opt2parg_bSet(evStepOptions[ev_class],NPA,pa))
 687                     evStepList[ev_class]=
 688                         scan_vecparams(opt2parg_str(evStepOptions[ev_class],NPA,pa),evStepOptions[ev_class],nvecs);
 689                 else { /*if list is not given fill with zeros */
 690                     snew(evStepList[ev_class],nvecs);
 691                     for (i=0; i<nvecs; i++)
 692                         evStepList[ev_class][i]=0.0;
 693                 }
 694             } else if (ev_class == evRADFIX && opt2parg_bSet(evStepOptions[ev_class],NPA,pa)) {
 695                 snew(evStepList[ev_class],nvecs);
 696                 for (i=0; i<nvecs; i++)
 697                     evStepList[ev_class][i]=radfix;
 698             } else if (ev_class == evFLOOD) {
 699                 snew(evStepList[ev_class],nvecs);
 700
 701                 /* Are we doing constant force flooding? In that case, we read in
 702                  * the fproj values from the command line */
 703                 if (opt2parg_bSet("-constF", NPA, pa))
 704                 {
 705                     evStepList[ev_class] = scan_vecparams(opt2parg_str("-constF",NPA,pa),"-constF",nvecs);
 706                 }
 707             } else {}; /*to avoid ambiguity   */
 708         } else { /* if there are no eigenvectors for this option set list to zero */
 709             listen[ev_class]=NULL;
 710             snew(listen[ev_class],1);
 711             listen[ev_class][0]=0;
 712         }
 713     }
 714
 715     /* print the interpreted list of eigenvectors - to give some feedback*/
 716     for (ev_class=0; ev_class<evNr; ++ev_class) {
 717         printf("Eigenvector list %7s consists of the indices: ",evOptions[ev_class]);
 718         i=0;
 719         while(listen[ev_class][i])
 720             printf("%d ",listen[ev_class][i++]);
 721         printf("\n");
 722     }
 723
 724     EigvecFile=NULL;
 725     EigvecFile=opt2fn("-f",NFILE,fnm);
 726
 727     /*read eigenvectors from eigvec.trr*/
 728     read_eigenvectors(EigvecFile,&nav,&bFit1,
 729                       &xref1,&edi_params.fitmas,&xav1,&edi_params.pcamas,&nvec1,&eignr1,&eigvec1,&eigval1);
 730
 731     bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
 732                        title,&top,&ePBC,&xtop,NULL,topbox,0);
 733     atoms=&top.atoms;
 734
 735
 736     printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",nav);
 737     get_index(atoms,indexfile,1,&i,&index,&grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
 738     if (i!=nav) {
 739         gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",
 740                   i,nav);
 741     }
 742     printf("\n");
 743
 744
 745     if (xref1==NULL)
 746     {
 747         if (bFit1)
 748         {
 749             /* if g_covar used different coordinate groups to fit and to do the PCA */
 750             printf("\nNote: the structure in %s should be the same\n"
 751                      "      as the one used for the fit in g_covar\n",ftp2fn(efTPS,NFILE,fnm));
 752             printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
 753         }
 754         else
 755         {
 756             printf("\nNote: Apparently no fitting was done in g_covar.\n"
 757                      "      However, you need to select a reference group for fitting in mdrun\n");
 758         }
 759         get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
 760         snew(xref1,nfit);
 761         for (i=0;i<nfit;i++)
 762             copy_rvec(xtop[ifit[i]],xref1[i]);
 763     }
 764     else
 765     {
 766         nfit=nav;
 767         ifit=index;
 768     }
 769
 770     if (opt2parg_bSet("-constF", NPA, pa))
 771     {
 772         /* Constant force flooding is special: Most of the normal flooding
 773          * options are not needed. */
 774         edi_params.flood.bConstForce = TRUE;
 775     }
 776     else
 777     {
 778         /* For normal flooding read eigenvalues and store them in evSteplist[evFLOOD] */
 779
 780         if ( listen[evFLOOD][0]!=0)
 781             read_eigenvalues(listen[evFLOOD],opt2fn("-eig",NFILE,fnm),evStepList[evFLOOD],bHesse,kB*T);
 782
 783         edi_params.flood.tau=tau;
 784         edi_params.flood.deltaF0=deltaF0;
 785         edi_params.flood.deltaF=deltaF;
 786         edi_params.presteps=eqSteps;
 787         edi_params.flood.kT=kB*T;
 788         edi_params.flood.bHarmonic=bHarmonic;
 789         if (bRestrain)
 790         {
 791             /* Trick: invert sign of Efl and alpha2 then this will give the same sign in the exponential and inverted sign outside */
 792             edi_params.flood.constEfl=-constEfl;
 793             edi_params.flood.alpha2=-sqr(alpha);
 794         }
 795         else
 796         {
 797             edi_params.flood.constEfl=constEfl;
 798             edi_params.flood.alpha2=sqr(alpha);
 799         }
 800     }
 801
 802     edi_params.ned=nav;
 803
 804   /*number of system atoms  */
 805   edi_params.nini=atoms->nr;
 806
 807
 808   /*store reference and average structure in edi_params*/
 809   make_t_edx(&edi_params.sref,nfit,xref1,ifit );
 810   make_t_edx(&edi_params.sav ,nav ,xav1 ,index);
 811
 812
 813   /* Store target positions in edi_params */
 814   if (opt2bSet("-tar",NFILE,fnm))
 815   {
 816       if (0 != listen[evFLOOD][0])
 817       {
 818           fprintf(stderr, "\nNote: Providing a TARGET structure has no effect when using flooding.\n"
 819                           "      You may want to use -ori to define the flooding potential center.\n\n");
 820       }
 821       get_structure(atoms,indexfile,TargetFile,&edi_params.star,nfit,ifit,nav,index);
 822   }
 823   else
 824   {
 825       make_t_edx(&edi_params.star,0,NULL,index);
 826   }
 827
 828   /* Store origin positions */
 829   if (opt2bSet("-ori",NFILE,fnm))
 830   {
 831       get_structure(atoms,indexfile,OriginFile,&edi_params.sori,nfit,ifit,nav,index);
 832   }
 833   else
 834   {
 835       make_t_edx(&edi_params.sori,0,NULL,index);
 836   }
 837
 838   /* Write edi-file */
 839   write_the_whole_thing(ffopen(EdiFile,"w"), &edi_params, eigvec1,nvec1, listen, evStepList);
 840   thanx(stderr);
 841
 842   return 0;
 843 }
 844