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