2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
5 * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
36 /* The make_edi program was generously contributed by Oliver Lange, based
37 * on the code from g_anaeig. You can reach him as olange@gwdg.de. He
38 * probably also holds copyright to the following code.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/fileio/pdbio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxana/eigio.h"
52 #include "gromacs/gmxana/gmx_ana.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/topology/index.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/utility/arraysize.h"
58 #include "gromacs/utility/cstringutil.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/futil.h"
61 #include "gromacs/utility/smalloc.h"
67 gmx_bool bConstForce; /* Do constant force flooding instead of
68 evaluating a flooding potential */
77 /* This type is for the average, reference, target, and origin structure */
80 int nr; /* number of atoms this structure contains */
81 int* anrs; /* atom index numbers */
82 rvec* x; /* positions */
83 real* sqrtm; /* sqrt of the masses used for mass-
84 * weighting of analysis */
90 int nini; /* total Nr of atoms */
91 gmx_bool fitmas; /* true if trans fit with cm */
92 gmx_bool pcamas; /* true if mass-weighted PCA */
93 int presteps; /* number of steps to run without any
94 * perturbations ... just monitoring */
95 int outfrq; /* freq (in steps) of writing to edo */
96 int maxedsteps; /* max nr of steps per cycle */
97 struct edix sref; /* reference positions, to these fitting
99 struct edix sav; /* average positions */
100 struct edix star; /* target positions */
101 struct edix sori; /* origin positions */
102 real slope; /* minimal slope in acceptance radexp */
103 int ned; /* Nr of atoms in essdyn buffer */
104 t_edflood flood; /* parameters especially for flooding */
108 static void make_t_edx(struct edix* edx, int natoms, rvec* pos, int index[])
115 static void write_t_edx(FILE* fp, struct edix edx, const char* comment)
117 /*here we copy only the pointers into the t_edx struct
118 no data is copied and edx.box is ignored */
120 fprintf(fp, "#%s \n %d \n", comment, edx.nr);
121 for (i = 0; i < edx.nr; i++)
123 fprintf(fp, "%d %f %f %f\n", (edx.anrs)[i] + 1, (edx.x)[i][XX], (edx.x)[i][YY], (edx.x)[i][ZZ]);
127 static int sscan_list(int* list[], const char* str, const char* listname)
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
134 listname is a string used in the errormessage*/
139 char *pos, *startpos, *step;
140 int n = std::strlen(str);
142 /*enums to define the different lexical stati */
155 int status = sBefore; /*status of the deterministic automat to scan str */
159 char* start = nullptr; /*holds the string of the number behind a ','*/
160 char* end = nullptr; /*holds the string of the number behind a '-' */
162 int nvecs = 0; /* counts the number of vectors in the list*/
167 std::strcpy(pos, str);
174 while ((c = *pos) != 0)
178 /* expect a number */
192 /* have read a number, expect ',' or '-' */
197 srenew(*list, nvecs + 1);
198 (*list)[nvecs++] = number = std::strtol(start, nullptr, 10);
211 else if (std::isdigit(c))
221 /* have read a '-' -> expect a number */
256 /* have read the number after a minus, expect ',' or ':' */
261 end_number = std::strtol(end, nullptr, 10);
262 number = std::strtol(start, nullptr, 10);
269 if (end_number <= number)
274 srenew(*list, nvecs + end_number - number + 1);
277 istep = strtol(step, nullptr, 10);
284 for (i = number; i <= end_number; i += istep)
286 (*list)[nvecs++] = i;
292 status = sSteppedRange;
295 else if (std::isdigit(c))
305 /* format error occured */
307 gmx_fatal(FARGS, "Error in the list of eigenvectors for %s at pos %td with char %c",
308 listname, pos - startpos, *(pos - 1));
309 /* logical error occured */
312 "Error in the list of eigenvectors for %s at pos %td: eigenvector 0 is "
314 listname, pos - startpos);
317 "Error in the list of eigenvectors for %s at pos %td: second index %d is "
318 "not bigger than %d",
319 listname, pos - startpos, end_number, number);
321 ++pos; /* read next character */
322 } /*scanner has finished */
324 /* append zero to list of eigenvectors */
325 srenew(*list, nvecs + 1);
332 write_eigvec(FILE* fp, int natoms, int eig_list[], rvec** eigvecs, int nvec, const char* grouptitle, real steps[])
334 /* eig_list is a zero-terminated list of indices into the eigvecs array.
335 eigvecs are coordinates of eigenvectors
336 grouptitle to write in the comment line
337 steps -- array with stepsizes for evLINFIX, evLINACC and evRADACC
342 while (eig_list[n++])
344 /*count selected eigenvecs*/
346 fprintf(fp, "# NUMBER OF EIGENVECTORS + %s\n %d\n", grouptitle, n - 1);
348 /* write list of eigenvector indicess */
349 for (n = 0; eig_list[n]; n++)
353 fprintf(fp, "%8d %g\n", eig_list[n], steps[n]);
357 fprintf(fp, "%8d %g\n", eig_list[n], 1.0);
362 /* dump coordinates of the selected eigenvectors */
365 for (i = 0; i < natoms; i++)
367 if (eig_list[n] > nvec)
370 "Selected eigenvector %d is higher than maximum number %d of available "
374 copy_rvec(eigvecs[eig_list[n] - 1][i], x);
375 fprintf(fp, "%8.5f %8.5f %8.5f\n", x[XX], x[YY], x[ZZ]);
382 /*enum referring to the different lists of eigenvectors*/
398 static void write_the_whole_thing(FILE* fp,
408 fprintf(fp, "#MAGIC\n %d \n#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n", MAGIC, edpars->nini,
409 int(edpars->fitmas), int(edpars->pcamas));
410 fprintf(fp, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n", edpars->outfrq, edpars->maxedsteps,
413 "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#INIT_DELTA_F\n %f\n#TAU\n %f\n#EFL_NULL\n "
414 "%f\n#ALPHA2\n %f\n#KT\n %f\n#HARMONIC\n %d\n#CONST_FORCE_FLOODING\n %d\n",
415 edpars->presteps, edpars->flood.deltaF0, edpars->flood.deltaF, edpars->flood.tau,
416 edpars->flood.constEfl, edpars->flood.alpha2, edpars->flood.kT,
417 int(edpars->flood.bHarmonic), int(edpars->flood.bConstForce));
419 /* Average and reference positions */
420 write_t_edx(fp, edpars->sref, "NREF, XREF");
421 write_t_edx(fp, edpars->sav, "NAV, XAV");
425 write_eigvec(fp, edpars->ned, eig_listen[evMON], eigvecs, nvec, "COMPONENTS GROUP 1", nullptr);
426 write_eigvec(fp, edpars->ned, eig_listen[evLINFIX], eigvecs, nvec, "COMPONENTS GROUP 2",
427 evStepList[evLINFIX]);
428 write_eigvec(fp, edpars->ned, eig_listen[evLINACC], eigvecs, nvec, "COMPONENTS GROUP 3",
429 evStepList[evLINACC]);
430 write_eigvec(fp, edpars->ned, eig_listen[evRADFIX], eigvecs, nvec, "COMPONENTS GROUP 4",
431 evStepList[evRADFIX]);
432 write_eigvec(fp, edpars->ned, eig_listen[evRADACC], eigvecs, nvec, "COMPONENTS GROUP 5", nullptr);
433 write_eigvec(fp, edpars->ned, eig_listen[evRADCON], eigvecs, nvec, "COMPONENTS GROUP 6", nullptr);
434 write_eigvec(fp, edpars->ned, eig_listen[evFLOOD], eigvecs, nvec, "COMPONENTS GROUP 7",
435 evStepList[evFLOOD]);
438 /*Target and Origin positions */
439 write_t_edx(fp, edpars->star, "NTARGET, XTARGET");
440 write_t_edx(fp, edpars->sori, "NORIGIN, XORIGIN");
443 static int read_conffile(const char* confin, rvec** x)
447 printf("read coordnumber from file %s\n", confin);
448 read_tps_conf(confin, &top, nullptr, x, nullptr, box, FALSE);
449 printf("number of coordinates in file %d\n", top.atoms.nr);
454 static void read_eigenvalues(const int vecs[], const char* eigfile, real values[], gmx_bool bHesse, real kT, int natoms_average_struct)
459 neig = read_xvg(eigfile, &eigval, &nrow);
461 fprintf(stderr, "Read %d eigenvalues\n", neig);
462 for (i = bHesse ? 6 : 0; i < neig; i++)
464 if (eigval[1][i] < -0.001 && bHesse)
467 "WARNING: The Hessian Matrix has negative eigenvalue %f, we set it to zero (no "
468 "flooding in this direction)\n\n",
472 if (eigval[1][i] < 0)
479 for (i = 0; vecs[i]; i++)
484 "ERROR: You have chosen one of the first 6 eigenvectors of the HESSE "
485 "Matrix. That does not make sense, since they correspond to the 6 "
486 "rotational and translational degrees of freedom.\n\n");
488 values[i] = eigval[1][vecs[i] - 1] / kT;
493 for (i = 0; vecs[i]; i++)
495 /* Make sure this eigenvalue does not correspond to one of the last 6 eigenvectors of the
496 * covariance matrix. These correspond to the rotational and translational degrees of
497 * freedom and will be zero within numerical accuracy.
499 * Note that the total number of eigenvectors produced by gmx covar depends on:
500 * 1) the total number of degrees of freedom of the system (3N, with N the number of atoms)
501 * 2) the number S of independent configurations fed into gmx covar.
502 * For long trajectories with lots of frames, usually S >= 3N + 1, so that one indeed gets
503 * 3N eigenvalues (of which the last 6 will have zero eigenvalues).
504 * For S < 3N + 1, however, the covariance matrix becomes rank deficient, and the number
505 * of possible eigenvalues is just S - 1. Since in make_edi we only know N but not S, we can
506 * only warn the user if he picked one of the last 6 of 3N.
508 if (vecs[i] > (3 * natoms_average_struct - 6))
511 "ERROR: You have chosen one of the last 6 eigenvectors of the COVARIANCE "
512 "Matrix. That does not make sense, since they correspond to the 6 "
513 "rotational and translational degrees of freedom.\n\n");
515 values[i] = 1 / eigval[1][vecs[i] - 1];
519 for (i = 0; i < nrow; i++)
527 static real* scan_vecparams(const char* str, const char* par, int nvecs)
529 char f0[256], f1[256]; /*format strings adapted every pass of the loop*/
534 snew(vec_params, nvecs);
538 for (i = 0; (i < nvecs); i++)
540 std::strcpy(f1, f0); /*f0 is the format string for the "to-be-ignored" numbers*/
541 std::strcat(f1, "%lf"); /*and f1 to read the actual number in this pass of the loop*/
542 if (sscanf(str, f1, &d) != 1)
544 gmx_fatal(FARGS, "Not enough elements for %s parameter (I need %d)", par, nvecs);
547 std::strcat(f0, "%*s");
554 static void init_edx(struct edix* edx)
562 filter2edx(struct edix* edx, int nindex, int index[], int ngro, const int igro[], const rvec* x, const char* structure)
564 /* filter2edx copies coordinates from x to edx which are given in index
570 srenew(edx->x, edx->nr);
571 srenew(edx->anrs, edx->nr);
572 for (i = 0; i < nindex; i++, ix++)
574 for (pos = 0; pos < ngro - 1 && igro[pos] != index[i]; ++pos) {} /*search element in igro*/
575 if (igro[pos] != index[i])
577 gmx_fatal(FARGS, "Couldn't find atom with index %d in structure %s", index[i], structure);
579 edx->anrs[ix] = index[i];
580 copy_rvec(x[pos], edx->x[ix]);
584 static void get_structure(const t_atoms* atoms,
585 const char* IndexFile,
586 const char* StructureFile,
593 int* igro; /*index corresponding to target or origin structure*/
600 ntar = read_conffile(StructureFile, &xtar);
601 printf("Select an index group of %d elements that corresponds to the atoms in the structure "
603 ntar, StructureFile);
604 get_index(atoms, IndexFile, 1, &ngro, &igro, &grpname);
607 gmx_fatal(FARGS, "You selected an index group with %d elements instead of %d", ngro, ntar);
610 filter2edx(edx, nfit, ifit, ngro, igro, xtar, StructureFile);
612 /* If average and reference/fitting structure differ, append the average structure as well */
613 if (ifit != index) /*if fit structure is different append these coordinates, too -- don't mind duplicates*/
615 filter2edx(edx, nav, index, ngro, igro, xtar, StructureFile);
619 int gmx_make_edi(int argc, char* argv[])
622 static const char* desc[] = {
623 "[THISMODULE] generates an essential dynamics (ED) sampling input file to be used with ",
624 "[TT]mdrun[tt] based on eigenvectors of a covariance matrix ([gmx-covar]) or from a",
625 "normal modes analysis ([gmx-nmeig]).",
626 "ED sampling can be used to manipulate the position along collective coordinates",
627 "(eigenvectors) of (biological) macromolecules during a simulation. Particularly,",
628 "it may be used to enhance the sampling efficiency of MD simulations by stimulating",
629 "the system to explore new regions along these collective coordinates. A number",
630 "of different algorithms are implemented to drive the system along the eigenvectors",
631 "([TT]-linfix[tt], [TT]-linacc[tt], [TT]-radfix[tt], [TT]-radacc[tt], [TT]-radcon[tt]),",
632 "to keep the position along a certain (set of) coordinate(s) fixed ([TT]-linfix[tt]),",
633 "or to only monitor the projections of the positions onto",
634 "these coordinates ([TT]-mon[tt]).[PAR]",
636 "A. Amadei, A.B.M. Linssen, B.L. de Groot, D.M.F. van Aalten and ",
637 "H.J.C. Berendsen; An efficient method for sampling the essential subspace ",
638 "of proteins., J. Biomol. Struct. Dyn. 13:615-626 (1996)[PAR]",
639 "B.L. de Groot, A. Amadei, D.M.F. van Aalten and H.J.C. Berendsen; ",
640 "Towards an exhaustive sampling of the configurational spaces of the ",
641 "two forms of the peptide hormone guanylin,",
642 "J. Biomol. Struct. Dyn. 13 : 741-751 (1996)[PAR]",
643 "B.L. de Groot, A.Amadei, R.M. Scheek, N.A.J. van Nuland and H.J.C. Berendsen; ",
644 "An extended sampling of the configurational space of HPr from E. coli",
645 "Proteins: Struct. Funct. Gen. 26: 314-322 (1996)",
646 "[PAR]You will be prompted for one or more index groups that correspond to the ",
648 "reference structure, target positions, etc.[PAR]",
650 "[TT]-mon[tt]: monitor projections of the coordinates onto selected eigenvectors.[PAR]",
651 "[TT]-linfix[tt]: perform fixed-step linear expansion along selected eigenvectors.[PAR]",
652 "[TT]-linacc[tt]: perform acceptance linear expansion along selected eigenvectors.",
653 "(steps in the desired directions will be accepted, others will be rejected).[PAR]",
654 "[TT]-radfix[tt]: perform fixed-step radius expansion along selected eigenvectors.[PAR]",
655 "[TT]-radacc[tt]: perform acceptance radius expansion along selected eigenvectors.",
656 "(steps in the desired direction will be accepted, others will be rejected).",
657 "[BB]Note:[bb] by default the starting MD structure will be taken as origin of the first",
658 "expansion cycle for radius expansion. If [TT]-ori[tt] is specified, you will be able",
659 "to read in a structure file that defines an external origin.[PAR]",
660 "[TT]-radcon[tt]: perform acceptance radius contraction along selected eigenvectors",
661 "towards a target structure specified with [TT]-tar[tt].[PAR]",
662 "NOTE: each eigenvector can be selected only once. [PAR]",
663 "[TT]-outfrq[tt]: frequency (in steps) of writing out projections etc. to [REF].xvg[ref] ",
666 "[TT]-slope[tt]: minimal slope in acceptance radius expansion. A new expansion",
667 "cycle will be started if the spontaneous increase of the radius (in nm/step)",
668 "is less than the value specified.[PAR]",
669 "[TT]-maxedsteps[tt]: maximum number of steps per cycle in radius expansion",
670 "before a new cycle is started.[PAR]",
671 "Note on the parallel implementation: since ED sampling is a 'global' thing",
672 "(collective coordinates etc.), at least on the 'protein' side, ED sampling",
673 "is not very parallel-friendly from an implementation point of view. Because",
674 "parallel ED requires some extra communication, expect the performance to be",
675 "lower as in a free MD simulation, especially on a large number of ranks and/or",
676 "when the ED group contains a lot of atoms. [PAR]",
677 "Please also note that if your ED group contains more than a single protein,",
678 "then the [REF].tpr[ref] file must contain the correct PBC representation of the ED group.",
679 "Take a look on the initial RMSD from the reference structure, which is printed",
680 "out at the start of the simulation; if this is much higher than expected, one",
681 "of the ED molecules might be shifted by a box vector. [PAR]",
682 "All ED-related output of [TT]mdrun[tt] (specify with [TT]-eo[tt]) is written to a ",
683 "[REF].xvg[ref] file as a function of time in intervals of OUTFRQ steps.[PAR]",
684 "[BB]Note[bb] that you can impose multiple ED constraints and flooding potentials in",
685 "a single simulation (on different molecules) if several [REF].edi[ref] files were ",
686 "concatenated first. The constraints are applied in the order they appear in ",
687 "the [REF].edi[ref] file. Depending on what was specified in the [REF].edi[ref] ",
688 "input file, the output file contains for each ED dataset",
690 " * the RMSD of the fitted molecule to the reference structure (for atoms involved in ",
691 " fitting prior to calculating the ED constraints)",
692 " * projections of the positions onto selected eigenvectors",
695 "with [TT]-flood[tt], you can specify which eigenvectors are used to compute a flooding ",
697 "which will lead to extra forces expelling the structure out of the region described",
698 "by the covariance matrix. If you switch -restrain the potential is inverted and the ",
699 "structure is kept in that region.",
701 "The origin is normally the average structure stored in the [TT]eigvec.trr[tt] file.",
702 "It can be changed with [TT]-ori[tt] to an arbitrary position in configuration space.",
703 "With [TT]-tau[tt], [TT]-deltaF0[tt], and [TT]-Eflnull[tt] you control the flooding ",
704 "behaviour. Efl is the flooding strength, it is updated according to the rule of ",
705 "adaptive flooding. Tau is the time constant of adaptive flooding, high ",
706 "[GRK]tau[grk] means slow adaption (i.e. growth). ",
707 "DeltaF0 is the flooding strength you want to reach after tau ps of simulation.",
708 "To use constant Efl set [TT]-tau[tt] to zero.",
710 "[TT]-alpha[tt] is a fudge parameter to control the width of the flooding potential. A ",
711 "value of 2 has been found",
712 "to give good results for most standard cases in flooding of proteins.",
713 "[GRK]alpha[grk] basically accounts for incomplete sampling, if you sampled further the ",
714 "width of the ensemble would",
715 "increase, this is mimicked by [GRK]alpha[grk] > 1.",
716 "For restraining, [GRK]alpha[grk] < 1 can give you smaller width in the restraining ",
719 "RESTART and FLOODING:",
720 "If you want to restart a crashed flooding simulation please find the values deltaF and ",
722 "the output file and manually put them into the [REF].edi[ref] file under DELTA_F0 and ",
726 /* Save all the params in this struct and then save it in an edi file.
727 * ignoring fields nmass,massnrs,mass,tmass,nfit,fitnrs,edo
729 static t_edipar edi_params;
733 evStepNr = evRADFIX + 1
735 static const char* evSelections[evNr] = { nullptr, nullptr, nullptr, nullptr, nullptr, nullptr };
736 static const char* evOptions[evNr] = { "-linfix", "-linacc", "-flood", "-radfix",
737 "-radacc", "-radcon", "-mon" };
738 static const char* evParams[evStepNr] = { nullptr, nullptr };
739 static const char* evStepOptions[evStepNr] = { "-linstep", "-accdir", "-not_used", "-radstep" };
740 static const char* ConstForceStr;
741 static real* evStepList[evStepNr];
742 static real radstep = 0.0;
743 static real deltaF0 = 150;
744 static real deltaF = 0;
745 static real tau = .1;
746 static real constEfl = 0.0;
747 static real alpha = 1;
748 static int eqSteps = 0;
749 static int* listen[evNr];
750 static real T = 300.0;
751 const real kB = 2.5 / 300.0; /* k_boltzmann in MD units */
752 static gmx_bool bRestrain = FALSE;
753 static gmx_bool bHesse = FALSE;
754 static gmx_bool bHarmonic = FALSE;
759 { &evSelections[evMON] },
760 "Indices of eigenvectors for projections of x (e.g. 1,2-5,9) or 1-100:10 means 1 11 21 "
765 { &evSelections[0] },
766 "Indices of eigenvectors for fixed increment linear sampling" },
770 { &evSelections[1] },
771 "Indices of eigenvectors for acceptance linear sampling" },
775 { &evSelections[3] },
776 "Indices of eigenvectors for fixed increment radius expansion" },
780 { &evSelections[4] },
781 "Indices of eigenvectors for acceptance radius expansion" },
785 { &evSelections[5] },
786 "Indices of eigenvectors for acceptance radius contraction" },
787 { "-flood", FALSE, etSTR, { &evSelections[2] }, "Indices of eigenvectors for flooding" },
791 { &edi_params.outfrq },
792 "Frequency (in steps) of writing output in [REF].xvg[ref] file" },
796 { &edi_params.slope },
797 "Minimal slope in acceptance radius expansion" },
802 "Stepsizes (nm/step) for fixed increment linear sampling (put in quotes! \"1.0 2.3 5.1 "
808 "Directions for acceptance linear sampling - only sign counts! (put in quotes! \"-1 +1 "
814 "Stepsize (nm/step) for fixed increment radius expansion" },
818 { &edi_params.maxedsteps },
819 "Maximum number of steps per cycle" },
824 "Number of steps to run without any perturbations " },
825 { "-deltaF0", FALSE, etREAL, { &deltaF0 }, "Target destabilization energy for flooding" },
830 "Start deltaF with this parameter - default 0, nonzero values only needed for restart" },
835 "Coupling constant for adaption of flooding strength according to deltaF0, 0 = infinity "
836 "i.e. constant flooding strength" },
841 "The starting value of the flooding strength. The flooding strength is updated "
842 "according to the adaptive flooding scheme. For a constant flooding strength use "
843 "[TT]-tau[tt] 0. " },
848 "T is temperature, the value is needed if you want to do flooding " },
853 "Scale width of gaussian flooding potential with alpha^2 " },
858 "Use the flooding potential with inverted sign -> effects as quasiharmonic restraining "
864 "The eigenvectors and eigenvalues are from a Hessian matrix" },
869 "The eigenvalues are interpreted as spring constant" },
874 "Constant force flooding: manually set the forces for the eigenvectors selected with "
876 "(put in quotes! \"1.0 2.3 5.1 -3.1\"). No other flooding parameters are needed when "
877 "specifying the forces directly." }
879 #define NPA asize(pa)
882 int nvec1, *eignr1 = nullptr;
883 rvec * xav1, **eigvec1 = nullptr;
884 t_atoms* atoms = nullptr;
885 int nav; /* Number of atoms in the average structure */
887 const char* indexfile;
890 int nfit; /* Number of atoms in the reference/fit structure */
891 int ev_class; /* parameter _class i.e. evMON, evRADFIX etc. */
893 real* eigval1 = nullptr; /* in V3.3 this is parameter of read_eigenvectors */
896 const char* TargetFile;
897 const char* OriginFile;
898 const char* EigvecFile;
900 gmx_output_env_t* oenv;
902 /*to read topology file*/
909 t_filenm fnm[] = { { efTRN, "-f", "eigenvec", ffREAD }, { efXVG, "-eig", "eigenval", ffOPTRD },
910 { efTPS, nullptr, nullptr, ffREAD }, { efNDX, nullptr, nullptr, ffOPTRD },
911 { efSTX, "-tar", "target", ffOPTRD }, { efSTX, "-ori", "origin", ffOPTRD },
912 { efEDI, "-o", "sam", ffWRITE } };
913 #define NFILE asize(fnm)
914 edi_params.outfrq = 100;
915 edi_params.slope = 0.0;
916 edi_params.maxedsteps = 0;
917 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
922 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
923 EdiFile = ftp2fn(efEDI, NFILE, fnm);
924 TargetFile = opt2fn_null("-tar", NFILE, fnm);
925 OriginFile = opt2fn_null("-ori", NFILE, fnm);
928 for (ev_class = 0; ev_class < evNr; ++ev_class)
930 if (opt2parg_bSet(evOptions[ev_class], NPA, pa))
932 /*get list of eigenvectors*/
933 nvecs = sscan_list(&(listen[ev_class]), opt2parg_str(evOptions[ev_class], NPA, pa),
934 evOptions[ev_class]);
935 if (ev_class < evStepNr - 2)
937 /*if apropriate get list of stepsizes for these eigenvectors*/
938 if (opt2parg_bSet(evStepOptions[ev_class], NPA, pa))
940 evStepList[ev_class] = scan_vecparams(opt2parg_str(evStepOptions[ev_class], NPA, pa),
941 evStepOptions[ev_class], nvecs);
943 else /*if list is not given fill with zeros */
945 snew(evStepList[ev_class], nvecs);
946 for (i = 0; i < nvecs; i++)
948 evStepList[ev_class][i] = 0.0;
952 else if (ev_class == evRADFIX)
954 snew(evStepList[ev_class], nvecs);
955 for (i = 0; i < nvecs; i++)
957 evStepList[ev_class][i] = radstep;
960 else if (ev_class == evFLOOD)
962 snew(evStepList[ev_class], nvecs);
964 /* Are we doing constant force flooding? In that case, we read in
965 * the fproj values from the command line */
966 if (opt2parg_bSet("-constF", NPA, pa))
968 evStepList[ev_class] =
969 scan_vecparams(opt2parg_str("-constF", NPA, pa), "-constF", nvecs);
974 } /*to avoid ambiguity */
976 else /* if there are no eigenvectors for this option set list to zero */
978 listen[ev_class] = nullptr;
979 snew(listen[ev_class], 1);
980 listen[ev_class][0] = 0;
984 /* print the interpreted list of eigenvectors - to give some feedback*/
985 for (ev_class = 0; ev_class < evNr; ++ev_class)
987 printf("Eigenvector list %7s consists of the indices: ", evOptions[ev_class]);
989 while (listen[ev_class][i])
991 printf("%d ", listen[ev_class][i++]);
996 EigvecFile = opt2fn("-f", NFILE, fnm);
998 /*read eigenvectors from eigvec.trr*/
999 read_eigenvectors(EigvecFile, &nav, &bFit1, &xref1, &edi_params.fitmas, &xav1,
1000 &edi_params.pcamas, &nvec1, &eignr1, &eigvec1, &eigval1);
1002 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtop, nullptr, topbox, false);
1006 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", nav);
1007 get_index(atoms, indexfile, 1, &i, &index, &grpname); /*if indexfile != NULL parameter 'atoms' is ignored */
1010 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, nav);
1015 if (xref1 == nullptr)
1019 /* if g_covar used different coordinate groups to fit and to do the PCA */
1020 printf("\nNote: the structure in %s should be the same\n"
1021 " as the one used for the fit in g_covar\n",
1022 ftp2fn(efTPS, NFILE, fnm));
1023 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1027 printf("\nNote: Apparently no fitting was done in g_covar.\n"
1028 " However, you need to select a reference group for fitting in mdrun\n");
1030 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1032 for (i = 0; i < nfit; i++)
1034 copy_rvec(xtop[ifit[i]], xref1[i]);
1043 if (opt2parg_bSet("-constF", NPA, pa))
1045 /* Constant force flooding is special: Most of the normal flooding
1046 * options are not needed. */
1047 edi_params.flood.bConstForce = TRUE;
1051 /* For normal flooding read eigenvalues and store them in evSteplist[evFLOOD] */
1053 if (listen[evFLOOD][0] != 0)
1055 read_eigenvalues(listen[evFLOOD], opt2fn("-eig", NFILE, fnm), evStepList[evFLOOD],
1056 bHesse, kB * T, nav);
1059 edi_params.flood.tau = tau;
1060 edi_params.flood.deltaF0 = deltaF0;
1061 edi_params.flood.deltaF = deltaF;
1062 edi_params.presteps = eqSteps;
1063 edi_params.flood.kT = kB * T;
1064 edi_params.flood.bHarmonic = bHarmonic;
1067 /* Trick: invert sign of Efl and alpha2 then this will give the same sign in the exponential and inverted sign outside */
1068 edi_params.flood.constEfl = -constEfl;
1069 edi_params.flood.alpha2 = -gmx::square(alpha);
1073 edi_params.flood.constEfl = constEfl;
1074 edi_params.flood.alpha2 = gmx::square(alpha);
1078 edi_params.ned = nav;
1080 /*number of system atoms */
1081 edi_params.nini = atoms->nr;
1084 /*store reference and average structure in edi_params*/
1085 make_t_edx(&edi_params.sref, nfit, xref1, ifit);
1086 make_t_edx(&edi_params.sav, nav, xav1, index);
1089 /* Store target positions in edi_params */
1090 if (opt2bSet("-tar", NFILE, fnm))
1092 if (0 != listen[evFLOOD][0])
1095 "\nNote: Providing a TARGET structure has no effect when using flooding.\n"
1096 " You may want to use -ori to define the flooding potential center.\n\n");
1098 get_structure(atoms, indexfile, TargetFile, &edi_params.star, nfit, ifit, nav, index);
1102 make_t_edx(&edi_params.star, 0, nullptr, index);
1105 /* Store origin positions */
1106 if (opt2bSet("-ori", NFILE, fnm))
1108 get_structure(atoms, indexfile, OriginFile, &edi_params.sori, nfit, ifit, nav, index);
1112 make_t_edx(&edi_params.sori, 0, nullptr, index);
1115 /* Write edi-file */
1116 write_the_whole_thing(gmx_ffopen(EdiFile, "w"), &edi_params, eigvec1, nvec1, listen, evStepList);