2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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.
48 #include "gromacs/commandline/pargs.h"
50 #include "gromacs/utility/smalloc.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/math/vec.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/fileio/pdbio.h"
57 #include "gromacs/fileio/confio.h"
58 #include "gromacs/fileio/tpxio.h"
59 #include "gromacs/fileio/matio.h"
61 #include "gromacs/fileio/xvgr.h"
66 #include "gromacs/utility/cstringutil.h"
72 gmx_bool bConstForce; /* Do constant force flooding instead of
73 evaluating a flooding potential */
82 /* This type is for the average, reference, target, and origin structure */
85 int nr; /* number of atoms this structure contains */
86 int *anrs; /* atom index numbers */
87 rvec *x; /* positions */
88 real *sqrtm; /* sqrt of the masses used for mass-
89 * weighting of analysis */
95 int nini; /* total Nr of atoms */
96 gmx_bool fitmas; /* true if trans fit with cm */
97 gmx_bool pcamas; /* true if mass-weighted PCA */
98 int presteps; /* number of steps to run without any
99 * perturbations ... just monitoring */
100 int outfrq; /* freq (in steps) of writing to edo */
101 int maxedsteps; /* max nr of steps per cycle */
102 struct edix sref; /* reference positions, to these fitting
104 struct edix sav; /* average positions */
105 struct edix star; /* target positions */
106 struct edix sori; /* origin positions */
107 real slope; /* minimal slope in acceptance radexp */
108 int ned; /* Nr of atoms in essdyn buffer */
109 t_edflood flood; /* parameters especially for flooding */
114 void make_t_edx(struct edix *edx, int natoms, rvec *pos, atom_id index[])
121 void write_t_edx(FILE *fp, struct edix edx, const char *comment)
123 /*here we copy only the pointers into the t_edx struct
124 no data is copied and edx.box is ignored */
126 fprintf(fp, "#%s \n %d \n", comment, edx.nr);
127 for (i = 0; i < edx.nr; i++)
129 fprintf(fp, "%d %f %f %f\n", (edx.anrs)[i]+1, (edx.x)[i][XX], (edx.x)[i][YY], (edx.x)[i][ZZ]);
133 int sscan_list(int *list[], const char *str, const char *listname)
135 /*this routine scans a string of the form 1,3-6,9 and returns the
136 selected numbers (in this case 1 3 4 5 6 9) in NULL-terminated array of integers.
137 memory for this list will be allocated in this routine -- sscan_list expects *list to
140 listname is a string used in the errormessage*/
145 char *pos, *startpos, *step;
148 /*enums to define the different lexical stati */
150 sBefore, sNumber, sMinus, sRange, sZero, sSmaller, sError, sSteppedRange
153 int status = sBefore; /*status of the deterministic automat to scan str */
157 char *start = NULL; /*holds the string of the number behind a ','*/
158 char *end = NULL; /*holds the string of the number behind a '-' */
160 int nvecs = 0; /* counts the number of vectors in the list*/
172 while ((c = *pos) != 0)
176 /* expect a number */
177 case sBefore: if (isdigit(c))
188 /* have read a number, expect ',' or '-' */
189 case sNumber: if (c == ',')
192 srenew(*list, nvecs+1);
193 (*list)[nvecs++] = number = strtol(start, NULL, 10);
203 status = sMinus; break;
214 /* have read a '-' -> expect a number */
219 status = sRange; break;
231 status = sError; break;
245 /* have read the number after a minus, expect ',' or ':' */
250 end_number = strtol(end, NULL, 10);
251 number = strtol(start, NULL, 10);
255 status = sZero; break;
257 if (end_number <= number)
259 status = sSmaller; break;
261 srenew(*list, nvecs+end_number-number+1);
264 istep = strtol(step, NULL, 10);
271 for (i = number; i <= end_number; i += istep)
273 (*list)[nvecs++] = i;
279 status = sSteppedRange;
291 /* format error occured */
293 gmx_fatal(FARGS, "Error in the list of eigenvectors for %s at pos %d with char %c", listname, pos-startpos, *(pos-1));
295 /* logical error occured */
297 gmx_fatal(FARGS, "Error in the list of eigenvectors for %s at pos %d: eigenvector 0 is not valid", listname, pos-startpos);
300 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);
303 ++pos; /* read next character */
304 } /*scanner has finished */
306 /* append zero to list of eigenvectors */
307 srenew(*list, nvecs+1);
313 void write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs, int nvec, const char *grouptitle, real steps[])
315 /* eig_list is a zero-terminated list of indices into the eigvecs array.
316 eigvecs are coordinates of eigenvectors
317 grouptitle to write in the comment line
318 steps -- array with stepsizes for evLINFIX, evLINACC and evRADACC
321 int n = 0, i; rvec x;
323 while (eig_list[n++])
325 ; /*count selected eigenvecs*/
328 fprintf(fp, "# NUMBER OF EIGENVECTORS + %s\n %d\n", grouptitle, n-1);
330 /* write list of eigenvector indicess */
331 for (n = 0; eig_list[n]; n++)
335 fprintf(fp, "%8d %g\n", eig_list[n], steps[n]);
339 fprintf(fp, "%8d %g\n", eig_list[n], 1.0);
344 /* dump coordinates of the selected eigenvectors */
348 for (i = 0; i < natoms; i++)
350 if (eig_list[n] > nvec)
352 gmx_fatal(FARGS, "Selected eigenvector %d is higher than maximum number %d of available eigenvectors", eig_list[n], nvec);
354 copy_rvec(eigvecs[eig_list[n]-1][i], x);
356 fprintf(fp, "%8.5f %8.5f %8.5f\n", x[XX], x[YY], x[ZZ]);
363 /*enum referring to the different lists of eigenvectors*/
365 evLINFIX, evLINACC, evFLOOD, evRADFIX, evRADACC, evRADCON, evMON, evNr
371 void write_the_whole_thing(FILE* fp, t_edipar *edpars, rvec** eigvecs,
372 int nvec, int *eig_listen[], real* evStepList[])
377 fprintf(fp, "#MAGIC\n %d \n#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
378 MAGIC, edpars->nini, edpars->fitmas, edpars->pcamas);
379 fprintf(fp, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
380 edpars->outfrq, edpars->maxedsteps, edpars->slope);
381 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",
382 edpars->presteps, edpars->flood.deltaF0, edpars->flood.deltaF, edpars->flood.tau, edpars->flood.constEfl,
383 edpars->flood.alpha2, edpars->flood.kT, edpars->flood.bHarmonic, edpars->flood.bConstForce);
385 /* Average and reference positions */
386 write_t_edx(fp, edpars->sref, "NREF, XREF");
387 write_t_edx(fp, edpars->sav, "NAV, XAV");
391 write_eigvec(fp, edpars->ned, eig_listen[evMON], eigvecs, nvec, "COMPONENTS GROUP 1", NULL);
392 write_eigvec(fp, edpars->ned, eig_listen[evLINFIX], eigvecs, nvec, "COMPONENTS GROUP 2", evStepList[evLINFIX]);
393 write_eigvec(fp, edpars->ned, eig_listen[evLINACC], eigvecs, nvec, "COMPONENTS GROUP 3", evStepList[evLINACC]);
394 write_eigvec(fp, edpars->ned, eig_listen[evRADFIX], eigvecs, nvec, "COMPONENTS GROUP 4", evStepList[evRADFIX]);
395 write_eigvec(fp, edpars->ned, eig_listen[evRADACC], eigvecs, nvec, "COMPONENTS GROUP 5", NULL);
396 write_eigvec(fp, edpars->ned, eig_listen[evRADCON], eigvecs, nvec, "COMPONENTS GROUP 6", NULL);
397 write_eigvec(fp, edpars->ned, eig_listen[evFLOOD], eigvecs, nvec, "COMPONENTS GROUP 7", evStepList[evFLOOD]);
400 /*Target and Origin positions */
401 write_t_edx(fp, edpars->star, "NTARGET, XTARGET");
402 write_t_edx(fp, edpars->sori, "NORIGIN, XORIGIN");
405 int read_conffile(const char *confin, char *title, rvec *x[])
407 /* read coordinates out of STX file */
411 printf("read coordnumber from file %s\n", confin);
412 get_stx_coordnum(confin, &natoms);
413 printf("number of coordinates in file %d\n", natoms);
414 /* if (natoms != ncoords)
415 gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n"
416 " does not match topology (= %d)",
417 confin,natoms,ncoords);
419 /* make space for coordinates and velocities */
420 init_t_atoms(&confat, natoms, FALSE);
422 read_stx_conf(confin, title, &confat, *x, NULL, NULL, box);
427 void read_eigenvalues(int vecs[], const char *eigfile, real values[],
428 gmx_bool bHesse, real kT)
433 neig = read_xvg(eigfile, &eigval, &nrow);
435 fprintf(stderr, "Read %d eigenvalues\n", neig);
436 for (i = bHesse ? 6 : 0; i < neig; i++)
438 if (eigval[1][i] < -0.001 && bHesse)
441 "WARNING: The Hessian Matrix has negative eigenvalue %f, we set it to zero (no flooding in this direction)\n\n", eigval[1][i]);
444 if (eigval[1][i] < 0)
451 for (i = 0; vecs[i]; i++)
455 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 values[i] = eigval[1][vecs[i]-1]/kT;
462 for (i = 0; vecs[i]; i++)
464 if (vecs[i] > (neig-6))
466 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 values[i] = 1/eigval[1][vecs[i]-1];
472 for (i = 0; i < nrow; i++)
480 static real *scan_vecparams(const char *str, const char * par, int nvecs)
482 char f0[256], f1[256]; /*format strings adapted every pass of the loop*/
487 snew(vec_params, nvecs);
491 for (i = 0; (i < nvecs); i++)
493 strcpy(f1, f0); /*f0 is the format string for the "to-be-ignored" numbers*/
494 strcat(f1, "%lf"); /*and f1 to read the actual number in this pass of the loop*/
495 if (sscanf(str, f1, &d) != 1)
497 gmx_fatal(FARGS, "Not enough elements for %s parameter (I need %d)", par, nvecs);
508 void init_edx(struct edix *edx)
515 void filter2edx(struct edix *edx, int nindex, atom_id index[], int ngro,
516 atom_id igro[], rvec *x, const char* structure)
518 /* filter2edx copies coordinates from x to edx which are given in index
524 srenew(edx->x, edx->nr);
525 srenew(edx->anrs, edx->nr);
526 for (i = 0; i < nindex; i++, ix++)
528 for (pos = 0; pos < ngro-1 && igro[pos] != index[i]; ++pos)
531 ; /*search element in igro*/
532 if (igro[pos] != index[i])
534 gmx_fatal(FARGS, "Couldn't find atom with index %d in structure %s", index[i], structure);
536 edx->anrs[ix] = index[i];
537 copy_rvec(x[pos], edx->x[ix]);
541 void get_structure(t_atoms *atoms, const char *IndexFile,
542 const char *StructureFile, struct edix *edx, int nfit,
543 atom_id ifit[], int nav, atom_id index[])
545 atom_id *igro; /*index corresponding to target or origin structure*/
553 ntar = read_conffile(StructureFile, title, &xtar);
554 printf("Select an index group of %d elements that corresponds to the atoms in the structure file %s\n",
555 ntar, StructureFile);
556 get_index(atoms, IndexFile, 1, &ngro, &igro, &grpname);
559 gmx_fatal(FARGS, "You selected an index group with %d elements instead of %d", ngro, ntar);
562 filter2edx(edx, nfit, ifit, ngro, igro, xtar, StructureFile);
564 /* If average and reference/fitting structure differ, append the average structure as well */
565 if (ifit != index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/
567 filter2edx(edx, nav, index, ngro, igro, xtar, StructureFile);
571 int gmx_make_edi(int argc, char *argv[])
574 static const char *desc[] = {
575 "[THISMODULE] generates an essential dynamics (ED) sampling input file to be used with [TT]mdrun[tt]",
576 "based on eigenvectors of a covariance matrix ([gmx-covar]) or from a",
577 "normal modes analysis ([gmx-nmeig]).",
578 "ED sampling can be used to manipulate the position along collective coordinates",
579 "(eigenvectors) of (biological) macromolecules during a simulation. Particularly,",
580 "it may be used to enhance the sampling efficiency of MD simulations by stimulating",
581 "the system to explore new regions along these collective coordinates. A number",
582 "of different algorithms are implemented to drive the system along the eigenvectors",
583 "([TT]-linfix[tt], [TT]-linacc[tt], [TT]-radfix[tt], [TT]-radacc[tt], [TT]-radcon[tt]),",
584 "to keep the position along a certain (set of) coordinate(s) fixed ([TT]-linfix[tt]),",
585 "or to only monitor the projections of the positions onto",
586 "these coordinates ([TT]-mon[tt]).[PAR]",
588 "A. Amadei, A.B.M. Linssen, B.L. de Groot, D.M.F. van Aalten and ",
589 "H.J.C. Berendsen; An efficient method for sampling the essential subspace ",
590 "of proteins., J. Biomol. Struct. Dyn. 13:615-626 (1996)[BR]",
591 "B.L. de Groot, A. Amadei, D.M.F. van Aalten and H.J.C. Berendsen; ",
592 "Towards an exhaustive sampling of the configurational spaces of the ",
593 "two forms of the peptide hormone guanylin,",
594 "J. Biomol. Struct. Dyn. 13 : 741-751 (1996)[BR]",
595 "B.L. de Groot, A.Amadei, R.M. Scheek, N.A.J. van Nuland and H.J.C. Berendsen; ",
596 "An extended sampling of the configurational space of HPr from E. coli",
597 "Proteins: Struct. Funct. Gen. 26: 314-322 (1996)",
598 "[PAR]You will be prompted for one or more index groups that correspond to the eigenvectors,",
599 "reference structure, target positions, etc.[PAR]",
601 "[TT]-mon[tt]: monitor projections of the coordinates onto selected eigenvectors.[PAR]",
602 "[TT]-linfix[tt]: perform fixed-step linear expansion along selected eigenvectors.[PAR]",
603 "[TT]-linacc[tt]: perform acceptance linear expansion along selected eigenvectors.",
604 "(steps in the desired directions will be accepted, others will be rejected).[PAR]",
605 "[TT]-radfix[tt]: perform fixed-step radius expansion along selected eigenvectors.[PAR]",
606 "[TT]-radacc[tt]: perform acceptance radius expansion along selected eigenvectors.",
607 "(steps in the desired direction will be accepted, others will be rejected).",
608 "[BB]Note:[bb] by default the starting MD structure will be taken as origin of the first",
609 "expansion cycle for radius expansion. If [TT]-ori[tt] is specified, you will be able",
610 "to read in a structure file that defines an external origin.[PAR]",
611 "[TT]-radcon[tt]: perform acceptance radius contraction along selected eigenvectors",
612 "towards a target structure specified with [TT]-tar[tt].[PAR]",
613 "NOTE: each eigenvector can be selected only once. [PAR]",
614 "[TT]-outfrq[tt]: frequency (in steps) of writing out projections etc. to [TT].xvg[tt] file[PAR]",
615 "[TT]-slope[tt]: minimal slope in acceptance radius expansion. A new expansion",
616 "cycle will be started if the spontaneous increase of the radius (in nm/step)",
617 "is less than the value specified.[PAR]",
618 "[TT]-maxedsteps[tt]: maximum number of steps per cycle in radius expansion",
619 "before a new cycle is started.[PAR]",
620 "Note on the parallel implementation: since ED sampling is a 'global' thing",
621 "(collective coordinates etc.), at least on the 'protein' side, ED sampling",
622 "is not very parallel-friendly from an implementation point of view. Because",
623 "parallel ED requires some extra communication, expect the performance to be",
624 "lower as in a free MD simulation, especially on a large number of nodes and/or",
625 "when the ED group contains a lot of atoms. [PAR]",
626 "Please also note that if your ED group contains more than a single protein,",
627 "then the [TT].tpr[tt] file must contain the correct PBC representation of the ED group.",
628 "Take a look on the initial RMSD from the reference structure, which is printed",
629 "out at the start of the simulation; if this is much higher than expected, one",
630 "of the ED molecules might be shifted by a box vector. [PAR]",
631 "All ED-related output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a [TT].xvg[tt] file",
632 "as a function of time in intervals of OUTFRQ steps.[PAR]",
633 "[BB]Note[bb] that you can impose multiple ED constraints and flooding potentials in",
634 "a single simulation (on different molecules) if several [TT].edi[tt] files were concatenated",
635 "first. The constraints are applied in the order they appear in the [TT].edi[tt] file. ",
636 "Depending on what was specified in the [TT].edi[tt] input file, the output file contains for each ED dataset[PAR]",
637 "[TT]*[tt] the RMSD of the fitted molecule to the reference structure (for atoms involved in fitting prior to calculating the ED constraints)[BR]",
638 "[TT]*[tt] projections of the positions onto selected eigenvectors[BR]",
641 "with [TT]-flood[tt], you can specify which eigenvectors are used to compute a flooding potential,",
642 "which will lead to extra forces expelling the structure out of the region described",
643 "by the covariance matrix. If you switch -restrain the potential is inverted and the structure",
644 "is kept in that region.",
646 "The origin is normally the average structure stored in the [TT]eigvec.trr[tt] file.",
647 "It can be changed with [TT]-ori[tt] to an arbitrary position in configuration space.",
648 "With [TT]-tau[tt], [TT]-deltaF0[tt], and [TT]-Eflnull[tt] you control the flooding behaviour.",
649 "Efl is the flooding strength, it is updated according to the rule of adaptive flooding.",
650 "Tau is the time constant of adaptive flooding, high [GRK]tau[grk] means slow adaption (i.e. growth). ",
651 "DeltaF0 is the flooding strength you want to reach after tau ps of simulation.",
652 "To use constant Efl set [TT]-tau[tt] to zero.",
654 "[TT]-alpha[tt] is a fudge parameter to control the width of the flooding potential. A value of 2 has been found",
655 "to give good results for most standard cases in flooding of proteins.",
656 "[GRK]alpha[grk] basically accounts for incomplete sampling, if you sampled further the width of the ensemble would",
657 "increase, this is mimicked by [GRK]alpha[grk] > 1.",
658 "For restraining, [GRK]alpha[grk] < 1 can give you smaller width in the restraining potential.",
660 "RESTART and FLOODING:",
661 "If you want to restart a crashed flooding simulation please find the values deltaF and Efl in",
662 "the output file and manually put them into the [TT].edi[tt] file under DELTA_F0 and EFL_NULL."
665 /* Save all the params in this struct and then save it in an edi file.
666 * ignoring fields nmass,massnrs,mass,tmass,nfit,fitnrs,edo
668 static t_edipar edi_params;
671 evStepNr = evRADFIX + 1
673 static const char* evSelections[evNr] = {NULL, NULL, NULL, NULL, NULL, NULL};
674 static const char* evOptions[evNr] = {"-linfix", "-linacc", "-flood", "-radfix", "-radacc", "-radcon", "-mon"};
675 static const char* evParams[evStepNr] = {NULL, NULL};
676 static const char* evStepOptions[evStepNr] = {"-linstep", "-accdir", "-not_used", "-radstep"};
677 static const char* ConstForceStr;
678 static real * evStepList[evStepNr];
679 static real radstep = 0.0;
680 static real deltaF0 = 150;
681 static real deltaF = 0;
682 static real tau = .1;
683 static real constEfl = 0.0;
684 static real alpha = 1;
685 static int eqSteps = 0;
686 static int * listen[evNr];
687 static real T = 300.0;
688 const real kB = 2.5 / 300.0; /* k_boltzmann in MD units */
689 static gmx_bool bRestrain = FALSE;
690 static gmx_bool bHesse = FALSE;
691 static gmx_bool bHarmonic = FALSE;
693 { "-mon", FALSE, etSTR, {&evSelections[evMON]},
694 "Indices of eigenvectors for projections of x (e.g. 1,2-5,9) or 1-100:10 means 1 11 21 31 ... 91" },
695 { "-linfix", FALSE, etSTR, {&evSelections[0]},
696 "Indices of eigenvectors for fixed increment linear sampling" },
697 { "-linacc", FALSE, etSTR, {&evSelections[1]},
698 "Indices of eigenvectors for acceptance linear sampling" },
699 { "-radfix", FALSE, etSTR, {&evSelections[3]},
700 "Indices of eigenvectors for fixed increment radius expansion" },
701 { "-radacc", FALSE, etSTR, {&evSelections[4]},
702 "Indices of eigenvectors for acceptance radius expansion" },
703 { "-radcon", FALSE, etSTR, {&evSelections[5]},
704 "Indices of eigenvectors for acceptance radius contraction" },
705 { "-flood", FALSE, etSTR, {&evSelections[2]},
706 "Indices of eigenvectors for flooding"},
707 { "-outfrq", FALSE, etINT, {&edi_params.outfrq},
708 "Freqency (in steps) of writing output in [TT].xvg[tt] file" },
709 { "-slope", FALSE, etREAL, { &edi_params.slope},
710 "Minimal slope in acceptance radius expansion"},
711 { "-linstep", FALSE, etSTR, {&evParams[0]},
712 "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 -3.1\")"},
713 { "-accdir", FALSE, etSTR, {&evParams[1]},
714 "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 -1.1\")"},
715 { "-radstep", FALSE, etREAL, {&radstep},
716 "Stepsize (nm/step) for fixed increment radius expansion"},
717 { "-maxedsteps", FALSE, etINT, {&edi_params.maxedsteps},
718 "Maximum number of steps per cycle" },
719 { "-eqsteps", FALSE, etINT, {&eqSteps},
720 "Number of steps to run without any perturbations "},
721 { "-deltaF0", FALSE, etREAL, {&deltaF0},
722 "Target destabilization energy for flooding"},
723 { "-deltaF", FALSE, etREAL, {&deltaF},
724 "Start deltaF with this parameter - default 0, nonzero values only needed for restart"},
725 { "-tau", FALSE, etREAL, {&tau},
726 "Coupling constant for adaption of flooding strength according to deltaF0, 0 = infinity i.e. constant flooding strength"},
727 { "-Eflnull", FALSE, etREAL, {&constEfl},
728 "The starting value of the flooding strength. The flooding strength is updated "
729 "according to the adaptive flooding scheme. For a constant flooding strength use [TT]-tau[tt] 0. "},
730 { "-T", FALSE, etREAL, {&T},
731 "T is temperature, the value is needed if you want to do flooding "},
732 { "-alpha", FALSE, etREAL, {&alpha},
733 "Scale width of gaussian flooding potential with alpha^2 "},
734 { "-restrain", FALSE, etBOOL, {&bRestrain},
735 "Use the flooding potential with inverted sign -> effects as quasiharmonic restraining potential"},
736 { "-hessian", FALSE, etBOOL, {&bHesse},
737 "The eigenvectors and eigenvalues are from a Hessian matrix"},
738 { "-harmonic", FALSE, etBOOL, {&bHarmonic},
739 "The eigenvalues are interpreted as spring constant"},
740 { "-constF", FALSE, etSTR, {&ConstForceStr},
741 "Constant force flooding: manually set the forces for the eigenvectors selected with -flood "
742 "(put in quotes! \"1.0 2.3 5.1 -3.1\"). No other flooding parameters are needed when specifying the forces directly."}
744 #define NPA asize(pa)
747 int nvec1, *eignr1 = NULL;
748 rvec *xav1, **eigvec1 = NULL;
749 t_atoms *atoms = NULL;
750 int nav; /* Number of atoms in the average structure */
752 const char *indexfile;
754 atom_id *index, *ifit;
755 int nfit; /* Number of atoms in the reference/fit structure */
756 int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */
758 real *eigval1 = NULL; /* in V3.3 this is parameter of read_eigenvectors */
761 const char *TargetFile;
762 const char *OriginFile;
763 const char *EigvecFile;
767 /*to read topology file*/
773 gmx_bool bTop, bFit1;
776 { efTRN, "-f", "eigenvec", ffREAD },
777 { efXVG, "-eig", "eigenval", ffOPTRD },
778 { efTPS, NULL, NULL, ffREAD },
779 { efNDX, NULL, NULL, ffOPTRD },
780 { efSTX, "-tar", "target", ffOPTRD},
781 { efSTX, "-ori", "origin", ffOPTRD},
782 { efEDI, "-o", "sam", ffWRITE }
784 #define NFILE asize(fnm)
785 edi_params.outfrq = 100; edi_params.slope = 0.0; edi_params.maxedsteps = 0;
786 if (!parse_common_args(&argc, argv, 0,
787 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
792 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
793 EdiFile = ftp2fn(efEDI, NFILE, fnm);
794 TargetFile = opt2fn_null("-tar", NFILE, fnm);
795 OriginFile = opt2fn_null("-ori", NFILE, fnm);
798 for (ev_class = 0; ev_class < evNr; ++ev_class)
800 if (opt2parg_bSet(evOptions[ev_class], NPA, pa))
802 /*get list of eigenvectors*/
803 nvecs = sscan_list(&(listen[ev_class]), opt2parg_str(evOptions[ev_class], NPA, pa), evOptions[ev_class]);
804 if (ev_class < evStepNr-2)
806 /*if apropriate get list of stepsizes for these eigenvectors*/
807 if (opt2parg_bSet(evStepOptions[ev_class], NPA, pa))
809 evStepList[ev_class] =
810 scan_vecparams(opt2parg_str(evStepOptions[ev_class], NPA, pa), evStepOptions[ev_class], nvecs);
812 else /*if list is not given fill with zeros */
814 snew(evStepList[ev_class], nvecs);
815 for (i = 0; i < nvecs; i++)
817 evStepList[ev_class][i] = 0.0;
821 else if (ev_class == evRADFIX)
823 snew(evStepList[ev_class], nvecs);
824 for (i = 0; i < nvecs; i++)
826 evStepList[ev_class][i] = radstep;
829 else if (ev_class == evFLOOD)
831 snew(evStepList[ev_class], nvecs);
833 /* Are we doing constant force flooding? In that case, we read in
834 * the fproj values from the command line */
835 if (opt2parg_bSet("-constF", NPA, pa))
837 evStepList[ev_class] = scan_vecparams(opt2parg_str("-constF", NPA, pa), "-constF", nvecs);
842 }; /*to avoid ambiguity */
844 else /* if there are no eigenvectors for this option set list to zero */
846 listen[ev_class] = NULL;
847 snew(listen[ev_class], 1);
848 listen[ev_class][0] = 0;
852 /* print the interpreted list of eigenvectors - to give some feedback*/
853 for (ev_class = 0; ev_class < evNr; ++ev_class)
855 printf("Eigenvector list %7s consists of the indices: ", evOptions[ev_class]);
857 while (listen[ev_class][i])
859 printf("%d ", listen[ev_class][i++]);
864 EigvecFile = opt2fn("-f", NFILE, fnm);
866 /*read eigenvectors from eigvec.trr*/
867 read_eigenvectors(EigvecFile, &nav, &bFit1,
868 &xref1, &edi_params.fitmas, &xav1, &edi_params.pcamas, &nvec1, &eignr1, &eigvec1, &eigval1);
870 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
871 title, &top, &ePBC, &xtop, NULL, topbox, 0);
875 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", nav);
876 get_index(atoms, indexfile, 1, &i, &index, &grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
879 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d",
889 /* if g_covar used different coordinate groups to fit and to do the PCA */
890 printf("\nNote: the structure in %s should be the same\n"
891 " as the one used for the fit in g_covar\n", ftp2fn(efTPS, NFILE, fnm));
892 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
896 printf("\nNote: Apparently no fitting was done in g_covar.\n"
897 " However, you need to select a reference group for fitting in mdrun\n");
899 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
901 for (i = 0; i < nfit; i++)
903 copy_rvec(xtop[ifit[i]], xref1[i]);
912 if (opt2parg_bSet("-constF", NPA, pa))
914 /* Constant force flooding is special: Most of the normal flooding
915 * options are not needed. */
916 edi_params.flood.bConstForce = TRUE;
920 /* For normal flooding read eigenvalues and store them in evSteplist[evFLOOD] */
922 if (listen[evFLOOD][0] != 0)
924 read_eigenvalues(listen[evFLOOD], opt2fn("-eig", NFILE, fnm), evStepList[evFLOOD], bHesse, kB*T);
927 edi_params.flood.tau = tau;
928 edi_params.flood.deltaF0 = deltaF0;
929 edi_params.flood.deltaF = deltaF;
930 edi_params.presteps = eqSteps;
931 edi_params.flood.kT = kB*T;
932 edi_params.flood.bHarmonic = bHarmonic;
935 /* Trick: invert sign of Efl and alpha2 then this will give the same sign in the exponential and inverted sign outside */
936 edi_params.flood.constEfl = -constEfl;
937 edi_params.flood.alpha2 = -sqr(alpha);
941 edi_params.flood.constEfl = constEfl;
942 edi_params.flood.alpha2 = sqr(alpha);
946 edi_params.ned = nav;
948 /*number of system atoms */
949 edi_params.nini = atoms->nr;
952 /*store reference and average structure in edi_params*/
953 make_t_edx(&edi_params.sref, nfit, xref1, ifit );
954 make_t_edx(&edi_params.sav, nav, xav1, index);
957 /* Store target positions in edi_params */
958 if (opt2bSet("-tar", NFILE, fnm))
960 if (0 != listen[evFLOOD][0])
962 fprintf(stderr, "\nNote: Providing a TARGET structure has no effect when using flooding.\n"
963 " You may want to use -ori to define the flooding potential center.\n\n");
965 get_structure(atoms, indexfile, TargetFile, &edi_params.star, nfit, ifit, nav, index);
969 make_t_edx(&edi_params.star, 0, NULL, index);
972 /* Store origin positions */
973 if (opt2bSet("-ori", NFILE, fnm))
975 get_structure(atoms, indexfile, OriginFile, &edi_params.sori, nfit, ifit, nav, index);
979 make_t_edx(&edi_params.sori, 0, NULL, index);
983 write_the_whole_thing(gmx_ffopen(EdiFile, "w"), &edi_params, eigvec1, nvec1, listen, evStepList);