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