src/gromacs/gmxana/gmx_make_edi.c

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