2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gromacs/commandline/filenm.h"
49 #include "gromacs/domdec/domdec_struct.h"
50 #include "gromacs/fileio/gmxfio.h"
51 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/gmxlib/network.h"
53 #include "gromacs/gmxlib/nrnb.h"
54 #include "gromacs/linearalgebra/nrjac.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/utilities.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/math/vectypes.h"
59 #include "gromacs/mdlib/broadcaststructs.h"
60 #include "gromacs/mdlib/constr.h"
61 #include "gromacs/mdlib/groupcoord.h"
62 #include "gromacs/mdlib/stat.h"
63 #include "gromacs/mdlib/update.h"
64 #include "gromacs/mdrunutility/handlerestart.h"
65 #include "gromacs/mdtypes/commrec.h"
66 #include "gromacs/mdtypes/edsamhistory.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/md_enums.h"
69 #include "gromacs/mdtypes/observableshistory.h"
70 #include "gromacs/mdtypes/state.h"
71 #include "gromacs/pbcutil/pbc.h"
72 #include "gromacs/topology/mtop_lookup.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/gmxassert.h"
77 #include "gromacs/utility/logger.h"
78 #include "gromacs/utility/smalloc.h"
79 #include "gromacs/utility/strconvert.h"
83 /*! \brief Identifies the type of ED: none, normal ED, flooding. */
84 enum class EssentialDynamicsType
86 //! No essential dynamics
88 //! Essential dynamics sampling
90 //! Essential dynamics flooding
94 /*! \brief Identify on which structure to operate. */
95 enum class EssentialDynamicsStructure
97 //! Reference structure
107 /*! \brief Essential dynamics vector.
108 * TODO: split into working data and input data
109 * NOTE: called eigvec, because traditionally eigenvectors from PCA analysis
110 * were used as essential dynamics vectors, however, vectors used for ED need
111 * not have any special properties
115 //! nr of eigenvectors
117 //! index nrs of eigenvectors
119 //! stepsizes (per eigenvector)
120 real* stpsz = nullptr;
121 //! eigenvector components
122 rvec** vec = nullptr;
123 //! instantaneous x projections
124 real* xproj = nullptr;
125 //! instantaneous f projections
126 real* fproj = nullptr;
127 //! instantaneous radius
129 //! starting or target projections
130 real* refproj = nullptr;
133 /*! \brief Essential dynamics vectors per method implementation.
137 //! only monitored, no constraints
139 //! fixed linear constraints
140 t_eigvec linfix = {};
141 //! acceptance linear constraints
142 t_eigvec linacc = {};
143 //! fixed radial constraints (exp)
144 t_eigvec radfix = {};
145 //! acceptance radial constraints (exp)
146 t_eigvec radacc = {};
147 //! acceptance rad. contraction constr.
148 t_eigvec radcon = {};
151 /*! \brief Essential dynamics flooding parameters and work data.
152 * TODO: split into working data and input parameters
153 * NOTE: The implementation here follows:
154 * O.E. Lange, L.V. Schafer, and H. Grubmuller,
155 * “Flooding in GROMACS: Accelerated barrier crossings in molecular dynamics,”
156 * J. Comp. Chem., 27 1693–1702 (2006)
160 /*! \brief Target destabilisation free energy.
161 * Controls the flooding potential strength.
162 * Linked to the expected speed-up of mean first passage time out of flooded minimum */
164 //! Do not calculate a flooding potential, instead flood with a constant force
165 bool bConstForce = false;
166 //! Relaxation time scale for slowest flooded degree of freedom
168 //! Current estimated destabilisation free energy
170 //! Flooding energy, acting as a proportionality factor for the flooding potential
172 //! Boltzmann constant times temperature, provided by user
174 //! The flooding potential
176 //! Integration time step
178 //! Inital flooding strenth
180 //! Empirical global scaling parameter of the essential dynamics vectors.
182 //! The forces from flooding in atom coordinate space (in contrast to projections onto essential dynamics vectors)
183 rvec* forces_cartesian = nullptr;
184 //! The vectors along which to flood
186 //! Use flooding for harmonic restraint on the eigenvector
187 bool bHarmonic = false;
188 //! The initial reference projection of the flooding vectors. Only with harmonic restraint.
189 real* initialReferenceProjection = nullptr;
190 //! The current reference projection is the initialReferenceProjection + step * slope. Only with harmonic restraint.
191 real* referenceProjectionSlope = nullptr;
196 /* This type is for the average, reference, target, and origin structure */
199 int nr = 0; /* number of atoms this structure contains */
200 int nr_loc = 0; /* number of atoms on local node */
201 int* anrs = nullptr; /* atom index numbers */
202 int* anrs_loc = nullptr; /* local atom index numbers */
203 int nalloc_loc = 0; /* allocation size of anrs_loc */
204 int* c_ind = nullptr; /* at which position of the whole anrs
205 * array is a local atom?, i.e.
206 * c_ind[0...nr_loc-1] gives the atom index
207 * with respect to the collective
208 * anrs[0...nr-1] array */
209 rvec* x = nullptr; /* positions for this structure */
210 rvec* x_old = nullptr; /* Last positions which have the correct PBC
211 representation of the ED group. In
212 combination with keeping track of the
213 shift vectors, the ED group can always
215 real* m = nullptr; /* masses */
216 real mtot = 0.; /* total mass (only used in sref) */
217 real* sqrtm = nullptr; /* sqrt of the masses used for mass-
218 * weighting of analysis (only used in sav) */
224 int nini = 0; /* total Nr of atoms */
225 gmx_bool fitmas = false; /* true if trans fit with cm */
226 gmx_bool pcamas = false; /* true if mass-weighted PCA */
227 int presteps = 0; /* number of steps to run without any
228 * perturbations ... just monitoring */
229 int outfrq = 0; /* freq (in steps) of writing to edo */
230 int maxedsteps = 0; /* max nr of steps per cycle */
232 /* all gmx_edx datasets are copied to all nodes in the parallel case */
233 gmx_edx sref = {}; /* reference positions, to these fitting
235 gmx_bool bRefEqAv = false; /* If true, reference & average indices
236 * are the same. Used for optimization */
237 gmx_edx sav = {}; /* average positions */
238 gmx_edx star = {}; /* target positions */
239 gmx_edx sori = {}; /* origin positions */
241 t_edvecs vecs = {}; /* eigenvectors */
242 real slope = 0; /* minimal slope in acceptance radexp */
244 t_edflood flood = {}; /* parameters especially for flooding */
245 struct t_ed_buffer* buf = nullptr; /* handle to local buffers */
253 EssentialDynamicsType eEDtype = EssentialDynamicsType::None;
254 //! output file pointer
256 std::vector<t_edpar> edpar;
257 gmx_bool bFirst = false;
259 gmx_edsam::~gmx_edsam()
263 /* edo is opened sometimes with xvgropen, sometimes with
264 * gmx_fio_fopen, so we use the least common denominator for
274 rvec old_transvec, older_transvec, transvec_compact;
275 rvec* xcoll; /* Positions from all nodes, this is the
276 collective set we work on.
277 These are the positions of atoms with
278 average structure indices */
279 rvec* xc_ref; /* same but with reference structure indices */
280 ivec* shifts_xcoll; /* Shifts for xcoll */
281 ivec* extra_shifts_xcoll; /* xcoll shift changes since last NS step */
282 ivec* shifts_xc_ref; /* Shifts for xc_ref */
283 ivec* extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
284 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
285 ED shifts for this ED group need to
290 /* definition of ED buffer structure */
293 struct t_fit_to_ref* fit_to_ref;
294 struct t_do_edfit* do_edfit;
295 struct t_do_edsam* do_edsam;
296 struct t_do_radcon* do_radcon;
301 class EssentialDynamics::Impl
304 // TODO: move all fields from gmx_edsam here and remove gmx_edsam
305 gmx_edsam essentialDynamics_;
307 EssentialDynamics::EssentialDynamics() : impl_(new Impl) {}
308 EssentialDynamics::~EssentialDynamics() = default;
310 gmx_edsam* EssentialDynamics::getLegacyED()
312 return &impl_->essentialDynamics_;
316 /* Function declarations */
317 static void fit_to_reference(rvec* xcoll, rvec transvec, matrix rotmat, t_edpar* edi);
318 static void translate_and_rotate(rvec* x, int nat, rvec transvec, matrix rotmat);
319 static real rmsd_from_structure(rvec* x, struct gmx_edx* s);
322 /*! \brief Read in the essential dynamics input file and return its contents.
323 * \param[in] fn the file name of the edi file to be read
324 * \param[in] nr_mdatoms the number of atoms in the simulation
325 * \returns A vector of containing the essentiail dyanmics parameters.
326 * NOTE: edi files may that it may contain several ED data sets from concatenated edi files.
327 * The standard case would be a single ED data set, though. */
328 std::vector<t_edpar> read_edi_file(const char* fn, int nr_mdatoms);
330 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam& ed, edsamhistory_t* EDstate);
331 static void init_edsamstate(const gmx_edsam& ed, edsamhistory_t* EDstate);
332 static void write_edo_legend(gmx_edsam* ed, int nED, const gmx_output_env_t* oenv);
333 /* End function declarations */
335 /* Define formats for the column width in the output file */
336 const char EDcol_sfmt[] = "%17s";
337 const char EDcol_efmt[] = "%17.5e";
338 const char EDcol_ffmt[] = "%17f";
340 /* Define a formats for reading, otherwise cppcheck complains for scanf without width limits */
341 const char max_ev_fmt_d[] = "%7d"; /* Number of eigenvectors. 9,999,999 should be enough */
342 const char max_ev_fmt_lf[] = "%12lf";
343 const char max_ev_fmt_dlf[] = "%7d%12lf";
344 const char max_ev_fmt_dlflflf[] = "%7d%12lf%12lf%12lf";
345 const char max_ev_fmt_lelele[] = "%12le%12le%12le";
349 /*!\brief Determines whether to perform essential dynamics constraints other than flooding.
350 * \param[in] edi the essential dynamics parameters
351 * \returns true if essential dyanmics constraints need to be performed
353 bool bNeedDoEdsam(const t_edpar& edi)
355 return (edi.vecs.mon.neig != 0) || (edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0)
356 || (edi.vecs.radfix.neig != 0) || (edi.vecs.radacc.neig != 0) || (edi.vecs.radcon.neig != 0);
361 /* Multiple ED groups will be labeled with letters instead of numbers
362 * to avoid confusion with eigenvector indices */
363 static char get_EDgroupChar(int nr_edi, int nED)
371 * nr_edi = 2 -> B ...
373 return 'A' + nr_edi - 1;
378 /*! \brief The mass-weighted inner product of two coordinate vectors.
379 * Does not subtract average positions, projection on single eigenvector is returned
380 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
381 * Average position is subtracted in ed_apply_constraints prior to calling projectx
382 * \param[in] edi Essential dynamics parameters
383 * \param[in] xcoll vector of atom coordinates
384 * \param[in] vec vector of coordinates to project onto
385 * \return mass-weighted projection.
387 real projectx(const t_edpar& edi, rvec* xcoll, rvec* vec)
393 for (i = 0; i < edi.sav.nr; i++)
395 proj += edi.sav.sqrtm[i] * iprod(vec[i], xcoll[i]);
400 /*!\brief Project coordinates onto vector after substracting average position.
401 * projection is stored in vec->refproj which is used for radacc, radfix,
402 * radcon and center of flooding potential.
403 * Average positions are first substracted from x, then added back again.
404 * \param[in] edi essential dynamics parameters with average position
405 * \param[in] x Coordinates to be projected
406 * \param[out] vec eigenvector, radius and refproj are overwritten here
408 void rad_project(const t_edpar& edi, rvec* x, t_eigvec* vec)
413 /* Subtract average positions */
414 for (i = 0; i < edi.sav.nr; i++)
416 rvec_dec(x[i], edi.sav.x[i]);
419 for (i = 0; i < vec->neig; i++)
421 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
422 rad += gmx::square((vec->refproj[i] - vec->xproj[i]));
424 vec->radius = sqrt(rad);
426 /* Add average positions */
427 for (i = 0; i < edi.sav.nr; i++)
429 rvec_inc(x[i], edi.sav.x[i]);
433 /*!\brief Projects coordinates onto eigenvectors and stores result in vec->xproj.
434 * Mass-weighting is applied. Subtracts average positions prior to projection and add
435 * them afterwards to retain the unchanged vector.
436 * \param[in] x The coordinates to project to an eigenvector
437 * \param[in,out] vec The eigenvectors
438 * \param[in] edi essential dynamics parameters holding average structure and masses
440 void project_to_eigvectors(rvec* x, t_eigvec* vec, const t_edpar& edi)
447 /* Subtract average positions */
448 for (int i = 0; i < edi.sav.nr; i++)
450 rvec_dec(x[i], edi.sav.x[i]);
453 for (int i = 0; i < vec->neig; i++)
455 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
458 /* Add average positions */
459 for (int i = 0; i < edi.sav.nr; i++)
461 rvec_inc(x[i], edi.sav.x[i]);
466 /* Project vector x onto all edi->vecs (mon, linfix,...) */
467 static void project(rvec* x, /* positions to project */
468 t_edpar* edi) /* edi data set */
470 /* It is not more work to subtract the average position in every
471 * subroutine again, because these routines are rarely used simultaneously */
472 project_to_eigvectors(x, &edi->vecs.mon, *edi);
473 project_to_eigvectors(x, &edi->vecs.linfix, *edi);
474 project_to_eigvectors(x, &edi->vecs.linacc, *edi);
475 project_to_eigvectors(x, &edi->vecs.radfix, *edi);
476 project_to_eigvectors(x, &edi->vecs.radacc, *edi);
477 project_to_eigvectors(x, &edi->vecs.radcon, *edi);
482 /*!\brief Evaluates the distance from reference to current eigenvector projection.
483 * \param[in] vec eigenvector
486 real calc_radius(const t_eigvec& vec)
490 for (int i = 0; i < vec.neig; i++)
492 rad += gmx::square((vec.refproj[i] - vec.xproj[i]));
495 return rad = sqrt(rad);
505 static void do_edfit(int natoms, rvec* xp, rvec* x, matrix R, t_edpar* edi)
507 /* this is a copy of do_fit with some modifications */
508 int c, r, n, j, i, irot;
509 double d[6], xnr, xpc;
514 struct t_do_edfit* loc;
517 if (edi->buf->do_edfit != nullptr)
524 snew(edi->buf->do_edfit, 1);
526 loc = edi->buf->do_edfit;
530 snew(loc->omega, 2 * DIM);
531 snew(loc->om, 2 * DIM);
532 for (i = 0; i < 2 * DIM; i++)
534 snew(loc->omega[i], 2 * DIM);
535 snew(loc->om[i], 2 * DIM);
539 for (i = 0; (i < 6); i++)
542 for (j = 0; (j < 6); j++)
544 loc->omega[i][j] = 0;
549 /* calculate the matrix U */
551 for (n = 0; (n < natoms); n++)
553 for (c = 0; (c < DIM); c++)
556 for (r = 0; (r < DIM); r++)
559 u[c][r] += xnr * xpc;
564 /* construct loc->omega */
565 /* loc->omega is symmetric -> loc->omega==loc->omega' */
566 for (r = 0; (r < 6); r++)
568 for (c = 0; (c <= r); c++)
570 if ((r >= 3) && (c < 3))
572 loc->omega[r][c] = u[r - 3][c];
573 loc->omega[c][r] = u[r - 3][c];
577 loc->omega[r][c] = 0;
578 loc->omega[c][r] = 0;
583 /* determine h and k */
584 jacobi(loc->omega, 6, d, loc->om, &irot);
588 fprintf(stderr, "IROT=0\n");
591 index = 0; /* For the compiler only */
593 for (j = 0; (j < 3); j++)
596 for (i = 0; (i < 6); i++)
605 for (i = 0; (i < 3); i++)
607 vh[j][i] = M_SQRT2 * loc->om[i][index];
608 vk[j][i] = M_SQRT2 * loc->om[i + DIM][index];
613 for (c = 0; (c < 3); c++)
615 for (r = 0; (r < 3); r++)
617 R[c][r] = vk[0][r] * vh[0][c] + vk[1][r] * vh[1][c] + vk[2][r] * vh[2][c];
622 for (c = 0; (c < 3); c++)
624 for (r = 0; (r < 3); r++)
626 R[c][r] = vk[0][r] * vh[0][c] + vk[1][r] * vh[1][c] - vk[2][r] * vh[2][c];
633 static void rmfit(int nat, rvec* xcoll, const rvec transvec, matrix rotmat)
640 * The inverse rotation is described by the transposed rotation matrix */
641 transpose(rotmat, tmat);
642 rotate_x(xcoll, nat, tmat);
644 /* Remove translation */
645 vec[XX] = -transvec[XX];
646 vec[YY] = -transvec[YY];
647 vec[ZZ] = -transvec[ZZ];
648 translate_x(xcoll, nat, vec);
652 /**********************************************************************************
653 ******************** FLOODING ****************************************************
654 **********************************************************************************
656 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
657 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
658 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
660 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
661 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
663 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
664 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
666 do_flood makes a copy of the positions,
667 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
668 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
669 fit. Then do_flood adds these forces to the forcefield-forces
670 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
672 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
673 structure is projected to the system of eigenvectors and then this position in the subspace is used as
674 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
675 i.e. the average structure as given in the make_edi file.
677 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
678 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
679 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
680 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
681 further adaption is applied, Efl will stay constant at zero.
683 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
684 used as spring constants for the harmonic potential.
685 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
686 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
688 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
689 the routine read_edi_file reads all of theses flooding files.
690 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
691 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
692 edi there is no interdependence whatsoever. The forces are added together.
694 To write energies into the .edr file, call the function
695 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
697 get_flood_energies(real Vfl[],int nnames);
700 - one could program the whole thing such that Efl, Vfl and deltaF is written to the .edr file. -- i dont know how to do that, yet.
702 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
703 two edsam files from two peptide chains
706 // TODO split this into multiple files
709 /*!\brief Output flooding simulation settings and results to file.
710 * \param[in] edi Essential dynamics input parameters
711 * \param[in] fp output file
712 * \param[in] rmsd rmsd to reference structure
714 void write_edo_flood(const t_edpar& edi, FILE* fp, real rmsd)
716 /* Output how well we fit to the reference structure */
717 fprintf(fp, EDcol_ffmt, rmsd);
719 for (int i = 0; i < edi.flood.vecs.neig; i++)
721 fprintf(fp, EDcol_efmt, edi.flood.vecs.xproj[i]);
723 /* Check whether the reference projection changes with time (this can happen
724 * in case flooding is used as harmonic restraint). If so, output the
725 * current reference projection */
726 if (edi.flood.bHarmonic && edi.flood.referenceProjectionSlope[i] != 0.0)
728 fprintf(fp, EDcol_efmt, edi.flood.vecs.refproj[i]);
731 /* Output Efl if we are doing adaptive flooding */
732 if (0 != edi.flood.tau)
734 fprintf(fp, EDcol_efmt, edi.flood.Efl);
736 fprintf(fp, EDcol_efmt, edi.flood.Vfl);
738 /* Output deltaF if we are doing adaptive flooding */
739 if (0 != edi.flood.tau)
741 fprintf(fp, EDcol_efmt, edi.flood.deltaF);
743 fprintf(fp, EDcol_efmt, edi.flood.vecs.fproj[i]);
749 /* From flood.xproj compute the Vfl(x) at this point */
750 static real flood_energy(t_edpar* edi, int64_t step)
752 /* compute flooding energy Vfl
753 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
754 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
755 it is already computed by make_edi and stored in stpsz[i]
757 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
764 /* Each time this routine is called (i.e. each time step), we add a small
765 * value to the reference projection. This way a harmonic restraint towards
766 * a moving reference is realized. If no value for the additive constant
767 * is provided in the edi file, the reference will not change. */
768 if (edi->flood.bHarmonic)
770 for (i = 0; i < edi->flood.vecs.neig; i++)
772 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i]
773 + step * edi->flood.referenceProjectionSlope[i];
778 /* Compute sum which will be the exponent of the exponential */
779 for (i = 0; i < edi->flood.vecs.neig; i++)
781 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
782 sum += edi->flood.vecs.stpsz[i] * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i])
783 * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i]);
786 /* Compute the Gauss function*/
787 if (edi->flood.bHarmonic)
789 Vfl = -0.5 * edi->flood.Efl * sum; /* minus sign because Efl is negative, if restrain is on. */
793 Vfl = edi->flood.Efl != 0
795 * std::exp(-edi->flood.kT / 2 / edi->flood.Efl / edi->flood.alpha2 * sum)
803 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
804 static void flood_forces(t_edpar* edi)
806 /* compute the forces in the subspace of the flooding eigenvectors
807 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
810 real energy = edi->flood.Vfl;
813 if (edi->flood.bHarmonic)
815 for (i = 0; i < edi->flood.vecs.neig; i++)
817 edi->flood.vecs.fproj[i] = edi->flood.Efl * edi->flood.vecs.stpsz[i]
818 * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i]);
823 for (i = 0; i < edi->flood.vecs.neig; i++)
825 /* if Efl is zero the forces are zero if not use the formula */
826 edi->flood.vecs.fproj[i] =
828 ? edi->flood.kT / edi->flood.Efl / edi->flood.alpha2 * energy
829 * edi->flood.vecs.stpsz[i]
830 * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i])
838 /*!\brief Raise forces from subspace into cartesian space.
839 * This function lifts the forces from the subspace to the cartesian space
840 * all the values not contained in the subspace are assumed to be zero and then
841 * a coordinate transformation from eigenvector to cartesian vectors is performed
842 * The nonexistent values don't have to be set to zero explicitly, they would occur
843 * as zero valued summands, hence we just stop to compute this part of the sum.
844 * For every atom we add all the contributions to this atom from all the different eigenvectors.
845 * NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
846 * field forces_cart prior the computation, but we compute the forces separately
847 * to have them accessible for diagnostics
849 * \param[in] edi Essential dynamics input parameters
850 * \param[out] forces_cart The cartesian forces
853 void flood_blowup(const t_edpar& edi, rvec* forces_cart)
855 const real* forces_sub = edi.flood.vecs.fproj;
856 /* Calculate the cartesian forces for the local atoms */
858 /* Clear forces first */
859 for (int j = 0; j < edi.sav.nr_loc; j++)
861 clear_rvec(forces_cart[j]);
864 /* Now compute atomwise */
865 for (int j = 0; j < edi.sav.nr_loc; j++)
867 /* Compute forces_cart[edi.sav.anrs[j]] */
868 for (int eig = 0; eig < edi.flood.vecs.neig; eig++)
871 /* Force vector is force * eigenvector (compute only atom j) */
872 svmul(forces_sub[eig], edi.flood.vecs.vec[eig][edi.sav.c_ind[j]], addedForce);
873 /* Add this vector to the cartesian forces */
874 rvec_inc(forces_cart[j], addedForce);
882 /* Update the values of Efl, deltaF depending on tau and Vfl */
883 static void update_adaption(t_edpar* edi)
885 /* this function updates the parameter Efl and deltaF according to the rules given in
886 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
889 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau) > 0.00000001)
891 edi->flood.Efl = edi->flood.Efl
892 + edi->flood.dt / edi->flood.tau * (edi->flood.deltaF0 - edi->flood.deltaF);
893 /* check if restrain (inverted flooding) -> don't let EFL become positive */
894 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
899 edi->flood.deltaF = (1 - edi->flood.dt / edi->flood.tau) * edi->flood.deltaF
900 + edi->flood.dt / edi->flood.tau * edi->flood.Vfl;
905 static void do_single_flood(FILE* edo,
912 gmx_bool bNS) /* Are we in a neighbor searching step? */
915 matrix rotmat; /* rotation matrix */
916 matrix tmat; /* inverse rotation */
917 rvec transvec; /* translation vector */
919 struct t_do_edsam* buf;
922 buf = edi->buf->do_edsam;
925 /* Broadcast the positions of the AVERAGE structure such that they are known on
926 * every processor. Each node contributes its local positions x and stores them in
927 * the collective ED array buf->xcoll */
928 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
929 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind,
930 edi->sav.x_old, box);
932 /* Only assembly REFERENCE positions if their indices differ from the average ones */
935 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref,
936 bNS, x, edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc,
937 edi->sref.c_ind, edi->sref.x_old, box);
940 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
941 * We do not need to update the shifts until the next NS step */
942 buf->bUpdateShifts = FALSE;
944 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
945 * as well as the indices in edi->sav.anrs */
947 /* Fit the reference indices to the reference structure */
950 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
954 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
957 /* Now apply the translation and rotation to the ED structure */
958 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
960 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
961 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
963 if (!edi->flood.bConstForce)
965 /* Compute Vfl(x) from flood.xproj */
966 edi->flood.Vfl = flood_energy(edi, step);
968 update_adaption(edi);
970 /* Compute the flooding forces */
974 /* Translate them into cartesian positions */
975 flood_blowup(*edi, edi->flood.forces_cartesian);
977 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
978 /* Each node rotates back its local forces */
979 transpose(rotmat, tmat);
980 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
982 /* Finally add forces to the main force variable */
983 for (i = 0; i < edi->sav.nr_loc; i++)
985 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
988 /* Output is written by the master process */
989 if (do_per_step(step, edi->outfrq) && MASTER(cr))
991 /* Output how well we fit to the reference */
994 /* Indices of reference and average structures are identical,
995 * thus we can calculate the rmsd to SREF using xcoll */
996 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
1000 /* We have to translate & rotate the reference atoms first */
1001 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1002 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1005 write_edo_flood(*edi, edo, rmsdev);
1010 /* Main flooding routine, called from do_force */
1011 extern void do_flood(const t_commrec* cr,
1012 const t_inputrec* ir,
1020 /* Write time to edo, when required. Output the time anyhow since we need
1021 * it in the output file for ED constraints. */
1022 if (MASTER(cr) && do_per_step(step, ed->edpar.begin()->outfrq))
1024 fprintf(ed->edo, "\n%12f", ir->init_t + step * ir->delta_t);
1027 if (ed->eEDtype != EssentialDynamicsType::Flooding)
1032 for (auto& edi : ed->edpar)
1034 /* Call flooding for one matrix */
1035 if (edi.flood.vecs.neig)
1037 do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS);
1043 /* Called by init_edi, configure some flooding related variables and structures,
1044 * print headers to output files */
1045 static void init_flood(t_edpar* edi, gmx_edsam* ed, real dt)
1050 edi->flood.Efl = edi->flood.constEfl;
1054 if (edi->flood.vecs.neig)
1056 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1057 ed->eEDtype = EssentialDynamicsType::Flooding;
1059 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig,
1060 edi->flood.vecs.neig > 1 ? "s" : "");
1062 if (edi->flood.bConstForce)
1064 /* We have used stpsz as a vehicle to carry the fproj values for constant
1065 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1066 * in const force flooding, fproj is never changed. */
1067 for (i = 0; i < edi->flood.vecs.neig; i++)
1069 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1071 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1072 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1080 /*********** Energy book keeping ******/
1081 static void get_flood_enx_names(t_edpar* edi, char** names, int* nnames) /* get header of energies */
1090 srenew(names, count);
1091 sprintf(buf, "Vfl_%d", count);
1092 names[count - 1] = gmx_strdup(buf);
1093 actual = actual->next_edi;
1096 *nnames = count - 1;
1100 static void get_flood_energies(t_edpar* edi, real Vfl[], int nnames)
1102 /*fl has to be big enough to capture nnames-many entries*/
1111 Vfl[count - 1] = actual->flood.Vfl;
1112 actual = actual->next_edi;
1115 if (nnames != count - 1)
1117 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1120 /************* END of FLOODING IMPLEMENTATION ****************************/
1124 /* This function opens the ED input and output files, reads in all datasets it finds
1125 * in the input file, and cross-checks whether the .edi file information is consistent
1126 * with the essential dynamics data found in the checkpoint file (if present).
1127 * gmx make_edi can be used to create an .edi input file.
1129 static std::unique_ptr<gmx::EssentialDynamics> ed_open(int natoms,
1130 ObservablesHistory* oh,
1131 const char* ediFileName,
1132 const char* edoFileName,
1133 const gmx::StartingBehavior startingBehavior,
1134 const gmx_output_env_t* oenv,
1135 const t_commrec* cr)
1137 auto edHandle = std::make_unique<gmx::EssentialDynamics>();
1138 auto ed = edHandle->getLegacyED();
1139 /* We want to perform ED (this switch might later be upgraded to EssentialDynamicsType::Flooding) */
1140 ed->eEDtype = EssentialDynamicsType::EDSampling;
1144 // If we start from a checkpoint file, we already have an edsamHistory struct
1145 if (oh->edsamHistory == nullptr)
1147 oh->edsamHistory = std::make_unique<edsamhistory_t>(edsamhistory_t{});
1149 edsamhistory_t* EDstate = oh->edsamHistory.get();
1151 /* Read the edi input file: */
1152 ed->edpar = read_edi_file(ediFileName, natoms);
1154 /* Make sure the checkpoint was produced in a run using this .edi file */
1155 if (EDstate->bFromCpt)
1157 crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1161 EDstate->nED = ed->edpar.size();
1163 init_edsamstate(*ed, EDstate);
1165 /* The master opens the ED output file */
1166 if (startingBehavior == gmx::StartingBehavior::RestartWithAppending)
1168 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1172 ed->edo = xvgropen(edoFileName, "Essential dynamics / flooding output", "Time (ps)",
1173 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1175 /* Make a descriptive legend */
1176 write_edo_legend(ed, EDstate->nED, oenv);
1183 /* Broadcasts the structure data */
1184 static void bc_ed_positions(const t_commrec* cr, struct gmx_edx* s, EssentialDynamicsStructure stype)
1186 snew_bc(cr, s->anrs, s->nr); /* Index numbers */
1187 snew_bc(cr, s->x, s->nr); /* Positions */
1188 nblock_bc(cr, s->nr, s->anrs);
1189 nblock_bc(cr, s->nr, s->x);
1191 /* For the average & reference structures we need an array for the collective indices,
1192 * and we need to broadcast the masses as well */
1193 if (stype == EssentialDynamicsStructure::Average || stype == EssentialDynamicsStructure::Reference)
1195 /* We need these additional variables in the parallel case: */
1196 snew(s->c_ind, s->nr); /* Collective indices */
1197 /* Local atom indices get assigned in dd_make_local_group_indices.
1198 * There, also memory is allocated */
1199 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1200 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1201 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1204 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1205 if (stype == EssentialDynamicsStructure::Reference)
1207 snew_bc(cr, s->m, s->nr);
1208 nblock_bc(cr, s->nr, s->m);
1211 /* For the average structure we might need the masses for mass-weighting */
1212 if (stype == EssentialDynamicsStructure::Average)
1214 snew_bc(cr, s->sqrtm, s->nr);
1215 nblock_bc(cr, s->nr, s->sqrtm);
1216 snew_bc(cr, s->m, s->nr);
1217 nblock_bc(cr, s->nr, s->m);
1222 /* Broadcasts the eigenvector data */
1223 static void bc_ed_vecs(const t_commrec* cr, t_eigvec* ev, int length)
1227 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1228 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1229 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1230 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1231 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1233 nblock_bc(cr, ev->neig, ev->ieig);
1234 nblock_bc(cr, ev->neig, ev->stpsz);
1235 nblock_bc(cr, ev->neig, ev->xproj);
1236 nblock_bc(cr, ev->neig, ev->fproj);
1237 nblock_bc(cr, ev->neig, ev->refproj);
1239 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1240 for (i = 0; i < ev->neig; i++)
1242 snew_bc(cr, ev->vec[i], length);
1243 nblock_bc(cr, length, ev->vec[i]);
1248 /* Broadcasts the ED / flooding data to other nodes
1249 * and allocates memory where needed */
1250 static void broadcast_ed_data(const t_commrec* cr, gmx_edsam* ed)
1252 /* Master lets the other nodes know if its ED only or also flooding */
1253 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1255 int numedis = ed->edpar.size();
1256 /* First let everybody know how many ED data sets to expect */
1257 gmx_bcast(sizeof(numedis), &numedis, cr);
1258 nblock_abc(cr, numedis, &(ed->edpar));
1259 /* Now transfer the ED data set(s) */
1260 for (auto& edi : ed->edpar)
1262 /* Broadcast a single ED data set */
1265 /* Broadcast positions */
1266 bc_ed_positions(cr, &(edi.sref),
1267 EssentialDynamicsStructure::Reference); /* reference positions (don't broadcast masses) */
1268 bc_ed_positions(cr, &(edi.sav),
1269 EssentialDynamicsStructure::Average); /* average positions (do broadcast masses as well) */
1270 bc_ed_positions(cr, &(edi.star), EssentialDynamicsStructure::Target); /* target positions */
1271 bc_ed_positions(cr, &(edi.sori), EssentialDynamicsStructure::Origin); /* origin positions */
1273 /* Broadcast eigenvectors */
1274 bc_ed_vecs(cr, &edi.vecs.mon, edi.sav.nr);
1275 bc_ed_vecs(cr, &edi.vecs.linfix, edi.sav.nr);
1276 bc_ed_vecs(cr, &edi.vecs.linacc, edi.sav.nr);
1277 bc_ed_vecs(cr, &edi.vecs.radfix, edi.sav.nr);
1278 bc_ed_vecs(cr, &edi.vecs.radacc, edi.sav.nr);
1279 bc_ed_vecs(cr, &edi.vecs.radcon, edi.sav.nr);
1280 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1281 bc_ed_vecs(cr, &edi.flood.vecs, edi.sav.nr);
1283 /* For harmonic restraints the reference projections can change with time */
1284 if (edi.flood.bHarmonic)
1286 snew_bc(cr, edi.flood.initialReferenceProjection, edi.flood.vecs.neig);
1287 snew_bc(cr, edi.flood.referenceProjectionSlope, edi.flood.vecs.neig);
1288 nblock_bc(cr, edi.flood.vecs.neig, edi.flood.initialReferenceProjection);
1289 nblock_bc(cr, edi.flood.vecs.neig, edi.flood.referenceProjectionSlope);
1295 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1296 static void init_edi(const gmx_mtop_t* mtop, t_edpar* edi)
1299 real totalmass = 0.0;
1302 /* NOTE Init_edi is executed on the master process only
1303 * The initialized data sets are then transmitted to the
1304 * other nodes in broadcast_ed_data */
1306 /* evaluate masses (reference structure) */
1307 snew(edi->sref.m, edi->sref.nr);
1309 for (i = 0; i < edi->sref.nr; i++)
1313 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1317 edi->sref.m[i] = 1.0;
1320 /* Check that every m > 0. Bad things will happen otherwise. */
1321 if (edi->sref.m[i] <= 0.0)
1324 "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1325 "For a mass-weighted fit, all reference structure atoms need to have a mass "
1327 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1328 "atoms from the reference structure by creating a proper index group.\n",
1329 i, edi->sref.anrs[i] + 1, edi->sref.m[i]);
1332 totalmass += edi->sref.m[i];
1334 edi->sref.mtot = totalmass;
1336 /* Masses m and sqrt(m) for the average structure. Note that m
1337 * is needed if forces have to be evaluated in do_edsam */
1338 snew(edi->sav.sqrtm, edi->sav.nr);
1339 snew(edi->sav.m, edi->sav.nr);
1340 for (i = 0; i < edi->sav.nr; i++)
1342 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1345 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1349 edi->sav.sqrtm[i] = 1.0;
1352 /* Check that every m > 0. Bad things will happen otherwise. */
1353 if (edi->sav.sqrtm[i] <= 0.0)
1356 "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1357 "For ED with mass-weighting, all average structure atoms need to have a mass "
1359 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1360 "atoms from the average structure by creating a proper index group.\n",
1361 i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1365 /* put reference structure in origin */
1366 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1370 translate_x(edi->sref.x, edi->sref.nr, com);
1372 /* Init ED buffer */
1377 static void check(const char* line, const char* label)
1379 if (!strstr(line, label))
1382 "Could not find input parameter %s at expected position in edsam input-file "
1383 "(.edi)\nline read instead is %s",
1389 static int read_checked_edint(FILE* file, const char* label)
1391 char line[STRLEN + 1];
1394 fgets2(line, STRLEN, file);
1396 fgets2(line, STRLEN, file);
1397 sscanf(line, max_ev_fmt_d, &idum);
1401 static bool read_checked_edbool(FILE* file, const char* label)
1403 return read_checked_edint(file, label) != 0;
1406 static int read_edint(FILE* file, bool* bEOF)
1408 char line[STRLEN + 1];
1413 eof = fgets2(line, STRLEN, file);
1419 eof = fgets2(line, STRLEN, file);
1425 sscanf(line, max_ev_fmt_d, &idum);
1431 static real read_checked_edreal(FILE* file, const char* label)
1433 char line[STRLEN + 1];
1437 fgets2(line, STRLEN, file);
1439 fgets2(line, STRLEN, file);
1440 sscanf(line, max_ev_fmt_lf, &rdum);
1441 return static_cast<real>(rdum); /* always read as double and convert to single */
1445 static void read_edx(FILE* file, int number, int* anrs, rvec* x)
1448 char line[STRLEN + 1];
1452 for (i = 0; i < number; i++)
1454 fgets2(line, STRLEN, file);
1455 sscanf(line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1456 anrs[i]--; /* we are reading FORTRAN indices */
1457 for (j = 0; j < 3; j++)
1459 x[i][j] = d[j]; /* always read as double and convert to single */
1466 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1467 * \param[in] in the file to read from
1468 * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1469 * \param[out] vec the eigen-vectors
1470 * \param[in] nEig the number of eigenvectors
1472 void scan_edvec(FILE* in, int numAtoms, rvec*** vec, int nEig)
1475 for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1477 snew((*vec)[iEigenvector], numAtoms);
1478 for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1480 char line[STRLEN + 1];
1482 fgets2(line, STRLEN, in);
1483 sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1484 (*vec)[iEigenvector][iAtom][XX] = x;
1485 (*vec)[iEigenvector][iAtom][YY] = y;
1486 (*vec)[iEigenvector][iAtom][ZZ] = z;
1490 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1491 * \param[in] in input file
1492 * \param[in] nrAtoms number of atoms/coordinates
1493 * \param[out] tvec the eigenvector
1495 void read_edvec(FILE* in, int nrAtoms, t_eigvec* tvec)
1497 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1498 if (tvec->neig <= 0)
1503 snew(tvec->ieig, tvec->neig);
1504 snew(tvec->stpsz, tvec->neig);
1505 for (int i = 0; i < tvec->neig; i++)
1507 char line[STRLEN + 1];
1508 fgets2(line, STRLEN, in);
1509 int numEigenVectors;
1511 const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1514 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1516 tvec->ieig[i] = numEigenVectors;
1517 tvec->stpsz[i] = stepSize;
1518 } /* end of loop over eigenvectors */
1520 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1522 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1524 * Eigenvector and its intial reference and slope are stored on the
1525 * same line in the .edi format. To avoid re-winding the .edi file,
1526 * the reference eigen-vector and reference data are read in one go.
1528 * \param[in] in input file
1529 * \param[in] nrAtoms number of atoms/coordinates
1530 * \param[out] tvec the eigenvector
1531 * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1532 * \param[out] referenceProjectionSlope The slope of the reference projections.
1534 bool readEdVecWithReferenceProjection(FILE* in,
1537 real** initialReferenceProjection,
1538 real** referenceProjectionSlope)
1540 bool providesReference = false;
1541 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1542 if (tvec->neig <= 0)
1547 snew(tvec->ieig, tvec->neig);
1548 snew(tvec->stpsz, tvec->neig);
1549 snew(*initialReferenceProjection, tvec->neig);
1550 snew(*referenceProjectionSlope, tvec->neig);
1551 for (int i = 0; i < tvec->neig; i++)
1553 char line[STRLEN + 1];
1554 fgets2(line, STRLEN, in);
1555 int numEigenVectors;
1556 double stepSize = 0.;
1557 double referenceProjection = 0.;
1558 double referenceSlope = 0.;
1559 const int nscan = sscanf(line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize,
1560 &referenceProjection, &referenceSlope);
1561 /* Zero out values which were not scanned */
1565 /* Every 4 values read, including reference position */
1566 providesReference = true;
1569 /* A reference position is provided */
1570 providesReference = true;
1571 /* No value for slope, set to 0 */
1572 referenceSlope = 0.0;
1575 /* No values for reference projection and slope, set to 0 */
1576 referenceProjection = 0.0;
1577 referenceSlope = 0.0;
1581 "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> "
1582 "<refproj> <refproj-slope>\n",
1585 (*initialReferenceProjection)[i] = referenceProjection;
1586 (*referenceProjectionSlope)[i] = referenceSlope;
1588 tvec->ieig[i] = numEigenVectors;
1589 tvec->stpsz[i] = stepSize;
1590 } /* end of loop over eigenvectors */
1592 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1593 return providesReference;
1596 /*!\brief Allocate working memory for the eigenvectors.
1597 * \param[in,out] tvec eigenvector for which memory will be allocated
1599 void setup_edvec(t_eigvec* tvec)
1601 snew(tvec->xproj, tvec->neig);
1602 snew(tvec->fproj, tvec->neig);
1603 snew(tvec->refproj, tvec->neig);
1608 /* Check if the same atom indices are used for reference and average positions */
1609 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1614 /* If the number of atoms differs between the two structures,
1615 * they cannot be identical */
1616 if (sref.nr != sav.nr)
1621 /* Now that we know that both stuctures have the same number of atoms,
1622 * check if also the indices are identical */
1623 for (i = 0; i < sav.nr; i++)
1625 if (sref.anrs[i] != sav.anrs[i])
1631 "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1639 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1640 constexpr int c_oldestUnsupportedMagicNumber = 666;
1641 //! Oldest (lowest) magic number indicating supported essential dynamics input
1642 constexpr int c_oldestSupportedMagicNumber = 669;
1643 //! Latest (highest) magic number indicating supported essential dynamics input
1644 constexpr int c_latestSupportedMagicNumber = 670;
1646 /*!\brief Set up essential dynamics work parameters.
1647 * \param[in] edi Essential dynamics working structure.
1649 void setup_edi(t_edpar* edi)
1651 snew(edi->sref.x_old, edi->sref.nr);
1652 edi->sref.sqrtm = nullptr;
1653 snew(edi->sav.x_old, edi->sav.nr);
1654 if (edi->star.nr > 0)
1656 edi->star.sqrtm = nullptr;
1658 if (edi->sori.nr > 0)
1660 edi->sori.sqrtm = nullptr;
1662 setup_edvec(&edi->vecs.linacc);
1663 setup_edvec(&edi->vecs.mon);
1664 setup_edvec(&edi->vecs.linfix);
1665 setup_edvec(&edi->vecs.linacc);
1666 setup_edvec(&edi->vecs.radfix);
1667 setup_edvec(&edi->vecs.radacc);
1668 setup_edvec(&edi->vecs.radcon);
1669 setup_edvec(&edi->flood.vecs);
1672 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1673 * \param[in] magicNumber indicates the essential dynamics input file version
1674 * \returns true if CONST_FORCE_FLOODING is supported
1676 bool ediFileHasConstForceFlooding(int magicNumber)
1678 return magicNumber > c_oldestSupportedMagicNumber;
1681 /*!\brief Reports reading success of the essential dynamics input file magic number.
1682 * \param[in] in pointer to input file
1683 * \param[out] magicNumber determines the edi file version
1684 * \returns true if setting file version from input file was successful.
1686 bool ediFileHasValidDataSet(FILE* in, int* magicNumber)
1688 /* the edi file is not free format, so expect problems if the input is corrupt. */
1689 bool reachedEndOfFile = true;
1690 /* check the magic number */
1691 *magicNumber = read_edint(in, &reachedEndOfFile);
1692 /* Check whether we have reached the end of the input file */
1693 return !reachedEndOfFile;
1696 /*!\brief Trigger fatal error if magic number is unsupported.
1697 * \param[in] magicNumber A magic number read from the edi file
1698 * \param[in] fn name of input file for error reporting
1700 void exitWhenMagicNumberIsInvalid(int magicNumber, const char* fn)
1703 if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1705 if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1708 "Wrong magic number: Use newest version of make_edi to produce edi file");
1712 gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1717 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1719 * \param[in] in input file
1720 * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1721 * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1722 * \param[in] fn the file name of the input file for error reporting
1723 * \returns edi essential dynamics parameters
1725 t_edpar read_edi(FILE* in, int nr_mdatoms, bool hasConstForceFlooding, const char* fn)
1729 /* check the number of atoms */
1730 edi.nini = read_edint(in, &bEOF);
1731 if (edi.nini != nr_mdatoms)
1733 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi.nini,
1737 /* Done checking. For the rest we blindly trust the input */
1738 edi.fitmas = read_checked_edbool(in, "FITMAS");
1739 edi.pcamas = read_checked_edbool(in, "ANALYSIS_MAS");
1740 edi.outfrq = read_checked_edint(in, "OUTFRQ");
1741 edi.maxedsteps = read_checked_edint(in, "MAXLEN");
1742 edi.slope = read_checked_edreal(in, "SLOPECRIT");
1744 edi.presteps = read_checked_edint(in, "PRESTEPS");
1745 edi.flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1746 edi.flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1747 edi.flood.tau = read_checked_edreal(in, "TAU");
1748 edi.flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1749 edi.flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1750 edi.flood.kT = read_checked_edreal(in, "KT");
1751 edi.flood.bHarmonic = read_checked_edbool(in, "HARMONIC");
1752 if (hasConstForceFlooding)
1754 edi.flood.bConstForce = read_checked_edbool(in, "CONST_FORCE_FLOODING");
1758 edi.flood.bConstForce = FALSE;
1760 edi.sref.nr = read_checked_edint(in, "NREF");
1762 /* allocate space for reference positions and read them */
1763 snew(edi.sref.anrs, edi.sref.nr);
1764 snew(edi.sref.x, edi.sref.nr);
1765 read_edx(in, edi.sref.nr, edi.sref.anrs, edi.sref.x);
1767 /* average positions. they define which atoms will be used for ED sampling */
1768 edi.sav.nr = read_checked_edint(in, "NAV");
1769 snew(edi.sav.anrs, edi.sav.nr);
1770 snew(edi.sav.x, edi.sav.nr);
1771 read_edx(in, edi.sav.nr, edi.sav.anrs, edi.sav.x);
1773 /* Check if the same atom indices are used for reference and average positions */
1774 edi.bRefEqAv = check_if_same(edi.sref, edi.sav);
1778 read_edvec(in, edi.sav.nr, &edi.vecs.mon);
1779 read_edvec(in, edi.sav.nr, &edi.vecs.linfix);
1780 read_edvec(in, edi.sav.nr, &edi.vecs.linacc);
1781 read_edvec(in, edi.sav.nr, &edi.vecs.radfix);
1782 read_edvec(in, edi.sav.nr, &edi.vecs.radacc);
1783 read_edvec(in, edi.sav.nr, &edi.vecs.radcon);
1785 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1786 bool bHaveReference = false;
1787 if (edi.flood.bHarmonic)
1789 bHaveReference = readEdVecWithReferenceProjection(in, edi.sav.nr, &edi.flood.vecs,
1790 &edi.flood.initialReferenceProjection,
1791 &edi.flood.referenceProjectionSlope);
1795 read_edvec(in, edi.sav.nr, &edi.flood.vecs);
1798 /* target positions */
1799 edi.star.nr = read_edint(in, &bEOF);
1800 if (edi.star.nr > 0)
1802 snew(edi.star.anrs, edi.star.nr);
1803 snew(edi.star.x, edi.star.nr);
1804 read_edx(in, edi.star.nr, edi.star.anrs, edi.star.x);
1807 /* positions defining origin of expansion circle */
1808 edi.sori.nr = read_edint(in, &bEOF);
1809 if (edi.sori.nr > 0)
1813 /* Both an -ori structure and a at least one manual reference point have been
1814 * specified. That's ambiguous and probably not intentional. */
1816 "ED: An origin structure has been provided and a at least one (moving) "
1818 " point was manually specified in the edi file. That is ambiguous. "
1821 snew(edi.sori.anrs, edi.sori.nr);
1822 snew(edi.sori.x, edi.sori.nr);
1823 read_edx(in, edi.sori.nr, edi.sori.anrs, edi.sori.x);
1828 std::vector<t_edpar> read_edi_file(const char* fn, int nr_mdatoms)
1830 std::vector<t_edpar> essentialDynamicsParameters;
1832 /* This routine is executed on the master only */
1834 /* Open the .edi parameter input file */
1835 in = gmx_fio_fopen(fn, "r");
1836 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1838 /* Now read a sequence of ED input parameter sets from the edi file */
1839 int ediFileMagicNumber;
1840 while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1842 exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1843 bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1844 auto edi = read_edi(in, nr_mdatoms, hasConstForceFlooding, fn);
1846 essentialDynamicsParameters.emplace_back(edi);
1848 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1849 /* We need to allocate space for the data: */
1851 if (essentialDynamicsParameters.empty())
1853 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1856 fprintf(stderr, "ED: Found %zu ED group%s.\n", essentialDynamicsParameters.size(),
1857 essentialDynamicsParameters.size() > 1 ? "s" : "");
1859 /* Close the .edi file again */
1862 return essentialDynamicsParameters;
1864 } // anonymous namespace
1869 rvec* xcopy; /* Working copy of the positions in fit_to_reference */
1872 /* Fit the current positions to the reference positions
1873 * Do not actually do the fit, just return rotation and translation.
1874 * Note that the COM of the reference structure was already put into
1875 * the origin by init_edi. */
1876 static void fit_to_reference(rvec* xcoll, /* The positions to be fitted */
1877 rvec transvec, /* The translation vector */
1878 matrix rotmat, /* The rotation matrix */
1879 t_edpar* edi) /* Just needed for do_edfit */
1881 rvec com; /* center of mass */
1883 struct t_fit_to_ref* loc;
1886 /* Allocate memory the first time this routine is called for each edi group */
1887 if (nullptr == edi->buf->fit_to_ref)
1889 snew(edi->buf->fit_to_ref, 1);
1890 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1892 loc = edi->buf->fit_to_ref;
1894 /* We do not touch the original positions but work on a copy. */
1895 for (i = 0; i < edi->sref.nr; i++)
1897 copy_rvec(xcoll[i], loc->xcopy[i]);
1900 /* Calculate the center of mass */
1901 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1903 transvec[XX] = -com[XX];
1904 transvec[YY] = -com[YY];
1905 transvec[ZZ] = -com[ZZ];
1907 /* Subtract the center of mass from the copy */
1908 translate_x(loc->xcopy, edi->sref.nr, transvec);
1910 /* Determine the rotation matrix */
1911 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1915 static void translate_and_rotate(rvec* x, /* The positions to be translated and rotated */
1916 int nat, /* How many positions are there? */
1917 rvec transvec, /* The translation vector */
1918 matrix rotmat) /* The rotation matrix */
1921 translate_x(x, nat, transvec);
1924 rotate_x(x, nat, rotmat);
1928 /* Gets the rms deviation of the positions to the structure s */
1929 /* fit_to_structure has to be called before calling this routine! */
1930 static real rmsd_from_structure(rvec* x, /* The positions under consideration */
1931 struct gmx_edx* s) /* The structure from which the rmsd shall be computed */
1937 for (i = 0; i < s->nr; i++)
1939 rmsd += distance2(s->x[i], x[i]);
1942 rmsd /= static_cast<real>(s->nr);
1949 void dd_make_local_ed_indices(gmx_domdec_t* dd, struct gmx_edsam* ed)
1951 if (ed->eEDtype != EssentialDynamicsType::None)
1953 /* Loop over ED groups */
1955 for (auto& edi : ed->edpar)
1957 /* Local atoms of the reference structure (for fitting), need only be assembled
1958 * if their indices differ from the average ones */
1961 dd_make_local_group_indices(dd->ga2la, edi.sref.nr, edi.sref.anrs, &edi.sref.nr_loc,
1962 &edi.sref.anrs_loc, &edi.sref.nalloc_loc, edi.sref.c_ind);
1965 /* Local atoms of the average structure (on these ED will be performed) */
1966 dd_make_local_group_indices(dd->ga2la, edi.sav.nr, edi.sav.anrs, &edi.sav.nr_loc,
1967 &edi.sav.anrs_loc, &edi.sav.nalloc_loc, edi.sav.c_ind);
1969 /* Indicate that the ED shift vectors for this structure need to be updated
1970 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1971 edi.buf->do_edsam->bUpdateShifts = TRUE;
1977 static inline void ed_unshift_single_coord(const matrix box, const rvec x, const ivec is, rvec xu)
1988 xu[XX] = x[XX] - tx * box[XX][XX] - ty * box[YY][XX] - tz * box[ZZ][XX];
1989 xu[YY] = x[YY] - ty * box[YY][YY] - tz * box[ZZ][YY];
1990 xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
1994 xu[XX] = x[XX] - tx * box[XX][XX];
1995 xu[YY] = x[YY] - ty * box[YY][YY];
1996 xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
2002 /*!\brief Apply fixed linear constraints to essential dynamics variable.
2003 * \param[in,out] xcoll The collected coordinates.
2004 * \param[in] edi the essential dynamics parameters
2005 * \param[in] step the current simulation step
2007 void do_linfix(rvec* xcoll, const t_edpar& edi, int64_t step)
2009 /* loop over linfix vectors */
2010 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2012 /* calculate the projection */
2013 real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
2015 /* calculate the correction */
2016 real preFactor = edi.vecs.linfix.refproj[i] + step * edi.vecs.linfix.stpsz[i] - proj;
2018 /* apply the correction */
2019 preFactor /= edi.sav.sqrtm[i];
2020 for (int j = 0; j < edi.sav.nr; j++)
2022 rvec differenceVector;
2023 svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
2024 rvec_inc(xcoll[j], differenceVector);
2029 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
2030 * \param[in,out] xcoll The collected coordinates.
2031 * \param[in] edi the essential dynamics parameters
2033 void do_linacc(rvec* xcoll, t_edpar* edi)
2035 /* loop over linacc vectors */
2036 for (int i = 0; i < edi->vecs.linacc.neig; i++)
2038 /* calculate the projection */
2039 real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
2041 /* calculate the correction */
2042 real preFactor = 0.0;
2043 if (edi->vecs.linacc.stpsz[i] > 0.0)
2045 if ((proj - edi->vecs.linacc.refproj[i]) < 0.0)
2047 preFactor = edi->vecs.linacc.refproj[i] - proj;
2050 if (edi->vecs.linacc.stpsz[i] < 0.0)
2052 if ((proj - edi->vecs.linacc.refproj[i]) > 0.0)
2054 preFactor = edi->vecs.linacc.refproj[i] - proj;
2058 /* apply the correction */
2059 preFactor /= edi->sav.sqrtm[i];
2060 for (int j = 0; j < edi->sav.nr; j++)
2062 rvec differenceVector;
2063 svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2064 rvec_inc(xcoll[j], differenceVector);
2067 /* new positions will act as reference */
2068 edi->vecs.linacc.refproj[i] = proj + preFactor;
2074 static void do_radfix(rvec* xcoll, t_edpar* edi)
2077 real *proj, rad = 0.0, ratio;
2081 if (edi->vecs.radfix.neig == 0)
2086 snew(proj, edi->vecs.radfix.neig);
2088 /* loop over radfix vectors */
2089 for (i = 0; i < edi->vecs.radfix.neig; i++)
2091 /* calculate the projections, radius */
2092 proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2093 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2097 ratio = (edi->vecs.radfix.stpsz[0] + edi->vecs.radfix.radius) / rad - 1.0;
2098 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2100 /* loop over radfix vectors */
2101 for (i = 0; i < edi->vecs.radfix.neig; i++)
2103 proj[i] -= edi->vecs.radfix.refproj[i];
2105 /* apply the correction */
2106 proj[i] /= edi->sav.sqrtm[i];
2108 for (j = 0; j < edi->sav.nr; j++)
2110 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2111 rvec_inc(xcoll[j], vec_dum);
2119 static void do_radacc(rvec* xcoll, t_edpar* edi)
2122 real *proj, rad = 0.0, ratio = 0.0;
2126 if (edi->vecs.radacc.neig == 0)
2131 snew(proj, edi->vecs.radacc.neig);
2133 /* loop over radacc vectors */
2134 for (i = 0; i < edi->vecs.radacc.neig; i++)
2136 /* calculate the projections, radius */
2137 proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2138 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2142 /* only correct when radius decreased */
2143 if (rad < edi->vecs.radacc.radius)
2145 ratio = edi->vecs.radacc.radius / rad - 1.0;
2149 edi->vecs.radacc.radius = rad;
2152 /* loop over radacc vectors */
2153 for (i = 0; i < edi->vecs.radacc.neig; i++)
2155 proj[i] -= edi->vecs.radacc.refproj[i];
2157 /* apply the correction */
2158 proj[i] /= edi->sav.sqrtm[i];
2160 for (j = 0; j < edi->sav.nr; j++)
2162 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2163 rvec_inc(xcoll[j], vec_dum);
2175 static void do_radcon(rvec* xcoll, t_edpar* edi)
2178 real rad = 0.0, ratio = 0.0;
2179 struct t_do_radcon* loc;
2184 if (edi->buf->do_radcon != nullptr)
2191 snew(edi->buf->do_radcon, 1);
2193 loc = edi->buf->do_radcon;
2195 if (edi->vecs.radcon.neig == 0)
2202 snew(loc->proj, edi->vecs.radcon.neig);
2205 /* loop over radcon vectors */
2206 for (i = 0; i < edi->vecs.radcon.neig; i++)
2208 /* calculate the projections, radius */
2209 loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2210 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2213 /* only correct when radius increased */
2214 if (rad > edi->vecs.radcon.radius)
2216 ratio = edi->vecs.radcon.radius / rad - 1.0;
2218 /* loop over radcon vectors */
2219 for (i = 0; i < edi->vecs.radcon.neig; i++)
2221 /* apply the correction */
2222 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2223 loc->proj[i] /= edi->sav.sqrtm[i];
2224 loc->proj[i] *= ratio;
2226 for (j = 0; j < edi->sav.nr; j++)
2228 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2229 rvec_inc(xcoll[j], vec_dum);
2235 edi->vecs.radcon.radius = rad;
2240 static void ed_apply_constraints(rvec* xcoll, t_edpar* edi, int64_t step)
2245 /* subtract the average positions */
2246 for (i = 0; i < edi->sav.nr; i++)
2248 rvec_dec(xcoll[i], edi->sav.x[i]);
2251 /* apply the constraints */
2254 do_linfix(xcoll, *edi, step);
2256 do_linacc(xcoll, edi);
2259 do_radfix(xcoll, edi);
2261 do_radacc(xcoll, edi);
2262 do_radcon(xcoll, edi);
2264 /* add back the average positions */
2265 for (i = 0; i < edi->sav.nr; i++)
2267 rvec_inc(xcoll[i], edi->sav.x[i]);
2274 /*!\brief Write out the projections onto the eigenvectors.
2275 * The order of output corresponds to ed_output_legend().
2276 * \param[in] edi The essential dyanmics parameters
2277 * \param[in] fp The output file
2278 * \param[in] rmsd the rmsd to the reference structure
2280 void write_edo(const t_edpar& edi, FILE* fp, real rmsd)
2282 /* Output how well we fit to the reference structure */
2283 fprintf(fp, EDcol_ffmt, rmsd);
2285 for (int i = 0; i < edi.vecs.mon.neig; i++)
2287 fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2290 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2292 fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2295 for (int i = 0; i < edi.vecs.linacc.neig; i++)
2297 fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2300 for (int i = 0; i < edi.vecs.radfix.neig; i++)
2302 fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2304 if (edi.vecs.radfix.neig)
2306 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2309 for (int i = 0; i < edi.vecs.radacc.neig; i++)
2311 fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2313 if (edi.vecs.radacc.neig)
2315 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2318 for (int i = 0; i < edi.vecs.radcon.neig; i++)
2320 fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2322 if (edi.vecs.radcon.neig)
2324 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2330 /* Returns if any constraints are switched on */
2331 static bool ed_constraints(EssentialDynamicsType edtype, const t_edpar& edi)
2333 if (edtype == EssentialDynamicsType::EDSampling || edtype == EssentialDynamicsType::Flooding)
2335 return ((edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0) || (edi.vecs.radfix.neig != 0)
2336 || (edi.vecs.radacc.neig != 0) || (edi.vecs.radcon.neig != 0));
2342 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for
2343 * flooding/ umbrella sampling simulations. */
2344 static void copyEvecReference(t_eigvec* floodvecs, real* initialReferenceProjection)
2349 if (nullptr == initialReferenceProjection)
2351 snew(initialReferenceProjection, floodvecs->neig);
2354 for (i = 0; i < floodvecs->neig; i++)
2356 initialReferenceProjection[i] = floodvecs->refproj[i];
2361 /* Call on MASTER only. Check whether the essential dynamics / flooding
2362 * groups of the checkpoint file are consistent with the provided .edi file. */
2363 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam& ed, edsamhistory_t* EDstate)
2365 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2368 "Essential dynamics and flooding can only be switched on (or off) at the\n"
2369 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2370 "it must also continue with/without ED constraints when checkpointing.\n"
2371 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2372 "from without a checkpoint.\n");
2375 for (size_t edinum = 0; edinum < ed.edpar.size(); ++edinum)
2377 /* Check number of atoms in the reference and average structures */
2378 if (EDstate->nref[edinum] != ed.edpar[edinum].sref.nr)
2381 "The number of reference structure atoms in ED group %c is\n"
2382 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2383 get_EDgroupChar(edinum + 1, 0), EDstate->nref[edinum], ed.edpar[edinum].sref.nr);
2385 if (EDstate->nav[edinum] != ed.edpar[edinum].sav.nr)
2388 "The number of average structure atoms in ED group %c is\n"
2389 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2390 get_EDgroupChar(edinum + 1, 0), EDstate->nav[edinum], ed.edpar[edinum].sav.nr);
2394 if (gmx::ssize(ed.edpar) != EDstate->nED)
2397 "The number of essential dynamics / flooding groups is not consistent.\n"
2398 "There are %d ED groups in the .cpt file, but %zu in the .edi file!\n"
2399 "Are you sure this is the correct .edi file?\n",
2400 EDstate->nED, ed.edpar.size());
2405 /* The edsamstate struct stores the information we need to make the ED group
2406 * whole again after restarts from a checkpoint file. Here we do the following:
2407 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2408 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2409 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2410 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2411 * all ED structures at the last time step. */
2412 static void init_edsamstate(const gmx_edsam& ed, edsamhistory_t* EDstate)
2414 snew(EDstate->old_sref_p, EDstate->nED);
2415 snew(EDstate->old_sav_p, EDstate->nED);
2417 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2418 if (!EDstate->bFromCpt)
2420 snew(EDstate->nref, EDstate->nED);
2421 snew(EDstate->nav, EDstate->nED);
2424 /* Loop over all ED/flooding data sets (usually only one, though) */
2425 for (int nr_edi = 0; nr_edi < EDstate->nED; nr_edi++)
2427 const auto& edi = ed.edpar[nr_edi];
2428 /* We always need the last reference and average positions such that
2429 * in the next time step we can make the ED group whole again
2430 * if the atoms do not have the correct PBC representation */
2431 if (EDstate->bFromCpt)
2433 /* Copy the last whole positions of reference and average group from .cpt */
2434 for (int i = 0; i < edi.sref.nr; i++)
2436 copy_rvec(EDstate->old_sref[nr_edi][i], edi.sref.x_old[i]);
2438 for (int i = 0; i < edi.sav.nr; i++)
2440 copy_rvec(EDstate->old_sav[nr_edi][i], edi.sav.x_old[i]);
2445 EDstate->nref[nr_edi] = edi.sref.nr;
2446 EDstate->nav[nr_edi] = edi.sav.nr;
2449 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2450 EDstate->old_sref_p[nr_edi] = edi.sref.x_old;
2451 EDstate->old_sav_p[nr_edi] = edi.sav.x_old;
2456 /* Adds 'buf' to 'str' */
2457 static void add_to_string(char** str, const char* buf)
2462 len = strlen(*str) + strlen(buf) + 1;
2468 static void add_to_string_aligned(char** str, const char* buf)
2470 char buf_aligned[STRLEN];
2472 sprintf(buf_aligned, EDcol_sfmt, buf);
2473 add_to_string(str, buf_aligned);
2477 static void nice_legend(const char*** setname,
2484 auto tmp = gmx::formatString("%c %s", EDgroupchar, value);
2485 add_to_string_aligned(LegendStr, tmp.c_str());
2486 tmp += gmx::formatString(" (%s)", unit);
2487 (*setname)[*nsets] = gmx_strdup(tmp.c_str());
2492 static void nice_legend_evec(const char*** setname,
2503 for (i = 0; i < evec->neig; i++)
2505 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2506 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2511 /* Makes a legend for the xvg output file. Call on MASTER only! */
2512 static void write_edo_legend(gmx_edsam* ed, int nED, const gmx_output_env_t* oenv)
2515 int nr_edi, nsets, n_flood, n_edsam;
2516 const char** setname;
2518 char* LegendStr = nullptr;
2521 auto edi = ed->edpar.begin();
2523 fprintf(ed->edo, "# Output will be written every %d step%s\n", edi->outfrq,
2524 edi->outfrq != 1 ? "s" : "");
2526 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2528 fprintf(ed->edo, "#\n");
2529 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n",
2530 get_EDgroupChar(nr_edi, nED));
2531 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2532 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig,
2533 edi->vecs.mon.neig != 1 ? "s" : "");
2534 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig,
2535 edi->vecs.linfix.neig != 1 ? "s" : "");
2536 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig,
2537 edi->vecs.linacc.neig != 1 ? "s" : "");
2538 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig,
2539 edi->vecs.radfix.neig != 1 ? "s" : "");
2540 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig,
2541 edi->vecs.radacc.neig != 1 ? "s" : "");
2542 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig,
2543 edi->vecs.radcon.neig != 1 ? "s" : "");
2544 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig,
2545 edi->flood.vecs.neig != 1 ? "s" : "");
2547 if (edi->flood.vecs.neig)
2549 /* If in any of the groups we find a flooding vector, flooding is turned on */
2550 ed->eEDtype = EssentialDynamicsType::Flooding;
2552 /* Print what flavor of flooding we will do */
2553 if (0 == edi->flood.tau) /* constant flooding strength */
2555 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2556 if (edi->flood.bHarmonic)
2558 fprintf(ed->edo, ", harmonic");
2561 else /* adaptive flooding */
2563 fprintf(ed->edo, ", adaptive");
2566 fprintf(ed->edo, "\n");
2571 /* Print a nice legend */
2573 LegendStr[0] = '\0';
2574 sprintf(buf, "# %6s", "time");
2575 add_to_string(&LegendStr, buf);
2577 /* Calculate the maximum number of columns we could end up with */
2578 edi = ed->edpar.begin();
2580 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2582 nsets += 5 + edi->vecs.mon.neig + edi->vecs.linfix.neig + edi->vecs.linacc.neig
2583 + edi->vecs.radfix.neig + edi->vecs.radacc.neig + edi->vecs.radcon.neig
2584 + 6 * edi->flood.vecs.neig;
2587 snew(setname, nsets);
2589 /* In the mdrun time step in a first function call (do_flood()) the flooding
2590 * forces are calculated and in a second function call (do_edsam()) the
2591 * ED constraints. To get a corresponding legend, we need to loop twice
2592 * over the edi groups and output first the flooding, then the ED part */
2594 /* The flooding-related legend entries, if flooding is done */
2596 if (EssentialDynamicsType::Flooding == ed->eEDtype)
2598 edi = ed->edpar.begin();
2599 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2601 /* Always write out the projection on the flooding EVs. Of course, this can also
2602 * be achieved with the monitoring option in do_edsam() (if switched on by the
2603 * user), but in that case the positions need to be communicated in do_edsam(),
2604 * which is not necessary when doing flooding only. */
2605 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2607 for (i = 0; i < edi->flood.vecs.neig; i++)
2609 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2610 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2612 /* Output the current reference projection if it changes with time;
2613 * this can happen when flooding is used as harmonic restraint */
2614 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2616 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2617 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2620 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2621 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2623 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2624 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2627 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2628 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2630 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2632 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2633 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2636 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2637 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2641 } /* End of flooding-related legend entries */
2645 /* Now the ED-related entries, if essential dynamics is done */
2646 edi = ed->edpar.begin();
2647 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2649 if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2651 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2653 /* Essential dynamics, projections on eigenvectors */
2654 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon,
2655 get_EDgroupChar(nr_edi, nED), "MON");
2656 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix,
2657 get_EDgroupChar(nr_edi, nED), "LINFIX");
2658 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc,
2659 get_EDgroupChar(nr_edi, nED), "LINACC");
2660 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix,
2661 get_EDgroupChar(nr_edi, nED), "RADFIX");
2662 if (edi->vecs.radfix.neig)
2664 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm",
2665 get_EDgroupChar(nr_edi, nED));
2667 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc,
2668 get_EDgroupChar(nr_edi, nED), "RADACC");
2669 if (edi->vecs.radacc.neig)
2671 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm",
2672 get_EDgroupChar(nr_edi, nED));
2674 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon,
2675 get_EDgroupChar(nr_edi, nED), "RADCON");
2676 if (edi->vecs.radcon.neig)
2678 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm",
2679 get_EDgroupChar(nr_edi, nED));
2683 } /* end of 'pure' essential dynamics legend entries */
2684 n_edsam = nsets - n_flood;
2686 xvgr_legend(ed->edo, nsets, setname, oenv);
2691 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2692 n_flood, 1 == n_flood ? "" : "s", n_edsam, 1 == n_edsam ? "" : "s");
2693 fprintf(ed->edo, "%s", LegendStr);
2700 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2701 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2702 std::unique_ptr<gmx::EssentialDynamics> init_edsam(const gmx::MDLogger& mdlog,
2703 const char* ediFileName,
2704 const char* edoFileName,
2705 const gmx_mtop_t* mtop,
2706 const t_inputrec* ir,
2707 const t_commrec* cr,
2708 gmx::Constraints* constr,
2709 const t_state* globalState,
2710 ObservablesHistory* oh,
2711 const gmx_output_env_t* oenv,
2712 const gmx::StartingBehavior startingBehavior)
2715 rvec* x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2716 rvec * xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2717 rvec fit_transvec; /* translation ... */
2718 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2719 rvec* ref_x_old = nullptr; /* helper pointer */
2724 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2726 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName))
2729 "The checkpoint file you provided is from an essential dynamics or flooding\n"
2730 "simulation. Please also set the .edi file on the command line with -ei.\n");
2736 "Activating essential dynamics via files passed to "
2737 "gmx mdrun -edi will change in future to be files passed to "
2738 "gmx grompp and the related .mdp options may change also.");
2740 /* Open input and output files, allocate space for ED data structure */
2741 auto edHandle = ed_open(mtop->natoms, oh, ediFileName, edoFileName, startingBehavior, oenv, cr);
2742 auto ed = edHandle->getLegacyED();
2743 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2744 constr->saveEdsamPointer(ed);
2746 /* Needed for initializing radacc radius in do_edsam */
2749 /* The input file is read by the master and the edi structures are
2750 * initialized here. Input is stored in ed->edpar. Then the edi
2751 * structures are transferred to the other nodes */
2754 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2755 * flooding vector, Essential dynamics can be applied to more than one structure
2756 * as well, but will be done in the order given in the edi file, so
2757 * expect different results for different order of edi file concatenation! */
2758 for (auto& edi : ed->edpar)
2760 init_edi(mtop, &edi);
2761 init_flood(&edi, ed, ir->delta_t);
2765 /* The master does the work here. The other nodes get the positions
2766 * not before dd_partition_system which is called after init_edsam */
2769 edsamhistory_t* EDstate = oh->edsamHistory.get();
2771 if (!EDstate->bFromCpt)
2773 /* Remove PBC, make molecule(s) subject to ED whole. */
2774 snew(x_pbc, mtop->natoms);
2775 copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop->natoms);
2776 do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
2778 /* Reset pointer to first ED data set which contains the actual ED data */
2779 auto edi = ed->edpar.begin();
2780 GMX_ASSERT(!ed->edpar.empty(), "Essential dynamics parameters is undefined");
2782 /* Loop over all ED/flooding data sets (usually only one, though) */
2783 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2785 /* For multiple ED groups we use the output frequency that was specified
2786 * in the first set */
2789 edi->outfrq = ed->edpar.begin()->outfrq;
2792 /* Extract the initial reference and average positions. When starting
2793 * from .cpt, these have already been read into sref.x_old
2794 * in init_edsamstate() */
2795 if (!EDstate->bFromCpt)
2797 /* If this is the first run (i.e. no checkpoint present) we assume
2798 * that the starting positions give us the correct PBC representation */
2799 for (i = 0; i < edi->sref.nr; i++)
2801 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2804 for (i = 0; i < edi->sav.nr; i++)
2806 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2810 /* Now we have the PBC-correct start positions of the reference and
2811 average structure. We copy that over to dummy arrays on which we
2812 can apply fitting to print out the RMSD. We srenew the memory since
2813 the size of the buffers is likely different for every ED group */
2814 srenew(xfit, edi->sref.nr);
2815 srenew(xstart, edi->sav.nr);
2818 /* Reference indices are the same as average indices */
2819 ref_x_old = edi->sav.x_old;
2823 ref_x_old = edi->sref.x_old;
2825 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2826 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2828 /* Make the fit to the REFERENCE structure, get translation and rotation */
2829 fit_to_reference(xfit, fit_transvec, fit_rotmat, &(*edi));
2831 /* Output how well we fit to the reference at the start */
2832 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2833 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2834 rmsd_from_structure(xfit, &edi->sref));
2835 if (EDstate->nED > 1)
2837 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2839 fprintf(stderr, "\n");
2841 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2842 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2844 /* calculate initial projections */
2845 project(xstart, &(*edi));
2847 /* For the target and origin structure both a reference (fit) and an
2848 * average structure can be provided in make_edi. If both structures
2849 * are the same, make_edi only stores one of them in the .edi file.
2850 * If they differ, first the fit and then the average structure is stored
2851 * in star (or sor), thus the number of entries in star/sor is
2852 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2853 * the size of the average group. */
2855 /* process target structure, if required */
2856 if (edi->star.nr > 0)
2858 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2860 /* get translation & rotation for fit of target structure to reference structure */
2861 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, &(*edi));
2863 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2864 if (edi->star.nr == edi->sav.nr)
2868 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2870 /* The last sav.nr indices of the target structure correspond to
2871 * the average structure, which must be projected */
2872 avindex = edi->star.nr - edi->sav.nr;
2874 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2878 rad_project(*edi, xstart, &edi->vecs.radcon);
2881 /* process structure that will serve as origin of expansion circle */
2882 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2885 "ED: Setting center of flooding potential (0 = average structure)\n");
2888 if (edi->sori.nr > 0)
2890 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2892 /* fit this structure to reference structure */
2893 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, &(*edi));
2895 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2896 if (edi->sori.nr == edi->sav.nr)
2900 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2902 /* For the projection, we need the last sav.nr indices of sori */
2903 avindex = edi->sori.nr - edi->sav.nr;
2906 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2907 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2908 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2911 "ED: The ORIGIN structure will define the flooding potential "
2913 /* Set center of flooding potential to the ORIGIN structure */
2914 rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2915 /* We already know that no (moving) reference position was provided,
2916 * therefore we can overwrite refproj[0]*/
2917 copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
2920 else /* No origin structure given */
2922 rad_project(*edi, xstart, &edi->vecs.radacc);
2923 rad_project(*edi, xstart, &edi->vecs.radfix);
2924 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2926 if (edi->flood.bHarmonic)
2929 "ED: A (possibly changing) ref. projection will define the "
2930 "flooding potential center.\n");
2931 for (i = 0; i < edi->flood.vecs.neig; i++)
2933 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
2939 "ED: The AVERAGE structure will define the flooding potential "
2941 /* Set center of flooding potential to the center of the covariance matrix,
2942 * i.e. the average structure, i.e. zero in the projected system */
2943 for (i = 0; i < edi->flood.vecs.neig; i++)
2945 edi->flood.vecs.refproj[i] = 0.0;
2950 /* For convenience, output the center of the flooding potential for the eigenvectors */
2951 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2953 for (i = 0; i < edi->flood.vecs.neig; i++)
2955 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e",
2956 edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2957 if (edi->flood.bHarmonic)
2959 fprintf(stdout, " (adding %11.4e/timestep)",
2960 edi->flood.referenceProjectionSlope[i]);
2962 fprintf(stdout, "\n");
2966 /* set starting projections for linsam */
2967 rad_project(*edi, xstart, &edi->vecs.linacc);
2968 rad_project(*edi, xstart, &edi->vecs.linfix);
2970 /* Prepare for the next edi data set: */
2973 /* Cleaning up on the master node: */
2974 if (!EDstate->bFromCpt)
2981 } /* end of MASTER only section */
2985 /* Broadcast the essential dynamics / flooding data to all nodes */
2986 broadcast_ed_data(cr, ed);
2990 /* In the single-CPU case, point the local atom numbers pointers to the global
2991 * one, so that we can use the same notation in serial and parallel case: */
2992 /* Loop over all ED data sets (usually only one, though) */
2993 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
2995 edi->sref.anrs_loc = edi->sref.anrs;
2996 edi->sav.anrs_loc = edi->sav.anrs;
2997 edi->star.anrs_loc = edi->star.anrs;
2998 edi->sori.anrs_loc = edi->sori.anrs;
2999 /* For the same reason as above, make a dummy c_ind array: */
3000 snew(edi->sav.c_ind, edi->sav.nr);
3001 /* Initialize the array */
3002 for (i = 0; i < edi->sav.nr; i++)
3004 edi->sav.c_ind[i] = i;
3006 /* In the general case we will need a different-sized array for the reference indices: */
3009 snew(edi->sref.c_ind, edi->sref.nr);
3010 for (i = 0; i < edi->sref.nr; i++)
3012 edi->sref.c_ind[i] = i;
3015 /* Point to the very same array in case of other structures: */
3016 edi->star.c_ind = edi->sav.c_ind;
3017 edi->sori.c_ind = edi->sav.c_ind;
3018 /* In the serial case, the local number of atoms is the global one: */
3019 edi->sref.nr_loc = edi->sref.nr;
3020 edi->sav.nr_loc = edi->sav.nr;
3021 edi->star.nr_loc = edi->star.nr;
3022 edi->sori.nr_loc = edi->sori.nr;
3026 /* Allocate space for ED buffer variables */
3027 /* Again, loop over ED data sets */
3028 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
3030 /* Allocate space for ED buffer variables */
3031 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
3032 snew(edi->buf->do_edsam, 1);
3034 /* Space for collective ED buffer variables */
3036 /* Collective positions of atoms with the average indices */
3037 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
3038 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
3039 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
3040 /* Collective positions of atoms with the reference indices */
3043 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
3044 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
3045 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
3048 /* Get memory for flooding forces */
3049 snew(edi->flood.forces_cartesian, edi->sav.nr);
3052 /* Flush the edo file so that the user can check some things
3053 * when the simulation has started */
3063 void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[], rvec v[], const matrix box, gmx_edsam* ed)
3065 int i, edinr, iupdate = 500;
3066 matrix rotmat; /* rotation matrix */
3067 rvec transvec; /* translation vector */
3068 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3069 real dt_1; /* 1/dt */
3070 struct t_do_edsam* buf;
3071 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3072 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3075 /* Check if ED sampling has to be performed */
3076 if (ed->eEDtype == EssentialDynamicsType::None)
3081 dt_1 = 1.0 / ir->delta_t;
3083 /* Loop over all ED groups (usually one) */
3085 for (auto& edi : ed->edpar)
3088 if (bNeedDoEdsam(edi))
3091 buf = edi.buf->do_edsam;
3095 /* initialize radacc radius for slope criterion */
3096 buf->oldrad = calc_radius(edi.vecs.radacc);
3099 /* Copy the positions into buf->xc* arrays and after ED
3100 * feed back corrections to the official positions */
3102 /* Broadcast the ED positions such that every node has all of them
3103 * Every node contributes its local positions xs and stores it in
3104 * the collective buf->xcoll array. Note that for edinr > 1
3105 * xs could already have been modified by an earlier ED */
3107 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll,
3108 PAR(cr) ? buf->bUpdateShifts : TRUE, xs, edi.sav.nr, edi.sav.nr_loc,
3109 edi.sav.anrs_loc, edi.sav.c_ind, edi.sav.x_old, box);
3111 /* Only assembly reference positions if their indices differ from the average ones */
3114 communicate_group_positions(
3115 cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref,
3116 PAR(cr) ? buf->bUpdateShifts : TRUE, xs, edi.sref.nr, edi.sref.nr_loc,
3117 edi.sref.anrs_loc, edi.sref.c_ind, edi.sref.x_old, box);
3120 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3121 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3122 * set bUpdateShifts=TRUE in the parallel case. */
3123 buf->bUpdateShifts = FALSE;
3125 /* Now all nodes have all of the ED positions in edi.sav->xcoll,
3126 * as well as the indices in edi.sav.anrs */
3128 /* Fit the reference indices to the reference structure */
3131 fit_to_reference(buf->xcoll, transvec, rotmat, &edi);
3135 fit_to_reference(buf->xc_ref, transvec, rotmat, &edi);
3138 /* Now apply the translation and rotation to the ED structure */
3139 translate_and_rotate(buf->xcoll, edi.sav.nr, transvec, rotmat);
3141 /* Find out how well we fit to the reference (just for output steps) */
3142 if (do_per_step(step, edi.outfrq) && MASTER(cr))
3146 /* Indices of reference and average structures are identical,
3147 * thus we can calculate the rmsd to SREF using xcoll */
3148 rmsdev = rmsd_from_structure(buf->xcoll, &edi.sref);
3152 /* We have to translate & rotate the reference atoms first */
3153 translate_and_rotate(buf->xc_ref, edi.sref.nr, transvec, rotmat);
3154 rmsdev = rmsd_from_structure(buf->xc_ref, &edi.sref);
3158 /* update radsam references, when required */
3159 if (do_per_step(step, edi.maxedsteps) && step >= edi.presteps)
3161 project(buf->xcoll, &edi);
3162 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3163 rad_project(edi, buf->xcoll, &edi.vecs.radfix);
3164 buf->oldrad = -1.e5;
3167 /* update radacc references, when required */
3168 if (do_per_step(step, iupdate) && step >= edi.presteps)
3170 edi.vecs.radacc.radius = calc_radius(edi.vecs.radacc);
3171 if (edi.vecs.radacc.radius - buf->oldrad < edi.slope)
3173 project(buf->xcoll, &edi);
3174 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3179 buf->oldrad = edi.vecs.radacc.radius;
3183 /* apply the constraints */
3184 if (step >= edi.presteps && ed_constraints(ed->eEDtype, edi))
3186 /* ED constraints should be applied already in the first MD step
3187 * (which is step 0), therefore we pass step+1 to the routine */
3188 ed_apply_constraints(buf->xcoll, &edi, step + 1 - ir->init_step);
3191 /* write to edo, when required */
3192 if (do_per_step(step, edi.outfrq))
3194 project(buf->xcoll, &edi);
3195 if (MASTER(cr) && !bSuppress)
3197 write_edo(edi, ed->edo, rmsdev);
3201 /* Copy back the positions unless monitoring only */
3202 if (ed_constraints(ed->eEDtype, edi))
3204 /* remove fitting */
3205 rmfit(edi.sav.nr, buf->xcoll, transvec, rotmat);
3207 /* Copy the ED corrected positions into the coordinate array */
3208 /* Each node copies its local part. In the serial case, nat_loc is the
3209 * total number of ED atoms */
3210 for (i = 0; i < edi.sav.nr_loc; i++)
3212 /* Unshift local ED coordinate and store in x_unsh */
3213 ed_unshift_single_coord(box, buf->xcoll[edi.sav.c_ind[i]],
3214 buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh);
3216 /* dx is the ED correction to the positions: */
3217 rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx);
3221 /* dv is the ED correction to the velocity: */
3222 svmul(dt_1, dx, dv);
3223 /* apply the velocity correction: */
3224 rvec_inc(v[edi.sav.anrs_loc[i]], dv);
3226 /* Finally apply the position correction due to ED: */
3227 copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]);
3230 } /* END of if ( bNeedDoEdsam(edi) ) */
3232 /* Prepare for the next ED group */
3234 } /* END of loop over ED groups */