src/gromacs/gmxana/gmx_make_edi.cpp

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