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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/commandline/filenm.h"
50 #include "gromacs/domdec/domdec_struct.h"
51 #include "gromacs/fileio/gmxfio.h"
52 #include "gromacs/fileio/xvgr.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/gmxlib/nrnb.h"
55 #include "gromacs/linearalgebra/nrjac.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/utilities.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/math/vectypes.h"
60 #include "gromacs/mdlib/broadcaststructs.h"
61 #include "gromacs/mdlib/constr.h"
62 #include "gromacs/mdlib/groupcoord.h"
63 #include "gromacs/mdlib/stat.h"
64 #include "gromacs/mdlib/update.h"
65 #include "gromacs/mdrunutility/handlerestart.h"
66 #include "gromacs/mdtypes/commrec.h"
67 #include "gromacs/mdtypes/edsamhistory.h"
68 #include "gromacs/mdtypes/inputrec.h"
69 #include "gromacs/mdtypes/md_enums.h"
70 #include "gromacs/mdtypes/observableshistory.h"
71 #include "gromacs/mdtypes/state.h"
72 #include "gromacs/pbcutil/pbc.h"
73 #include "gromacs/topology/mtop_lookup.h"
74 #include "gromacs/topology/topology.h"
75 #include "gromacs/utility/cstringutil.h"
76 #include "gromacs/utility/fatalerror.h"
77 #include "gromacs/utility/gmxassert.h"
78 #include "gromacs/utility/logger.h"
79 #include "gromacs/utility/smalloc.h"
80 #include "gromacs/utility/strconvert.h"
84 /*! \brief Identifies the type of ED: none, normal ED, flooding. */
85 enum class EssentialDynamicsType
87 //! No essential dynamics
89 //! Essential dynamics sampling
91 //! Essential dynamics flooding
95 /*! \brief Identify on which structure to operate. */
96 enum class EssentialDynamicsStructure
98 //! Reference structure
100 //! Average Structure
108 /*! \brief Essential dynamics vector.
109 * TODO: split into working data and input data
110 * NOTE: called eigvec, because traditionally eigenvectors from PCA analysis
111 * were used as essential dynamics vectors, however, vectors used for ED need
112 * not have any special properties
116 //! nr of eigenvectors
118 //! index nrs of eigenvectors
120 //! stepsizes (per eigenvector)
121 real* stpsz = nullptr;
122 //! eigenvector components
123 rvec** vec = nullptr;
124 //! instantaneous x projections
125 real* xproj = nullptr;
126 //! instantaneous f projections
127 real* fproj = nullptr;
128 //! instantaneous radius
130 //! starting or target projections
131 real* refproj = nullptr;
134 /*! \brief Essential dynamics vectors per method implementation.
138 //! only monitored, no constraints
140 //! fixed linear constraints
141 t_eigvec linfix = {};
142 //! acceptance linear constraints
143 t_eigvec linacc = {};
144 //! fixed radial constraints (exp)
145 t_eigvec radfix = {};
146 //! acceptance radial constraints (exp)
147 t_eigvec radacc = {};
148 //! acceptance rad. contraction constr.
149 t_eigvec radcon = {};
152 /*! \brief Essential dynamics flooding parameters and work data.
153 * TODO: split into working data and input parameters
154 * NOTE: The implementation here follows:
155 * O.E. Lange, L.V. Schafer, and H. Grubmuller,
156 * “Flooding in GROMACS: Accelerated barrier crossings in molecular dynamics,”
157 * J. Comp. Chem., 27 1693–1702 (2006)
161 /*! \brief Target destabilisation free energy.
162 * Controls the flooding potential strength.
163 * Linked to the expected speed-up of mean first passage time out of flooded minimum */
165 //! Do not calculate a flooding potential, instead flood with a constant force
166 bool bConstForce = false;
167 //! Relaxation time scale for slowest flooded degree of freedom
169 //! Current estimated destabilisation free energy
171 //! Flooding energy, acting as a proportionality factor for the flooding potential
173 //! Boltzmann constant times temperature, provided by user
175 //! The flooding potential
177 //! Integration time step
179 //! Inital flooding strenth
181 //! Empirical global scaling parameter of the essential dynamics vectors.
183 //! The forces from flooding in atom coordinate space (in contrast to projections onto essential dynamics vectors)
184 rvec* forces_cartesian = nullptr;
185 //! The vectors along which to flood
187 //! Use flooding for harmonic restraint on the eigenvector
188 bool bHarmonic = false;
189 //! The initial reference projection of the flooding vectors. Only with harmonic restraint.
190 real* initialReferenceProjection = nullptr;
191 //! The current reference projection is the initialReferenceProjection + step * slope. Only with harmonic restraint.
192 real* referenceProjectionSlope = nullptr;
197 /* This type is for the average, reference, target, and origin structure */
200 int nr = 0; /* number of atoms this structure contains */
201 int nr_loc = 0; /* number of atoms on local node */
202 int* anrs = nullptr; /* atom index numbers */
203 int* anrs_loc = nullptr; /* local atom index numbers */
204 int nalloc_loc = 0; /* allocation size of anrs_loc */
205 int* c_ind = nullptr; /* at which position of the whole anrs
206 * array is a local atom?, i.e.
207 * c_ind[0...nr_loc-1] gives the atom index
208 * with respect to the collective
209 * anrs[0...nr-1] array */
210 rvec* x = nullptr; /* positions for this structure */
211 rvec* x_old = nullptr; /* Last positions which have the correct PBC
212 representation of the ED group. In
213 combination with keeping track of the
214 shift vectors, the ED group can always
216 real* m = nullptr; /* masses */
217 real mtot = 0.; /* total mass (only used in sref) */
218 real* sqrtm = nullptr; /* sqrt of the masses used for mass-
219 * weighting of analysis (only used in sav) */
225 int nini = 0; /* total Nr of atoms */
226 gmx_bool fitmas = false; /* true if trans fit with cm */
227 gmx_bool pcamas = false; /* true if mass-weighted PCA */
228 int presteps = 0; /* number of steps to run without any
229 * perturbations ... just monitoring */
230 int outfrq = 0; /* freq (in steps) of writing to edo */
231 int maxedsteps = 0; /* max nr of steps per cycle */
233 /* all gmx_edx datasets are copied to all nodes in the parallel case */
234 gmx_edx sref = {}; /* reference positions, to these fitting
236 gmx_bool bRefEqAv = false; /* If true, reference & average indices
237 * are the same. Used for optimization */
238 gmx_edx sav = {}; /* average positions */
239 gmx_edx star = {}; /* target positions */
240 gmx_edx sori = {}; /* origin positions */
242 t_edvecs vecs = {}; /* eigenvectors */
243 real slope = 0; /* minimal slope in acceptance radexp */
245 t_edflood flood = {}; /* parameters especially for flooding */
246 struct t_ed_buffer* buf = nullptr; /* handle to local buffers */
254 EssentialDynamicsType eEDtype = EssentialDynamicsType::None;
255 //! output file pointer
257 std::vector<t_edpar> edpar;
258 gmx_bool bFirst = false;
260 gmx_edsam::~gmx_edsam()
264 /* edo is opened sometimes with xvgropen, sometimes with
265 * gmx_fio_fopen, so we use the least common denominator for
275 rvec old_transvec, older_transvec, transvec_compact;
276 rvec* xcoll; /* Positions from all nodes, this is the
277 collective set we work on.
278 These are the positions of atoms with
279 average structure indices */
280 rvec* xc_ref; /* same but with reference structure indices */
281 ivec* shifts_xcoll; /* Shifts for xcoll */
282 ivec* extra_shifts_xcoll; /* xcoll shift changes since last NS step */
283 ivec* shifts_xc_ref; /* Shifts for xc_ref */
284 ivec* extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
285 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
286 ED shifts for this ED group need to
291 /* definition of ED buffer structure */
294 struct t_fit_to_ref* fit_to_ref;
295 struct t_do_edfit* do_edfit;
296 struct t_do_edsam* do_edsam;
297 struct t_do_radcon* do_radcon;
302 class EssentialDynamics::Impl
305 // TODO: move all fields from gmx_edsam here and remove gmx_edsam
306 gmx_edsam essentialDynamics_;
308 EssentialDynamics::EssentialDynamics() : impl_(new Impl) {}
309 EssentialDynamics::~EssentialDynamics() = default;
311 gmx_edsam* EssentialDynamics::getLegacyED()
313 return &impl_->essentialDynamics_;
317 /* Function declarations */
318 static void fit_to_reference(rvec* xcoll, rvec transvec, matrix rotmat, t_edpar* edi);
319 static void translate_and_rotate(rvec* x, int nat, rvec transvec, matrix rotmat);
320 static real rmsd_from_structure(rvec* x, struct gmx_edx* s);
323 /*! \brief Read in the essential dynamics input file and return its contents.
324 * \param[in] fn the file name of the edi file to be read
325 * \param[in] nr_mdatoms the number of atoms in the simulation
326 * \returns A vector of containing the essentiail dyanmics parameters.
327 * NOTE: edi files may that it may contain several ED data sets from concatenated edi files.
328 * The standard case would be a single ED data set, though. */
329 std::vector<t_edpar> read_edi_file(const char* fn, int nr_mdatoms);
331 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam& ed, edsamhistory_t* EDstate);
332 static void init_edsamstate(const gmx_edsam& ed, edsamhistory_t* EDstate);
333 static void write_edo_legend(gmx_edsam* ed, int nED, const gmx_output_env_t* oenv);
334 /* End function declarations */
336 /* Define formats for the column width in the output file */
337 const char EDcol_sfmt[] = "%17s";
338 const char EDcol_efmt[] = "%17.5e";
339 const char EDcol_ffmt[] = "%17f";
341 /* Define a formats for reading, otherwise cppcheck complains for scanf without width limits */
342 const char max_ev_fmt_d[] = "%7d"; /* Number of eigenvectors. 9,999,999 should be enough */
343 const char max_ev_fmt_lf[] = "%12lf";
344 const char max_ev_fmt_dlf[] = "%7d%12lf";
345 const char max_ev_fmt_dlflflf[] = "%7d%12lf%12lf%12lf";
346 const char max_ev_fmt_lelele[] = "%12le%12le%12le";
350 /*!\brief Determines whether to perform essential dynamics constraints other than flooding.
351 * \param[in] edi the essential dynamics parameters
352 * \returns true if essential dyanmics constraints need to be performed
354 bool bNeedDoEdsam(const t_edpar& edi)
356 return (edi.vecs.mon.neig != 0) || (edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0)
357 || (edi.vecs.radfix.neig != 0) || (edi.vecs.radacc.neig != 0) || (edi.vecs.radcon.neig != 0);
362 /* Multiple ED groups will be labeled with letters instead of numbers
363 * to avoid confusion with eigenvector indices */
364 static char get_EDgroupChar(int nr_edi, int nED)
372 * nr_edi = 2 -> B ...
374 return 'A' + nr_edi - 1;
379 /*! \brief The mass-weighted inner product of two coordinate vectors.
380 * Does not subtract average positions, projection on single eigenvector is returned
381 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
382 * Average position is subtracted in ed_apply_constraints prior to calling projectx
383 * \param[in] edi Essential dynamics parameters
384 * \param[in] xcoll vector of atom coordinates
385 * \param[in] vec vector of coordinates to project onto
386 * \return mass-weighted projection.
388 real projectx(const t_edpar& edi, rvec* xcoll, rvec* vec)
394 for (i = 0; i < edi.sav.nr; i++)
396 proj += edi.sav.sqrtm[i] * iprod(vec[i], xcoll[i]);
401 /*!\brief Project coordinates onto vector after substracting average position.
402 * projection is stored in vec->refproj which is used for radacc, radfix,
403 * radcon and center of flooding potential.
404 * Average positions are first substracted from x, then added back again.
405 * \param[in] edi essential dynamics parameters with average position
406 * \param[in] x Coordinates to be projected
407 * \param[out] vec eigenvector, radius and refproj are overwritten here
409 void rad_project(const t_edpar& edi, rvec* x, t_eigvec* vec)
414 /* Subtract average positions */
415 for (i = 0; i < edi.sav.nr; i++)
417 rvec_dec(x[i], edi.sav.x[i]);
420 for (i = 0; i < vec->neig; i++)
422 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
423 rad += gmx::square((vec->refproj[i] - vec->xproj[i]));
425 vec->radius = sqrt(rad);
427 /* Add average positions */
428 for (i = 0; i < edi.sav.nr; i++)
430 rvec_inc(x[i], edi.sav.x[i]);
434 /*!\brief Projects coordinates onto eigenvectors and stores result in vec->xproj.
435 * Mass-weighting is applied. Subtracts average positions prior to projection and add
436 * them afterwards to retain the unchanged vector.
437 * \param[in] x The coordinates to project to an eigenvector
438 * \param[in,out] vec The eigenvectors
439 * \param[in] edi essential dynamics parameters holding average structure and masses
441 void project_to_eigvectors(rvec* x, t_eigvec* vec, const t_edpar& edi)
448 /* Subtract average positions */
449 for (int i = 0; i < edi.sav.nr; i++)
451 rvec_dec(x[i], edi.sav.x[i]);
454 for (int i = 0; i < vec->neig; i++)
456 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
459 /* Add average positions */
460 for (int i = 0; i < edi.sav.nr; i++)
462 rvec_inc(x[i], edi.sav.x[i]);
467 /* Project vector x onto all edi->vecs (mon, linfix,...) */
468 static void project(rvec* x, /* positions to project */
469 t_edpar* edi) /* edi data set */
471 /* It is not more work to subtract the average position in every
472 * subroutine again, because these routines are rarely used simultaneously */
473 project_to_eigvectors(x, &edi->vecs.mon, *edi);
474 project_to_eigvectors(x, &edi->vecs.linfix, *edi);
475 project_to_eigvectors(x, &edi->vecs.linacc, *edi);
476 project_to_eigvectors(x, &edi->vecs.radfix, *edi);
477 project_to_eigvectors(x, &edi->vecs.radacc, *edi);
478 project_to_eigvectors(x, &edi->vecs.radcon, *edi);
483 /*!\brief Evaluates the distance from reference to current eigenvector projection.
484 * \param[in] vec eigenvector
487 real calc_radius(const t_eigvec& vec)
491 for (int i = 0; i < vec.neig; i++)
493 rad += gmx::square((vec.refproj[i] - vec.xproj[i]));
506 static void do_edfit(int natoms, rvec* xp, rvec* x, matrix R, t_edpar* edi)
508 /* this is a copy of do_fit with some modifications */
509 int c, r, n, j, i, irot;
510 double d[6], xnr, xpc;
515 struct t_do_edfit* loc;
518 if (edi->buf->do_edfit != nullptr)
525 snew(edi->buf->do_edfit, 1);
527 loc = edi->buf->do_edfit;
531 snew(loc->omega, 2 * DIM);
532 snew(loc->om, 2 * DIM);
533 for (i = 0; i < 2 * DIM; i++)
535 snew(loc->omega[i], 2 * DIM);
536 snew(loc->om[i], 2 * DIM);
540 for (i = 0; (i < 6); i++)
543 for (j = 0; (j < 6); j++)
545 loc->omega[i][j] = 0;
550 /* calculate the matrix U */
552 for (n = 0; (n < natoms); n++)
554 for (c = 0; (c < DIM); c++)
557 for (r = 0; (r < DIM); r++)
560 u[c][r] += xnr * xpc;
565 /* construct loc->omega */
566 /* loc->omega is symmetric -> loc->omega==loc->omega' */
567 for (r = 0; (r < 6); r++)
569 for (c = 0; (c <= r); c++)
571 if ((r >= 3) && (c < 3))
573 loc->omega[r][c] = u[r - 3][c];
574 loc->omega[c][r] = u[r - 3][c];
578 loc->omega[r][c] = 0;
579 loc->omega[c][r] = 0;
584 /* determine h and k */
585 jacobi(loc->omega, 6, d, loc->om, &irot);
589 fprintf(stderr, "IROT=0\n");
592 index = 0; /* For the compiler only */
594 for (j = 0; (j < 3); j++)
597 for (i = 0; (i < 6); i++)
606 for (i = 0; (i < 3); i++)
608 vh[j][i] = M_SQRT2 * loc->om[i][index];
609 vk[j][i] = M_SQRT2 * loc->om[i + DIM][index];
614 for (c = 0; (c < 3); c++)
616 for (r = 0; (r < 3); r++)
618 R[c][r] = vk[0][r] * vh[0][c] + vk[1][r] * vh[1][c] + vk[2][r] * vh[2][c];
623 for (c = 0; (c < 3); c++)
625 for (r = 0; (r < 3); r++)
627 R[c][r] = vk[0][r] * vh[0][c] + vk[1][r] * vh[1][c] - vk[2][r] * vh[2][c];
634 static void rmfit(int nat, rvec* xcoll, const rvec transvec, matrix rotmat)
641 * The inverse rotation is described by the transposed rotation matrix */
642 transpose(rotmat, tmat);
643 rotate_x(xcoll, nat, tmat);
645 /* Remove translation */
646 vec[XX] = -transvec[XX];
647 vec[YY] = -transvec[YY];
648 vec[ZZ] = -transvec[ZZ];
649 translate_x(xcoll, nat, vec);
653 /**********************************************************************************
654 ******************** FLOODING ****************************************************
655 **********************************************************************************
657 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
658 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
659 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
661 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
662 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
664 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
665 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
667 do_flood makes a copy of the positions,
668 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
669 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
670 fit. Then do_flood adds these forces to the forcefield-forces
671 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
673 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
674 structure is projected to the system of eigenvectors and then this position in the subspace is used as
675 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
676 i.e. the average structure as given in the make_edi file.
678 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
679 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
680 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
681 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
682 further adaption is applied, Efl will stay constant at zero.
684 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
685 used as spring constants for the harmonic potential.
686 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
687 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
689 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
690 the routine read_edi_file reads all of theses flooding files.
691 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
692 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
693 edi there is no interdependence whatsoever. The forces are added together.
695 To write energies into the .edr file, call the function
696 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
698 get_flood_energies(real Vfl[],int nnames);
701 - 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.
703 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
704 two edsam files from two peptide chains
707 // TODO split this into multiple files
710 /*!\brief Output flooding simulation settings and results to file.
711 * \param[in] edi Essential dynamics input parameters
712 * \param[in] fp output file
713 * \param[in] rmsd rmsd to reference structure
715 void write_edo_flood(const t_edpar& edi, FILE* fp, real rmsd)
717 /* Output how well we fit to the reference structure */
718 fprintf(fp, EDcol_ffmt, rmsd);
720 for (int i = 0; i < edi.flood.vecs.neig; i++)
722 fprintf(fp, EDcol_efmt, edi.flood.vecs.xproj[i]);
724 /* Check whether the reference projection changes with time (this can happen
725 * in case flooding is used as harmonic restraint). If so, output the
726 * current reference projection */
727 if (edi.flood.bHarmonic && edi.flood.referenceProjectionSlope[i] != 0.0)
729 fprintf(fp, EDcol_efmt, edi.flood.vecs.refproj[i]);
732 /* Output Efl if we are doing adaptive flooding */
733 if (0 != edi.flood.tau)
735 fprintf(fp, EDcol_efmt, edi.flood.Efl);
737 fprintf(fp, EDcol_efmt, edi.flood.Vfl);
739 /* Output deltaF if we are doing adaptive flooding */
740 if (0 != edi.flood.tau)
742 fprintf(fp, EDcol_efmt, edi.flood.deltaF);
744 fprintf(fp, EDcol_efmt, edi.flood.vecs.fproj[i]);
750 /* From flood.xproj compute the Vfl(x) at this point */
751 static real flood_energy(t_edpar* edi, int64_t step)
753 /* compute flooding energy Vfl
754 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
755 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
756 it is already computed by make_edi and stored in stpsz[i]
758 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
765 /* Each time this routine is called (i.e. each time step), we add a small
766 * value to the reference projection. This way a harmonic restraint towards
767 * a moving reference is realized. If no value for the additive constant
768 * is provided in the edi file, the reference will not change. */
769 if (edi->flood.bHarmonic)
771 for (i = 0; i < edi->flood.vecs.neig; i++)
773 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i]
774 + step * edi->flood.referenceProjectionSlope[i];
779 /* Compute sum which will be the exponent of the exponential */
780 for (i = 0; i < edi->flood.vecs.neig; i++)
782 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
783 sum += edi->flood.vecs.stpsz[i] * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i])
784 * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i]);
787 /* Compute the Gauss function*/
788 if (edi->flood.bHarmonic)
790 Vfl = -0.5 * edi->flood.Efl * sum; /* minus sign because Efl is negative, if restrain is on. */
794 Vfl = edi->flood.Efl != 0
796 * std::exp(-edi->flood.kT / 2 / edi->flood.Efl / edi->flood.alpha2 * sum)
804 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
805 static void flood_forces(t_edpar* edi)
807 /* compute the forces in the subspace of the flooding eigenvectors
808 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
811 real energy = edi->flood.Vfl;
814 if (edi->flood.bHarmonic)
816 for (i = 0; i < edi->flood.vecs.neig; i++)
818 edi->flood.vecs.fproj[i] = edi->flood.Efl * edi->flood.vecs.stpsz[i]
819 * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i]);
824 for (i = 0; i < edi->flood.vecs.neig; i++)
826 /* if Efl is zero the forces are zero if not use the formula */
827 edi->flood.vecs.fproj[i] =
829 ? edi->flood.kT / edi->flood.Efl / edi->flood.alpha2 * energy
830 * edi->flood.vecs.stpsz[i]
831 * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i])
839 /*!\brief Raise forces from subspace into cartesian space.
840 * This function lifts the forces from the subspace to the cartesian space
841 * all the values not contained in the subspace are assumed to be zero and then
842 * a coordinate transformation from eigenvector to cartesian vectors is performed
843 * The nonexistent values don't have to be set to zero explicitly, they would occur
844 * as zero valued summands, hence we just stop to compute this part of the sum.
845 * For every atom we add all the contributions to this atom from all the different eigenvectors.
846 * NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
847 * field forces_cart prior the computation, but we compute the forces separately
848 * to have them accessible for diagnostics
850 * \param[in] edi Essential dynamics input parameters
851 * \param[out] forces_cart The cartesian forces
854 void flood_blowup(const t_edpar& edi, rvec* forces_cart)
856 const real* forces_sub = edi.flood.vecs.fproj;
857 /* Calculate the cartesian forces for the local atoms */
859 /* Clear forces first */
860 for (int j = 0; j < edi.sav.nr_loc; j++)
862 clear_rvec(forces_cart[j]);
865 /* Now compute atomwise */
866 for (int j = 0; j < edi.sav.nr_loc; j++)
868 /* Compute forces_cart[edi.sav.anrs[j]] */
869 for (int eig = 0; eig < edi.flood.vecs.neig; eig++)
872 /* Force vector is force * eigenvector (compute only atom j) */
873 svmul(forces_sub[eig], edi.flood.vecs.vec[eig][edi.sav.c_ind[j]], addedForce);
874 /* Add this vector to the cartesian forces */
875 rvec_inc(forces_cart[j], addedForce);
883 /* Update the values of Efl, deltaF depending on tau and Vfl */
884 static void update_adaption(t_edpar* edi)
886 /* this function updates the parameter Efl and deltaF according to the rules given in
887 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
890 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau) > 0.00000001)
892 edi->flood.Efl = edi->flood.Efl
893 + edi->flood.dt / edi->flood.tau * (edi->flood.deltaF0 - edi->flood.deltaF);
894 /* check if restrain (inverted flooding) -> don't let EFL become positive */
895 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
900 edi->flood.deltaF = (1 - edi->flood.dt / edi->flood.tau) * edi->flood.deltaF
901 + edi->flood.dt / edi->flood.tau * edi->flood.Vfl;
906 static void do_single_flood(FILE* edo,
913 gmx_bool bNS) /* Are we in a neighbor searching step? */
916 matrix rotmat; /* rotation matrix */
917 matrix tmat; /* inverse rotation */
918 rvec transvec; /* translation vector */
920 struct t_do_edsam* buf;
923 buf = edi->buf->do_edsam;
926 /* Broadcast the positions of the AVERAGE structure such that they are known on
927 * every processor. Each node contributes its local positions x and stores them in
928 * the collective ED array buf->xcoll */
929 communicate_group_positions(cr,
932 buf->extra_shifts_xcoll,
942 /* Only assembly REFERENCE positions if their indices differ from the average ones */
945 communicate_group_positions(cr,
948 buf->extra_shifts_xc_ref,
959 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
960 * We do not need to update the shifts until the next NS step */
961 buf->bUpdateShifts = FALSE;
963 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
964 * as well as the indices in edi->sav.anrs */
966 /* Fit the reference indices to the reference structure */
969 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
973 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
976 /* Now apply the translation and rotation to the ED structure */
977 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
979 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
980 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
982 if (!edi->flood.bConstForce)
984 /* Compute Vfl(x) from flood.xproj */
985 edi->flood.Vfl = flood_energy(edi, step);
987 update_adaption(edi);
989 /* Compute the flooding forces */
993 /* Translate them into cartesian positions */
994 flood_blowup(*edi, edi->flood.forces_cartesian);
996 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
997 /* Each node rotates back its local forces */
998 transpose(rotmat, tmat);
999 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
1001 /* Finally add forces to the main force variable */
1002 for (i = 0; i < edi->sav.nr_loc; i++)
1004 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
1007 /* Output is written by the master process */
1008 if (do_per_step(step, edi->outfrq) && MASTER(cr))
1010 /* Output how well we fit to the reference */
1013 /* Indices of reference and average structures are identical,
1014 * thus we can calculate the rmsd to SREF using xcoll */
1015 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
1019 /* We have to translate & rotate the reference atoms first */
1020 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1021 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1024 write_edo_flood(*edi, edo, rmsdev);
1029 /* Main flooding routine, called from do_force */
1030 extern void do_flood(const t_commrec* cr,
1031 const t_inputrec* ir,
1039 /* Write time to edo, when required. Output the time anyhow since we need
1040 * it in the output file for ED constraints. */
1041 if (MASTER(cr) && do_per_step(step, ed->edpar.begin()->outfrq))
1043 fprintf(ed->edo, "\n%12f", ir->init_t + step * ir->delta_t);
1046 if (ed->eEDtype != EssentialDynamicsType::Flooding)
1051 for (auto& edi : ed->edpar)
1053 /* Call flooding for one matrix */
1054 if (edi.flood.vecs.neig)
1056 do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS);
1062 /* Called by init_edi, configure some flooding related variables and structures,
1063 * print headers to output files */
1064 static void init_flood(t_edpar* edi, gmx_edsam* ed, real dt)
1069 edi->flood.Efl = edi->flood.constEfl;
1073 if (edi->flood.vecs.neig)
1075 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1076 ed->eEDtype = EssentialDynamicsType::Flooding;
1079 "ED: Flooding %d eigenvector%s.\n",
1080 edi->flood.vecs.neig,
1081 edi->flood.vecs.neig > 1 ? "s" : "");
1083 if (edi->flood.bConstForce)
1085 /* We have used stpsz as a vehicle to carry the fproj values for constant
1086 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1087 * in const force flooding, fproj is never changed. */
1088 for (i = 0; i < edi->flood.vecs.neig; i++)
1090 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1093 "ED: applying on eigenvector %d a constant force of %g\n",
1094 edi->flood.vecs.ieig[i],
1095 edi->flood.vecs.fproj[i]);
1103 /*********** Energy book keeping ******/
1104 static void get_flood_enx_names(t_edpar* edi, char** names, int* nnames) /* get header of energies */
1113 srenew(names, count);
1114 sprintf(buf, "Vfl_%d", count);
1115 names[count - 1] = gmx_strdup(buf);
1116 actual = actual->next_edi;
1119 *nnames = count - 1;
1123 static void get_flood_energies(t_edpar* edi, real Vfl[], int nnames)
1125 /*fl has to be big enough to capture nnames-many entries*/
1134 Vfl[count - 1] = actual->flood.Vfl;
1135 actual = actual->next_edi;
1138 if (nnames != count - 1)
1140 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1143 /************* END of FLOODING IMPLEMENTATION ****************************/
1147 /* This function opens the ED input and output files, reads in all datasets it finds
1148 * in the input file, and cross-checks whether the .edi file information is consistent
1149 * with the essential dynamics data found in the checkpoint file (if present).
1150 * gmx make_edi can be used to create an .edi input file.
1152 static std::unique_ptr<gmx::EssentialDynamics> ed_open(int natoms,
1153 ObservablesHistory* oh,
1154 const char* ediFileName,
1155 const char* edoFileName,
1156 const gmx::StartingBehavior startingBehavior,
1157 const gmx_output_env_t* oenv,
1158 const t_commrec* cr)
1160 auto edHandle = std::make_unique<gmx::EssentialDynamics>();
1161 auto ed = edHandle->getLegacyED();
1162 /* We want to perform ED (this switch might later be upgraded to EssentialDynamicsType::Flooding) */
1163 ed->eEDtype = EssentialDynamicsType::EDSampling;
1167 // If we start from a checkpoint file, we already have an edsamHistory struct
1168 if (oh->edsamHistory == nullptr)
1170 oh->edsamHistory = std::make_unique<edsamhistory_t>(edsamhistory_t{});
1172 edsamhistory_t* EDstate = oh->edsamHistory.get();
1174 /* Read the edi input file: */
1175 ed->edpar = read_edi_file(ediFileName, natoms);
1177 /* Make sure the checkpoint was produced in a run using this .edi file */
1178 if (EDstate->bFromCpt)
1180 crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1184 EDstate->nED = ed->edpar.size();
1186 init_edsamstate(*ed, EDstate);
1188 /* The master opens the ED output file */
1189 if (startingBehavior == gmx::StartingBehavior::RestartWithAppending)
1191 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1195 ed->edo = xvgropen(edoFileName,
1196 "Essential dynamics / flooding output",
1198 "RMSDs (nm), projections on EVs (nm), ...",
1201 /* Make a descriptive legend */
1202 write_edo_legend(ed, EDstate->nED, oenv);
1209 /* Broadcasts the structure data */
1210 static void bc_ed_positions(const t_commrec* cr, struct gmx_edx* s, EssentialDynamicsStructure stype)
1212 snew_bc(MASTER(cr), s->anrs, s->nr); /* Index numbers */
1213 snew_bc(MASTER(cr), s->x, s->nr); /* Positions */
1214 nblock_bc(cr->mpi_comm_mygroup, s->nr, s->anrs);
1215 nblock_bc(cr->mpi_comm_mygroup, s->nr, s->x);
1217 /* For the average & reference structures we need an array for the collective indices,
1218 * and we need to broadcast the masses as well */
1219 if (stype == EssentialDynamicsStructure::Average || stype == EssentialDynamicsStructure::Reference)
1221 /* We need these additional variables in the parallel case: */
1222 snew(s->c_ind, s->nr); /* Collective indices */
1223 /* Local atom indices get assigned in dd_make_local_group_indices.
1224 * There, also memory is allocated */
1225 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1226 snew_bc(MASTER(cr), s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1227 nblock_bc(cr->mpi_comm_mygroup, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1230 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1231 if (stype == EssentialDynamicsStructure::Reference)
1233 snew_bc(MASTER(cr), s->m, s->nr);
1234 nblock_bc(cr->mpi_comm_mygroup, s->nr, s->m);
1237 /* For the average structure we might need the masses for mass-weighting */
1238 if (stype == EssentialDynamicsStructure::Average)
1240 snew_bc(MASTER(cr), s->sqrtm, s->nr);
1241 nblock_bc(cr->mpi_comm_mygroup, s->nr, s->sqrtm);
1242 snew_bc(MASTER(cr), s->m, s->nr);
1243 nblock_bc(cr->mpi_comm_mygroup, s->nr, s->m);
1248 /* Broadcasts the eigenvector data */
1249 static void bc_ed_vecs(const t_commrec* cr, t_eigvec* ev, int length)
1253 snew_bc(MASTER(cr), ev->ieig, ev->neig); /* index numbers of eigenvector */
1254 snew_bc(MASTER(cr), ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1255 snew_bc(MASTER(cr), ev->xproj, ev->neig); /* instantaneous x projection */
1256 snew_bc(MASTER(cr), ev->fproj, ev->neig); /* instantaneous f projection */
1257 snew_bc(MASTER(cr), ev->refproj, ev->neig); /* starting or target projection */
1259 nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->ieig);
1260 nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->stpsz);
1261 nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->xproj);
1262 nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->fproj);
1263 nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->refproj);
1265 snew_bc(MASTER(cr), ev->vec, ev->neig); /* Eigenvector components */
1266 for (i = 0; i < ev->neig; i++)
1268 snew_bc(MASTER(cr), ev->vec[i], length);
1269 nblock_bc(cr->mpi_comm_mygroup, length, ev->vec[i]);
1274 /* Broadcasts the ED / flooding data to other nodes
1275 * and allocates memory where needed */
1276 static void broadcast_ed_data(const t_commrec* cr, gmx_edsam* ed)
1278 /* Master lets the other nodes know if its ED only or also flooding */
1279 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr->mpi_comm_mygroup);
1281 int numedis = ed->edpar.size();
1282 /* First let everybody know how many ED data sets to expect */
1283 gmx_bcast(sizeof(numedis), &numedis, cr->mpi_comm_mygroup);
1284 nblock_abc(MASTER(cr), cr->mpi_comm_mygroup, numedis, &(ed->edpar));
1285 /* Now transfer the ED data set(s) */
1286 for (auto& edi : ed->edpar)
1288 /* Broadcast a single ED data set */
1289 block_bc(cr->mpi_comm_mygroup, edi);
1291 /* Broadcast positions */
1292 bc_ed_positions(cr, &(edi.sref), EssentialDynamicsStructure::Reference); /* reference positions (don't broadcast masses) */
1293 bc_ed_positions(cr, &(edi.sav), EssentialDynamicsStructure::Average); /* average positions (do broadcast masses as well) */
1294 bc_ed_positions(cr, &(edi.star), EssentialDynamicsStructure::Target); /* target positions */
1295 bc_ed_positions(cr, &(edi.sori), EssentialDynamicsStructure::Origin); /* origin positions */
1297 /* Broadcast eigenvectors */
1298 bc_ed_vecs(cr, &edi.vecs.mon, edi.sav.nr);
1299 bc_ed_vecs(cr, &edi.vecs.linfix, edi.sav.nr);
1300 bc_ed_vecs(cr, &edi.vecs.linacc, edi.sav.nr);
1301 bc_ed_vecs(cr, &edi.vecs.radfix, edi.sav.nr);
1302 bc_ed_vecs(cr, &edi.vecs.radacc, edi.sav.nr);
1303 bc_ed_vecs(cr, &edi.vecs.radcon, edi.sav.nr);
1304 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1305 bc_ed_vecs(cr, &edi.flood.vecs, edi.sav.nr);
1307 /* For harmonic restraints the reference projections can change with time */
1308 if (edi.flood.bHarmonic)
1310 snew_bc(MASTER(cr), edi.flood.initialReferenceProjection, edi.flood.vecs.neig);
1311 snew_bc(MASTER(cr), edi.flood.referenceProjectionSlope, edi.flood.vecs.neig);
1312 nblock_bc(cr->mpi_comm_mygroup, edi.flood.vecs.neig, edi.flood.initialReferenceProjection);
1313 nblock_bc(cr->mpi_comm_mygroup, edi.flood.vecs.neig, edi.flood.referenceProjectionSlope);
1319 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1320 static void init_edi(const gmx_mtop_t* mtop, t_edpar* edi)
1323 real totalmass = 0.0;
1326 /* NOTE Init_edi is executed on the master process only
1327 * The initialized data sets are then transmitted to the
1328 * other nodes in broadcast_ed_data */
1330 /* evaluate masses (reference structure) */
1331 snew(edi->sref.m, edi->sref.nr);
1333 for (i = 0; i < edi->sref.nr; i++)
1337 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1341 edi->sref.m[i] = 1.0;
1344 /* Check that every m > 0. Bad things will happen otherwise. */
1345 if (edi->sref.m[i] <= 0.0)
1348 "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1349 "For a mass-weighted fit, all reference structure atoms need to have a mass "
1351 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1352 "atoms from the reference structure by creating a proper index group.\n",
1354 edi->sref.anrs[i] + 1,
1358 totalmass += edi->sref.m[i];
1360 edi->sref.mtot = totalmass;
1362 /* Masses m and sqrt(m) for the average structure. Note that m
1363 * is needed if forces have to be evaluated in do_edsam */
1364 snew(edi->sav.sqrtm, edi->sav.nr);
1365 snew(edi->sav.m, edi->sav.nr);
1366 for (i = 0; i < edi->sav.nr; i++)
1368 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1371 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1375 edi->sav.sqrtm[i] = 1.0;
1378 /* Check that every m > 0. Bad things will happen otherwise. */
1379 if (edi->sav.sqrtm[i] <= 0.0)
1382 "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1383 "For ED with mass-weighting, all average structure atoms need to have a mass "
1385 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1386 "atoms from the average structure by creating a proper index group.\n",
1388 edi->sav.anrs[i] + 1,
1393 /* put reference structure in origin */
1394 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1398 translate_x(edi->sref.x, edi->sref.nr, com);
1400 /* Init ED buffer */
1405 static void check(const char* line, const char* label)
1407 if (!strstr(line, label))
1410 "Could not find input parameter %s at expected position in edsam input-file "
1411 "(.edi)\nline read instead is %s",
1418 static int read_checked_edint(FILE* file, const char* label)
1420 char line[STRLEN + 1];
1423 fgets2(line, STRLEN, file);
1425 fgets2(line, STRLEN, file);
1426 sscanf(line, max_ev_fmt_d, &idum);
1430 static bool read_checked_edbool(FILE* file, const char* label)
1432 return read_checked_edint(file, label) != 0;
1435 static int read_edint(FILE* file, bool* bEOF)
1437 char line[STRLEN + 1];
1442 eof = fgets2(line, STRLEN, file);
1448 eof = fgets2(line, STRLEN, file);
1454 sscanf(line, max_ev_fmt_d, &idum);
1460 static real read_checked_edreal(FILE* file, const char* label)
1462 char line[STRLEN + 1];
1466 fgets2(line, STRLEN, file);
1468 fgets2(line, STRLEN, file);
1469 sscanf(line, max_ev_fmt_lf, &rdum);
1470 return static_cast<real>(rdum); /* always read as double and convert to single */
1474 static void read_edx(FILE* file, int number, int* anrs, rvec* x)
1477 char line[STRLEN + 1];
1481 for (i = 0; i < number; i++)
1483 fgets2(line, STRLEN, file);
1484 sscanf(line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1485 anrs[i]--; /* we are reading FORTRAN indices */
1486 for (j = 0; j < 3; j++)
1488 x[i][j] = d[j]; /* always read as double and convert to single */
1495 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1496 * \param[in] in the file to read from
1497 * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1498 * \param[out] vec the eigen-vectors
1499 * \param[in] nEig the number of eigenvectors
1501 void scan_edvec(FILE* in, int numAtoms, rvec*** vec, int nEig)
1504 for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1506 snew((*vec)[iEigenvector], numAtoms);
1507 for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1509 char line[STRLEN + 1];
1511 fgets2(line, STRLEN, in);
1512 sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1513 (*vec)[iEigenvector][iAtom][XX] = x;
1514 (*vec)[iEigenvector][iAtom][YY] = y;
1515 (*vec)[iEigenvector][iAtom][ZZ] = z;
1519 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1520 * \param[in] in input file
1521 * \param[in] nrAtoms number of atoms/coordinates
1522 * \param[out] tvec the eigenvector
1524 void read_edvec(FILE* in, int nrAtoms, t_eigvec* tvec)
1526 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1527 if (tvec->neig <= 0)
1532 snew(tvec->ieig, tvec->neig);
1533 snew(tvec->stpsz, tvec->neig);
1534 for (int i = 0; i < tvec->neig; i++)
1536 char line[STRLEN + 1];
1537 fgets2(line, STRLEN, in);
1538 int numEigenVectors;
1540 const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1543 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1545 tvec->ieig[i] = numEigenVectors;
1546 tvec->stpsz[i] = stepSize;
1547 } /* end of loop over eigenvectors */
1549 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1551 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1553 * Eigenvector and its intial reference and slope are stored on the
1554 * same line in the .edi format. To avoid re-winding the .edi file,
1555 * the reference eigen-vector and reference data are read in one go.
1557 * \param[in] in input file
1558 * \param[in] nrAtoms number of atoms/coordinates
1559 * \param[out] tvec the eigenvector
1560 * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1561 * \param[out] referenceProjectionSlope The slope of the reference projections.
1563 bool readEdVecWithReferenceProjection(FILE* in,
1566 real** initialReferenceProjection,
1567 real** referenceProjectionSlope)
1569 bool providesReference = false;
1570 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1571 if (tvec->neig <= 0)
1576 snew(tvec->ieig, tvec->neig);
1577 snew(tvec->stpsz, tvec->neig);
1578 snew(*initialReferenceProjection, tvec->neig);
1579 snew(*referenceProjectionSlope, tvec->neig);
1580 for (int i = 0; i < tvec->neig; i++)
1582 char line[STRLEN + 1];
1583 fgets2(line, STRLEN, in);
1584 int numEigenVectors;
1585 double stepSize = 0.;
1586 double referenceProjection = 0.;
1587 double referenceSlope = 0.;
1588 const int nscan = sscanf(
1589 line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize, &referenceProjection, &referenceSlope);
1590 /* Zero out values which were not scanned */
1594 /* Every 4 values read, including reference position */
1595 providesReference = true;
1598 /* A reference position is provided */
1599 providesReference = true;
1600 /* No value for slope, set to 0 */
1601 referenceSlope = 0.0;
1604 /* No values for reference projection and slope, set to 0 */
1605 referenceProjection = 0.0;
1606 referenceSlope = 0.0;
1610 "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> "
1611 "<refproj> <refproj-slope>\n",
1614 (*initialReferenceProjection)[i] = referenceProjection;
1615 (*referenceProjectionSlope)[i] = referenceSlope;
1617 tvec->ieig[i] = numEigenVectors;
1618 tvec->stpsz[i] = stepSize;
1619 } /* end of loop over eigenvectors */
1621 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1622 return providesReference;
1625 /*!\brief Allocate working memory for the eigenvectors.
1626 * \param[in,out] tvec eigenvector for which memory will be allocated
1628 void setup_edvec(t_eigvec* tvec)
1630 snew(tvec->xproj, tvec->neig);
1631 snew(tvec->fproj, tvec->neig);
1632 snew(tvec->refproj, tvec->neig);
1637 /* Check if the same atom indices are used for reference and average positions */
1638 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1643 /* If the number of atoms differs between the two structures,
1644 * they cannot be identical */
1645 if (sref.nr != sav.nr)
1650 /* Now that we know that both stuctures have the same number of atoms,
1651 * check if also the indices are identical */
1652 for (i = 0; i < sav.nr; i++)
1654 if (sref.anrs[i] != sav.anrs[i])
1660 "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1668 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1669 constexpr int c_oldestUnsupportedMagicNumber = 666;
1670 //! Oldest (lowest) magic number indicating supported essential dynamics input
1671 constexpr int c_oldestSupportedMagicNumber = 669;
1672 //! Latest (highest) magic number indicating supported essential dynamics input
1673 constexpr int c_latestSupportedMagicNumber = 670;
1675 /*!\brief Set up essential dynamics work parameters.
1676 * \param[in] edi Essential dynamics working structure.
1678 void setup_edi(t_edpar* edi)
1680 snew(edi->sref.x_old, edi->sref.nr);
1681 edi->sref.sqrtm = nullptr;
1682 snew(edi->sav.x_old, edi->sav.nr);
1683 if (edi->star.nr > 0)
1685 edi->star.sqrtm = nullptr;
1687 if (edi->sori.nr > 0)
1689 edi->sori.sqrtm = nullptr;
1691 setup_edvec(&edi->vecs.linacc);
1692 setup_edvec(&edi->vecs.mon);
1693 setup_edvec(&edi->vecs.linfix);
1694 setup_edvec(&edi->vecs.linacc);
1695 setup_edvec(&edi->vecs.radfix);
1696 setup_edvec(&edi->vecs.radacc);
1697 setup_edvec(&edi->vecs.radcon);
1698 setup_edvec(&edi->flood.vecs);
1701 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1702 * \param[in] magicNumber indicates the essential dynamics input file version
1703 * \returns true if CONST_FORCE_FLOODING is supported
1705 bool ediFileHasConstForceFlooding(int magicNumber)
1707 return magicNumber > c_oldestSupportedMagicNumber;
1710 /*!\brief Reports reading success of the essential dynamics input file magic number.
1711 * \param[in] in pointer to input file
1712 * \param[out] magicNumber determines the edi file version
1713 * \returns true if setting file version from input file was successful.
1715 bool ediFileHasValidDataSet(FILE* in, int* magicNumber)
1717 /* the edi file is not free format, so expect problems if the input is corrupt. */
1718 bool reachedEndOfFile = true;
1719 /* check the magic number */
1720 *magicNumber = read_edint(in, &reachedEndOfFile);
1721 /* Check whether we have reached the end of the input file */
1722 return !reachedEndOfFile;
1725 /*!\brief Trigger fatal error if magic number is unsupported.
1726 * \param[in] magicNumber A magic number read from the edi file
1727 * \param[in] fn name of input file for error reporting
1729 void exitWhenMagicNumberIsInvalid(int magicNumber, const char* fn)
1732 if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1734 if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1737 "Wrong magic number: Use newest version of make_edi to produce edi file");
1741 gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1746 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1748 * \param[in] in input file
1749 * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1750 * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1751 * \param[in] fn the file name of the input file for error reporting
1752 * \returns edi essential dynamics parameters
1754 t_edpar read_edi(FILE* in, int nr_mdatoms, bool hasConstForceFlooding, const char* fn)
1758 /* check the number of atoms */
1759 edi.nini = read_edint(in, &bEOF);
1760 if (edi.nini != nr_mdatoms)
1762 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi.nini, nr_mdatoms);
1765 /* Done checking. For the rest we blindly trust the input */
1766 edi.fitmas = read_checked_edbool(in, "FITMAS");
1767 edi.pcamas = read_checked_edbool(in, "ANALYSIS_MAS");
1768 edi.outfrq = read_checked_edint(in, "OUTFRQ");
1769 edi.maxedsteps = read_checked_edint(in, "MAXLEN");
1770 edi.slope = read_checked_edreal(in, "SLOPECRIT");
1772 edi.presteps = read_checked_edint(in, "PRESTEPS");
1773 edi.flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1774 edi.flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1775 edi.flood.tau = read_checked_edreal(in, "TAU");
1776 edi.flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1777 edi.flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1778 edi.flood.kT = read_checked_edreal(in, "KT");
1779 edi.flood.bHarmonic = read_checked_edbool(in, "HARMONIC");
1780 if (hasConstForceFlooding)
1782 edi.flood.bConstForce = read_checked_edbool(in, "CONST_FORCE_FLOODING");
1786 edi.flood.bConstForce = FALSE;
1788 edi.sref.nr = read_checked_edint(in, "NREF");
1790 /* allocate space for reference positions and read them */
1791 snew(edi.sref.anrs, edi.sref.nr);
1792 snew(edi.sref.x, edi.sref.nr);
1793 read_edx(in, edi.sref.nr, edi.sref.anrs, edi.sref.x);
1795 /* average positions. they define which atoms will be used for ED sampling */
1796 edi.sav.nr = read_checked_edint(in, "NAV");
1797 snew(edi.sav.anrs, edi.sav.nr);
1798 snew(edi.sav.x, edi.sav.nr);
1799 read_edx(in, edi.sav.nr, edi.sav.anrs, edi.sav.x);
1801 /* Check if the same atom indices are used for reference and average positions */
1802 edi.bRefEqAv = check_if_same(edi.sref, edi.sav);
1806 read_edvec(in, edi.sav.nr, &edi.vecs.mon);
1807 read_edvec(in, edi.sav.nr, &edi.vecs.linfix);
1808 read_edvec(in, edi.sav.nr, &edi.vecs.linacc);
1809 read_edvec(in, edi.sav.nr, &edi.vecs.radfix);
1810 read_edvec(in, edi.sav.nr, &edi.vecs.radacc);
1811 read_edvec(in, edi.sav.nr, &edi.vecs.radcon);
1813 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1814 bool bHaveReference = false;
1815 if (edi.flood.bHarmonic)
1817 bHaveReference = readEdVecWithReferenceProjection(in,
1820 &edi.flood.initialReferenceProjection,
1821 &edi.flood.referenceProjectionSlope);
1825 read_edvec(in, edi.sav.nr, &edi.flood.vecs);
1828 /* target positions */
1829 edi.star.nr = read_edint(in, &bEOF);
1830 if (edi.star.nr > 0)
1832 snew(edi.star.anrs, edi.star.nr);
1833 snew(edi.star.x, edi.star.nr);
1834 read_edx(in, edi.star.nr, edi.star.anrs, edi.star.x);
1837 /* positions defining origin of expansion circle */
1838 edi.sori.nr = read_edint(in, &bEOF);
1839 if (edi.sori.nr > 0)
1843 /* Both an -ori structure and a at least one manual reference point have been
1844 * specified. That's ambiguous and probably not intentional. */
1846 "ED: An origin structure has been provided and a at least one (moving) "
1848 " point was manually specified in the edi file. That is ambiguous. "
1851 snew(edi.sori.anrs, edi.sori.nr);
1852 snew(edi.sori.x, edi.sori.nr);
1853 read_edx(in, edi.sori.nr, edi.sori.anrs, edi.sori.x);
1858 std::vector<t_edpar> read_edi_file(const char* fn, int nr_mdatoms)
1860 std::vector<t_edpar> essentialDynamicsParameters;
1862 /* This routine is executed on the master only */
1864 /* Open the .edi parameter input file */
1865 in = gmx_fio_fopen(fn, "r");
1866 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1868 /* Now read a sequence of ED input parameter sets from the edi file */
1869 int ediFileMagicNumber;
1870 while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1872 exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1873 bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1874 auto edi = read_edi(in, nr_mdatoms, hasConstForceFlooding, fn);
1876 essentialDynamicsParameters.emplace_back(edi);
1878 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1879 /* We need to allocate space for the data: */
1881 if (essentialDynamicsParameters.empty())
1883 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1887 "ED: Found %zu ED group%s.\n",
1888 essentialDynamicsParameters.size(),
1889 essentialDynamicsParameters.size() > 1 ? "s" : "");
1891 /* Close the .edi file again */
1894 return essentialDynamicsParameters;
1896 } // anonymous namespace
1901 rvec* xcopy; /* Working copy of the positions in fit_to_reference */
1904 /* Fit the current positions to the reference positions
1905 * Do not actually do the fit, just return rotation and translation.
1906 * Note that the COM of the reference structure was already put into
1907 * the origin by init_edi. */
1908 static void fit_to_reference(rvec* xcoll, /* The positions to be fitted */
1909 rvec transvec, /* The translation vector */
1910 matrix rotmat, /* The rotation matrix */
1911 t_edpar* edi) /* Just needed for do_edfit */
1913 rvec com; /* center of mass */
1915 struct t_fit_to_ref* loc;
1918 /* Allocate memory the first time this routine is called for each edi group */
1919 if (nullptr == edi->buf->fit_to_ref)
1921 snew(edi->buf->fit_to_ref, 1);
1922 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1924 loc = edi->buf->fit_to_ref;
1926 /* We do not touch the original positions but work on a copy. */
1927 for (i = 0; i < edi->sref.nr; i++)
1929 copy_rvec(xcoll[i], loc->xcopy[i]);
1932 /* Calculate the center of mass */
1933 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1935 transvec[XX] = -com[XX];
1936 transvec[YY] = -com[YY];
1937 transvec[ZZ] = -com[ZZ];
1939 /* Subtract the center of mass from the copy */
1940 translate_x(loc->xcopy, edi->sref.nr, transvec);
1942 /* Determine the rotation matrix */
1943 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1947 static void translate_and_rotate(rvec* x, /* The positions to be translated and rotated */
1948 int nat, /* How many positions are there? */
1949 rvec transvec, /* The translation vector */
1950 matrix rotmat) /* The rotation matrix */
1953 translate_x(x, nat, transvec);
1956 rotate_x(x, nat, rotmat);
1960 /* Gets the rms deviation of the positions to the structure s */
1961 /* fit_to_structure has to be called before calling this routine! */
1962 static real rmsd_from_structure(rvec* x, /* The positions under consideration */
1963 struct gmx_edx* s) /* The structure from which the rmsd shall be computed */
1969 for (i = 0; i < s->nr; i++)
1971 rmsd += distance2(s->x[i], x[i]);
1974 rmsd /= static_cast<real>(s->nr);
1981 void dd_make_local_ed_indices(gmx_domdec_t* dd, struct gmx_edsam* ed)
1983 if (ed->eEDtype != EssentialDynamicsType::None)
1985 /* Loop over ED groups */
1987 for (auto& edi : ed->edpar)
1989 /* Local atoms of the reference structure (for fitting), need only be assembled
1990 * if their indices differ from the average ones */
1993 dd_make_local_group_indices(dd->ga2la,
1998 &edi.sref.nalloc_loc,
2002 /* Local atoms of the average structure (on these ED will be performed) */
2003 dd_make_local_group_indices(dd->ga2la,
2008 &edi.sav.nalloc_loc,
2011 /* Indicate that the ED shift vectors for this structure need to be updated
2012 * at the next call to communicate_group_positions, since obviously we are in a NS step */
2013 edi.buf->do_edsam->bUpdateShifts = TRUE;
2019 static inline void ed_unshift_single_coord(const matrix box, const rvec x, const ivec is, rvec xu)
2030 xu[XX] = x[XX] - tx * box[XX][XX] - ty * box[YY][XX] - tz * box[ZZ][XX];
2031 xu[YY] = x[YY] - ty * box[YY][YY] - tz * box[ZZ][YY];
2032 xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
2036 xu[XX] = x[XX] - tx * box[XX][XX];
2037 xu[YY] = x[YY] - ty * box[YY][YY];
2038 xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
2044 /*!\brief Apply fixed linear constraints to essential dynamics variable.
2045 * \param[in,out] xcoll The collected coordinates.
2046 * \param[in] edi the essential dynamics parameters
2047 * \param[in] step the current simulation step
2049 void do_linfix(rvec* xcoll, const t_edpar& edi, int64_t step)
2051 /* loop over linfix vectors */
2052 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2054 /* calculate the projection */
2055 real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
2057 /* calculate the correction */
2058 real preFactor = edi.vecs.linfix.refproj[i] + step * edi.vecs.linfix.stpsz[i] - proj;
2060 /* apply the correction */
2061 preFactor /= edi.sav.sqrtm[i];
2062 for (int j = 0; j < edi.sav.nr; j++)
2064 rvec differenceVector;
2065 svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
2066 rvec_inc(xcoll[j], differenceVector);
2071 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
2072 * \param[in,out] xcoll The collected coordinates.
2073 * \param[in] edi the essential dynamics parameters
2075 void do_linacc(rvec* xcoll, t_edpar* edi)
2077 /* loop over linacc vectors */
2078 for (int i = 0; i < edi->vecs.linacc.neig; i++)
2080 /* calculate the projection */
2081 real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
2083 /* calculate the correction */
2084 real preFactor = 0.0;
2085 if (edi->vecs.linacc.stpsz[i] > 0.0)
2087 if ((proj - edi->vecs.linacc.refproj[i]) < 0.0)
2089 preFactor = edi->vecs.linacc.refproj[i] - proj;
2092 if (edi->vecs.linacc.stpsz[i] < 0.0)
2094 if ((proj - edi->vecs.linacc.refproj[i]) > 0.0)
2096 preFactor = edi->vecs.linacc.refproj[i] - proj;
2100 /* apply the correction */
2101 preFactor /= edi->sav.sqrtm[i];
2102 for (int j = 0; j < edi->sav.nr; j++)
2104 rvec differenceVector;
2105 svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2106 rvec_inc(xcoll[j], differenceVector);
2109 /* new positions will act as reference */
2110 edi->vecs.linacc.refproj[i] = proj + preFactor;
2116 static void do_radfix(rvec* xcoll, t_edpar* edi)
2119 real *proj, rad = 0.0, ratio;
2123 if (edi->vecs.radfix.neig == 0)
2128 snew(proj, edi->vecs.radfix.neig);
2130 /* loop over radfix vectors */
2131 for (i = 0; i < edi->vecs.radfix.neig; i++)
2133 /* calculate the projections, radius */
2134 proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2135 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2139 ratio = (edi->vecs.radfix.stpsz[0] + edi->vecs.radfix.radius) / rad - 1.0;
2140 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2142 /* loop over radfix vectors */
2143 for (i = 0; i < edi->vecs.radfix.neig; i++)
2145 proj[i] -= edi->vecs.radfix.refproj[i];
2147 /* apply the correction */
2148 proj[i] /= edi->sav.sqrtm[i];
2150 for (j = 0; j < edi->sav.nr; j++)
2152 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2153 rvec_inc(xcoll[j], vec_dum);
2161 static void do_radacc(rvec* xcoll, t_edpar* edi)
2164 real *proj, rad = 0.0, ratio = 0.0;
2168 if (edi->vecs.radacc.neig == 0)
2173 snew(proj, edi->vecs.radacc.neig);
2175 /* loop over radacc vectors */
2176 for (i = 0; i < edi->vecs.radacc.neig; i++)
2178 /* calculate the projections, radius */
2179 proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2180 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2184 /* only correct when radius decreased */
2185 if (rad < edi->vecs.radacc.radius)
2187 ratio = edi->vecs.radacc.radius / rad - 1.0;
2191 edi->vecs.radacc.radius = rad;
2194 /* loop over radacc vectors */
2195 for (i = 0; i < edi->vecs.radacc.neig; i++)
2197 proj[i] -= edi->vecs.radacc.refproj[i];
2199 /* apply the correction */
2200 proj[i] /= edi->sav.sqrtm[i];
2202 for (j = 0; j < edi->sav.nr; j++)
2204 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2205 rvec_inc(xcoll[j], vec_dum);
2217 static void do_radcon(rvec* xcoll, t_edpar* edi)
2220 real rad = 0.0, ratio = 0.0;
2221 struct t_do_radcon* loc;
2226 if (edi->buf->do_radcon != nullptr)
2233 snew(edi->buf->do_radcon, 1);
2235 loc = edi->buf->do_radcon;
2237 if (edi->vecs.radcon.neig == 0)
2244 snew(loc->proj, edi->vecs.radcon.neig);
2247 /* loop over radcon vectors */
2248 for (i = 0; i < edi->vecs.radcon.neig; i++)
2250 /* calculate the projections, radius */
2251 loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2252 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2255 /* only correct when radius increased */
2256 if (rad > edi->vecs.radcon.radius)
2258 ratio = edi->vecs.radcon.radius / rad - 1.0;
2260 /* loop over radcon vectors */
2261 for (i = 0; i < edi->vecs.radcon.neig; i++)
2263 /* apply the correction */
2264 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2265 loc->proj[i] /= edi->sav.sqrtm[i];
2266 loc->proj[i] *= ratio;
2268 for (j = 0; j < edi->sav.nr; j++)
2270 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2271 rvec_inc(xcoll[j], vec_dum);
2277 edi->vecs.radcon.radius = rad;
2282 static void ed_apply_constraints(rvec* xcoll, t_edpar* edi, int64_t step)
2287 /* subtract the average positions */
2288 for (i = 0; i < edi->sav.nr; i++)
2290 rvec_dec(xcoll[i], edi->sav.x[i]);
2293 /* apply the constraints */
2296 do_linfix(xcoll, *edi, step);
2298 do_linacc(xcoll, edi);
2301 do_radfix(xcoll, edi);
2303 do_radacc(xcoll, edi);
2304 do_radcon(xcoll, edi);
2306 /* add back the average positions */
2307 for (i = 0; i < edi->sav.nr; i++)
2309 rvec_inc(xcoll[i], edi->sav.x[i]);
2316 /*!\brief Write out the projections onto the eigenvectors.
2317 * The order of output corresponds to ed_output_legend().
2318 * \param[in] edi The essential dyanmics parameters
2319 * \param[in] fp The output file
2320 * \param[in] rmsd the rmsd to the reference structure
2322 void write_edo(const t_edpar& edi, FILE* fp, real rmsd)
2324 /* Output how well we fit to the reference structure */
2325 fprintf(fp, EDcol_ffmt, rmsd);
2327 for (int i = 0; i < edi.vecs.mon.neig; i++)
2329 fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2332 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2334 fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2337 for (int i = 0; i < edi.vecs.linacc.neig; i++)
2339 fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2342 for (int i = 0; i < edi.vecs.radfix.neig; i++)
2344 fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2346 if (edi.vecs.radfix.neig)
2348 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2351 for (int i = 0; i < edi.vecs.radacc.neig; i++)
2353 fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2355 if (edi.vecs.radacc.neig)
2357 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2360 for (int i = 0; i < edi.vecs.radcon.neig; i++)
2362 fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2364 if (edi.vecs.radcon.neig)
2366 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2372 /* Returns if any constraints are switched on */
2373 static bool ed_constraints(EssentialDynamicsType edtype, const t_edpar& edi)
2375 if (edtype == EssentialDynamicsType::EDSampling || edtype == EssentialDynamicsType::Flooding)
2377 return ((edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0) || (edi.vecs.radfix.neig != 0)
2378 || (edi.vecs.radacc.neig != 0) || (edi.vecs.radcon.neig != 0));
2384 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for
2385 * flooding/ umbrella sampling simulations. */
2386 static void copyEvecReference(t_eigvec* floodvecs, real* initialReferenceProjection)
2391 if (nullptr == initialReferenceProjection)
2393 snew(initialReferenceProjection, floodvecs->neig);
2396 for (i = 0; i < floodvecs->neig; i++)
2398 initialReferenceProjection[i] = floodvecs->refproj[i];
2403 /* Call on MASTER only. Check whether the essential dynamics / flooding
2404 * groups of the checkpoint file are consistent with the provided .edi file. */
2405 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam& ed, edsamhistory_t* EDstate)
2407 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2410 "Essential dynamics and flooding can only be switched on (or off) at the\n"
2411 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2412 "it must also continue with/without ED constraints when checkpointing.\n"
2413 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2414 "from without a checkpoint.\n");
2417 for (size_t edinum = 0; edinum < ed.edpar.size(); ++edinum)
2419 /* Check number of atoms in the reference and average structures */
2420 if (EDstate->nref[edinum] != ed.edpar[edinum].sref.nr)
2423 "The number of reference structure atoms in ED group %c is\n"
2424 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2425 get_EDgroupChar(edinum + 1, 0),
2426 EDstate->nref[edinum],
2427 ed.edpar[edinum].sref.nr);
2429 if (EDstate->nav[edinum] != ed.edpar[edinum].sav.nr)
2432 "The number of average structure atoms in ED group %c is\n"
2433 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2434 get_EDgroupChar(edinum + 1, 0),
2435 EDstate->nav[edinum],
2436 ed.edpar[edinum].sav.nr);
2440 if (gmx::ssize(ed.edpar) != EDstate->nED)
2443 "The number of essential dynamics / flooding groups is not consistent.\n"
2444 "There are %d ED groups in the .cpt file, but %zu in the .edi file!\n"
2445 "Are you sure this is the correct .edi file?\n",
2452 /* The edsamstate struct stores the information we need to make the ED group
2453 * whole again after restarts from a checkpoint file. Here we do the following:
2454 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2455 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2456 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2457 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2458 * all ED structures at the last time step. */
2459 static void init_edsamstate(const gmx_edsam& ed, edsamhistory_t* EDstate)
2461 snew(EDstate->old_sref_p, EDstate->nED);
2462 snew(EDstate->old_sav_p, EDstate->nED);
2464 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2465 if (!EDstate->bFromCpt)
2467 snew(EDstate->nref, EDstate->nED);
2468 snew(EDstate->nav, EDstate->nED);
2471 /* Loop over all ED/flooding data sets (usually only one, though) */
2472 for (int nr_edi = 0; nr_edi < EDstate->nED; nr_edi++)
2474 const auto& edi = ed.edpar[nr_edi];
2475 /* We always need the last reference and average positions such that
2476 * in the next time step we can make the ED group whole again
2477 * if the atoms do not have the correct PBC representation */
2478 if (EDstate->bFromCpt)
2480 /* Copy the last whole positions of reference and average group from .cpt */
2481 for (int i = 0; i < edi.sref.nr; i++)
2483 copy_rvec(EDstate->old_sref[nr_edi][i], edi.sref.x_old[i]);
2485 for (int i = 0; i < edi.sav.nr; i++)
2487 copy_rvec(EDstate->old_sav[nr_edi][i], edi.sav.x_old[i]);
2492 EDstate->nref[nr_edi] = edi.sref.nr;
2493 EDstate->nav[nr_edi] = edi.sav.nr;
2496 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2497 EDstate->old_sref_p[nr_edi] = edi.sref.x_old;
2498 EDstate->old_sav_p[nr_edi] = edi.sav.x_old;
2503 /* Adds 'buf' to 'str' */
2504 static void add_to_string(char** str, const char* buf)
2509 len = strlen(*str) + strlen(buf) + 1;
2515 static void add_to_string_aligned(char** str, const char* buf)
2517 char buf_aligned[STRLEN];
2519 sprintf(buf_aligned, EDcol_sfmt, buf);
2520 add_to_string(str, buf_aligned);
2524 static void nice_legend(const char*** setname,
2531 auto tmp = gmx::formatString("%c %s", EDgroupchar, value);
2532 add_to_string_aligned(LegendStr, tmp.c_str());
2533 tmp += gmx::formatString(" (%s)", unit);
2534 (*setname)[*nsets] = gmx_strdup(tmp.c_str());
2539 static void nice_legend_evec(const char*** setname,
2550 for (i = 0; i < evec->neig; i++)
2552 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2553 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2558 /* Makes a legend for the xvg output file. Call on MASTER only! */
2559 static void write_edo_legend(gmx_edsam* ed, int nED, const gmx_output_env_t* oenv)
2562 int nr_edi, nsets, n_flood, n_edsam;
2563 const char** setname;
2565 char* LegendStr = nullptr;
2568 auto edi = ed->edpar.begin();
2570 fprintf(ed->edo, "# Output will be written every %d step%s\n", edi->outfrq, edi->outfrq != 1 ? "s" : "");
2572 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2574 fprintf(ed->edo, "#\n");
2576 "# Summary of applied con/restraints for the ED group %c\n",
2577 get_EDgroupChar(nr_edi, nED));
2578 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2579 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2581 "# LINFIX : %d vec%s\n",
2582 edi->vecs.linfix.neig,
2583 edi->vecs.linfix.neig != 1 ? "s" : "");
2585 "# LINACC : %d vec%s\n",
2586 edi->vecs.linacc.neig,
2587 edi->vecs.linacc.neig != 1 ? "s" : "");
2589 "# RADFIX : %d vec%s\n",
2590 edi->vecs.radfix.neig,
2591 edi->vecs.radfix.neig != 1 ? "s" : "");
2593 "# RADACC : %d vec%s\n",
2594 edi->vecs.radacc.neig,
2595 edi->vecs.radacc.neig != 1 ? "s" : "");
2597 "# RADCON : %d vec%s\n",
2598 edi->vecs.radcon.neig,
2599 edi->vecs.radcon.neig != 1 ? "s" : "");
2600 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2602 if (edi->flood.vecs.neig)
2604 /* If in any of the groups we find a flooding vector, flooding is turned on */
2605 ed->eEDtype = EssentialDynamicsType::Flooding;
2607 /* Print what flavor of flooding we will do */
2608 if (0 == edi->flood.tau) /* constant flooding strength */
2610 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2611 if (edi->flood.bHarmonic)
2613 fprintf(ed->edo, ", harmonic");
2616 else /* adaptive flooding */
2618 fprintf(ed->edo, ", adaptive");
2621 fprintf(ed->edo, "\n");
2626 /* Print a nice legend */
2628 LegendStr[0] = '\0';
2629 sprintf(buf, "# %6s", "time");
2630 add_to_string(&LegendStr, buf);
2632 /* Calculate the maximum number of columns we could end up with */
2633 edi = ed->edpar.begin();
2635 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2637 nsets += 5 + edi->vecs.mon.neig + edi->vecs.linfix.neig + edi->vecs.linacc.neig
2638 + edi->vecs.radfix.neig + edi->vecs.radacc.neig + edi->vecs.radcon.neig
2639 + 6 * edi->flood.vecs.neig;
2642 snew(setname, nsets);
2644 /* In the mdrun time step in a first function call (do_flood()) the flooding
2645 * forces are calculated and in a second function call (do_edsam()) the
2646 * ED constraints. To get a corresponding legend, we need to loop twice
2647 * over the edi groups and output first the flooding, then the ED part */
2649 /* The flooding-related legend entries, if flooding is done */
2651 if (EssentialDynamicsType::Flooding == ed->eEDtype)
2653 edi = ed->edpar.begin();
2654 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2656 /* Always write out the projection on the flooding EVs. Of course, this can also
2657 * be achieved with the monitoring option in do_edsam() (if switched on by the
2658 * user), but in that case the positions need to be communicated in do_edsam(),
2659 * which is not necessary when doing flooding only. */
2660 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2662 for (i = 0; i < edi->flood.vecs.neig; i++)
2664 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2665 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2667 /* Output the current reference projection if it changes with time;
2668 * this can happen when flooding is used as harmonic restraint */
2669 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2671 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2672 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2675 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2676 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2678 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2679 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2682 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2683 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2685 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2687 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2688 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2691 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2692 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2696 } /* End of flooding-related legend entries */
2700 /* Now the ED-related entries, if essential dynamics is done */
2701 edi = ed->edpar.begin();
2702 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2704 if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2706 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2708 /* Essential dynamics, projections on eigenvectors */
2709 nice_legend_evec(&setname,
2713 get_EDgroupChar(nr_edi, nED),
2715 nice_legend_evec(&setname,
2719 get_EDgroupChar(nr_edi, nED),
2721 nice_legend_evec(&setname,
2725 get_EDgroupChar(nr_edi, nED),
2727 nice_legend_evec(&setname,
2731 get_EDgroupChar(nr_edi, nED),
2733 if (edi->vecs.radfix.neig)
2735 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2737 nice_legend_evec(&setname,
2741 get_EDgroupChar(nr_edi, nED),
2743 if (edi->vecs.radacc.neig)
2745 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2747 nice_legend_evec(&setname,
2751 get_EDgroupChar(nr_edi, nED),
2753 if (edi->vecs.radcon.neig)
2755 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2759 } /* end of 'pure' essential dynamics legend entries */
2760 n_edsam = nsets - n_flood;
2762 xvgr_legend(ed->edo, nsets, setname, oenv);
2767 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2769 1 == n_flood ? "" : "s",
2771 1 == n_edsam ? "" : "s");
2772 fprintf(ed->edo, "%s", LegendStr);
2779 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2780 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2781 std::unique_ptr<gmx::EssentialDynamics> init_edsam(const gmx::MDLogger& mdlog,
2782 const char* ediFileName,
2783 const char* edoFileName,
2784 const gmx_mtop_t* mtop,
2785 const t_inputrec* ir,
2786 const t_commrec* cr,
2787 gmx::Constraints* constr,
2788 const t_state* globalState,
2789 ObservablesHistory* oh,
2790 const gmx_output_env_t* oenv,
2791 const gmx::StartingBehavior startingBehavior)
2794 rvec* x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2795 rvec * xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2796 rvec fit_transvec; /* translation ... */
2797 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2798 rvec* ref_x_old = nullptr; /* helper pointer */
2803 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2805 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName))
2808 "The checkpoint file you provided is from an essential dynamics or flooding\n"
2809 "simulation. Please also set the .edi file on the command line with -ei.\n");
2815 "Activating essential dynamics via files passed to "
2816 "gmx mdrun -edi will change in future to be files passed to "
2817 "gmx grompp and the related .mdp options may change also.");
2819 /* Open input and output files, allocate space for ED data structure */
2820 auto edHandle = ed_open(mtop->natoms, oh, ediFileName, edoFileName, startingBehavior, oenv, cr);
2821 auto ed = edHandle->getLegacyED();
2822 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2823 constr->saveEdsamPointer(ed);
2825 /* Needed for initializing radacc radius in do_edsam */
2828 /* The input file is read by the master and the edi structures are
2829 * initialized here. Input is stored in ed->edpar. Then the edi
2830 * structures are transferred to the other nodes */
2833 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2834 * flooding vector, Essential dynamics can be applied to more than one structure
2835 * as well, but will be done in the order given in the edi file, so
2836 * expect different results for different order of edi file concatenation! */
2837 for (auto& edi : ed->edpar)
2839 init_edi(mtop, &edi);
2840 init_flood(&edi, ed, ir->delta_t);
2844 /* The master does the work here. The other nodes get the positions
2845 * not before dd_partition_system which is called after init_edsam */
2848 edsamhistory_t* EDstate = oh->edsamHistory.get();
2850 if (!EDstate->bFromCpt)
2852 /* Remove PBC, make molecule(s) subject to ED whole. */
2853 snew(x_pbc, mtop->natoms);
2854 copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop->natoms);
2855 do_pbc_first_mtop(nullptr, ir->pbcType, globalState->box, mtop, x_pbc);
2857 /* Reset pointer to first ED data set which contains the actual ED data */
2858 auto edi = ed->edpar.begin();
2859 GMX_ASSERT(!ed->edpar.empty(), "Essential dynamics parameters is undefined");
2861 /* Loop over all ED/flooding data sets (usually only one, though) */
2862 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2864 /* For multiple ED groups we use the output frequency that was specified
2865 * in the first set */
2868 edi->outfrq = ed->edpar.begin()->outfrq;
2871 /* Extract the initial reference and average positions. When starting
2872 * from .cpt, these have already been read into sref.x_old
2873 * in init_edsamstate() */
2874 if (!EDstate->bFromCpt)
2876 /* If this is the first run (i.e. no checkpoint present) we assume
2877 * that the starting positions give us the correct PBC representation */
2878 for (i = 0; i < edi->sref.nr; i++)
2880 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2883 for (i = 0; i < edi->sav.nr; i++)
2885 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2889 /* Now we have the PBC-correct start positions of the reference and
2890 average structure. We copy that over to dummy arrays on which we
2891 can apply fitting to print out the RMSD. We srenew the memory since
2892 the size of the buffers is likely different for every ED group */
2893 srenew(xfit, edi->sref.nr);
2894 srenew(xstart, edi->sav.nr);
2897 /* Reference indices are the same as average indices */
2898 ref_x_old = edi->sav.x_old;
2902 ref_x_old = edi->sref.x_old;
2904 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2905 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2907 /* Make the fit to the REFERENCE structure, get translation and rotation */
2908 fit_to_reference(xfit, fit_transvec, fit_rotmat, &(*edi));
2910 /* Output how well we fit to the reference at the start */
2911 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2913 "ED: Initial RMSD from reference after fit = %f nm",
2914 rmsd_from_structure(xfit, &edi->sref));
2915 if (EDstate->nED > 1)
2917 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2919 fprintf(stderr, "\n");
2921 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2922 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2924 /* calculate initial projections */
2925 project(xstart, &(*edi));
2927 /* For the target and origin structure both a reference (fit) and an
2928 * average structure can be provided in make_edi. If both structures
2929 * are the same, make_edi only stores one of them in the .edi file.
2930 * If they differ, first the fit and then the average structure is stored
2931 * in star (or sor), thus the number of entries in star/sor is
2932 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2933 * the size of the average group. */
2935 /* process target structure, if required */
2936 if (edi->star.nr > 0)
2938 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2940 /* get translation & rotation for fit of target structure to reference structure */
2941 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, &(*edi));
2943 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2944 if (edi->star.nr == edi->sav.nr)
2948 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2950 /* The last sav.nr indices of the target structure correspond to
2951 * the average structure, which must be projected */
2952 avindex = edi->star.nr - edi->sav.nr;
2954 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2958 rad_project(*edi, xstart, &edi->vecs.radcon);
2961 /* process structure that will serve as origin of expansion circle */
2962 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2965 "ED: Setting center of flooding potential (0 = average structure)\n");
2968 if (edi->sori.nr > 0)
2970 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2972 /* fit this structure to reference structure */
2973 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, &(*edi));
2975 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2976 if (edi->sori.nr == edi->sav.nr)
2980 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2982 /* For the projection, we need the last sav.nr indices of sori */
2983 avindex = edi->sori.nr - edi->sav.nr;
2986 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2987 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2988 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2991 "ED: The ORIGIN structure will define the flooding potential "
2993 /* Set center of flooding potential to the ORIGIN structure */
2994 rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2995 /* We already know that no (moving) reference position was provided,
2996 * therefore we can overwrite refproj[0]*/
2997 copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
3000 else /* No origin structure given */
3002 rad_project(*edi, xstart, &edi->vecs.radacc);
3003 rad_project(*edi, xstart, &edi->vecs.radfix);
3004 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
3006 if (edi->flood.bHarmonic)
3009 "ED: A (possibly changing) ref. projection will define the "
3010 "flooding potential center.\n");
3011 for (i = 0; i < edi->flood.vecs.neig; i++)
3013 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
3019 "ED: The AVERAGE structure will define the flooding potential "
3021 /* Set center of flooding potential to the center of the covariance matrix,
3022 * i.e. the average structure, i.e. zero in the projected system */
3023 for (i = 0; i < edi->flood.vecs.neig; i++)
3025 edi->flood.vecs.refproj[i] = 0.0;
3030 /* For convenience, output the center of the flooding potential for the eigenvectors */
3031 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
3033 for (i = 0; i < edi->flood.vecs.neig; i++)
3036 "ED: EV %d flooding potential center: %11.4e",
3037 edi->flood.vecs.ieig[i],
3038 edi->flood.vecs.refproj[i]);
3039 if (edi->flood.bHarmonic)
3041 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.referenceProjectionSlope[i]);
3043 fprintf(stdout, "\n");
3047 /* set starting projections for linsam */
3048 rad_project(*edi, xstart, &edi->vecs.linacc);
3049 rad_project(*edi, xstart, &edi->vecs.linfix);
3051 /* Prepare for the next edi data set: */
3054 /* Cleaning up on the master node: */
3055 if (!EDstate->bFromCpt)
3062 } /* end of MASTER only section */
3066 /* Broadcast the essential dynamics / flooding data to all nodes */
3067 broadcast_ed_data(cr, ed);
3071 /* In the single-CPU case, point the local atom numbers pointers to the global
3072 * one, so that we can use the same notation in serial and parallel case: */
3073 /* Loop over all ED data sets (usually only one, though) */
3074 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
3076 edi->sref.anrs_loc = edi->sref.anrs;
3077 edi->sav.anrs_loc = edi->sav.anrs;
3078 edi->star.anrs_loc = edi->star.anrs;
3079 edi->sori.anrs_loc = edi->sori.anrs;
3080 /* For the same reason as above, make a dummy c_ind array: */
3081 snew(edi->sav.c_ind, edi->sav.nr);
3082 /* Initialize the array */
3083 for (i = 0; i < edi->sav.nr; i++)
3085 edi->sav.c_ind[i] = i;
3087 /* In the general case we will need a different-sized array for the reference indices: */
3090 snew(edi->sref.c_ind, edi->sref.nr);
3091 for (i = 0; i < edi->sref.nr; i++)
3093 edi->sref.c_ind[i] = i;
3096 /* Point to the very same array in case of other structures: */
3097 edi->star.c_ind = edi->sav.c_ind;
3098 edi->sori.c_ind = edi->sav.c_ind;
3099 /* In the serial case, the local number of atoms is the global one: */
3100 edi->sref.nr_loc = edi->sref.nr;
3101 edi->sav.nr_loc = edi->sav.nr;
3102 edi->star.nr_loc = edi->star.nr;
3103 edi->sori.nr_loc = edi->sori.nr;
3107 /* Allocate space for ED buffer variables */
3108 /* Again, loop over ED data sets */
3109 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
3111 /* Allocate space for ED buffer variables */
3112 snew_bc(MASTER(cr), edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
3113 snew(edi->buf->do_edsam, 1);
3115 /* Space for collective ED buffer variables */
3117 /* Collective positions of atoms with the average indices */
3118 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
3119 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
3120 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
3121 /* Collective positions of atoms with the reference indices */
3124 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
3125 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
3126 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
3129 /* Get memory for flooding forces */
3130 snew(edi->flood.forces_cartesian, edi->sav.nr);
3133 /* Flush the edo file so that the user can check some things
3134 * when the simulation has started */
3144 void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[], rvec v[], const matrix box, gmx_edsam* ed)
3146 int i, edinr, iupdate = 500;
3147 matrix rotmat; /* rotation matrix */
3148 rvec transvec; /* translation vector */
3149 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3150 real dt_1; /* 1/dt */
3151 struct t_do_edsam* buf;
3152 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3153 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3156 /* Check if ED sampling has to be performed */
3157 if (ed->eEDtype == EssentialDynamicsType::None)
3162 dt_1 = 1.0 / ir->delta_t;
3164 /* Loop over all ED groups (usually one) */
3166 for (auto& edi : ed->edpar)
3169 if (bNeedDoEdsam(edi))
3172 buf = edi.buf->do_edsam;
3176 /* initialize radacc radius for slope criterion */
3177 buf->oldrad = calc_radius(edi.vecs.radacc);
3180 /* Copy the positions into buf->xc* arrays and after ED
3181 * feed back corrections to the official positions */
3183 /* Broadcast the ED positions such that every node has all of them
3184 * Every node contributes its local positions xs and stores it in
3185 * the collective buf->xcoll array. Note that for edinr > 1
3186 * xs could already have been modified by an earlier ED */
3188 communicate_group_positions(cr,
3191 buf->extra_shifts_xcoll,
3192 PAR(cr) ? buf->bUpdateShifts : TRUE,
3201 /* Only assembly reference positions if their indices differ from the average ones */
3204 communicate_group_positions(cr,
3207 buf->extra_shifts_xc_ref,
3208 PAR(cr) ? buf->bUpdateShifts : TRUE,
3218 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3219 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3220 * set bUpdateShifts=TRUE in the parallel case. */
3221 buf->bUpdateShifts = FALSE;
3223 /* Now all nodes have all of the ED positions in edi.sav->xcoll,
3224 * as well as the indices in edi.sav.anrs */
3226 /* Fit the reference indices to the reference structure */
3229 fit_to_reference(buf->xcoll, transvec, rotmat, &edi);
3233 fit_to_reference(buf->xc_ref, transvec, rotmat, &edi);
3236 /* Now apply the translation and rotation to the ED structure */
3237 translate_and_rotate(buf->xcoll, edi.sav.nr, transvec, rotmat);
3239 /* Find out how well we fit to the reference (just for output steps) */
3240 if (do_per_step(step, edi.outfrq) && MASTER(cr))
3244 /* Indices of reference and average structures are identical,
3245 * thus we can calculate the rmsd to SREF using xcoll */
3246 rmsdev = rmsd_from_structure(buf->xcoll, &edi.sref);
3250 /* We have to translate & rotate the reference atoms first */
3251 translate_and_rotate(buf->xc_ref, edi.sref.nr, transvec, rotmat);
3252 rmsdev = rmsd_from_structure(buf->xc_ref, &edi.sref);
3256 /* update radsam references, when required */
3257 if (do_per_step(step, edi.maxedsteps) && step >= edi.presteps)
3259 project(buf->xcoll, &edi);
3260 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3261 rad_project(edi, buf->xcoll, &edi.vecs.radfix);
3262 buf->oldrad = -1.e5;
3265 /* update radacc references, when required */
3266 if (do_per_step(step, iupdate) && step >= edi.presteps)
3268 edi.vecs.radacc.radius = calc_radius(edi.vecs.radacc);
3269 if (edi.vecs.radacc.radius - buf->oldrad < edi.slope)
3271 project(buf->xcoll, &edi);
3272 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3277 buf->oldrad = edi.vecs.radacc.radius;
3281 /* apply the constraints */
3282 if (step >= edi.presteps && ed_constraints(ed->eEDtype, edi))
3284 /* ED constraints should be applied already in the first MD step
3285 * (which is step 0), therefore we pass step+1 to the routine */
3286 ed_apply_constraints(buf->xcoll, &edi, step + 1 - ir->init_step);
3289 /* write to edo, when required */
3290 if (do_per_step(step, edi.outfrq))
3292 project(buf->xcoll, &edi);
3293 if (MASTER(cr) && !bSuppress)
3295 write_edo(edi, ed->edo, rmsdev);
3299 /* Copy back the positions unless monitoring only */
3300 if (ed_constraints(ed->eEDtype, edi))
3302 /* remove fitting */
3303 rmfit(edi.sav.nr, buf->xcoll, transvec, rotmat);
3305 /* Copy the ED corrected positions into the coordinate array */
3306 /* Each node copies its local part. In the serial case, nat_loc is the
3307 * total number of ED atoms */
3308 for (i = 0; i < edi.sav.nr_loc; i++)
3310 /* Unshift local ED coordinate and store in x_unsh */
3311 ed_unshift_single_coord(
3312 box, buf->xcoll[edi.sav.c_ind[i]], buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh);
3314 /* dx is the ED correction to the positions: */
3315 rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx);
3319 /* dv is the ED correction to the velocity: */
3320 svmul(dt_1, dx, dv);
3321 /* apply the velocity correction: */
3322 rvec_inc(v[edi.sav.anrs_loc[i]], dv);
3324 /* Finally apply the position correction due to ED: */
3325 copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]);
3328 } /* END of if ( bNeedDoEdsam(edi) ) */
3330 /* Prepare for the next ED group */
3332 } /* END of loop over ED groups */