Checks for zero-mass atoms in ED, better ordering of options in make_edi
[alexxy/gromacs.git] / src / tools / make_edi.c
index 263a557d7e23163e0d1faef84c72bee364bd21c4..6ff8af6bbaeb37bd39273c37c27e7b4e7b66ebe7 100644 (file)
@@ -5,7 +5,7 @@
  *                 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
@@ -15,7 +15,7 @@
  * 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
 #include "rmpbc.h"
 #include "txtdump.h"
 #include "eigio.h"
-#include "index.h" 
+#include "index.h"
 
 
 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        kT;
     real        constEfl;
-    real        alpha2; 
+    real        alpha2;
 } t_edflood;
 
 
@@ -93,7 +95,7 @@ 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   
+    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            */
@@ -141,26 +143,26 @@ int sscan_list(int *list[], const char *str, const char *listname) {
 
    /*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]=',';
    pos[n+1]='1';
    pos[n+2]='\0';
 
    *list = NULL;
-     
+
    while ((c=*pos)!=0) {
      switch(status) {
         /* expect a number */
@@ -184,22 +186,22 @@ int sscan_list(int *list[], const char *str, const char *listname) {
              else status=sError; break;
 
         /* have read a '-' -> expect a number */
-     case sMinus: 
+     case sMinus:
        if (isdigit(c)) {
         end=pos;
         status=sRange; break;
        } else status=sError; break;
-       
+
      case sSteppedRange:
        if (isdigit(c)) {
         if (step) {
-          status=sError; break; 
+          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==',') {
@@ -228,7 +230,7 @@ int sscan_list(int *list[], const char *str, const char *listname) {
 
        /* 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)); 
+       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:
@@ -246,7 +248,7 @@ int sscan_list(int *list[], const char *str, const char *listname) {
    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
@@ -257,18 +259,18 @@ void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs,int nvec,
   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],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;
@@ -277,21 +279,21 @@ void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs,int 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]);      
-    }    
+      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,  evEND };
+enum { evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON , evMON,  evNr };
 #define oldMAGIC 666
-#define MAGIC 669
+#define MAGIC 670
 
 
-void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs, 
-                           int nvec, int *eig_listen[], real* evStepList[]) 
+void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
+                           int nvec, int *eig_listen[], real* evStepList[])
 {
 /* write edi-file */
 
@@ -300,15 +302,16 @@ void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
         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",
-        edpars->presteps,edpars->flood.deltaF0,edpars->flood.deltaF,edpars->flood.tau,edpars->flood.constEfl,edpars->flood.alpha2,edpars->flood.kT,edpars->flood.bHarmonic);
-    
+    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]);
@@ -321,9 +324,9 @@ void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
     /*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[]) 
+int read_conffile(const char *confin,char *title,rvec *x[])
 {
 /* read coordinates out of STX file  */
   int natoms;
@@ -342,23 +345,23 @@ int read_conffile(const char *confin,char *title,rvec *x[])
     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) 
+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++) { 
+  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;
   }
@@ -400,9 +403,9 @@ static real *scan_vecparams(const char *str,const char * par, int nvecs)
       tcap += d;
       strcat(f0,"%*s");
     }
-  }  
+  }
   return vec_params;
-}    
+}
 
 
 void init_edx(struct edix *edx) {
@@ -411,12 +414,12 @@ void init_edx(struct edix *edx) {
   snew(edx->anrs,1);
 }
 
-void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro, 
-                atom_id igro[],rvec *x,const char* structure) 
+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;
@@ -433,7 +436,7 @@ void filter2edx(struct edix *edx,int nindex, atom_id index[],int ngro,
 
 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 ifit[],int natoms, atom_id index[])
 {
   atom_id *igro;  /*index corresponding to target or origin structure*/
   int ngro;
@@ -441,7 +444,7 @@ void get_structure(t_atoms *atoms,const char *IndexFile,
   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",
@@ -460,7 +463,7 @@ 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", 
+      "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,",
@@ -484,7 +487,7 @@ int main(int argc,char *argv[])
       "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.",
@@ -507,7 +510,7 @@ int main(int argc,char *argv[])
       "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 much extra communication, expect the performance to be",
+      "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]",
@@ -542,18 +545,19 @@ int main(int argc,char *argv[])
       "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  { evSTEPEND = evRADFIX + 1}; 
-    static const char* evSelections[evEND]= {NULL,NULL,NULL,NULL,NULL,NULL};
-    static const char* evOptions[evEND] = {"-linfix","-linacc","-flood","-radfix","-radacc","-radcon","-mon"};
-    static const char* evParams[evSTEPEND] ={NULL,NULL};
-    static const char* evStepOptions[evSTEPEND] = {"-linstep","-accdir","-not_used","-radstep"};
-    static real* evStepList[evSTEPEND];
+    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;
@@ -561,7 +565,7 @@ int main(int argc,char *argv[])
     static real constEfl=0.0;
     static real alpha=1;
     static int eqSteps=0;
-    static int* listen[evEND];
+    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;
@@ -574,49 +578,53 @@ int main(int argc,char *argv[])
         "Indices of eigenvectors for fixed increment linear sampling" },
     { "-linacc", FALSE, etSTR, {&evSelections[1]},
         "Indices of eigenvectors for acceptance linear sampling" },
-    { "-flood",  FALSE, etSTR, {&evSelections[2]},
-        "Indices of eigenvectors for flooding"},
     { "-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},
-        "Max nr of steps per cycle" },
+        "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  - used for flooding"},
+        "Target destabilization energy for flooding"},
     { "-deltaF", FALSE, etREAL, {&deltaF},
-        "Start deltaF with this parameter - default 0, i.e. nonzero values only needed for restart"},
-    { "-tau", FALSE, etREAL, {&tau}, 
+        "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"},
-    { "-eqsteps", FALSE, etINT, {&eqSteps},
-        "Number of steps to run without any perturbations "},
     { "-Eflnull", FALSE, etREAL, {&constEfl},
-        "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. "},
+        "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 "},
-    { "-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"},
     { "-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}, 
+    { "-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;
@@ -630,14 +638,14 @@ int main(int argc,char *argv[])
     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;
@@ -645,10 +653,10 @@ int main(int argc,char *argv[])
     matrix     topbox;
     rvec       *xtop;
     gmx_bool bTop, bFit1;
-    
+
     t_filenm fnm[] = {
     { efTRN, "-f",    "eigenvec",    ffREAD  },
-    { efXVG, "-eig",  "eigenval",    ffOPTRD }, 
+    { efXVG, "-eig",  "eigenval",    ffOPTRD },
     { efTPS, NULL,    NULL,          ffREAD },
     { efNDX, NULL,    NULL,  ffOPTRD },
     { efSTX, "-tar", "target", ffOPTRD},
@@ -660,64 +668,69 @@ int main(int argc,char *argv[])
     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<evEND; ++ev_class) {
+
+
+    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<evSTEPEND-2) {
+            if (ev_class<evStepNr-2) {
                 /*if apropriate get list of stepsizes for these eigenvectors*/
-                if (opt2parg_bSet(evStepOptions[ev_class],NPA,pa)) 
+                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++) 
+                    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<evEND; ++ev_class) {
+    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); 
-    
-    
+
+    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) {
@@ -725,8 +738,8 @@ int main(int argc,char *argv[])
                   i,natoms);
     }
     printf("\n");
-    
-    
+
+
     if (xref1==NULL)
     {
         if (bFit1)
@@ -752,39 +765,49 @@ int main(int argc,char *argv[])
         ifit=index;
     }
 
-    /* if -flood 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);
-  /*number of eigenvectors*/
-  edi_params.ned=natoms;
-  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);
-  }
-  
+    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))
   {