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