src/gromacs/gmxana/gmx_make_edi.c

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