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