src/gromacs/gmxana/gmx_make_edi.c

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