[alexxy/gromacs.git] / src / tools / make_edi.c

/*
 *
 *                This source code is part of
 *
 *                 G   R   O   M   A   C   S
 *
 *          GROningen MAchine for Chemical Simulations
 *
 *                        VERSION 3.2.0
 *
 * The make_edi program was generously contributed by Oliver Lange, based
 * on the code from g_anaeig. You can reach him as olange@gwdg.de.
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * If you want to redistribute modifications, please consider that
 * scientific software is very special. Version control is crucial -
 * bugs must be traceable. We will be happy to consider code for
 * inclusion in the official distribution, but derived work must not
 * be called official GROMACS. Details are found in the README & COPYING
 * files - if they are missing, get the official version at www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the papers on the package - you can find them in the top README file.
 *
 * For more info, check our website at http://www.gromacs.org
 *
 * And Hey:
 * Gromacs Runs One Microsecond At Cannonball Speeds
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <math.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include "readinp.h"
#include "statutil.h"
#include "sysstuff.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
#include "gmx_fatal.h"
#include "vec.h"
#include "pbc.h"
#include "copyrite.h"
#include "futil.h"
#include "statutil.h"
#include "rdgroup.h"
#include "pdbio.h"
#include "confio.h"
#include "tpxio.h"
#include "matio.h"
#include "mshift.h"
#include "xvgr.h"
#include "do_fit.h"
#include "rmpbc.h"
#include "txtdump.h"
#include "eigio.h"
#include "index.h"

/* Suppress Cygwin compiler warnings from using newlib version of
 * ctype.h */
#ifdef GMX_CYGWIN
#undef isdigit
#endif

typedef struct
{
    real        deltaF0;
    gmx_bool    bHarmonic;
    gmx_bool    bConstForce;   /* Do constant force flooding instead of
                                  evaluating a flooding potential             */
    real        tau;
    real        deltaF;
    real        kT;
    real        constEfl;
    real        alpha2;
} t_edflood;


/* This type is for the average, reference, target, and origin structure   */
typedef struct edix
{
    int         nr;             /* number of atoms this structure contains */
    int         *anrs;          /* atom index numbers                      */
    rvec        *x;             /* positions                               */
    real        *sqrtm;         /* sqrt of the masses used for mass-
                                 * weighting of analysis                   */
} t_edix;


typedef struct edipar
{
    int         nini;           /* total Nr of atoms                    */
    gmx_bool    fitmas;         /* true if trans fit with cm            */
    gmx_bool    pcamas;         /* true if mass-weighted PCA            */
    int         presteps;       /* number of steps to run without any
                                 *    perturbations ... just monitoring */
    int         outfrq;         /* freq (in steps) of writing to edo    */
    int         maxedsteps;     /* max nr of steps per cycle            */
    struct edix sref;           /* reference positions, to these fitting
                                 * will be done                         */
    struct edix sav;            /* average positions                    */
    struct edix star;           /* target positions                     */
    struct edix sori;           /* origin positions                     */
    real        slope;          /* minimal slope in acceptance radexp   */
    int         ned;            /* Nr of atoms in essdyn buffer         */
    t_edflood   flood;          /* parameters especially for flooding   */
} t_edipar;



void make_t_edx(struct edix *edx, int natoms, rvec *pos, atom_id index[]) {
  edx->nr=natoms;
  edx->anrs=index;
  edx->x=pos;
}

void write_t_edx(FILE *fp, struct edix edx, const char *comment) {
 /*here we copy only the pointers into the t_edx struct
  no data is copied and edx.box is ignored  */
 int i;
  fprintf(fp,"#%s \n %d \n",comment,edx.nr);
  for (i=0;i<edx.nr;i++) {
    fprintf(fp,"%d  %f  %f  %f\n",(edx.anrs)[i]+1,(edx.x)[i][XX],(edx.x)[i][YY],(edx.x)[i][ZZ]);
  }
}

int sscan_list(int *list[], const char *str, const char *listname) {
 /*this routine scans a string of the form 1,3-6,9 and returns the
    selected numbers (in this case 1 3 4 5 6 9) in NULL-terminated array of integers.
    memory for this list will be allocated  in this routine -- sscan_list expects *list to
    be a NULL-Pointer

    listname is a string used in the errormessage*/


   int i,istep;
   char c;
   char *pos,*startpos,*step;
   int n=strlen(str);

   /*enums to define the different lexical stati */
   enum { sBefore, sNumber, sMinus, sRange, sZero, sSmaller, sError, sSteppedRange };

   int status=sBefore; /*status of the deterministic automat to scan str   */
   int number=0;
   int end_number=0;

   char *start=NULL; /*holds the string of the number behind a ','*/
   char *end=NULL; /*holds the string of the number behind a '-' */

   int nvecs=0; /* counts the number of vectors in the list*/

   step=NULL;
   snew(pos,n+4);
   startpos=pos;
   strcpy(pos,str);
   pos[n]=',';
   pos[n+1]='1';
   pos[n+2]='\0';

   *list = NULL;

   while ((c=*pos)!=0) {
     switch(status) {
        /* expect a number */
        case sBefore: if (isdigit(c)) {
              start=pos;
              status=sNumber;
              break;
           } else status=sError; break;

        /* have read a number, expect ',' or '-' */
        case sNumber: if (c==',') {
             /*store number*/
             srenew(*list,nvecs+1);
             (*list)[nvecs++]=number=strtol(start,NULL,10);
             status=sBefore;
             if (number==0)
                 status=sZero;
              break;
             } else if (c=='-') { status=sMinus; break; }
             else if (isdigit(c)) break;
             else status=sError; break;

        /* have read a '-' -> expect a number */
     case sMinus:
       if (isdigit(c)) {
         end=pos;
         status=sRange; break;
       } else status=sError; break;

     case sSteppedRange:
       if (isdigit(c)) {
         if (step) {
           status=sError; break;
         } else
           step=pos;
         status=sRange;
         break;
       } else status=sError; break;

        /* have read the number after a minus, expect ',' or ':' */
        case sRange:
            if (c==',') {
               /*store numbers*/
               end_number=strtol(end,NULL,10);
               number=strtol(start,NULL,10);
               status=sBefore;
               if (number==0) {
                  status=sZero; break;
               }
               if (end_number<=number) {
                  status=sSmaller; break;
               }
               srenew(*list,nvecs+end_number-number+1);
               if (step) {
                 istep=strtol(step,NULL,10);
                 step=NULL;
               } else istep=1;
               for (i=number;i<=end_number;i+=istep)
                    (*list)[nvecs++]=i;
               break;
            } else if (c==':') {
              status = sSteppedRange;
              break;
            } else if (isdigit(c)) break; else status=sError;  break;

       /* format error occured */
       case sError:
       gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d with char %c",listname,pos-startpos,*(pos-1));

       /* logical error occured */
       case sZero:
               gmx_fatal(FARGS,"Error in the list of eigenvectors for %s at pos %d: eigenvector 0 is not valid",listname,pos-startpos);
       case sSmaller:
               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);

     }
   ++pos; /* read next character */
   } /*scanner has finished */

   /* append zero to list of eigenvectors */
   srenew(*list,nvecs+1);
   (*list)[nvecs]=0;
   sfree(startpos);
   return nvecs;
} /*sscan_list*/

void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs,int nvec, const char *grouptitle,real steps[]) {
/* eig_list is a zero-terminated list of indices into the eigvecs array.
   eigvecs are coordinates of eigenvectors
   grouptitle to write in the comment line
   steps  -- array with stepsizes for evLINFIX, evLINACC and evRADACC
*/

  int n=0,i; rvec x;
  real sum;
  while (eig_list[n++]);  /*count selected eigenvecs*/

  fprintf(fp,"# NUMBER OF EIGENVECTORS + %s\n %d\n",grouptitle,n-1);

  /* write list of eigenvector indicess */
  for(n=0;eig_list[n];n++) {
    if (steps)
      fprintf(fp,"%8d   %g\n",eig_list[n],steps[n]);
    else
      fprintf(fp,"%8d   %g\n",eig_list[n],1.0);
  }
  n=0;

  /* dump coordinates of the selected eigenvectors */
  while (eig_list[n]) {
    sum=0;
    for (i=0; i<natoms; i++) {
      if (eig_list[n]>nvec)
        gmx_fatal(FARGS,"Selected eigenvector %d is higher than maximum number %d of available eigenvectors",eig_list[n],nvec);
      copy_rvec(eigvecs[eig_list[n]-1][i],x);
      sum+=norm2(x);
      fprintf(fp,"%8.5f %8.5f %8.5f\n",x[XX],x[YY],x[ZZ]);
    }
    n++;
  }
}


/*enum referring to the different lists of eigenvectors*/
enum { evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON , evMON,  evNr };
#define oldMAGIC 666
#define MAGIC 670


void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
                           int nvec, int *eig_listen[], real* evStepList[])
{
/* write edi-file */

    /*Header*/
    fprintf(fp,"#MAGIC\n %d \n#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
        MAGIC,edpars->nini,edpars->fitmas,edpars->pcamas);
    fprintf(fp,"#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
        edpars->outfrq,edpars->maxedsteps,edpars->slope);
    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",
        edpars->presteps,edpars->flood.deltaF0,edpars->flood.deltaF,edpars->flood.tau,edpars->flood.constEfl,
        edpars->flood.alpha2,edpars->flood.kT,edpars->flood.bHarmonic,edpars->flood.bConstForce);

    /* Average and reference positions */
    write_t_edx(fp,edpars->sref,"NREF, XREF");
    write_t_edx(fp,edpars->sav,"NAV, XAV");

    /*Eigenvectors */

    write_eigvec(fp, edpars->ned, eig_listen[evMON],eigvecs,nvec,"COMPONENTS GROUP 1",NULL);
    write_eigvec(fp, edpars->ned, eig_listen[evLINFIX],eigvecs,nvec,"COMPONENTS GROUP 2",evStepList[evLINFIX]);
    write_eigvec(fp, edpars->ned, eig_listen[evLINACC],eigvecs,nvec,"COMPONENTS GROUP 3",evStepList[evLINACC]);
    write_eigvec(fp, edpars->ned, eig_listen[evRADFIX],eigvecs,nvec,"COMPONENTS GROUP 4",evStepList[evRADFIX]);
    write_eigvec(fp, edpars->ned, eig_listen[evRADACC],eigvecs,nvec,"COMPONENTS GROUP 5",NULL);
    write_eigvec(fp, edpars->ned, eig_listen[evRADCON],eigvecs,nvec,"COMPONENTS GROUP 6",NULL);
    write_eigvec(fp, edpars->ned, eig_listen[evFLOOD],eigvecs,nvec,"COMPONENTS GROUP 7",evStepList[evFLOOD]);


    /*Target and Origin positions */
    write_t_edx(fp,edpars->star,"NTARGET, XTARGET");
    write_t_edx(fp,edpars->sori,"NORIGIN, XORIGIN");
}

int read_conffile(const char *confin,char *title,rvec *x[])
{
/* read coordinates out of STX file  */
  int natoms;
  t_atoms  confat;
  matrix box;
  printf("read coordnumber from file %s\n",confin);
  get_stx_coordnum(confin,&natoms);
  printf("number of coordinates in file %d\n",natoms);
/*  if (natoms != ncoords)
     gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n"
           "             does not match topology (= %d)",
           confin,natoms,ncoords);
  else {*/
    /* make space for coordinates and velocities */
    init_t_atoms(&confat,natoms,FALSE);
    snew(*x,natoms);
    read_stx_conf(confin,title,&confat,*x,NULL,NULL,box);
    return natoms;
}


void read_eigenvalues(int vecs[],const char *eigfile, real values[],
                      gmx_bool bHesse, real kT)
{
  int  neig,nrow,i;
  double **eigval;

  neig = read_xvg(eigfile,&eigval,&nrow);

  fprintf(stderr,"Read %d eigenvalues\n",neig);
  for(i=bHesse ? 6 : 0 ; i<neig; i++) {
    if (eigval[1][i] < -0.001 && bHesse)
      fprintf(stderr,
              "WARNING: The Hessian Matrix has negative eigenvalue %f, we set it to zero (no flooding in this direction)\n\n",eigval[1][i]);

    if (eigval[1][i] < 0)
      eigval[1][i] = 0;
  }
  if (bHesse)
    for (i=0; vecs[i]; i++) {
      if (vecs[i]<7)
        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");
      values[i]=eigval[1][vecs[i]-1]/kT;
    }
  else
    for (i=0; vecs[i]; i++) {
      if (vecs[i]>(neig-6))
        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");
      values[i]=1/eigval[1][vecs[i]-1];
    }
  /* free memory */
  for (i=0; i<nrow; i++)
    sfree(eigval[i]);
  sfree(eigval);
}


static real *scan_vecparams(const char *str,const char * par, int nvecs)
{
  char   f0[256],f1[256];             /*format strings adapted every pass of the loop*/
  double d,tcap=0;
  int    i;
  real   *vec_params;

  snew(vec_params,nvecs);
  if (str) {
    f0[0] = '\0';
    for(i=0; (i<nvecs); i++) {
      strcpy(f1,f0);  /*f0 is the format string for the "to-be-ignored" numbers*/
      strcat(f1,"%lf"); /*and f1 to read the actual number in this pass of the loop*/
      if (sscanf(str,f1,&d) != 1)
        gmx_fatal(FARGS,"Not enough elements for %s parameter (I need %d)",par,nvecs);
      vec_params[i] = d;
      tcap += d;
      strcat(f0,"%*s");
    }
  }
  return vec_params;
}


void init_edx(struct edix *edx) {
  edx->nr=0;
  snew(edx->x,1);
  snew(edx->anrs,1);
}

void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
                atom_id igro[],rvec *x,const char* structure)
{
/* filter2edx copies coordinates from x to edx which are given in index
*/

   int pos,i;
   int ix=edx->nr;
   edx->nr+=nindex;
   srenew(edx->x,edx->nr);
   srenew(edx->anrs,edx->nr);
   for (i=0;i<nindex;i++,ix++) {
         for (pos=0; pos<ngro-1 && igro[pos]!=index[i] ; ++pos) {};  /*search element in igro*/
         if (igro[pos]!=index[i])
              gmx_fatal(FARGS,"Couldn't find atom with index %d in structure %s",index[i],structure);
         edx->anrs[ix]=index[i];
         copy_rvec(x[pos],edx->x[ix]);
   }
}

void get_structure(t_atoms *atoms,const char *IndexFile,
                   const char *StructureFile,struct edix *edx,int nfit,
                   atom_id ifit[],int natoms, atom_id index[])
{
  atom_id *igro;  /*index corresponding to target or origin structure*/
  int ngro;
  int ntar;
  rvec *xtar;
  char  title[STRLEN];
  char* grpname;


  ntar=read_conffile(StructureFile,title,&xtar);
  printf("Select an index group of %d elements that corresponds to the atoms in the structure file %s\n",
            ntar,StructureFile);
  get_index(atoms,IndexFile,1,&ngro,&igro,&grpname);
  if (ngro!=ntar)
     gmx_fatal(FARGS,"You selected an index group with %d elements instead of %d",ngro,ntar);
  init_edx(edx);
  filter2edx(edx,nfit,ifit,ngro,igro,xtar,StructureFile);
  if (ifit!=index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/
     filter2edx(edx,natoms,index,ngro,igro,xtar,StructureFile);
}

int main(int argc,char *argv[])
{

  static const char *desc[] = {
      "[TT]make_edi[tt] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]",
      "based on eigenvectors of a covariance matrix ([TT]g_covar[tt]) or from a",
      "normal modes anaysis ([TT]g_nmeig[tt]).",
      "ED sampling can be used to manipulate the position along collective coordinates",
      "(eigenvectors) of (biological) macromolecules during a simulation. Particularly,",
      "it may be used to enhance the sampling efficiency of MD simulations by stimulating",
      "the system to explore new regions along these collective coordinates. A number",
      "of different algorithms are implemented to drive the system along the eigenvectors",
      "([TT]-linfix[tt], [TT]-linacc[tt], [TT]-radfix[tt], [TT]-radacc[tt], [TT]-radcon[tt]),",
      "to keep the position along a certain (set of) coordinate(s) fixed ([TT]-linfix[tt]),",
      "or to only monitor the projections of the positions onto",
      "these coordinates ([TT]-mon[tt]).[PAR]",
      "References:[BR]",
      "A. Amadei, A.B.M. Linssen, B.L. de Groot, D.M.F. van Aalten and ",
      "H.J.C. Berendsen; An efficient method for sampling the essential subspace ",
      "of proteins., J. Biomol. Struct. Dyn. 13:615-626 (1996)[BR]",
      "B.L. de Groot, A. Amadei, D.M.F. van Aalten and H.J.C. Berendsen; ",
      "Towards an exhaustive sampling of the configurational spaces of the ",
      "two forms of the peptide hormone guanylin,",
      "J. Biomol. Struct. Dyn. 13 : 741-751 (1996)[BR]",
      "B.L. de Groot, A.Amadei, R.M. Scheek, N.A.J. van Nuland and H.J.C. Berendsen; ",
      "An extended sampling of the configurational space of HPr from E. coli",
      "Proteins: Struct. Funct. Gen. 26: 314-322 (1996)",
      "[PAR]You will be prompted for one or more index groups that correspond to the eigenvectors,",
      "reference structure, target positions, etc.[PAR]",

      "[TT]-mon[tt]: monitor projections of the coordinates onto selected eigenvectors.[PAR]",
      "[TT]-linfix[tt]: perform fixed-step linear expansion along selected eigenvectors.[PAR]",
      "[TT]-linacc[tt]: perform acceptance linear expansion along selected eigenvectors.",
      "(steps in the desired directions will be accepted, others will be rejected).[PAR]",
      "[TT]-radfix[tt]: perform fixed-step radius expansion along selected eigenvectors.[PAR]",
      "[TT]-radacc[tt]: perform acceptance radius expansion along selected eigenvectors.",
      "(steps in the desired direction will be accepted, others will be rejected).",
      "[BB]Note:[bb] by default the starting MD structure will be taken as origin of the first",
      "expansion cycle for radius expansion. If [TT]-ori[tt] is specified, you will be able",
      "to read in a structure file that defines an external origin.[PAR]",
      "[TT]-radcon[tt]: perform acceptance radius contraction along selected eigenvectors",
      "towards a target structure specified with [TT]-tar[tt].[PAR]",
      "NOTE: each eigenvector can be selected only once. [PAR]",
      "[TT]-outfrq[tt]: frequency (in steps) of writing out projections etc. to [TT].edo[tt] file[PAR]",
      "[TT]-slope[tt]: minimal slope in acceptance radius expansion. A new expansion",
      "cycle will be started if the spontaneous increase of the radius (in nm/step)",
      "is less than the value specified.[PAR]",
      "[TT]-maxedsteps[tt]: maximum number of steps per cycle in radius expansion",
      "before a new cycle is started.[PAR]",
      "Note on the parallel implementation: since ED sampling is a 'global' thing",
      "(collective coordinates etc.), at least on the 'protein' side, ED sampling",
      "is not very parallel-friendly from an implentation point of view. Because",
      "parallel ED requires some extra communication, expect the performance to be",
      "lower as in a free MD simulation, especially on a large number of nodes. [PAR]",
      "All output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a .edo file. In the output",
      "file, per OUTFRQ step the following information is present: [PAR]",
      "[TT]*[tt] the step number[BR]",
      "[TT]*[tt] the number of the ED dataset. ([BB]Note[bb] that you can impose multiple ED constraints in",
      "a single simulation (on different molecules) if several [TT].edi[tt] files were concatenated",
      "first. The constraints are applied in the order they appear in the [TT].edi[tt] file.) [BR]",
      "[TT]*[tt] RMSD (for atoms involved in fitting prior to calculating the ED constraints)[BR]",
      "* projections of the positions onto selected eigenvectors[BR]",
      "[PAR][PAR]",
      "FLOODING:[PAR]",
      "with [TT]-flood[tt], you can specify which eigenvectors are used to compute a flooding potential,",
      "which will lead to extra forces expelling the structure out of the region described",
      "by the covariance matrix. If you switch -restrain the potential is inverted and the structure",
      "is kept in that region.",
      "[PAR]",
      "The origin is normally the average structure stored in the [TT]eigvec.trr[tt] file.",
      "It can be changed with [TT]-ori[tt] to an arbitrary position in configurational space.",
      "With [TT]-tau[tt], [TT]-deltaF0[tt], and [TT]-Eflnull[tt] you control the flooding behaviour.",
      "Efl is the flooding strength, it is updated according to the rule of adaptive flooding.",
      "Tau is the time constant of adaptive flooding, high [GRK]tau[grk] means slow adaption (i.e. growth). ",
      "DeltaF0 is the flooding strength you want to reach after tau ps of simulation.",
      "To use constant Efl set [TT]-tau[tt] to zero.",
      "[PAR]",
      "[TT]-alpha[tt] is a fudge parameter to control the width of the flooding potential. A value of 2 has been found",
      "to give good results for most standard cases in flooding of proteins.",
      "[GRK]alpha[grk] basically accounts for incomplete sampling, if you sampled further the width of the ensemble would",
      "increase, this is mimicked by [GRK]alpha[grk] > 1.",
      "For restraining, [GRK]alpha[grk] < 1 can give you smaller width in the restraining potential.",
      "[PAR]",
      "RESTART and FLOODING:",
      "If you want to restart a crashed flooding simulation please find the values deltaF and Efl in",
      "the output file and manually put them into the [TT].edi[tt] file under DELTA_F0 and EFL_NULL."
  };

    /* Save all the params in this struct and then save it in an edi file.
    * ignoring fields nmass,massnrs,mass,tmass,nfit,fitnrs,edo
    */
    static t_edipar edi_params;

    enum  { evStepNr = evRADFIX + 1};
    static const char* evSelections[evNr]= {NULL,NULL,NULL,NULL,NULL,NULL};
    static const char* evOptions[evNr] = {"-linfix","-linacc","-flood","-radfix","-radacc","-radcon","-mon"};
    static const char* evParams[evStepNr] ={NULL,NULL};
    static const char* evStepOptions[evStepNr] = {"-linstep","-accdir","-not_used","-radstep"};
    static const char* ConstForceStr;
    static real* evStepList[evStepNr];
    static real  radfix=0.0;
    static real deltaF0=150;
    static real deltaF=0;
    static real tau=.1;
    static real constEfl=0.0;
    static real alpha=1;
    static int eqSteps=0;
    static int* listen[evNr];
    static real T=300.0;
    const real kB = 2.5 / 300.0; /* k_boltzmann in MD units */
    static gmx_bool bRestrain = FALSE;
    static gmx_bool bHesse=FALSE;
    static gmx_bool bHarmonic=FALSE;
    t_pargs pa[] = {
    { "-mon", FALSE, etSTR, {&evSelections[evMON]},
        "Indices of eigenvectors for projections of x (e.g. 1,2-5,9) or 1-100:10 means 1 11 21 31 ... 91" },
    { "-linfix", FALSE, etSTR, {&evSelections[0]},
        "Indices of eigenvectors for fixed increment linear sampling" },
    { "-linacc", FALSE, etSTR, {&evSelections[1]},
        "Indices of eigenvectors for acceptance linear sampling" },
    { "-radfix", FALSE, etSTR, {&evSelections[3]},
        "Indices of eigenvectors for fixed increment radius expansion" },
    { "-radacc", FALSE, etSTR, {&evSelections[4]},
        "Indices of eigenvectors for acceptance radius expansion" },
    { "-radcon", FALSE, etSTR, {&evSelections[5]},
        "Indices of eigenvectors for acceptance radius contraction" },
    { "-flood",  FALSE, etSTR, {&evSelections[2]},
        "Indices of eigenvectors for flooding"},
    { "-outfrq", FALSE, etINT, {&edi_params.outfrq},
        "Freqency (in steps) of writing output in [TT].edo[tt] file" },
    { "-slope", FALSE, etREAL, { &edi_params.slope},
        "Minimal slope in acceptance radius expansion"},
    { "-linstep", FALSE, etSTR, {&evParams[0]},
        "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"},
    { "-accdir", FALSE, etSTR, {&evParams[1]},
        "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"},
    { "-radstep", FALSE, etREAL, {&radfix},
        "Stepsize (nm/step) for fixed increment radius expansion"},
    { "-maxedsteps", FALSE, etINT, {&edi_params.maxedsteps},
        "Maximum number of steps per cycle" },
    { "-eqsteps", FALSE, etINT, {&eqSteps},
        "Number of steps to run without any perturbations "},
    { "-deltaF0", FALSE,etREAL, {&deltaF0},
        "Target destabilization energy for flooding"},
    { "-deltaF", FALSE, etREAL, {&deltaF},
        "Start deltaF with this parameter - default 0, nonzero values only needed for restart"},
    { "-tau", FALSE, etREAL, {&tau},
        "Coupling constant for adaption of flooding strength according to deltaF0, 0 = infinity i.e. constant flooding strength"},
    { "-Eflnull", FALSE, etREAL, {&constEfl},
        "The starting value of the flooding strength. The flooding strength is updated "
        "according to the adaptive flooding scheme. For a constant flooding strength use [TT]-tau[tt] 0. "},
    { "-T", FALSE, etREAL, {&T},
        "T is temperature, the value is needed if you want to do flooding "},
    { "-alpha",FALSE,etREAL,{&alpha},
        "Scale width of gaussian flooding potential with alpha^2 "},
    { "-restrain",FALSE, etBOOL, {&bRestrain},
        "Use the flooding potential with inverted sign -> effects as quasiharmonic restraining potential"},
    { "-hessian",FALSE, etBOOL, {&bHesse},
        "The eigenvectors and eigenvalues are from a Hessian matrix"},
    { "-harmonic",FALSE, etBOOL, {&bHarmonic},
        "The eigenvalues are interpreted as spring constant"},
    { "-constF", FALSE, etSTR, {&ConstForceStr},
        "Constant force flooding: manually set the forces for the eigenvectors selected with -flood "
        "(put in quotes! \"1.0 2.3 5.1 -3.1\"). No other flooding parameters are needed when specifying the forces directly."}
    };
#define NPA asize(pa)

    rvec       *xref1;
    int        nvec1,*eignr1=NULL;
    rvec       *xav1,**eigvec1=NULL;
    t_atoms    *atoms=NULL;
    int natoms;
    char       *grpname;
    const char *indexfile;
    int        i;
    atom_id    *index,*ifit;
    int        nfit;
    int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */
    int nvecs;
    real *eigval1=NULL; /* in V3.3 this is parameter of read_eigenvectors */

    const char *EdiFile;
    const char *TargetFile;
    const char *OriginFile;
    const char *EigvecFile;

    output_env_t oenv;

    /*to read topology file*/
    t_topology top;
    int        ePBC;
    char       title[STRLEN];
    matrix     topbox;
    rvec       *xtop;
    gmx_bool bTop, bFit1;

    t_filenm fnm[] = {
    { efTRN, "-f",    "eigenvec",    ffREAD  },
    { efXVG, "-eig",  "eigenval",    ffOPTRD },
    { efTPS, NULL,    NULL,          ffREAD },
    { efNDX, NULL,    NULL,  ffOPTRD },
    { efSTX, "-tar", "target", ffOPTRD},
    { efSTX, "-ori", "origin", ffOPTRD},
    { efEDI, "-o", "sam", ffWRITE }
    };
#define NFILE asize(fnm)
    edi_params.outfrq=100; edi_params.slope=0.0; edi_params.maxedsteps=0;
    CopyRight(stderr,argv[0]);
    parse_common_args(&argc,argv, 0 ,
                      NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL,&oenv);

    indexfile=ftp2fn_null(efNDX,NFILE,fnm);
    EdiFile=ftp2fn(efEDI,NFILE,fnm);
    TargetFile      = opt2fn_null("-tar",NFILE,fnm);
    OriginFile      = opt2fn_null("-ori",NFILE,fnm);


    for (ev_class=0; ev_class<evNr; ++ev_class) {
        if (opt2parg_bSet(evOptions[ev_class],NPA,pa)) {
            /*get list of eigenvectors*/
            nvecs=sscan_list(&(listen[ev_class]),opt2parg_str(evOptions[ev_class],NPA,pa),evOptions[ev_class]);
            if (ev_class<evStepNr-2) {
                /*if apropriate get list of stepsizes for these eigenvectors*/
                if (opt2parg_bSet(evStepOptions[ev_class],NPA,pa))
                    evStepList[ev_class]=
                        scan_vecparams(opt2parg_str(evStepOptions[ev_class],NPA,pa),evStepOptions[ev_class],nvecs);
                else { /*if list is not given fill with zeros */
                    snew(evStepList[ev_class],nvecs);
                    for (i=0; i<nvecs; i++)
                        evStepList[ev_class][i]=0.0;
                }
            } else if (ev_class == evRADFIX && opt2parg_bSet(evStepOptions[ev_class],NPA,pa)) {
                snew(evStepList[ev_class],nvecs);
                for (i=0; i<nvecs; i++)
                    evStepList[ev_class][i]=radfix;
            } else if (ev_class == evFLOOD) {
                snew(evStepList[ev_class],nvecs);

                /* Are we doing constant force flooding? In that case, we read in
                 * the fproj values from the command line */
                if (opt2parg_bSet("-constF", NPA, pa))
                {
                    evStepList[ev_class] = scan_vecparams(opt2parg_str("-constF",NPA,pa),"-constF",nvecs);
                }
            } else {}; /*to avoid ambiguity   */
        } else { /* if there are no eigenvectors for this option set list to zero */
            listen[ev_class]=NULL;
            snew(listen[ev_class],1);
            listen[ev_class][0]=0;
        }
    }

    /* print the interpreted list of eigenvectors - to give some feedback*/
    for (ev_class=0; ev_class<evNr; ++ev_class) {
        printf("Eigenvector list %7s consists of the indices: ",evOptions[ev_class]);
        i=0;
        while(listen[ev_class][i])
            printf("%d ",listen[ev_class][i++]);
        printf("\n");
    }

    EigvecFile=NULL;
    EigvecFile=opt2fn("-f",NFILE,fnm);

    /*read eigenvectors from eigvec.trr*/
    read_eigenvectors(EigvecFile,&natoms,&bFit1,
                      &xref1,&edi_params.fitmas,&xav1,&edi_params.pcamas,&nvec1,&eignr1,&eigvec1,&eigval1);

    bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),
                       title,&top,&ePBC,&xtop,NULL,topbox,0);
    atoms=&top.atoms;


    printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
    get_index(atoms,indexfile,1,&i,&index,&grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
    if (i!=natoms) {
        gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",
                  i,natoms);
    }
    printf("\n");


    if (xref1==NULL)
    {
        if (bFit1)
        {
            /* if g_covar used different coordinate groups to fit and to do the PCA */
            printf("\nNote: the structure in %s should be the same\n"
                     "      as the one used for the fit in g_covar\n",ftp2fn(efTPS,NFILE,fnm));
            printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
        }
        else
        {
            printf("\nNote: Apparently no fitting was done in g_covar.\n"
                     "      However, you need to select a reference group for fitting in mdrun\n");
        }
        get_index(atoms,indexfile,1,&nfit,&ifit,&grpname);
        snew(xref1,nfit);
        for (i=0;i<nfit;i++)
            copy_rvec(xtop[ifit[i]],xref1[i]);
    }
    else
    {
        nfit=natoms;
        ifit=index;
    }

    if (opt2parg_bSet("-constF", NPA, pa))
    {
        /* Constant force flooding is special: Most of the normal flooding
         * options are not needed. */
        edi_params.flood.bConstForce = TRUE;
    }
    else
    {
        /* For normal flooding read eigenvalues and store them in evSteplist[evFLOOD] */

        if ( listen[evFLOOD][0]!=0)
            read_eigenvalues(listen[evFLOOD],opt2fn("-eig",NFILE,fnm),evStepList[evFLOOD],bHesse,kB*T);

        edi_params.flood.tau=tau;
        edi_params.flood.deltaF0=deltaF0;
        edi_params.flood.deltaF=deltaF;
        edi_params.presteps=eqSteps;
        edi_params.flood.kT=kB*T;
        edi_params.flood.bHarmonic=bHarmonic;
        if (bRestrain)
        {
            /* Trick: invert sign of Efl and alpha2 then this will give the same sign in the exponential and inverted sign outside */
            edi_params.flood.constEfl=-constEfl;
            edi_params.flood.alpha2=-sqr(alpha);
        }
        else
        {
            edi_params.flood.constEfl=constEfl;
            edi_params.flood.alpha2=sqr(alpha);
        }
    }

    edi_params.ned=natoms;

  /*number of system atoms  */
  edi_params.nini=atoms->nr;


  /*store reference and average structure in edi_params*/
  make_t_edx(&edi_params.sref,nfit,xref1,ifit);
  make_t_edx(&edi_params.sav,natoms,xav1,index);


  /* Store target positions in edi_params */
  if (opt2bSet("-tar",NFILE,fnm))
  {
      if (0 != listen[evFLOOD][0])
      {
          fprintf(stderr, "\nNote: Providing a TARGET structure has no effect when using flooding.\n"
                          "      You may want to use -ori to define the flooding potential center.\n\n");
      }
      get_structure(atoms,indexfile,TargetFile,&edi_params.star,nfit,ifit,natoms,index);
  }
  else
  {
      make_t_edx(&edi_params.star,0,NULL,index);
  }

  /* Store origin positions */
  if (opt2bSet("-ori",NFILE,fnm))
  {
      get_structure(atoms,indexfile,OriginFile,&edi_params.sori,nfit,ifit,natoms,index);
  }
  else
  {
      make_t_edx(&edi_params.sori,0,NULL,index);
  }

  /* Write edi-file */
  write_the_whole_thing(ffopen(EdiFile,"w"), &edi_params, eigvec1,nvec1, listen, evStepList);
  thanx(stderr);

  return 0;
}