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