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