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, 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]));
496 return rad = sqrt(rad);
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, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
930 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind,
931 edi->sav.x_old, box);
933 /* Only assembly REFERENCE positions if their indices differ from the average ones */
936 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref,
937 bNS, x, edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc,
938 edi->sref.c_ind, edi->sref.x_old, box);
941 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
942 * We do not need to update the shifts until the next NS step */
943 buf->bUpdateShifts = FALSE;
945 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
946 * as well as the indices in edi->sav.anrs */
948 /* Fit the reference indices to the reference structure */
951 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
955 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
958 /* Now apply the translation and rotation to the ED structure */
959 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
961 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
962 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
964 if (!edi->flood.bConstForce)
966 /* Compute Vfl(x) from flood.xproj */
967 edi->flood.Vfl = flood_energy(edi, step);
969 update_adaption(edi);
971 /* Compute the flooding forces */
975 /* Translate them into cartesian positions */
976 flood_blowup(*edi, edi->flood.forces_cartesian);
978 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
979 /* Each node rotates back its local forces */
980 transpose(rotmat, tmat);
981 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
983 /* Finally add forces to the main force variable */
984 for (i = 0; i < edi->sav.nr_loc; i++)
986 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
989 /* Output is written by the master process */
990 if (do_per_step(step, edi->outfrq) && MASTER(cr))
992 /* Output how well we fit to the reference */
995 /* Indices of reference and average structures are identical,
996 * thus we can calculate the rmsd to SREF using xcoll */
997 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
1001 /* We have to translate & rotate the reference atoms first */
1002 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1003 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1006 write_edo_flood(*edi, edo, rmsdev);
1011 /* Main flooding routine, called from do_force */
1012 extern void do_flood(const t_commrec* cr,
1013 const t_inputrec* ir,
1021 /* Write time to edo, when required. Output the time anyhow since we need
1022 * it in the output file for ED constraints. */
1023 if (MASTER(cr) && do_per_step(step, ed->edpar.begin()->outfrq))
1025 fprintf(ed->edo, "\n%12f", ir->init_t + step * ir->delta_t);
1028 if (ed->eEDtype != EssentialDynamicsType::Flooding)
1033 for (auto& edi : ed->edpar)
1035 /* Call flooding for one matrix */
1036 if (edi.flood.vecs.neig)
1038 do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS);
1044 /* Called by init_edi, configure some flooding related variables and structures,
1045 * print headers to output files */
1046 static void init_flood(t_edpar* edi, gmx_edsam* ed, real dt)
1051 edi->flood.Efl = edi->flood.constEfl;
1055 if (edi->flood.vecs.neig)
1057 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1058 ed->eEDtype = EssentialDynamicsType::Flooding;
1060 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig,
1061 edi->flood.vecs.neig > 1 ? "s" : "");
1063 if (edi->flood.bConstForce)
1065 /* We have used stpsz as a vehicle to carry the fproj values for constant
1066 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1067 * in const force flooding, fproj is never changed. */
1068 for (i = 0; i < edi->flood.vecs.neig; i++)
1070 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1072 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1073 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1081 /*********** Energy book keeping ******/
1082 static void get_flood_enx_names(t_edpar* edi, char** names, int* nnames) /* get header of energies */
1091 srenew(names, count);
1092 sprintf(buf, "Vfl_%d", count);
1093 names[count - 1] = gmx_strdup(buf);
1094 actual = actual->next_edi;
1097 *nnames = count - 1;
1101 static void get_flood_energies(t_edpar* edi, real Vfl[], int nnames)
1103 /*fl has to be big enough to capture nnames-many entries*/
1112 Vfl[count - 1] = actual->flood.Vfl;
1113 actual = actual->next_edi;
1116 if (nnames != count - 1)
1118 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1121 /************* END of FLOODING IMPLEMENTATION ****************************/
1125 /* This function opens the ED input and output files, reads in all datasets it finds
1126 * in the input file, and cross-checks whether the .edi file information is consistent
1127 * with the essential dynamics data found in the checkpoint file (if present).
1128 * gmx make_edi can be used to create an .edi input file.
1130 static std::unique_ptr<gmx::EssentialDynamics> ed_open(int natoms,
1131 ObservablesHistory* oh,
1132 const char* ediFileName,
1133 const char* edoFileName,
1134 const gmx::StartingBehavior startingBehavior,
1135 const gmx_output_env_t* oenv,
1136 const t_commrec* cr)
1138 auto edHandle = std::make_unique<gmx::EssentialDynamics>();
1139 auto ed = edHandle->getLegacyED();
1140 /* We want to perform ED (this switch might later be upgraded to EssentialDynamicsType::Flooding) */
1141 ed->eEDtype = EssentialDynamicsType::EDSampling;
1145 // If we start from a checkpoint file, we already have an edsamHistory struct
1146 if (oh->edsamHistory == nullptr)
1148 oh->edsamHistory = std::make_unique<edsamhistory_t>(edsamhistory_t{});
1150 edsamhistory_t* EDstate = oh->edsamHistory.get();
1152 /* Read the edi input file: */
1153 ed->edpar = read_edi_file(ediFileName, natoms);
1155 /* Make sure the checkpoint was produced in a run using this .edi file */
1156 if (EDstate->bFromCpt)
1158 crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1162 EDstate->nED = ed->edpar.size();
1164 init_edsamstate(*ed, EDstate);
1166 /* The master opens the ED output file */
1167 if (startingBehavior == gmx::StartingBehavior::RestartWithAppending)
1169 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1173 ed->edo = xvgropen(edoFileName, "Essential dynamics / flooding output", "Time (ps)",
1174 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1176 /* Make a descriptive legend */
1177 write_edo_legend(ed, EDstate->nED, oenv);
1184 /* Broadcasts the structure data */
1185 static void bc_ed_positions(const t_commrec* cr, struct gmx_edx* s, EssentialDynamicsStructure stype)
1187 snew_bc(cr, s->anrs, s->nr); /* Index numbers */
1188 snew_bc(cr, s->x, s->nr); /* Positions */
1189 nblock_bc(cr, s->nr, s->anrs);
1190 nblock_bc(cr, s->nr, s->x);
1192 /* For the average & reference structures we need an array for the collective indices,
1193 * and we need to broadcast the masses as well */
1194 if (stype == EssentialDynamicsStructure::Average || stype == EssentialDynamicsStructure::Reference)
1196 /* We need these additional variables in the parallel case: */
1197 snew(s->c_ind, s->nr); /* Collective indices */
1198 /* Local atom indices get assigned in dd_make_local_group_indices.
1199 * There, also memory is allocated */
1200 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1201 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1202 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1205 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1206 if (stype == EssentialDynamicsStructure::Reference)
1208 snew_bc(cr, s->m, s->nr);
1209 nblock_bc(cr, s->nr, s->m);
1212 /* For the average structure we might need the masses for mass-weighting */
1213 if (stype == EssentialDynamicsStructure::Average)
1215 snew_bc(cr, s->sqrtm, s->nr);
1216 nblock_bc(cr, s->nr, s->sqrtm);
1217 snew_bc(cr, s->m, s->nr);
1218 nblock_bc(cr, s->nr, s->m);
1223 /* Broadcasts the eigenvector data */
1224 static void bc_ed_vecs(const t_commrec* cr, t_eigvec* ev, int length)
1228 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1229 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1230 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1231 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1232 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1234 nblock_bc(cr, ev->neig, ev->ieig);
1235 nblock_bc(cr, ev->neig, ev->stpsz);
1236 nblock_bc(cr, ev->neig, ev->xproj);
1237 nblock_bc(cr, ev->neig, ev->fproj);
1238 nblock_bc(cr, ev->neig, ev->refproj);
1240 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1241 for (i = 0; i < ev->neig; i++)
1243 snew_bc(cr, ev->vec[i], length);
1244 nblock_bc(cr, length, ev->vec[i]);
1249 /* Broadcasts the ED / flooding data to other nodes
1250 * and allocates memory where needed */
1251 static void broadcast_ed_data(const t_commrec* cr, gmx_edsam* ed)
1253 /* Master lets the other nodes know if its ED only or also flooding */
1254 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1256 int numedis = ed->edpar.size();
1257 /* First let everybody know how many ED data sets to expect */
1258 gmx_bcast(sizeof(numedis), &numedis, cr);
1259 nblock_abc(cr, numedis, &(ed->edpar));
1260 /* Now transfer the ED data set(s) */
1261 for (auto& edi : ed->edpar)
1263 /* Broadcast a single ED data set */
1266 /* Broadcast positions */
1267 bc_ed_positions(cr, &(edi.sref),
1268 EssentialDynamicsStructure::Reference); /* reference positions (don't broadcast masses) */
1269 bc_ed_positions(cr, &(edi.sav),
1270 EssentialDynamicsStructure::Average); /* average positions (do broadcast masses as well) */
1271 bc_ed_positions(cr, &(edi.star), EssentialDynamicsStructure::Target); /* target positions */
1272 bc_ed_positions(cr, &(edi.sori), EssentialDynamicsStructure::Origin); /* origin positions */
1274 /* Broadcast eigenvectors */
1275 bc_ed_vecs(cr, &edi.vecs.mon, edi.sav.nr);
1276 bc_ed_vecs(cr, &edi.vecs.linfix, edi.sav.nr);
1277 bc_ed_vecs(cr, &edi.vecs.linacc, edi.sav.nr);
1278 bc_ed_vecs(cr, &edi.vecs.radfix, edi.sav.nr);
1279 bc_ed_vecs(cr, &edi.vecs.radacc, edi.sav.nr);
1280 bc_ed_vecs(cr, &edi.vecs.radcon, edi.sav.nr);
1281 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1282 bc_ed_vecs(cr, &edi.flood.vecs, edi.sav.nr);
1284 /* For harmonic restraints the reference projections can change with time */
1285 if (edi.flood.bHarmonic)
1287 snew_bc(cr, edi.flood.initialReferenceProjection, edi.flood.vecs.neig);
1288 snew_bc(cr, edi.flood.referenceProjectionSlope, edi.flood.vecs.neig);
1289 nblock_bc(cr, edi.flood.vecs.neig, edi.flood.initialReferenceProjection);
1290 nblock_bc(cr, edi.flood.vecs.neig, edi.flood.referenceProjectionSlope);
1296 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1297 static void init_edi(const gmx_mtop_t* mtop, t_edpar* edi)
1300 real totalmass = 0.0;
1303 /* NOTE Init_edi is executed on the master process only
1304 * The initialized data sets are then transmitted to the
1305 * other nodes in broadcast_ed_data */
1307 /* evaluate masses (reference structure) */
1308 snew(edi->sref.m, edi->sref.nr);
1310 for (i = 0; i < edi->sref.nr; i++)
1314 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1318 edi->sref.m[i] = 1.0;
1321 /* Check that every m > 0. Bad things will happen otherwise. */
1322 if (edi->sref.m[i] <= 0.0)
1325 "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1326 "For a mass-weighted fit, all reference structure atoms need to have a mass "
1328 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1329 "atoms from the reference structure by creating a proper index group.\n",
1330 i, edi->sref.anrs[i] + 1, edi->sref.m[i]);
1333 totalmass += edi->sref.m[i];
1335 edi->sref.mtot = totalmass;
1337 /* Masses m and sqrt(m) for the average structure. Note that m
1338 * is needed if forces have to be evaluated in do_edsam */
1339 snew(edi->sav.sqrtm, edi->sav.nr);
1340 snew(edi->sav.m, edi->sav.nr);
1341 for (i = 0; i < edi->sav.nr; i++)
1343 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1346 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1350 edi->sav.sqrtm[i] = 1.0;
1353 /* Check that every m > 0. Bad things will happen otherwise. */
1354 if (edi->sav.sqrtm[i] <= 0.0)
1357 "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1358 "For ED with mass-weighting, all average structure atoms need to have a mass "
1360 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1361 "atoms from the average structure by creating a proper index group.\n",
1362 i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1366 /* put reference structure in origin */
1367 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1371 translate_x(edi->sref.x, edi->sref.nr, com);
1373 /* Init ED buffer */
1378 static void check(const char* line, const char* label)
1380 if (!strstr(line, label))
1383 "Could not find input parameter %s at expected position in edsam input-file "
1384 "(.edi)\nline read instead is %s",
1390 static int read_checked_edint(FILE* file, const char* label)
1392 char line[STRLEN + 1];
1395 fgets2(line, STRLEN, file);
1397 fgets2(line, STRLEN, file);
1398 sscanf(line, max_ev_fmt_d, &idum);
1402 static bool read_checked_edbool(FILE* file, const char* label)
1404 return read_checked_edint(file, label) != 0;
1407 static int read_edint(FILE* file, bool* bEOF)
1409 char line[STRLEN + 1];
1414 eof = fgets2(line, STRLEN, file);
1420 eof = fgets2(line, STRLEN, file);
1426 sscanf(line, max_ev_fmt_d, &idum);
1432 static real read_checked_edreal(FILE* file, const char* label)
1434 char line[STRLEN + 1];
1438 fgets2(line, STRLEN, file);
1440 fgets2(line, STRLEN, file);
1441 sscanf(line, max_ev_fmt_lf, &rdum);
1442 return static_cast<real>(rdum); /* always read as double and convert to single */
1446 static void read_edx(FILE* file, int number, int* anrs, rvec* x)
1449 char line[STRLEN + 1];
1453 for (i = 0; i < number; i++)
1455 fgets2(line, STRLEN, file);
1456 sscanf(line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1457 anrs[i]--; /* we are reading FORTRAN indices */
1458 for (j = 0; j < 3; j++)
1460 x[i][j] = d[j]; /* always read as double and convert to single */
1467 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1468 * \param[in] in the file to read from
1469 * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1470 * \param[out] vec the eigen-vectors
1471 * \param[in] nEig the number of eigenvectors
1473 void scan_edvec(FILE* in, int numAtoms, rvec*** vec, int nEig)
1476 for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1478 snew((*vec)[iEigenvector], numAtoms);
1479 for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1481 char line[STRLEN + 1];
1483 fgets2(line, STRLEN, in);
1484 sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1485 (*vec)[iEigenvector][iAtom][XX] = x;
1486 (*vec)[iEigenvector][iAtom][YY] = y;
1487 (*vec)[iEigenvector][iAtom][ZZ] = z;
1491 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1492 * \param[in] in input file
1493 * \param[in] nrAtoms number of atoms/coordinates
1494 * \param[out] tvec the eigenvector
1496 void read_edvec(FILE* in, int nrAtoms, t_eigvec* tvec)
1498 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1499 if (tvec->neig <= 0)
1504 snew(tvec->ieig, tvec->neig);
1505 snew(tvec->stpsz, tvec->neig);
1506 for (int i = 0; i < tvec->neig; i++)
1508 char line[STRLEN + 1];
1509 fgets2(line, STRLEN, in);
1510 int numEigenVectors;
1512 const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1515 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1517 tvec->ieig[i] = numEigenVectors;
1518 tvec->stpsz[i] = stepSize;
1519 } /* end of loop over eigenvectors */
1521 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1523 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1525 * Eigenvector and its intial reference and slope are stored on the
1526 * same line in the .edi format. To avoid re-winding the .edi file,
1527 * the reference eigen-vector and reference data are read in one go.
1529 * \param[in] in input file
1530 * \param[in] nrAtoms number of atoms/coordinates
1531 * \param[out] tvec the eigenvector
1532 * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1533 * \param[out] referenceProjectionSlope The slope of the reference projections.
1535 bool readEdVecWithReferenceProjection(FILE* in,
1538 real** initialReferenceProjection,
1539 real** referenceProjectionSlope)
1541 bool providesReference = false;
1542 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1543 if (tvec->neig <= 0)
1548 snew(tvec->ieig, tvec->neig);
1549 snew(tvec->stpsz, tvec->neig);
1550 snew(*initialReferenceProjection, tvec->neig);
1551 snew(*referenceProjectionSlope, tvec->neig);
1552 for (int i = 0; i < tvec->neig; i++)
1554 char line[STRLEN + 1];
1555 fgets2(line, STRLEN, in);
1556 int numEigenVectors;
1557 double stepSize = 0.;
1558 double referenceProjection = 0.;
1559 double referenceSlope = 0.;
1560 const int nscan = sscanf(line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize,
1561 &referenceProjection, &referenceSlope);
1562 /* Zero out values which were not scanned */
1566 /* Every 4 values read, including reference position */
1567 providesReference = true;
1570 /* A reference position is provided */
1571 providesReference = true;
1572 /* No value for slope, set to 0 */
1573 referenceSlope = 0.0;
1576 /* No values for reference projection and slope, set to 0 */
1577 referenceProjection = 0.0;
1578 referenceSlope = 0.0;
1582 "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> "
1583 "<refproj> <refproj-slope>\n",
1586 (*initialReferenceProjection)[i] = referenceProjection;
1587 (*referenceProjectionSlope)[i] = referenceSlope;
1589 tvec->ieig[i] = numEigenVectors;
1590 tvec->stpsz[i] = stepSize;
1591 } /* end of loop over eigenvectors */
1593 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1594 return providesReference;
1597 /*!\brief Allocate working memory for the eigenvectors.
1598 * \param[in,out] tvec eigenvector for which memory will be allocated
1600 void setup_edvec(t_eigvec* tvec)
1602 snew(tvec->xproj, tvec->neig);
1603 snew(tvec->fproj, tvec->neig);
1604 snew(tvec->refproj, tvec->neig);
1609 /* Check if the same atom indices are used for reference and average positions */
1610 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1615 /* If the number of atoms differs between the two structures,
1616 * they cannot be identical */
1617 if (sref.nr != sav.nr)
1622 /* Now that we know that both stuctures have the same number of atoms,
1623 * check if also the indices are identical */
1624 for (i = 0; i < sav.nr; i++)
1626 if (sref.anrs[i] != sav.anrs[i])
1632 "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1640 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1641 constexpr int c_oldestUnsupportedMagicNumber = 666;
1642 //! Oldest (lowest) magic number indicating supported essential dynamics input
1643 constexpr int c_oldestSupportedMagicNumber = 669;
1644 //! Latest (highest) magic number indicating supported essential dynamics input
1645 constexpr int c_latestSupportedMagicNumber = 670;
1647 /*!\brief Set up essential dynamics work parameters.
1648 * \param[in] edi Essential dynamics working structure.
1650 void setup_edi(t_edpar* edi)
1652 snew(edi->sref.x_old, edi->sref.nr);
1653 edi->sref.sqrtm = nullptr;
1654 snew(edi->sav.x_old, edi->sav.nr);
1655 if (edi->star.nr > 0)
1657 edi->star.sqrtm = nullptr;
1659 if (edi->sori.nr > 0)
1661 edi->sori.sqrtm = nullptr;
1663 setup_edvec(&edi->vecs.linacc);
1664 setup_edvec(&edi->vecs.mon);
1665 setup_edvec(&edi->vecs.linfix);
1666 setup_edvec(&edi->vecs.linacc);
1667 setup_edvec(&edi->vecs.radfix);
1668 setup_edvec(&edi->vecs.radacc);
1669 setup_edvec(&edi->vecs.radcon);
1670 setup_edvec(&edi->flood.vecs);
1673 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1674 * \param[in] magicNumber indicates the essential dynamics input file version
1675 * \returns true if CONST_FORCE_FLOODING is supported
1677 bool ediFileHasConstForceFlooding(int magicNumber)
1679 return magicNumber > c_oldestSupportedMagicNumber;
1682 /*!\brief Reports reading success of the essential dynamics input file magic number.
1683 * \param[in] in pointer to input file
1684 * \param[out] magicNumber determines the edi file version
1685 * \returns true if setting file version from input file was successful.
1687 bool ediFileHasValidDataSet(FILE* in, int* magicNumber)
1689 /* the edi file is not free format, so expect problems if the input is corrupt. */
1690 bool reachedEndOfFile = true;
1691 /* check the magic number */
1692 *magicNumber = read_edint(in, &reachedEndOfFile);
1693 /* Check whether we have reached the end of the input file */
1694 return !reachedEndOfFile;
1697 /*!\brief Trigger fatal error if magic number is unsupported.
1698 * \param[in] magicNumber A magic number read from the edi file
1699 * \param[in] fn name of input file for error reporting
1701 void exitWhenMagicNumberIsInvalid(int magicNumber, const char* fn)
1704 if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1706 if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1709 "Wrong magic number: Use newest version of make_edi to produce edi file");
1713 gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1718 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1720 * \param[in] in input file
1721 * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1722 * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1723 * \param[in] fn the file name of the input file for error reporting
1724 * \returns edi essential dynamics parameters
1726 t_edpar read_edi(FILE* in, int nr_mdatoms, bool hasConstForceFlooding, const char* fn)
1730 /* check the number of atoms */
1731 edi.nini = read_edint(in, &bEOF);
1732 if (edi.nini != nr_mdatoms)
1734 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi.nini,
1738 /* Done checking. For the rest we blindly trust the input */
1739 edi.fitmas = read_checked_edbool(in, "FITMAS");
1740 edi.pcamas = read_checked_edbool(in, "ANALYSIS_MAS");
1741 edi.outfrq = read_checked_edint(in, "OUTFRQ");
1742 edi.maxedsteps = read_checked_edint(in, "MAXLEN");
1743 edi.slope = read_checked_edreal(in, "SLOPECRIT");
1745 edi.presteps = read_checked_edint(in, "PRESTEPS");
1746 edi.flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1747 edi.flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1748 edi.flood.tau = read_checked_edreal(in, "TAU");
1749 edi.flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1750 edi.flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1751 edi.flood.kT = read_checked_edreal(in, "KT");
1752 edi.flood.bHarmonic = read_checked_edbool(in, "HARMONIC");
1753 if (hasConstForceFlooding)
1755 edi.flood.bConstForce = read_checked_edbool(in, "CONST_FORCE_FLOODING");
1759 edi.flood.bConstForce = FALSE;
1761 edi.sref.nr = read_checked_edint(in, "NREF");
1763 /* allocate space for reference positions and read them */
1764 snew(edi.sref.anrs, edi.sref.nr);
1765 snew(edi.sref.x, edi.sref.nr);
1766 read_edx(in, edi.sref.nr, edi.sref.anrs, edi.sref.x);
1768 /* average positions. they define which atoms will be used for ED sampling */
1769 edi.sav.nr = read_checked_edint(in, "NAV");
1770 snew(edi.sav.anrs, edi.sav.nr);
1771 snew(edi.sav.x, edi.sav.nr);
1772 read_edx(in, edi.sav.nr, edi.sav.anrs, edi.sav.x);
1774 /* Check if the same atom indices are used for reference and average positions */
1775 edi.bRefEqAv = check_if_same(edi.sref, edi.sav);
1779 read_edvec(in, edi.sav.nr, &edi.vecs.mon);
1780 read_edvec(in, edi.sav.nr, &edi.vecs.linfix);
1781 read_edvec(in, edi.sav.nr, &edi.vecs.linacc);
1782 read_edvec(in, edi.sav.nr, &edi.vecs.radfix);
1783 read_edvec(in, edi.sav.nr, &edi.vecs.radacc);
1784 read_edvec(in, edi.sav.nr, &edi.vecs.radcon);
1786 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1787 bool bHaveReference = false;
1788 if (edi.flood.bHarmonic)
1790 bHaveReference = readEdVecWithReferenceProjection(in, edi.sav.nr, &edi.flood.vecs,
1791 &edi.flood.initialReferenceProjection,
1792 &edi.flood.referenceProjectionSlope);
1796 read_edvec(in, edi.sav.nr, &edi.flood.vecs);
1799 /* target positions */
1800 edi.star.nr = read_edint(in, &bEOF);
1801 if (edi.star.nr > 0)
1803 snew(edi.star.anrs, edi.star.nr);
1804 snew(edi.star.x, edi.star.nr);
1805 read_edx(in, edi.star.nr, edi.star.anrs, edi.star.x);
1808 /* positions defining origin of expansion circle */
1809 edi.sori.nr = read_edint(in, &bEOF);
1810 if (edi.sori.nr > 0)
1814 /* Both an -ori structure and a at least one manual reference point have been
1815 * specified. That's ambiguous and probably not intentional. */
1817 "ED: An origin structure has been provided and a at least one (moving) "
1819 " point was manually specified in the edi file. That is ambiguous. "
1822 snew(edi.sori.anrs, edi.sori.nr);
1823 snew(edi.sori.x, edi.sori.nr);
1824 read_edx(in, edi.sori.nr, edi.sori.anrs, edi.sori.x);
1829 std::vector<t_edpar> read_edi_file(const char* fn, int nr_mdatoms)
1831 std::vector<t_edpar> essentialDynamicsParameters;
1833 /* This routine is executed on the master only */
1835 /* Open the .edi parameter input file */
1836 in = gmx_fio_fopen(fn, "r");
1837 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1839 /* Now read a sequence of ED input parameter sets from the edi file */
1840 int ediFileMagicNumber;
1841 while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1843 exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1844 bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1845 auto edi = read_edi(in, nr_mdatoms, hasConstForceFlooding, fn);
1847 essentialDynamicsParameters.emplace_back(edi);
1849 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1850 /* We need to allocate space for the data: */
1852 if (essentialDynamicsParameters.empty())
1854 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1857 fprintf(stderr, "ED: Found %zu ED group%s.\n", essentialDynamicsParameters.size(),
1858 essentialDynamicsParameters.size() > 1 ? "s" : "");
1860 /* Close the .edi file again */
1863 return essentialDynamicsParameters;
1865 } // anonymous namespace
1870 rvec* xcopy; /* Working copy of the positions in fit_to_reference */
1873 /* Fit the current positions to the reference positions
1874 * Do not actually do the fit, just return rotation and translation.
1875 * Note that the COM of the reference structure was already put into
1876 * the origin by init_edi. */
1877 static void fit_to_reference(rvec* xcoll, /* The positions to be fitted */
1878 rvec transvec, /* The translation vector */
1879 matrix rotmat, /* The rotation matrix */
1880 t_edpar* edi) /* Just needed for do_edfit */
1882 rvec com; /* center of mass */
1884 struct t_fit_to_ref* loc;
1887 /* Allocate memory the first time this routine is called for each edi group */
1888 if (nullptr == edi->buf->fit_to_ref)
1890 snew(edi->buf->fit_to_ref, 1);
1891 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1893 loc = edi->buf->fit_to_ref;
1895 /* We do not touch the original positions but work on a copy. */
1896 for (i = 0; i < edi->sref.nr; i++)
1898 copy_rvec(xcoll[i], loc->xcopy[i]);
1901 /* Calculate the center of mass */
1902 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1904 transvec[XX] = -com[XX];
1905 transvec[YY] = -com[YY];
1906 transvec[ZZ] = -com[ZZ];
1908 /* Subtract the center of mass from the copy */
1909 translate_x(loc->xcopy, edi->sref.nr, transvec);
1911 /* Determine the rotation matrix */
1912 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1916 static void translate_and_rotate(rvec* x, /* The positions to be translated and rotated */
1917 int nat, /* How many positions are there? */
1918 rvec transvec, /* The translation vector */
1919 matrix rotmat) /* The rotation matrix */
1922 translate_x(x, nat, transvec);
1925 rotate_x(x, nat, rotmat);
1929 /* Gets the rms deviation of the positions to the structure s */
1930 /* fit_to_structure has to be called before calling this routine! */
1931 static real rmsd_from_structure(rvec* x, /* The positions under consideration */
1932 struct gmx_edx* s) /* The structure from which the rmsd shall be computed */
1938 for (i = 0; i < s->nr; i++)
1940 rmsd += distance2(s->x[i], x[i]);
1943 rmsd /= static_cast<real>(s->nr);
1950 void dd_make_local_ed_indices(gmx_domdec_t* dd, struct gmx_edsam* ed)
1952 if (ed->eEDtype != EssentialDynamicsType::None)
1954 /* Loop over ED groups */
1956 for (auto& edi : ed->edpar)
1958 /* Local atoms of the reference structure (for fitting), need only be assembled
1959 * if their indices differ from the average ones */
1962 dd_make_local_group_indices(dd->ga2la, edi.sref.nr, edi.sref.anrs, &edi.sref.nr_loc,
1963 &edi.sref.anrs_loc, &edi.sref.nalloc_loc, edi.sref.c_ind);
1966 /* Local atoms of the average structure (on these ED will be performed) */
1967 dd_make_local_group_indices(dd->ga2la, edi.sav.nr, edi.sav.anrs, &edi.sav.nr_loc,
1968 &edi.sav.anrs_loc, &edi.sav.nalloc_loc, edi.sav.c_ind);
1970 /* Indicate that the ED shift vectors for this structure need to be updated
1971 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1972 edi.buf->do_edsam->bUpdateShifts = TRUE;
1978 static inline void ed_unshift_single_coord(const matrix box, const rvec x, const ivec is, rvec xu)
1989 xu[XX] = x[XX] - tx * box[XX][XX] - ty * box[YY][XX] - tz * box[ZZ][XX];
1990 xu[YY] = x[YY] - ty * box[YY][YY] - tz * box[ZZ][YY];
1991 xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
1995 xu[XX] = x[XX] - tx * box[XX][XX];
1996 xu[YY] = x[YY] - ty * box[YY][YY];
1997 xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
2003 /*!\brief Apply fixed linear constraints to essential dynamics variable.
2004 * \param[in,out] xcoll The collected coordinates.
2005 * \param[in] edi the essential dynamics parameters
2006 * \param[in] step the current simulation step
2008 void do_linfix(rvec* xcoll, const t_edpar& edi, int64_t step)
2010 /* loop over linfix vectors */
2011 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2013 /* calculate the projection */
2014 real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
2016 /* calculate the correction */
2017 real preFactor = edi.vecs.linfix.refproj[i] + step * edi.vecs.linfix.stpsz[i] - proj;
2019 /* apply the correction */
2020 preFactor /= edi.sav.sqrtm[i];
2021 for (int j = 0; j < edi.sav.nr; j++)
2023 rvec differenceVector;
2024 svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
2025 rvec_inc(xcoll[j], differenceVector);
2030 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
2031 * \param[in,out] xcoll The collected coordinates.
2032 * \param[in] edi the essential dynamics parameters
2034 void do_linacc(rvec* xcoll, t_edpar* edi)
2036 /* loop over linacc vectors */
2037 for (int i = 0; i < edi->vecs.linacc.neig; i++)
2039 /* calculate the projection */
2040 real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
2042 /* calculate the correction */
2043 real preFactor = 0.0;
2044 if (edi->vecs.linacc.stpsz[i] > 0.0)
2046 if ((proj - edi->vecs.linacc.refproj[i]) < 0.0)
2048 preFactor = edi->vecs.linacc.refproj[i] - proj;
2051 if (edi->vecs.linacc.stpsz[i] < 0.0)
2053 if ((proj - edi->vecs.linacc.refproj[i]) > 0.0)
2055 preFactor = edi->vecs.linacc.refproj[i] - proj;
2059 /* apply the correction */
2060 preFactor /= edi->sav.sqrtm[i];
2061 for (int j = 0; j < edi->sav.nr; j++)
2063 rvec differenceVector;
2064 svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2065 rvec_inc(xcoll[j], differenceVector);
2068 /* new positions will act as reference */
2069 edi->vecs.linacc.refproj[i] = proj + preFactor;
2075 static void do_radfix(rvec* xcoll, t_edpar* edi)
2078 real *proj, rad = 0.0, ratio;
2082 if (edi->vecs.radfix.neig == 0)
2087 snew(proj, edi->vecs.radfix.neig);
2089 /* loop over radfix vectors */
2090 for (i = 0; i < edi->vecs.radfix.neig; i++)
2092 /* calculate the projections, radius */
2093 proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2094 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2098 ratio = (edi->vecs.radfix.stpsz[0] + edi->vecs.radfix.radius) / rad - 1.0;
2099 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2101 /* loop over radfix vectors */
2102 for (i = 0; i < edi->vecs.radfix.neig; i++)
2104 proj[i] -= edi->vecs.radfix.refproj[i];
2106 /* apply the correction */
2107 proj[i] /= edi->sav.sqrtm[i];
2109 for (j = 0; j < edi->sav.nr; j++)
2111 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2112 rvec_inc(xcoll[j], vec_dum);
2120 static void do_radacc(rvec* xcoll, t_edpar* edi)
2123 real *proj, rad = 0.0, ratio = 0.0;
2127 if (edi->vecs.radacc.neig == 0)
2132 snew(proj, edi->vecs.radacc.neig);
2134 /* loop over radacc vectors */
2135 for (i = 0; i < edi->vecs.radacc.neig; i++)
2137 /* calculate the projections, radius */
2138 proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2139 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2143 /* only correct when radius decreased */
2144 if (rad < edi->vecs.radacc.radius)
2146 ratio = edi->vecs.radacc.radius / rad - 1.0;
2150 edi->vecs.radacc.radius = rad;
2153 /* loop over radacc vectors */
2154 for (i = 0; i < edi->vecs.radacc.neig; i++)
2156 proj[i] -= edi->vecs.radacc.refproj[i];
2158 /* apply the correction */
2159 proj[i] /= edi->sav.sqrtm[i];
2161 for (j = 0; j < edi->sav.nr; j++)
2163 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2164 rvec_inc(xcoll[j], vec_dum);
2176 static void do_radcon(rvec* xcoll, t_edpar* edi)
2179 real rad = 0.0, ratio = 0.0;
2180 struct t_do_radcon* loc;
2185 if (edi->buf->do_radcon != nullptr)
2192 snew(edi->buf->do_radcon, 1);
2194 loc = edi->buf->do_radcon;
2196 if (edi->vecs.radcon.neig == 0)
2203 snew(loc->proj, edi->vecs.radcon.neig);
2206 /* loop over radcon vectors */
2207 for (i = 0; i < edi->vecs.radcon.neig; i++)
2209 /* calculate the projections, radius */
2210 loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2211 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2214 /* only correct when radius increased */
2215 if (rad > edi->vecs.radcon.radius)
2217 ratio = edi->vecs.radcon.radius / rad - 1.0;
2219 /* loop over radcon vectors */
2220 for (i = 0; i < edi->vecs.radcon.neig; i++)
2222 /* apply the correction */
2223 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2224 loc->proj[i] /= edi->sav.sqrtm[i];
2225 loc->proj[i] *= ratio;
2227 for (j = 0; j < edi->sav.nr; j++)
2229 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2230 rvec_inc(xcoll[j], vec_dum);
2236 edi->vecs.radcon.radius = rad;
2241 static void ed_apply_constraints(rvec* xcoll, t_edpar* edi, int64_t step)
2246 /* subtract the average positions */
2247 for (i = 0; i < edi->sav.nr; i++)
2249 rvec_dec(xcoll[i], edi->sav.x[i]);
2252 /* apply the constraints */
2255 do_linfix(xcoll, *edi, step);
2257 do_linacc(xcoll, edi);
2260 do_radfix(xcoll, edi);
2262 do_radacc(xcoll, edi);
2263 do_radcon(xcoll, edi);
2265 /* add back the average positions */
2266 for (i = 0; i < edi->sav.nr; i++)
2268 rvec_inc(xcoll[i], edi->sav.x[i]);
2275 /*!\brief Write out the projections onto the eigenvectors.
2276 * The order of output corresponds to ed_output_legend().
2277 * \param[in] edi The essential dyanmics parameters
2278 * \param[in] fp The output file
2279 * \param[in] rmsd the rmsd to the reference structure
2281 void write_edo(const t_edpar& edi, FILE* fp, real rmsd)
2283 /* Output how well we fit to the reference structure */
2284 fprintf(fp, EDcol_ffmt, rmsd);
2286 for (int i = 0; i < edi.vecs.mon.neig; i++)
2288 fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2291 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2293 fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2296 for (int i = 0; i < edi.vecs.linacc.neig; i++)
2298 fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2301 for (int i = 0; i < edi.vecs.radfix.neig; i++)
2303 fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2305 if (edi.vecs.radfix.neig)
2307 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2310 for (int i = 0; i < edi.vecs.radacc.neig; i++)
2312 fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2314 if (edi.vecs.radacc.neig)
2316 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2319 for (int i = 0; i < edi.vecs.radcon.neig; i++)
2321 fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2323 if (edi.vecs.radcon.neig)
2325 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2331 /* Returns if any constraints are switched on */
2332 static bool ed_constraints(EssentialDynamicsType edtype, const t_edpar& edi)
2334 if (edtype == EssentialDynamicsType::EDSampling || edtype == EssentialDynamicsType::Flooding)
2336 return ((edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0) || (edi.vecs.radfix.neig != 0)
2337 || (edi.vecs.radacc.neig != 0) || (edi.vecs.radcon.neig != 0));
2343 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for
2344 * flooding/ umbrella sampling simulations. */
2345 static void copyEvecReference(t_eigvec* floodvecs, real* initialReferenceProjection)
2350 if (nullptr == initialReferenceProjection)
2352 snew(initialReferenceProjection, floodvecs->neig);
2355 for (i = 0; i < floodvecs->neig; i++)
2357 initialReferenceProjection[i] = floodvecs->refproj[i];
2362 /* Call on MASTER only. Check whether the essential dynamics / flooding
2363 * groups of the checkpoint file are consistent with the provided .edi file. */
2364 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam& ed, edsamhistory_t* EDstate)
2366 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2369 "Essential dynamics and flooding can only be switched on (or off) at the\n"
2370 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2371 "it must also continue with/without ED constraints when checkpointing.\n"
2372 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2373 "from without a checkpoint.\n");
2376 for (size_t edinum = 0; edinum < ed.edpar.size(); ++edinum)
2378 /* Check number of atoms in the reference and average structures */
2379 if (EDstate->nref[edinum] != ed.edpar[edinum].sref.nr)
2382 "The number of reference structure atoms in ED group %c is\n"
2383 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2384 get_EDgroupChar(edinum + 1, 0), EDstate->nref[edinum], ed.edpar[edinum].sref.nr);
2386 if (EDstate->nav[edinum] != ed.edpar[edinum].sav.nr)
2389 "The number of average structure atoms in ED group %c is\n"
2390 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2391 get_EDgroupChar(edinum + 1, 0), EDstate->nav[edinum], ed.edpar[edinum].sav.nr);
2395 if (gmx::ssize(ed.edpar) != EDstate->nED)
2398 "The number of essential dynamics / flooding groups is not consistent.\n"
2399 "There are %d ED groups in the .cpt file, but %zu in the .edi file!\n"
2400 "Are you sure this is the correct .edi file?\n",
2401 EDstate->nED, ed.edpar.size());
2406 /* The edsamstate struct stores the information we need to make the ED group
2407 * whole again after restarts from a checkpoint file. Here we do the following:
2408 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2409 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2410 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2411 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2412 * all ED structures at the last time step. */
2413 static void init_edsamstate(const gmx_edsam& ed, edsamhistory_t* EDstate)
2415 snew(EDstate->old_sref_p, EDstate->nED);
2416 snew(EDstate->old_sav_p, EDstate->nED);
2418 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2419 if (!EDstate->bFromCpt)
2421 snew(EDstate->nref, EDstate->nED);
2422 snew(EDstate->nav, EDstate->nED);
2425 /* Loop over all ED/flooding data sets (usually only one, though) */
2426 for (int nr_edi = 0; nr_edi < EDstate->nED; nr_edi++)
2428 const auto& edi = ed.edpar[nr_edi];
2429 /* We always need the last reference and average positions such that
2430 * in the next time step we can make the ED group whole again
2431 * if the atoms do not have the correct PBC representation */
2432 if (EDstate->bFromCpt)
2434 /* Copy the last whole positions of reference and average group from .cpt */
2435 for (int i = 0; i < edi.sref.nr; i++)
2437 copy_rvec(EDstate->old_sref[nr_edi][i], edi.sref.x_old[i]);
2439 for (int i = 0; i < edi.sav.nr; i++)
2441 copy_rvec(EDstate->old_sav[nr_edi][i], edi.sav.x_old[i]);
2446 EDstate->nref[nr_edi] = edi.sref.nr;
2447 EDstate->nav[nr_edi] = edi.sav.nr;
2450 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2451 EDstate->old_sref_p[nr_edi] = edi.sref.x_old;
2452 EDstate->old_sav_p[nr_edi] = edi.sav.x_old;
2457 /* Adds 'buf' to 'str' */
2458 static void add_to_string(char** str, const char* buf)
2463 len = strlen(*str) + strlen(buf) + 1;
2469 static void add_to_string_aligned(char** str, const char* buf)
2471 char buf_aligned[STRLEN];
2473 sprintf(buf_aligned, EDcol_sfmt, buf);
2474 add_to_string(str, buf_aligned);
2478 static void nice_legend(const char*** setname,
2485 auto tmp = gmx::formatString("%c %s", EDgroupchar, value);
2486 add_to_string_aligned(LegendStr, tmp.c_str());
2487 tmp += gmx::formatString(" (%s)", unit);
2488 (*setname)[*nsets] = gmx_strdup(tmp.c_str());
2493 static void nice_legend_evec(const char*** setname,
2504 for (i = 0; i < evec->neig; i++)
2506 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2507 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2512 /* Makes a legend for the xvg output file. Call on MASTER only! */
2513 static void write_edo_legend(gmx_edsam* ed, int nED, const gmx_output_env_t* oenv)
2516 int nr_edi, nsets, n_flood, n_edsam;
2517 const char** setname;
2519 char* LegendStr = nullptr;
2522 auto edi = ed->edpar.begin();
2524 fprintf(ed->edo, "# Output will be written every %d step%s\n", edi->outfrq,
2525 edi->outfrq != 1 ? "s" : "");
2527 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2529 fprintf(ed->edo, "#\n");
2530 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n",
2531 get_EDgroupChar(nr_edi, nED));
2532 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2533 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig,
2534 edi->vecs.mon.neig != 1 ? "s" : "");
2535 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig,
2536 edi->vecs.linfix.neig != 1 ? "s" : "");
2537 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig,
2538 edi->vecs.linacc.neig != 1 ? "s" : "");
2539 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig,
2540 edi->vecs.radfix.neig != 1 ? "s" : "");
2541 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig,
2542 edi->vecs.radacc.neig != 1 ? "s" : "");
2543 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig,
2544 edi->vecs.radcon.neig != 1 ? "s" : "");
2545 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig,
2546 edi->flood.vecs.neig != 1 ? "s" : "");
2548 if (edi->flood.vecs.neig)
2550 /* If in any of the groups we find a flooding vector, flooding is turned on */
2551 ed->eEDtype = EssentialDynamicsType::Flooding;
2553 /* Print what flavor of flooding we will do */
2554 if (0 == edi->flood.tau) /* constant flooding strength */
2556 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2557 if (edi->flood.bHarmonic)
2559 fprintf(ed->edo, ", harmonic");
2562 else /* adaptive flooding */
2564 fprintf(ed->edo, ", adaptive");
2567 fprintf(ed->edo, "\n");
2572 /* Print a nice legend */
2574 LegendStr[0] = '\0';
2575 sprintf(buf, "# %6s", "time");
2576 add_to_string(&LegendStr, buf);
2578 /* Calculate the maximum number of columns we could end up with */
2579 edi = ed->edpar.begin();
2581 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2583 nsets += 5 + edi->vecs.mon.neig + edi->vecs.linfix.neig + edi->vecs.linacc.neig
2584 + edi->vecs.radfix.neig + edi->vecs.radacc.neig + edi->vecs.radcon.neig
2585 + 6 * edi->flood.vecs.neig;
2588 snew(setname, nsets);
2590 /* In the mdrun time step in a first function call (do_flood()) the flooding
2591 * forces are calculated and in a second function call (do_edsam()) the
2592 * ED constraints. To get a corresponding legend, we need to loop twice
2593 * over the edi groups and output first the flooding, then the ED part */
2595 /* The flooding-related legend entries, if flooding is done */
2597 if (EssentialDynamicsType::Flooding == ed->eEDtype)
2599 edi = ed->edpar.begin();
2600 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2602 /* Always write out the projection on the flooding EVs. Of course, this can also
2603 * be achieved with the monitoring option in do_edsam() (if switched on by the
2604 * user), but in that case the positions need to be communicated in do_edsam(),
2605 * which is not necessary when doing flooding only. */
2606 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2608 for (i = 0; i < edi->flood.vecs.neig; i++)
2610 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2611 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2613 /* Output the current reference projection if it changes with time;
2614 * this can happen when flooding is used as harmonic restraint */
2615 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2617 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2618 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2621 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2622 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2624 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2625 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2628 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2629 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2631 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2633 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2634 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2637 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2638 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2642 } /* End of flooding-related legend entries */
2646 /* Now the ED-related entries, if essential dynamics is done */
2647 edi = ed->edpar.begin();
2648 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2650 if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2652 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2654 /* Essential dynamics, projections on eigenvectors */
2655 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon,
2656 get_EDgroupChar(nr_edi, nED), "MON");
2657 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix,
2658 get_EDgroupChar(nr_edi, nED), "LINFIX");
2659 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc,
2660 get_EDgroupChar(nr_edi, nED), "LINACC");
2661 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix,
2662 get_EDgroupChar(nr_edi, nED), "RADFIX");
2663 if (edi->vecs.radfix.neig)
2665 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm",
2666 get_EDgroupChar(nr_edi, nED));
2668 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc,
2669 get_EDgroupChar(nr_edi, nED), "RADACC");
2670 if (edi->vecs.radacc.neig)
2672 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm",
2673 get_EDgroupChar(nr_edi, nED));
2675 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon,
2676 get_EDgroupChar(nr_edi, nED), "RADCON");
2677 if (edi->vecs.radcon.neig)
2679 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm",
2680 get_EDgroupChar(nr_edi, nED));
2684 } /* end of 'pure' essential dynamics legend entries */
2685 n_edsam = nsets - n_flood;
2687 xvgr_legend(ed->edo, nsets, setname, oenv);
2692 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2693 n_flood, 1 == n_flood ? "" : "s", n_edsam, 1 == n_edsam ? "" : "s");
2694 fprintf(ed->edo, "%s", LegendStr);
2701 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2702 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2703 std::unique_ptr<gmx::EssentialDynamics> init_edsam(const gmx::MDLogger& mdlog,
2704 const char* ediFileName,
2705 const char* edoFileName,
2706 const gmx_mtop_t* mtop,
2707 const t_inputrec* ir,
2708 const t_commrec* cr,
2709 gmx::Constraints* constr,
2710 const t_state* globalState,
2711 ObservablesHistory* oh,
2712 const gmx_output_env_t* oenv,
2713 const gmx::StartingBehavior startingBehavior)
2716 rvec* x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2717 rvec * xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2718 rvec fit_transvec; /* translation ... */
2719 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2720 rvec* ref_x_old = nullptr; /* helper pointer */
2725 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2727 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName))
2730 "The checkpoint file you provided is from an essential dynamics or flooding\n"
2731 "simulation. Please also set the .edi file on the command line with -ei.\n");
2737 "Activating essential dynamics via files passed to "
2738 "gmx mdrun -edi will change in future to be files passed to "
2739 "gmx grompp and the related .mdp options may change also.");
2741 /* Open input and output files, allocate space for ED data structure */
2742 auto edHandle = ed_open(mtop->natoms, oh, ediFileName, edoFileName, startingBehavior, oenv, cr);
2743 auto ed = edHandle->getLegacyED();
2744 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2745 constr->saveEdsamPointer(ed);
2747 /* Needed for initializing radacc radius in do_edsam */
2750 /* The input file is read by the master and the edi structures are
2751 * initialized here. Input is stored in ed->edpar. Then the edi
2752 * structures are transferred to the other nodes */
2755 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2756 * flooding vector, Essential dynamics can be applied to more than one structure
2757 * as well, but will be done in the order given in the edi file, so
2758 * expect different results for different order of edi file concatenation! */
2759 for (auto& edi : ed->edpar)
2761 init_edi(mtop, &edi);
2762 init_flood(&edi, ed, ir->delta_t);
2766 /* The master does the work here. The other nodes get the positions
2767 * not before dd_partition_system which is called after init_edsam */
2770 edsamhistory_t* EDstate = oh->edsamHistory.get();
2772 if (!EDstate->bFromCpt)
2774 /* Remove PBC, make molecule(s) subject to ED whole. */
2775 snew(x_pbc, mtop->natoms);
2776 copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop->natoms);
2777 do_pbc_first_mtop(nullptr, ir->pbcType, globalState->box, mtop, x_pbc);
2779 /* Reset pointer to first ED data set which contains the actual ED data */
2780 auto edi = ed->edpar.begin();
2781 GMX_ASSERT(!ed->edpar.empty(), "Essential dynamics parameters is undefined");
2783 /* Loop over all ED/flooding data sets (usually only one, though) */
2784 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2786 /* For multiple ED groups we use the output frequency that was specified
2787 * in the first set */
2790 edi->outfrq = ed->edpar.begin()->outfrq;
2793 /* Extract the initial reference and average positions. When starting
2794 * from .cpt, these have already been read into sref.x_old
2795 * in init_edsamstate() */
2796 if (!EDstate->bFromCpt)
2798 /* If this is the first run (i.e. no checkpoint present) we assume
2799 * that the starting positions give us the correct PBC representation */
2800 for (i = 0; i < edi->sref.nr; i++)
2802 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2805 for (i = 0; i < edi->sav.nr; i++)
2807 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2811 /* Now we have the PBC-correct start positions of the reference and
2812 average structure. We copy that over to dummy arrays on which we
2813 can apply fitting to print out the RMSD. We srenew the memory since
2814 the size of the buffers is likely different for every ED group */
2815 srenew(xfit, edi->sref.nr);
2816 srenew(xstart, edi->sav.nr);
2819 /* Reference indices are the same as average indices */
2820 ref_x_old = edi->sav.x_old;
2824 ref_x_old = edi->sref.x_old;
2826 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2827 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2829 /* Make the fit to the REFERENCE structure, get translation and rotation */
2830 fit_to_reference(xfit, fit_transvec, fit_rotmat, &(*edi));
2832 /* Output how well we fit to the reference at the start */
2833 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2834 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2835 rmsd_from_structure(xfit, &edi->sref));
2836 if (EDstate->nED > 1)
2838 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2840 fprintf(stderr, "\n");
2842 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2843 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2845 /* calculate initial projections */
2846 project(xstart, &(*edi));
2848 /* For the target and origin structure both a reference (fit) and an
2849 * average structure can be provided in make_edi. If both structures
2850 * are the same, make_edi only stores one of them in the .edi file.
2851 * If they differ, first the fit and then the average structure is stored
2852 * in star (or sor), thus the number of entries in star/sor is
2853 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2854 * the size of the average group. */
2856 /* process target structure, if required */
2857 if (edi->star.nr > 0)
2859 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2861 /* get translation & rotation for fit of target structure to reference structure */
2862 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, &(*edi));
2864 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2865 if (edi->star.nr == edi->sav.nr)
2869 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2871 /* The last sav.nr indices of the target structure correspond to
2872 * the average structure, which must be projected */
2873 avindex = edi->star.nr - edi->sav.nr;
2875 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2879 rad_project(*edi, xstart, &edi->vecs.radcon);
2882 /* process structure that will serve as origin of expansion circle */
2883 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2886 "ED: Setting center of flooding potential (0 = average structure)\n");
2889 if (edi->sori.nr > 0)
2891 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2893 /* fit this structure to reference structure */
2894 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, &(*edi));
2896 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2897 if (edi->sori.nr == edi->sav.nr)
2901 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2903 /* For the projection, we need the last sav.nr indices of sori */
2904 avindex = edi->sori.nr - edi->sav.nr;
2907 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2908 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2909 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2912 "ED: The ORIGIN structure will define the flooding potential "
2914 /* Set center of flooding potential to the ORIGIN structure */
2915 rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2916 /* We already know that no (moving) reference position was provided,
2917 * therefore we can overwrite refproj[0]*/
2918 copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
2921 else /* No origin structure given */
2923 rad_project(*edi, xstart, &edi->vecs.radacc);
2924 rad_project(*edi, xstart, &edi->vecs.radfix);
2925 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2927 if (edi->flood.bHarmonic)
2930 "ED: A (possibly changing) ref. projection will define the "
2931 "flooding potential center.\n");
2932 for (i = 0; i < edi->flood.vecs.neig; i++)
2934 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
2940 "ED: The AVERAGE structure will define the flooding potential "
2942 /* Set center of flooding potential to the center of the covariance matrix,
2943 * i.e. the average structure, i.e. zero in the projected system */
2944 for (i = 0; i < edi->flood.vecs.neig; i++)
2946 edi->flood.vecs.refproj[i] = 0.0;
2951 /* For convenience, output the center of the flooding potential for the eigenvectors */
2952 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2954 for (i = 0; i < edi->flood.vecs.neig; i++)
2956 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e",
2957 edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2958 if (edi->flood.bHarmonic)
2960 fprintf(stdout, " (adding %11.4e/timestep)",
2961 edi->flood.referenceProjectionSlope[i]);
2963 fprintf(stdout, "\n");
2967 /* set starting projections for linsam */
2968 rad_project(*edi, xstart, &edi->vecs.linacc);
2969 rad_project(*edi, xstart, &edi->vecs.linfix);
2971 /* Prepare for the next edi data set: */
2974 /* Cleaning up on the master node: */
2975 if (!EDstate->bFromCpt)
2982 } /* end of MASTER only section */
2986 /* Broadcast the essential dynamics / flooding data to all nodes */
2987 broadcast_ed_data(cr, ed);
2991 /* In the single-CPU case, point the local atom numbers pointers to the global
2992 * one, so that we can use the same notation in serial and parallel case: */
2993 /* Loop over all ED data sets (usually only one, though) */
2994 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
2996 edi->sref.anrs_loc = edi->sref.anrs;
2997 edi->sav.anrs_loc = edi->sav.anrs;
2998 edi->star.anrs_loc = edi->star.anrs;
2999 edi->sori.anrs_loc = edi->sori.anrs;
3000 /* For the same reason as above, make a dummy c_ind array: */
3001 snew(edi->sav.c_ind, edi->sav.nr);
3002 /* Initialize the array */
3003 for (i = 0; i < edi->sav.nr; i++)
3005 edi->sav.c_ind[i] = i;
3007 /* In the general case we will need a different-sized array for the reference indices: */
3010 snew(edi->sref.c_ind, edi->sref.nr);
3011 for (i = 0; i < edi->sref.nr; i++)
3013 edi->sref.c_ind[i] = i;
3016 /* Point to the very same array in case of other structures: */
3017 edi->star.c_ind = edi->sav.c_ind;
3018 edi->sori.c_ind = edi->sav.c_ind;
3019 /* In the serial case, the local number of atoms is the global one: */
3020 edi->sref.nr_loc = edi->sref.nr;
3021 edi->sav.nr_loc = edi->sav.nr;
3022 edi->star.nr_loc = edi->star.nr;
3023 edi->sori.nr_loc = edi->sori.nr;
3027 /* Allocate space for ED buffer variables */
3028 /* Again, loop over ED data sets */
3029 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
3031 /* Allocate space for ED buffer variables */
3032 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
3033 snew(edi->buf->do_edsam, 1);
3035 /* Space for collective ED buffer variables */
3037 /* Collective positions of atoms with the average indices */
3038 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
3039 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
3040 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
3041 /* Collective positions of atoms with the reference indices */
3044 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
3045 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
3046 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
3049 /* Get memory for flooding forces */
3050 snew(edi->flood.forces_cartesian, edi->sav.nr);
3053 /* Flush the edo file so that the user can check some things
3054 * when the simulation has started */
3064 void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[], rvec v[], const matrix box, gmx_edsam* ed)
3066 int i, edinr, iupdate = 500;
3067 matrix rotmat; /* rotation matrix */
3068 rvec transvec; /* translation vector */
3069 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3070 real dt_1; /* 1/dt */
3071 struct t_do_edsam* buf;
3072 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3073 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3076 /* Check if ED sampling has to be performed */
3077 if (ed->eEDtype == EssentialDynamicsType::None)
3082 dt_1 = 1.0 / ir->delta_t;
3084 /* Loop over all ED groups (usually one) */
3086 for (auto& edi : ed->edpar)
3089 if (bNeedDoEdsam(edi))
3092 buf = edi.buf->do_edsam;
3096 /* initialize radacc radius for slope criterion */
3097 buf->oldrad = calc_radius(edi.vecs.radacc);
3100 /* Copy the positions into buf->xc* arrays and after ED
3101 * feed back corrections to the official positions */
3103 /* Broadcast the ED positions such that every node has all of them
3104 * Every node contributes its local positions xs and stores it in
3105 * the collective buf->xcoll array. Note that for edinr > 1
3106 * xs could already have been modified by an earlier ED */
3108 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll,
3109 PAR(cr) ? buf->bUpdateShifts : TRUE, xs, edi.sav.nr, edi.sav.nr_loc,
3110 edi.sav.anrs_loc, edi.sav.c_ind, edi.sav.x_old, box);
3112 /* Only assembly reference positions if their indices differ from the average ones */
3115 communicate_group_positions(
3116 cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref,
3117 PAR(cr) ? buf->bUpdateShifts : TRUE, xs, edi.sref.nr, edi.sref.nr_loc,
3118 edi.sref.anrs_loc, edi.sref.c_ind, edi.sref.x_old, box);
3121 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3122 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3123 * set bUpdateShifts=TRUE in the parallel case. */
3124 buf->bUpdateShifts = FALSE;
3126 /* Now all nodes have all of the ED positions in edi.sav->xcoll,
3127 * as well as the indices in edi.sav.anrs */
3129 /* Fit the reference indices to the reference structure */
3132 fit_to_reference(buf->xcoll, transvec, rotmat, &edi);
3136 fit_to_reference(buf->xc_ref, transvec, rotmat, &edi);
3139 /* Now apply the translation and rotation to the ED structure */
3140 translate_and_rotate(buf->xcoll, edi.sav.nr, transvec, rotmat);
3142 /* Find out how well we fit to the reference (just for output steps) */
3143 if (do_per_step(step, edi.outfrq) && MASTER(cr))
3147 /* Indices of reference and average structures are identical,
3148 * thus we can calculate the rmsd to SREF using xcoll */
3149 rmsdev = rmsd_from_structure(buf->xcoll, &edi.sref);
3153 /* We have to translate & rotate the reference atoms first */
3154 translate_and_rotate(buf->xc_ref, edi.sref.nr, transvec, rotmat);
3155 rmsdev = rmsd_from_structure(buf->xc_ref, &edi.sref);
3159 /* update radsam references, when required */
3160 if (do_per_step(step, edi.maxedsteps) && step >= edi.presteps)
3162 project(buf->xcoll, &edi);
3163 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3164 rad_project(edi, buf->xcoll, &edi.vecs.radfix);
3165 buf->oldrad = -1.e5;
3168 /* update radacc references, when required */
3169 if (do_per_step(step, iupdate) && step >= edi.presteps)
3171 edi.vecs.radacc.radius = calc_radius(edi.vecs.radacc);
3172 if (edi.vecs.radacc.radius - buf->oldrad < edi.slope)
3174 project(buf->xcoll, &edi);
3175 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3180 buf->oldrad = edi.vecs.radacc.radius;
3184 /* apply the constraints */
3185 if (step >= edi.presteps && ed_constraints(ed->eEDtype, edi))
3187 /* ED constraints should be applied already in the first MD step
3188 * (which is step 0), therefore we pass step+1 to the routine */
3189 ed_apply_constraints(buf->xcoll, &edi, step + 1 - ir->init_step);
3192 /* write to edo, when required */
3193 if (do_per_step(step, edi.outfrq))
3195 project(buf->xcoll, &edi);
3196 if (MASTER(cr) && !bSuppress)
3198 write_edo(edi, ed->edo, rmsdev);
3202 /* Copy back the positions unless monitoring only */
3203 if (ed_constraints(ed->eEDtype, edi))
3205 /* remove fitting */
3206 rmfit(edi.sav.nr, buf->xcoll, transvec, rotmat);
3208 /* Copy the ED corrected positions into the coordinate array */
3209 /* Each node copies its local part. In the serial case, nat_loc is the
3210 * total number of ED atoms */
3211 for (i = 0; i < edi.sav.nr_loc; i++)
3213 /* Unshift local ED coordinate and store in x_unsh */
3214 ed_unshift_single_coord(box, buf->xcoll[edi.sav.c_ind[i]],
3215 buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh);
3217 /* dx is the ED correction to the positions: */
3218 rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx);
3222 /* dv is the ED correction to the velocity: */
3223 svmul(dt_1, dx, dv);
3224 /* apply the velocity correction: */
3225 rvec_inc(v[edi.sav.anrs_loc[i]], dv);
3227 /* Finally apply the position correction due to ED: */
3228 copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]);
3231 } /* END of if ( bNeedDoEdsam(edi) ) */
3233 /* Prepare for the next ED group */
3235 } /* END of loop over ED groups */