src/gromacs/gmxana/gmx_make_edi.c

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