2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
46 #include "gromacs/commandline/filenm.h"
47 #include "gromacs/compat/make_unique.h"
48 #include "gromacs/domdec/domdec_struct.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxlib/network.h"
52 #include "gromacs/gmxlib/nrnb.h"
53 #include "gromacs/linearalgebra/nrjac.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/math/vectypes.h"
58 #include "gromacs/mdlib/broadcaststructs.h"
59 #include "gromacs/mdlib/constr.h"
60 #include "gromacs/mdlib/groupcoord.h"
61 #include "gromacs/mdlib/mdrun.h"
62 #include "gromacs/mdlib/sim_util.h"
63 #include "gromacs/mdlib/update.h"
64 #include "gromacs/mdtypes/commrec.h"
65 #include "gromacs/mdtypes/edsamhistory.h"
66 #include "gromacs/mdtypes/inputrec.h"
67 #include "gromacs/mdtypes/md_enums.h"
68 #include "gromacs/mdtypes/observableshistory.h"
69 #include "gromacs/mdtypes/state.h"
70 #include "gromacs/pbcutil/pbc.h"
71 #include "gromacs/topology/mtop_lookup.h"
72 #include "gromacs/topology/topology.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/gmxassert.h"
76 #include "gromacs/utility/logger.h"
77 #include "gromacs/utility/smalloc.h"
78 #include "gromacs/utility/strconvert.h"
82 /*! \brief Identifies the type of ED: none, normal ED, flooding. */
83 enum class EssentialDynamicsType
85 //! No essential dynamics
87 //! Essential dynamics sampling
89 //! Essential dynamics flooding
93 /*! \brief Identify on which structure to operate. */
94 enum class EssentialDynamicsStructure
96 //! Reference structure
106 /*! \brief Essential dynamics vector.
107 * TODO: split into working data and input data
108 * NOTE: called eigvec, because traditionally eigenvectors from PCA analysis
109 * were used as essential dynamics vectors, however, vectors used for ED need
110 * not have any special properties
114 //! nr of eigenvectors
116 //! index nrs of eigenvectors
118 //! stepsizes (per eigenvector)
119 real *stpsz = nullptr;
120 //! eigenvector components
121 rvec **vec = nullptr;
122 //! instantaneous x projections
123 real *xproj = nullptr;
124 //! instantaneous f projections
125 real *fproj = nullptr;
126 //! instantaneous radius
128 //! starting or target projections
129 real *refproj = nullptr;
132 /*! \brief Essential dynamics vectors per method implementation.
136 //! only monitored, no constraints
138 //! fixed linear constraints
139 t_eigvec linfix = {};
140 //! acceptance linear constraints
141 t_eigvec linacc = {};
142 //! fixed radial constraints (exp)
143 t_eigvec radfix = {};
144 //! acceptance radial constraints (exp)
145 t_eigvec radacc = {};
146 //! acceptance rad. contraction constr.
147 t_eigvec radcon = {};
150 /*! \brief Essential dynamics flooding parameters and work data.
151 * TODO: split into working data and input parameters
152 * NOTE: The implementation here follows:
153 * O.E. Lange, L.V. Schafer, and H. Grubmuller,
154 * “Flooding in GROMACS: Accelerated barrier crossings in molecular dynamics,”
155 * J. Comp. Chem., 27 1693–1702 (2006)
159 /*! \brief Target destabilisation free energy.
160 * Controls the flooding potential strength.
161 * Linked to the expected speed-up of mean first passage time out of flooded minimum */
163 //! Do not calculate a flooding potential, instead flood with a constant force
164 bool bConstForce = false;
165 //! Relaxation time scale for slowest flooded degree of freedom
167 //! Current estimated destabilisation free energy
169 //! Flooding energy, acting as a proportionality factor for the flooding potential
171 //! Boltzmann constant times temperature, provided by user
173 //! The flooding potential
175 //! Integration time step
177 //! Inital flooding strenth
179 //! Empirical global scaling parameter of the essential dynamics vectors.
181 //! The forces from flooding in atom coordinate space (in contrast to projections onto essential dynamics vectors)
182 rvec *forces_cartesian = nullptr;
183 //! The vectors along which to flood
185 //! Use flooding for harmonic restraint on the eigenvector
186 bool bHarmonic = false;
187 //! The initial reference projection of the flooding vectors. Only with harmonic restraint.
188 real *initialReferenceProjection = nullptr;
189 //! The current reference projection is the initialReferenceProjection + step * slope. Only with harmonic restraint.
190 real *referenceProjectionSlope = nullptr;
195 /* This type is for the average, reference, target, and origin structure */
198 int nr = 0; /* number of atoms this structure contains */
199 int nr_loc = 0; /* number of atoms on local node */
200 int *anrs = nullptr; /* atom index numbers */
201 int *anrs_loc = nullptr; /* local atom index numbers */
202 int nalloc_loc = 0; /* allocation size of anrs_loc */
203 int *c_ind = nullptr; /* at which position of the whole anrs
204 * array is a local atom?, i.e.
205 * c_ind[0...nr_loc-1] gives the atom index
206 * with respect to the collective
207 * anrs[0...nr-1] array */
208 rvec *x = nullptr; /* positions for this structure */
209 rvec *x_old = nullptr; /* Last positions which have the correct PBC
210 representation of the ED group. In
211 combination with keeping track of the
212 shift vectors, the ED group can always
214 real *m = nullptr; /* masses */
215 real mtot = 0.; /* total mass (only used in sref) */
216 real *sqrtm = nullptr; /* sqrt of the masses used for mass-
217 * weighting of analysis (only used in sav) */
223 int nini = 0; /* total Nr of atoms */
224 gmx_bool fitmas = false; /* true if trans fit with cm */
225 gmx_bool pcamas = false; /* true if mass-weighted PCA */
226 int presteps = 0; /* number of steps to run without any
227 * perturbations ... just monitoring */
228 int outfrq = 0; /* freq (in steps) of writing to edo */
229 int maxedsteps = 0; /* max nr of steps per cycle */
231 /* all gmx_edx datasets are copied to all nodes in the parallel case */
232 gmx_edx sref = {}; /* reference positions, to these fitting
234 gmx_bool bRefEqAv = false; /* If true, reference & average indices
235 * are the same. Used for optimization */
236 gmx_edx sav = {}; /* average positions */
237 gmx_edx star = {}; /* target positions */
238 gmx_edx sori = {}; /* origin positions */
240 t_edvecs vecs = {}; /* eigenvectors */
241 real slope = 0; /* minimal slope in acceptance radexp */
243 t_edflood flood = {}; /* parameters especially for flooding */
244 struct t_ed_buffer *buf = nullptr; /* handle to local buffers */
252 EssentialDynamicsType eEDtype = EssentialDynamicsType::None;
253 //! output file pointer
255 std::vector<t_edpar> edpar;
256 gmx_bool bFirst = false;
258 gmx_edsam::~gmx_edsam()
262 /* edo is opened sometimes with xvgropen, sometimes with
263 * gmx_fio_fopen, so we use the least common denominator for
273 rvec old_transvec, older_transvec, transvec_compact;
274 rvec *xcoll; /* Positions from all nodes, this is the
275 collective set we work on.
276 These are the positions of atoms with
277 average structure indices */
278 rvec *xc_ref; /* same but with reference structure indices */
279 ivec *shifts_xcoll; /* Shifts for xcoll */
280 ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
281 ivec *shifts_xc_ref; /* Shifts for xc_ref */
282 ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
283 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
284 ED shifts for this ED group need to
289 /* definition of ED buffer structure */
292 struct t_fit_to_ref * fit_to_ref;
293 struct t_do_edfit * do_edfit;
294 struct t_do_edsam * do_edsam;
295 struct t_do_radcon * do_radcon;
300 class EssentialDynamics::Impl
303 // TODO: move all fields from gmx_edsam here and remove gmx_edsam
304 gmx_edsam essentialDynamics_;
306 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)
357 || (edi.vecs.linfix.neig != 0)
358 || (edi.vecs.linacc.neig != 0)
359 || (edi.vecs.radfix.neig != 0)
360 || (edi.vecs.radacc.neig != 0)
361 || (edi.vecs.radcon.neig != 0);
366 /* Multiple ED groups will be labeled with letters instead of numbers
367 * to avoid confusion with eigenvector indices */
368 static char get_EDgroupChar(int nr_edi, int nED)
376 * nr_edi = 2 -> B ...
378 return 'A' + nr_edi - 1;
383 /*! \brief The mass-weighted inner product of two coordinate vectors.
384 * Does not subtract average positions, projection on single eigenvector is returned
385 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
386 * Average position is subtracted in ed_apply_constraints prior to calling projectx
387 * \param[in] edi Essential dynamics parameters
388 * \param[in] xcoll vector of atom coordinates
389 * \param[in] vec vector of coordinates to project onto
390 * \return mass-weighted projection.
392 real projectx(const t_edpar &edi, rvec *xcoll, rvec *vec)
398 for (i = 0; i < edi.sav.nr; i++)
400 proj += edi.sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
405 /*!\brief Project coordinates onto vector after substracting average position.
406 * projection is stored in vec->refproj which is used for radacc, radfix,
407 * radcon and center of flooding potential.
408 * Average positions are first substracted from x, then added back again.
409 * \param[in] edi essential dynamics parameters with average position
410 * \param[in] x Coordinates to be projected
411 * \param[out] vec eigenvector, radius and refproj are overwritten here
413 void rad_project(const t_edpar &edi, rvec *x, t_eigvec *vec)
418 /* Subtract average positions */
419 for (i = 0; i < edi.sav.nr; i++)
421 rvec_dec(x[i], edi.sav.x[i]);
424 for (i = 0; i < vec->neig; i++)
426 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
427 rad += gmx::square((vec->refproj[i]-vec->xproj[i]));
429 vec->radius = sqrt(rad);
431 /* Add average positions */
432 for (i = 0; i < edi.sav.nr; i++)
434 rvec_inc(x[i], edi.sav.x[i]);
438 /*!\brief Projects coordinates onto eigenvectors and stores result in vec->xproj.
439 * Mass-weighting is applied. Subtracts average positions prior to projection and add
440 * them afterwards to retain the unchanged vector.
441 * \param[in] x The coordinates to project to an eigenvector
442 * \param[in,out] vec The eigenvectors
443 * \param[in] edi essential dynamics parameters holding average structure and masses
445 void project_to_eigvectors(rvec *x, t_eigvec *vec, const t_edpar &edi)
452 /* Subtract average positions */
453 for (int i = 0; i < edi.sav.nr; i++)
455 rvec_dec(x[i], edi.sav.x[i]);
458 for (int i = 0; i < vec->neig; i++)
460 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
463 /* Add average positions */
464 for (int i = 0; i < edi.sav.nr; i++)
466 rvec_inc(x[i], edi.sav.x[i]);
471 /* Project vector x onto all edi->vecs (mon, linfix,...) */
472 static void project(rvec *x, /* positions to project */
473 t_edpar *edi) /* edi data set */
475 /* It is not more work to subtract the average position in every
476 * subroutine again, because these routines are rarely used simultaneously */
477 project_to_eigvectors(x, &edi->vecs.mon, *edi);
478 project_to_eigvectors(x, &edi->vecs.linfix, *edi);
479 project_to_eigvectors(x, &edi->vecs.linacc, *edi);
480 project_to_eigvectors(x, &edi->vecs.radfix, *edi);
481 project_to_eigvectors(x, &edi->vecs.radacc, *edi);
482 project_to_eigvectors(x, &edi->vecs.radcon, *edi);
487 /*!\brief Evaluates the distance from reference to current eigenvector projection.
488 * \param[in] vec eigenvector
491 real calc_radius(const t_eigvec &vec)
495 for (int i = 0; i < vec.neig; i++)
497 rad += gmx::square((vec.refproj[i]-vec.xproj[i]));
500 return rad = sqrt(rad);
509 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
511 /* this is a copy of do_fit with some modifications */
512 int c, r, n, j, i, irot;
513 double d[6], xnr, xpc;
518 struct t_do_edfit *loc;
521 if (edi->buf->do_edfit != nullptr)
528 snew(edi->buf->do_edfit, 1);
530 loc = edi->buf->do_edfit;
534 snew(loc->omega, 2*DIM);
535 snew(loc->om, 2*DIM);
536 for (i = 0; i < 2*DIM; i++)
538 snew(loc->omega[i], 2*DIM);
539 snew(loc->om[i], 2*DIM);
543 for (i = 0; (i < 6); i++)
546 for (j = 0; (j < 6); j++)
548 loc->omega[i][j] = 0;
553 /* calculate the matrix U */
555 for (n = 0; (n < natoms); n++)
557 for (c = 0; (c < DIM); c++)
560 for (r = 0; (r < DIM); r++)
568 /* construct loc->omega */
569 /* loc->omega is symmetric -> loc->omega==loc->omega' */
570 for (r = 0; (r < 6); r++)
572 for (c = 0; (c <= r); c++)
574 if ((r >= 3) && (c < 3))
576 loc->omega[r][c] = u[r-3][c];
577 loc->omega[c][r] = u[r-3][c];
581 loc->omega[r][c] = 0;
582 loc->omega[c][r] = 0;
587 /* determine h and k */
588 jacobi(loc->omega, 6, d, loc->om, &irot);
592 fprintf(stderr, "IROT=0\n");
595 index = 0; /* For the compiler only */
597 for (j = 0; (j < 3); j++)
600 for (i = 0; (i < 6); i++)
609 for (i = 0; (i < 3); i++)
611 vh[j][i] = M_SQRT2*loc->om[i][index];
612 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
617 for (c = 0; (c < 3); c++)
619 for (r = 0; (r < 3); r++)
621 R[c][r] = vk[0][r]*vh[0][c]+
628 for (c = 0; (c < 3); c++)
630 for (r = 0; (r < 3); r++)
632 R[c][r] = vk[0][r]*vh[0][c]+
641 static void rmfit(int nat, rvec *xcoll, const rvec transvec, matrix rotmat)
648 * The inverse rotation is described by the transposed rotation matrix */
649 transpose(rotmat, tmat);
650 rotate_x(xcoll, nat, tmat);
652 /* Remove translation */
653 vec[XX] = -transvec[XX];
654 vec[YY] = -transvec[YY];
655 vec[ZZ] = -transvec[ZZ];
656 translate_x(xcoll, nat, vec);
660 /**********************************************************************************
661 ******************** FLOODING ****************************************************
662 **********************************************************************************
664 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
665 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
666 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
668 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
669 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
671 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
672 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
674 do_flood makes a copy of the positions,
675 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
676 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
677 fit. Then do_flood adds these forces to the forcefield-forces
678 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
680 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
681 structure is projected to the system of eigenvectors and then this position in the subspace is used as
682 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
683 i.e. the average structure as given in the make_edi file.
685 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
686 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
687 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
688 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
689 further adaption is applied, Efl will stay constant at zero.
691 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
692 used as spring constants for the harmonic potential.
693 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
694 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
696 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
697 the routine read_edi_file reads all of theses flooding files.
698 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
699 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
700 edi there is no interdependence whatsoever. The forces are added together.
702 To write energies into the .edr file, call the function
703 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
705 get_flood_energies(real Vfl[],int nnames);
708 - 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.
710 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
711 two edsam files from two peptide chains
714 // TODO split this into multiple files
717 /*!\brief Output flooding simulation settings and results to file.
718 * \param[in] edi Essential dynamics input parameters
719 * \param[in] fp output file
720 * \param[in] rmsd rmsd to reference structure
722 void write_edo_flood(const t_edpar &edi, FILE *fp, real rmsd)
724 /* Output how well we fit to the reference structure */
725 fprintf(fp, EDcol_ffmt, rmsd);
727 for (int i = 0; i < edi.flood.vecs.neig; i++)
729 fprintf(fp, EDcol_efmt, edi.flood.vecs.xproj[i]);
731 /* Check whether the reference projection changes with time (this can happen
732 * in case flooding is used as harmonic restraint). If so, output the
733 * current reference projection */
734 if (edi.flood.bHarmonic && edi.flood.referenceProjectionSlope[i] != 0.0)
736 fprintf(fp, EDcol_efmt, edi.flood.vecs.refproj[i]);
739 /* Output Efl if we are doing adaptive flooding */
740 if (0 != edi.flood.tau)
742 fprintf(fp, EDcol_efmt, edi.flood.Efl);
744 fprintf(fp, EDcol_efmt, edi.flood.Vfl);
746 /* Output deltaF if we are doing adaptive flooding */
747 if (0 != edi.flood.tau)
749 fprintf(fp, EDcol_efmt, edi.flood.deltaF);
751 fprintf(fp, EDcol_efmt, edi.flood.vecs.fproj[i]);
757 /* From flood.xproj compute the Vfl(x) at this point */
758 static real flood_energy(t_edpar *edi, int64_t step)
760 /* compute flooding energy Vfl
761 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
762 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
763 it is already computed by make_edi and stored in stpsz[i]
765 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
772 /* Each time this routine is called (i.e. each time step), we add a small
773 * value to the reference projection. This way a harmonic restraint towards
774 * a moving reference is realized. If no value for the additive constant
775 * is provided in the edi file, the reference will not change. */
776 if (edi->flood.bHarmonic)
778 for (i = 0; i < edi->flood.vecs.neig; i++)
780 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i] + step * edi->flood.referenceProjectionSlope[i];
785 /* Compute sum which will be the exponent of the exponential */
786 for (i = 0; i < edi->flood.vecs.neig; i++)
788 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
789 sum += edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i])*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
792 /* Compute the Gauss function*/
793 if (edi->flood.bHarmonic)
795 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
799 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*std::exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
806 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
807 static void flood_forces(t_edpar *edi)
809 /* compute the forces in the subspace of the flooding eigenvectors
810 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
813 real energy = edi->flood.Vfl;
816 if (edi->flood.bHarmonic)
818 for (i = 0; i < edi->flood.vecs.neig; i++)
820 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
825 for (i = 0; i < edi->flood.vecs.neig; i++)
827 /* if Efl is zero the forces are zero if not use the formula */
828 edi->flood.vecs.fproj[i] = edi->flood.Efl != 0 ? edi->flood.kT/edi->flood.Efl/edi->flood.alpha2*energy*edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]) : 0;
835 /*!\brief Raise forces from subspace into cartesian space.
836 * This function lifts the forces from the subspace to the cartesian space
837 * all the values not contained in the subspace are assumed to be zero and then
838 * a coordinate transformation from eigenvector to cartesian vectors is performed
839 * The nonexistent values don't have to be set to zero explicitly, they would occur
840 * as zero valued summands, hence we just stop to compute this part of the sum.
841 * For every atom we add all the contributions to this atom from all the different eigenvectors.
842 * NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
843 * field forces_cart prior the computation, but we compute the forces separately
844 * to have them accessible for diagnostics
846 * \param[in] edi Essential dynamics input parameters
847 * \param[out] forces_cart The cartesian forces
850 void flood_blowup(const t_edpar &edi, rvec *forces_cart)
852 const real * forces_sub = edi.flood.vecs.fproj;
853 /* Calculate the cartesian forces for the local atoms */
855 /* Clear forces first */
856 for (int j = 0; j < edi.sav.nr_loc; j++)
858 clear_rvec(forces_cart[j]);
861 /* Now compute atomwise */
862 for (int j = 0; j < edi.sav.nr_loc; j++)
864 /* Compute forces_cart[edi.sav.anrs[j]] */
865 for (int eig = 0; eig < edi.flood.vecs.neig; eig++)
868 /* Force vector is force * eigenvector (compute only atom j) */
869 svmul(forces_sub[eig], edi.flood.vecs.vec[eig][edi.sav.c_ind[j]], addedForce);
870 /* Add this vector to the cartesian forces */
871 rvec_inc(forces_cart[j], addedForce);
879 /* Update the values of Efl, deltaF depending on tau and Vfl */
880 static void update_adaption(t_edpar *edi)
882 /* this function updates the parameter Efl and deltaF according to the rules given in
883 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
886 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
888 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
889 /* check if restrain (inverted flooding) -> don't let EFL become positive */
890 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
895 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
900 static void do_single_flood(
908 gmx_bool bNS) /* Are we in a neighbor searching step? */
911 matrix rotmat; /* rotation matrix */
912 matrix tmat; /* inverse rotation */
913 rvec transvec; /* translation vector */
915 struct t_do_edsam *buf;
918 buf = edi->buf->do_edsam;
921 /* Broadcast the positions of the AVERAGE structure such that they are known on
922 * every processor. Each node contributes its local positions x and stores them in
923 * the collective ED array buf->xcoll */
924 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
925 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
927 /* Only assembly REFERENCE positions if their indices differ from the average ones */
930 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
931 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
934 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
935 * We do not need to update the shifts until the next NS step */
936 buf->bUpdateShifts = FALSE;
938 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
939 * as well as the indices in edi->sav.anrs */
941 /* Fit the reference indices to the reference structure */
944 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
948 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
951 /* Now apply the translation and rotation to the ED structure */
952 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
954 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
955 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
957 if (!edi->flood.bConstForce)
959 /* Compute Vfl(x) from flood.xproj */
960 edi->flood.Vfl = flood_energy(edi, step);
962 update_adaption(edi);
964 /* Compute the flooding forces */
968 /* Translate them into cartesian positions */
969 flood_blowup(*edi, edi->flood.forces_cartesian);
971 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
972 /* Each node rotates back its local forces */
973 transpose(rotmat, tmat);
974 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
976 /* Finally add forces to the main force variable */
977 for (i = 0; i < edi->sav.nr_loc; i++)
979 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
982 /* Output is written by the master process */
983 if (do_per_step(step, edi->outfrq) && MASTER(cr))
985 /* Output how well we fit to the reference */
988 /* Indices of reference and average structures are identical,
989 * thus we can calculate the rmsd to SREF using xcoll */
990 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
994 /* We have to translate & rotate the reference atoms first */
995 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
996 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
999 write_edo_flood(*edi, edo, rmsdev);
1004 /* Main flooding routine, called from do_force */
1005 extern void do_flood(const t_commrec *cr,
1006 const t_inputrec *ir,
1014 /* Write time to edo, when required. Output the time anyhow since we need
1015 * it in the output file for ED constraints. */
1016 if (MASTER(cr) && do_per_step(step, ed->edpar.begin()->outfrq))
1018 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
1021 if (ed->eEDtype != EssentialDynamicsType::Flooding)
1026 for (auto &edi : ed->edpar)
1028 /* Call flooding for one matrix */
1029 if (edi.flood.vecs.neig)
1031 do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS);
1037 /* Called by init_edi, configure some flooding related variables and structures,
1038 * print headers to output files */
1039 static void init_flood(t_edpar *edi, gmx_edsam * ed, real dt)
1044 edi->flood.Efl = edi->flood.constEfl;
1048 if (edi->flood.vecs.neig)
1050 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1051 ed->eEDtype = EssentialDynamicsType::Flooding;
1053 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
1055 if (edi->flood.bConstForce)
1057 /* We have used stpsz as a vehicle to carry the fproj values for constant
1058 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1059 * in const force flooding, fproj is never changed. */
1060 for (i = 0; i < edi->flood.vecs.neig; i++)
1062 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1064 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1065 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1073 /*********** Energy book keeping ******/
1074 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1083 srenew(names, count);
1084 sprintf(buf, "Vfl_%d", count);
1085 names[count-1] = gmx_strdup(buf);
1086 actual = actual->next_edi;
1093 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1095 /*fl has to be big enough to capture nnames-many entries*/
1104 Vfl[count-1] = actual->flood.Vfl;
1105 actual = actual->next_edi;
1108 if (nnames != count-1)
1110 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1113 /************* END of FLOODING IMPLEMENTATION ****************************/
1117 /* This function opens the ED input and output files, reads in all datasets it finds
1118 * in the input file, and cross-checks whether the .edi file information is consistent
1119 * with the essential dynamics data found in the checkpoint file (if present).
1120 * gmx make_edi can be used to create an .edi input file.
1122 static std::unique_ptr<gmx::EssentialDynamics> ed_open(
1124 ObservablesHistory *oh,
1125 const char *ediFileName,
1126 const char *edoFileName,
1128 const gmx_output_env_t *oenv,
1129 const t_commrec *cr)
1131 auto edHandle = gmx::compat::make_unique<gmx::EssentialDynamics>();
1132 auto ed = edHandle->getLegacyED();
1133 /* We want to perform ED (this switch might later be upgraded to EssentialDynamicsType::Flooding) */
1134 ed->eEDtype = EssentialDynamicsType::EDSampling;
1138 // If we start from a checkpoint file, we already have an edsamHistory struct
1139 if (oh->edsamHistory == nullptr)
1141 oh->edsamHistory = gmx::compat::make_unique<edsamhistory_t>(edsamhistory_t {});
1143 edsamhistory_t *EDstate = oh->edsamHistory.get();
1145 /* Read the edi input file: */
1146 ed->edpar = read_edi_file(ediFileName, natoms);
1148 /* Make sure the checkpoint was produced in a run using this .edi file */
1149 if (EDstate->bFromCpt)
1151 crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1155 EDstate->nED = ed->edpar.size();
1157 init_edsamstate(*ed, EDstate);
1159 /* The master opens the ED output file */
1162 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1166 ed->edo = xvgropen(edoFileName,
1167 "Essential dynamics / flooding output",
1169 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1171 /* Make a descriptive legend */
1172 write_edo_legend(ed, EDstate->nED, oenv);
1179 /* Broadcasts the structure data */
1180 static void bc_ed_positions(const t_commrec *cr, struct gmx_edx *s, EssentialDynamicsStructure stype)
1182 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1183 snew_bc(cr, s->x, s->nr ); /* Positions */
1184 nblock_bc(cr, s->nr, s->anrs );
1185 nblock_bc(cr, s->nr, s->x );
1187 /* For the average & reference structures we need an array for the collective indices,
1188 * and we need to broadcast the masses as well */
1189 if (stype == EssentialDynamicsStructure::Average || stype == EssentialDynamicsStructure::Reference)
1191 /* We need these additional variables in the parallel case: */
1192 snew(s->c_ind, s->nr ); /* Collective indices */
1193 /* Local atom indices get assigned in dd_make_local_group_indices.
1194 * There, also memory is allocated */
1195 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1196 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1197 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1200 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1201 if (stype == EssentialDynamicsStructure::Reference)
1203 snew_bc(cr, s->m, s->nr);
1204 nblock_bc(cr, s->nr, s->m);
1207 /* For the average structure we might need the masses for mass-weighting */
1208 if (stype == EssentialDynamicsStructure::Average)
1210 snew_bc(cr, s->sqrtm, s->nr);
1211 nblock_bc(cr, s->nr, s->sqrtm);
1212 snew_bc(cr, s->m, s->nr);
1213 nblock_bc(cr, s->nr, s->m);
1218 /* Broadcasts the eigenvector data */
1219 static void bc_ed_vecs(const t_commrec *cr, t_eigvec *ev, int length)
1223 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1224 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1225 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1226 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1227 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1229 nblock_bc(cr, ev->neig, ev->ieig );
1230 nblock_bc(cr, ev->neig, ev->stpsz );
1231 nblock_bc(cr, ev->neig, ev->xproj );
1232 nblock_bc(cr, ev->neig, ev->fproj );
1233 nblock_bc(cr, ev->neig, ev->refproj);
1235 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1236 for (i = 0; i < ev->neig; i++)
1238 snew_bc(cr, ev->vec[i], length);
1239 nblock_bc(cr, length, ev->vec[i]);
1245 /* Broadcasts the ED / flooding data to other nodes
1246 * and allocates memory where needed */
1247 static void broadcast_ed_data(const t_commrec *cr, gmx_edsam * ed)
1249 /* Master lets the other nodes know if its ED only or also flooding */
1250 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1252 int numedis = ed->edpar.size();
1253 /* First let everybody know how many ED data sets to expect */
1254 gmx_bcast(sizeof(numedis), &numedis, cr);
1255 nblock_abc(cr, numedis, &(ed->edpar));
1256 /* Now transfer the ED data set(s) */
1257 for (auto &edi : ed->edpar)
1259 /* Broadcast a single ED data set */
1262 /* Broadcast positions */
1263 bc_ed_positions(cr, &(edi.sref), EssentialDynamicsStructure::Reference); /* reference positions (don't broadcast masses) */
1264 bc_ed_positions(cr, &(edi.sav ), EssentialDynamicsStructure::Average ); /* average positions (do broadcast masses as well) */
1265 bc_ed_positions(cr, &(edi.star), EssentialDynamicsStructure::Target); /* target positions */
1266 bc_ed_positions(cr, &(edi.sori), EssentialDynamicsStructure::Origin); /* origin positions */
1268 /* Broadcast eigenvectors */
1269 bc_ed_vecs(cr, &edi.vecs.mon, edi.sav.nr);
1270 bc_ed_vecs(cr, &edi.vecs.linfix, edi.sav.nr);
1271 bc_ed_vecs(cr, &edi.vecs.linacc, edi.sav.nr);
1272 bc_ed_vecs(cr, &edi.vecs.radfix, edi.sav.nr);
1273 bc_ed_vecs(cr, &edi.vecs.radacc, edi.sav.nr);
1274 bc_ed_vecs(cr, &edi.vecs.radcon, edi.sav.nr);
1275 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1276 bc_ed_vecs(cr, &edi.flood.vecs, edi.sav.nr);
1278 /* For harmonic restraints the reference projections can change with time */
1279 if (edi.flood.bHarmonic)
1281 snew_bc(cr, edi.flood.initialReferenceProjection, edi.flood.vecs.neig);
1282 snew_bc(cr, edi.flood.referenceProjectionSlope, edi.flood.vecs.neig);
1283 nblock_bc(cr, edi.flood.vecs.neig, edi.flood.initialReferenceProjection);
1284 nblock_bc(cr, edi.flood.vecs.neig, edi.flood.referenceProjectionSlope);
1290 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1291 static void init_edi(const gmx_mtop_t *mtop, t_edpar *edi)
1294 real totalmass = 0.0;
1297 /* NOTE Init_edi is executed on the master process only
1298 * The initialized data sets are then transmitted to the
1299 * other nodes in broadcast_ed_data */
1301 /* evaluate masses (reference structure) */
1302 snew(edi->sref.m, edi->sref.nr);
1304 for (i = 0; i < edi->sref.nr; i++)
1308 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1312 edi->sref.m[i] = 1.0;
1315 /* Check that every m > 0. Bad things will happen otherwise. */
1316 if (edi->sref.m[i] <= 0.0)
1318 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1319 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1320 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1321 "atoms from the reference structure by creating a proper index group.\n",
1322 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1325 totalmass += edi->sref.m[i];
1327 edi->sref.mtot = totalmass;
1329 /* Masses m and sqrt(m) for the average structure. Note that m
1330 * is needed if forces have to be evaluated in do_edsam */
1331 snew(edi->sav.sqrtm, edi->sav.nr );
1332 snew(edi->sav.m, edi->sav.nr );
1333 for (i = 0; i < edi->sav.nr; i++)
1335 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1338 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1342 edi->sav.sqrtm[i] = 1.0;
1345 /* Check that every m > 0. Bad things will happen otherwise. */
1346 if (edi->sav.sqrtm[i] <= 0.0)
1348 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1349 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1350 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1351 "atoms from the average structure by creating a proper index group.\n",
1352 i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1356 /* put reference structure in origin */
1357 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1361 translate_x(edi->sref.x, edi->sref.nr, com);
1363 /* Init ED buffer */
1368 static void check(const char *line, const char *label)
1370 if (!strstr(line, label))
1372 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1377 static int read_checked_edint(FILE *file, const char *label)
1379 char line[STRLEN+1];
1382 fgets2 (line, STRLEN, file);
1384 fgets2 (line, STRLEN, file);
1385 sscanf (line, max_ev_fmt_d, &idum);
1389 static bool read_checked_edbool(FILE *file, const char *label)
1391 return read_checked_edint(file, label) != 0;
1394 static int read_edint(FILE *file, bool *bEOF)
1396 char line[STRLEN+1];
1401 eof = fgets2 (line, STRLEN, file);
1407 eof = fgets2 (line, STRLEN, file);
1413 sscanf (line, max_ev_fmt_d, &idum);
1419 static real read_checked_edreal(FILE *file, const char *label)
1421 char line[STRLEN+1];
1425 fgets2 (line, STRLEN, file);
1427 fgets2 (line, STRLEN, file);
1428 sscanf (line, max_ev_fmt_lf, &rdum);
1429 return static_cast<real>(rdum); /* always read as double and convert to single */
1433 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1436 char line[STRLEN+1];
1440 for (i = 0; i < number; i++)
1442 fgets2 (line, STRLEN, file);
1443 sscanf (line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1444 anrs[i]--; /* we are reading FORTRAN indices */
1445 for (j = 0; j < 3; j++)
1447 x[i][j] = d[j]; /* always read as double and convert to single */
1454 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1455 * \param[in] in the file to read from
1456 * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1457 * \param[out] vec the eigen-vectors
1458 * \param[in] nEig the number of eigenvectors
1460 void scan_edvec(FILE *in, int numAtoms, rvec ***vec, int nEig)
1463 for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1465 snew((*vec)[iEigenvector], numAtoms);
1466 for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1468 char line[STRLEN+1];
1470 fgets2(line, STRLEN, in);
1471 sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1472 (*vec)[iEigenvector][iAtom][XX] = x;
1473 (*vec)[iEigenvector][iAtom][YY] = y;
1474 (*vec)[iEigenvector][iAtom][ZZ] = z;
1479 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1480 * \param[in] in input file
1481 * \param[in] nrAtoms number of atoms/coordinates
1482 * \param[out] tvec the eigenvector
1484 void read_edvec(FILE *in, int nrAtoms, t_eigvec *tvec)
1486 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1487 if (tvec->neig <= 0)
1492 snew(tvec->ieig, tvec->neig);
1493 snew(tvec->stpsz, tvec->neig);
1494 for (int i = 0; i < tvec->neig; i++)
1496 char line[STRLEN+1];
1497 fgets2(line, STRLEN, in);
1498 int numEigenVectors;
1500 const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1503 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1505 tvec->ieig[i] = numEigenVectors;
1506 tvec->stpsz[i] = stepSize;
1507 } /* end of loop over eigenvectors */
1509 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1511 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1513 * Eigenvector and its intial reference and slope are stored on the
1514 * same line in the .edi format. To avoid re-winding the .edi file,
1515 * the reference eigen-vector and reference data are read in one go.
1517 * \param[in] in input file
1518 * \param[in] nrAtoms number of atoms/coordinates
1519 * \param[out] tvec the eigenvector
1520 * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1521 * \param[out] referenceProjectionSlope The slope of the reference projections.
1523 bool readEdVecWithReferenceProjection(FILE *in, int nrAtoms, t_eigvec *tvec, real ** initialReferenceProjection, real ** referenceProjectionSlope)
1525 bool providesReference = false;
1526 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1527 if (tvec->neig <= 0)
1532 snew(tvec->ieig, tvec->neig);
1533 snew(tvec->stpsz, tvec->neig);
1534 snew(*initialReferenceProjection, tvec->neig);
1535 snew(*referenceProjectionSlope, tvec->neig);
1536 for (int i = 0; i < tvec->neig; i++)
1538 char line[STRLEN+1];
1539 fgets2 (line, STRLEN, in);
1540 int numEigenVectors;
1541 double stepSize = 0.;
1542 double referenceProjection = 0.;
1543 double referenceSlope = 0.;
1544 const int nscan = sscanf(line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize, &referenceProjection, &referenceSlope);
1545 /* Zero out values which were not scanned */
1549 /* Every 4 values read, including reference position */
1550 providesReference = true;
1553 /* A reference position is provided */
1554 providesReference = true;
1555 /* No value for slope, set to 0 */
1556 referenceSlope = 0.0;
1559 /* No values for reference projection and slope, set to 0 */
1560 referenceProjection = 0.0;
1561 referenceSlope = 0.0;
1564 gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1566 (*initialReferenceProjection)[i] = referenceProjection;
1567 (*referenceProjectionSlope)[i] = referenceSlope;
1569 tvec->ieig[i] = numEigenVectors;
1570 tvec->stpsz[i] = stepSize;
1571 } /* end of loop over eigenvectors */
1573 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1574 return providesReference;
1577 /*!\brief Allocate working memory for the eigenvectors.
1578 * \param[in,out] tvec eigenvector for which memory will be allocated
1580 void setup_edvec(t_eigvec *tvec)
1582 snew(tvec->xproj, tvec->neig);
1583 snew(tvec->fproj, tvec->neig);
1584 snew(tvec->refproj, tvec->neig);
1589 /* Check if the same atom indices are used for reference and average positions */
1590 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1595 /* If the number of atoms differs between the two structures,
1596 * they cannot be identical */
1597 if (sref.nr != sav.nr)
1602 /* Now that we know that both stuctures have the same number of atoms,
1603 * check if also the indices are identical */
1604 for (i = 0; i < sav.nr; i++)
1606 if (sref.anrs[i] != sav.anrs[i])
1611 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1619 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1620 constexpr int c_oldestUnsupportedMagicNumber = 666;
1621 //! Oldest (lowest) magic number indicating supported essential dynamics input
1622 constexpr int c_oldestSupportedMagicNumber = 669;
1623 //! Latest (highest) magic number indicating supported essential dynamics input
1624 constexpr int c_latestSupportedMagicNumber = 670;
1626 /*!\brief Set up essential dynamics work parameters.
1627 * \param[in] edi Essential dynamics working structure.
1629 void setup_edi(t_edpar *edi)
1631 snew(edi->sref.x_old, edi->sref.nr);
1632 edi->sref.sqrtm = nullptr;
1633 snew(edi->sav.x_old, edi->sav.nr);
1634 if (edi->star.nr > 0)
1636 edi->star.sqrtm = nullptr;
1638 if (edi->sori.nr > 0)
1640 edi->sori.sqrtm = nullptr;
1642 setup_edvec(&edi->vecs.linacc);
1643 setup_edvec(&edi->vecs.mon);
1644 setup_edvec(&edi->vecs.linfix);
1645 setup_edvec(&edi->vecs.linacc);
1646 setup_edvec(&edi->vecs.radfix);
1647 setup_edvec(&edi->vecs.radacc);
1648 setup_edvec(&edi->vecs.radcon);
1649 setup_edvec(&edi->flood.vecs);
1652 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1653 * \param[in] magicNumber indicates the essential dynamics input file version
1654 * \returns true if CONST_FORCE_FLOODING is supported
1656 bool ediFileHasConstForceFlooding(int magicNumber)
1658 return magicNumber > c_oldestSupportedMagicNumber;
1661 /*!\brief Reports reading success of the essential dynamics input file magic number.
1662 * \param[in] in pointer to input file
1663 * \param[out] magicNumber determines the edi file version
1664 * \returns true if setting file version from input file was successful.
1666 bool ediFileHasValidDataSet(FILE *in, int * magicNumber)
1668 /* the edi file is not free format, so expect problems if the input is corrupt. */
1669 bool reachedEndOfFile = true;
1670 /* check the magic number */
1671 *magicNumber = read_edint(in, &reachedEndOfFile);
1672 /* Check whether we have reached the end of the input file */
1673 return !reachedEndOfFile;
1676 /*!\brief Trigger fatal error if magic number is unsupported.
1677 * \param[in] magicNumber A magic number read from the edi file
1678 * \param[in] fn name of input file for error reporting
1680 void exitWhenMagicNumberIsInvalid(int magicNumber, const char * fn)
1683 if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1685 if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1687 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1691 gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1696 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1698 * \param[in] in input file
1699 * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1700 * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1701 * \param[in] fn the file name of the input file for error reporting
1702 * \returns edi essential dynamics parameters
1704 t_edpar read_edi(FILE* in, int nr_mdatoms, bool hasConstForceFlooding, const char *fn)
1708 /* check the number of atoms */
1709 edi.nini = read_edint(in, &bEOF);
1710 if (edi.nini != nr_mdatoms)
1712 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi.nini, nr_mdatoms);
1715 /* Done checking. For the rest we blindly trust the input */
1716 edi.fitmas = read_checked_edbool(in, "FITMAS");
1717 edi.pcamas = read_checked_edbool(in, "ANALYSIS_MAS");
1718 edi.outfrq = read_checked_edint(in, "OUTFRQ");
1719 edi.maxedsteps = read_checked_edint(in, "MAXLEN");
1720 edi.slope = read_checked_edreal(in, "SLOPECRIT");
1722 edi.presteps = read_checked_edint(in, "PRESTEPS");
1723 edi.flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1724 edi.flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1725 edi.flood.tau = read_checked_edreal(in, "TAU");
1726 edi.flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1727 edi.flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1728 edi.flood.kT = read_checked_edreal(in, "KT");
1729 edi.flood.bHarmonic = read_checked_edbool(in, "HARMONIC");
1730 if (hasConstForceFlooding)
1732 edi.flood.bConstForce = read_checked_edbool(in, "CONST_FORCE_FLOODING");
1736 edi.flood.bConstForce = FALSE;
1738 edi.sref.nr = read_checked_edint(in, "NREF");
1740 /* allocate space for reference positions and read them */
1741 snew(edi.sref.anrs, edi.sref.nr);
1742 snew(edi.sref.x, edi.sref.nr);
1743 read_edx(in, edi.sref.nr, edi.sref.anrs, edi.sref.x);
1745 /* average positions. they define which atoms will be used for ED sampling */
1746 edi.sav.nr = read_checked_edint(in, "NAV");
1747 snew(edi.sav.anrs, edi.sav.nr);
1748 snew(edi.sav.x, edi.sav.nr);
1749 read_edx(in, edi.sav.nr, edi.sav.anrs, edi.sav.x);
1751 /* Check if the same atom indices are used for reference and average positions */
1752 edi.bRefEqAv = check_if_same(edi.sref, edi.sav);
1756 read_edvec(in, edi.sav.nr, &edi.vecs.mon);
1757 read_edvec(in, edi.sav.nr, &edi.vecs.linfix);
1758 read_edvec(in, edi.sav.nr, &edi.vecs.linacc);
1759 read_edvec(in, edi.sav.nr, &edi.vecs.radfix);
1760 read_edvec(in, edi.sav.nr, &edi.vecs.radacc);
1761 read_edvec(in, edi.sav.nr, &edi.vecs.radcon);
1763 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1764 bool bHaveReference = false;
1765 if (edi.flood.bHarmonic)
1767 bHaveReference = readEdVecWithReferenceProjection(in, edi.sav.nr, &edi.flood.vecs, &edi.flood.initialReferenceProjection, &edi.flood.referenceProjectionSlope);
1771 read_edvec(in, edi.sav.nr, &edi.flood.vecs);
1774 /* target positions */
1775 edi.star.nr = read_edint(in, &bEOF);
1776 if (edi.star.nr > 0)
1778 snew(edi.star.anrs, edi.star.nr);
1779 snew(edi.star.x, edi.star.nr);
1780 read_edx(in, edi.star.nr, edi.star.anrs, edi.star.x);
1783 /* positions defining origin of expansion circle */
1784 edi.sori.nr = read_edint(in, &bEOF);
1785 if (edi.sori.nr > 0)
1789 /* Both an -ori structure and a at least one manual reference point have been
1790 * specified. That's ambiguous and probably not intentional. */
1791 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1792 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1794 snew(edi.sori.anrs, edi.sori.nr);
1795 snew(edi.sori.x, edi.sori.nr);
1796 read_edx(in, edi.sori.nr, edi.sori.anrs, edi.sori.x);
1801 std::vector<t_edpar> read_edi_file(const char *fn, int nr_mdatoms)
1803 std::vector<t_edpar> essentialDynamicsParameters;
1805 /* This routine is executed on the master only */
1807 /* Open the .edi parameter input file */
1808 in = gmx_fio_fopen(fn, "r");
1809 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1811 /* Now read a sequence of ED input parameter sets from the edi file */
1812 int ediFileMagicNumber;
1813 while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1815 exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1816 bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1817 auto edi = read_edi(in, nr_mdatoms, hasConstForceFlooding, fn);
1819 essentialDynamicsParameters.emplace_back(edi);
1821 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1822 /* We need to allocate space for the data: */
1824 if (essentialDynamicsParameters.empty())
1826 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1829 fprintf(stderr, "ED: Found %zu ED group%s.\n", essentialDynamicsParameters.size(), essentialDynamicsParameters.size() > 1 ? "s" : "");
1831 /* Close the .edi file again */
1834 return essentialDynamicsParameters;
1836 } // anonymous namespace
1839 struct t_fit_to_ref {
1840 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1843 /* Fit the current positions to the reference positions
1844 * Do not actually do the fit, just return rotation and translation.
1845 * Note that the COM of the reference structure was already put into
1846 * the origin by init_edi. */
1847 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1848 rvec transvec, /* The translation vector */
1849 matrix rotmat, /* The rotation matrix */
1850 t_edpar *edi) /* Just needed for do_edfit */
1852 rvec com; /* center of mass */
1854 struct t_fit_to_ref *loc;
1857 /* Allocate memory the first time this routine is called for each edi group */
1858 if (nullptr == edi->buf->fit_to_ref)
1860 snew(edi->buf->fit_to_ref, 1);
1861 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1863 loc = edi->buf->fit_to_ref;
1865 /* We do not touch the original positions but work on a copy. */
1866 for (i = 0; i < edi->sref.nr; i++)
1868 copy_rvec(xcoll[i], loc->xcopy[i]);
1871 /* Calculate the center of mass */
1872 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1874 transvec[XX] = -com[XX];
1875 transvec[YY] = -com[YY];
1876 transvec[ZZ] = -com[ZZ];
1878 /* Subtract the center of mass from the copy */
1879 translate_x(loc->xcopy, edi->sref.nr, transvec);
1881 /* Determine the rotation matrix */
1882 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1886 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1887 int nat, /* How many positions are there? */
1888 rvec transvec, /* The translation vector */
1889 matrix rotmat) /* The rotation matrix */
1892 translate_x(x, nat, transvec);
1895 rotate_x(x, nat, rotmat);
1899 /* Gets the rms deviation of the positions to the structure s */
1900 /* fit_to_structure has to be called before calling this routine! */
1901 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1902 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1908 for (i = 0; i < s->nr; i++)
1910 rmsd += distance2(s->x[i], x[i]);
1913 rmsd /= static_cast<real>(s->nr);
1920 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1922 if (ed->eEDtype != EssentialDynamicsType::None)
1924 /* Loop over ED groups */
1926 for (auto &edi : ed->edpar)
1928 /* Local atoms of the reference structure (for fitting), need only be assembled
1929 * if their indices differ from the average ones */
1932 dd_make_local_group_indices(dd->ga2la, edi.sref.nr, edi.sref.anrs,
1933 &edi.sref.nr_loc, &edi.sref.anrs_loc, &edi.sref.nalloc_loc, edi.sref.c_ind);
1936 /* Local atoms of the average structure (on these ED will be performed) */
1937 dd_make_local_group_indices(dd->ga2la, edi.sav.nr, edi.sav.anrs,
1938 &edi.sav.nr_loc, &edi.sav.anrs_loc, &edi.sav.nalloc_loc, edi.sav.c_ind);
1940 /* Indicate that the ED shift vectors for this structure need to be updated
1941 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1942 edi.buf->do_edsam->bUpdateShifts = TRUE;
1948 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1959 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1960 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1961 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1965 xu[XX] = x[XX]-tx*box[XX][XX];
1966 xu[YY] = x[YY]-ty*box[YY][YY];
1967 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1973 /*!\brief Apply fixed linear constraints to essential dynamics variable.
1974 * \param[in,out] xcoll The collected coordinates.
1975 * \param[in] edi the essential dynamics parameters
1976 * \param[in] step the current simulation step
1978 void do_linfix(rvec *xcoll, const t_edpar &edi, int64_t step)
1980 /* loop over linfix vectors */
1981 for (int i = 0; i < edi.vecs.linfix.neig; i++)
1983 /* calculate the projection */
1984 real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
1986 /* calculate the correction */
1987 real preFactor = edi.vecs.linfix.refproj[i] + step*edi.vecs.linfix.stpsz[i] - proj;
1989 /* apply the correction */
1990 preFactor /= edi.sav.sqrtm[i];
1991 for (int j = 0; j < edi.sav.nr; j++)
1993 rvec differenceVector;
1994 svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
1995 rvec_inc(xcoll[j], differenceVector);
2000 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
2001 * \param[in,out] xcoll The collected coordinates.
2002 * \param[in] edi the essential dynamics parameters
2004 void do_linacc(rvec *xcoll, t_edpar *edi)
2006 /* loop over linacc vectors */
2007 for (int i = 0; i < edi->vecs.linacc.neig; i++)
2009 /* calculate the projection */
2010 real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
2012 /* calculate the correction */
2013 real preFactor = 0.0;
2014 if (edi->vecs.linacc.stpsz[i] > 0.0)
2016 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
2018 preFactor = edi->vecs.linacc.refproj[i] - proj;
2021 if (edi->vecs.linacc.stpsz[i] < 0.0)
2023 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
2025 preFactor = edi->vecs.linacc.refproj[i] - proj;
2029 /* apply the correction */
2030 preFactor /= edi->sav.sqrtm[i];
2031 for (int j = 0; j < edi->sav.nr; j++)
2033 rvec differenceVector;
2034 svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2035 rvec_inc(xcoll[j], differenceVector);
2038 /* new positions will act as reference */
2039 edi->vecs.linacc.refproj[i] = proj + preFactor;
2045 static void do_radfix(rvec *xcoll, t_edpar *edi)
2048 real *proj, rad = 0.0, ratio;
2052 if (edi->vecs.radfix.neig == 0)
2057 snew(proj, edi->vecs.radfix.neig);
2059 /* loop over radfix vectors */
2060 for (i = 0; i < edi->vecs.radfix.neig; i++)
2062 /* calculate the projections, radius */
2063 proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2064 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2068 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2069 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2071 /* loop over radfix vectors */
2072 for (i = 0; i < edi->vecs.radfix.neig; i++)
2074 proj[i] -= edi->vecs.radfix.refproj[i];
2076 /* apply the correction */
2077 proj[i] /= edi->sav.sqrtm[i];
2079 for (j = 0; j < edi->sav.nr; j++)
2081 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2082 rvec_inc(xcoll[j], vec_dum);
2090 static void do_radacc(rvec *xcoll, t_edpar *edi)
2093 real *proj, rad = 0.0, ratio = 0.0;
2097 if (edi->vecs.radacc.neig == 0)
2102 snew(proj, edi->vecs.radacc.neig);
2104 /* loop over radacc vectors */
2105 for (i = 0; i < edi->vecs.radacc.neig; i++)
2107 /* calculate the projections, radius */
2108 proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2109 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2113 /* only correct when radius decreased */
2114 if (rad < edi->vecs.radacc.radius)
2116 ratio = edi->vecs.radacc.radius/rad - 1.0;
2120 edi->vecs.radacc.radius = rad;
2123 /* loop over radacc vectors */
2124 for (i = 0; i < edi->vecs.radacc.neig; i++)
2126 proj[i] -= edi->vecs.radacc.refproj[i];
2128 /* apply the correction */
2129 proj[i] /= edi->sav.sqrtm[i];
2131 for (j = 0; j < edi->sav.nr; j++)
2133 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2134 rvec_inc(xcoll[j], vec_dum);
2141 struct t_do_radcon {
2145 static void do_radcon(rvec *xcoll, t_edpar *edi)
2148 real rad = 0.0, ratio = 0.0;
2149 struct t_do_radcon *loc;
2154 if (edi->buf->do_radcon != nullptr)
2161 snew(edi->buf->do_radcon, 1);
2163 loc = edi->buf->do_radcon;
2165 if (edi->vecs.radcon.neig == 0)
2172 snew(loc->proj, edi->vecs.radcon.neig);
2175 /* loop over radcon vectors */
2176 for (i = 0; i < edi->vecs.radcon.neig; i++)
2178 /* calculate the projections, radius */
2179 loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2180 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2183 /* only correct when radius increased */
2184 if (rad > edi->vecs.radcon.radius)
2186 ratio = edi->vecs.radcon.radius/rad - 1.0;
2188 /* loop over radcon vectors */
2189 for (i = 0; i < edi->vecs.radcon.neig; i++)
2191 /* apply the correction */
2192 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2193 loc->proj[i] /= edi->sav.sqrtm[i];
2194 loc->proj[i] *= ratio;
2196 for (j = 0; j < edi->sav.nr; j++)
2198 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2199 rvec_inc(xcoll[j], vec_dum);
2206 edi->vecs.radcon.radius = rad;
2212 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, int64_t step)
2217 /* subtract the average positions */
2218 for (i = 0; i < edi->sav.nr; i++)
2220 rvec_dec(xcoll[i], edi->sav.x[i]);
2223 /* apply the constraints */
2226 do_linfix(xcoll, *edi, step);
2228 do_linacc(xcoll, edi);
2231 do_radfix(xcoll, edi);
2233 do_radacc(xcoll, edi);
2234 do_radcon(xcoll, edi);
2236 /* add back the average positions */
2237 for (i = 0; i < edi->sav.nr; i++)
2239 rvec_inc(xcoll[i], edi->sav.x[i]);
2246 /*!\brief Write out the projections onto the eigenvectors.
2247 * The order of output corresponds to ed_output_legend().
2248 * \param[in] edi The essential dyanmics parameters
2249 * \param[in] fp The output file
2250 * \param[in] rmsd the rmsd to the reference structure
2252 void write_edo(const t_edpar &edi, FILE *fp, real rmsd)
2254 /* Output how well we fit to the reference structure */
2255 fprintf(fp, EDcol_ffmt, rmsd);
2257 for (int i = 0; i < edi.vecs.mon.neig; i++)
2259 fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2262 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2264 fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2267 for (int i = 0; i < edi.vecs.linacc.neig; i++)
2269 fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2272 for (int i = 0; i < edi.vecs.radfix.neig; i++)
2274 fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2276 if (edi.vecs.radfix.neig)
2278 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2281 for (int i = 0; i < edi.vecs.radacc.neig; i++)
2283 fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2285 if (edi.vecs.radacc.neig)
2287 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2290 for (int i = 0; i < edi.vecs.radcon.neig; i++)
2292 fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2294 if (edi.vecs.radcon.neig)
2296 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2302 /* Returns if any constraints are switched on */
2303 static bool ed_constraints(EssentialDynamicsType edtype, const t_edpar &edi)
2305 if (edtype == EssentialDynamicsType::EDSampling || edtype == EssentialDynamicsType::Flooding)
2307 return ((edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0) ||
2308 (edi.vecs.radfix.neig != 0) || (edi.vecs.radacc.neig != 0) ||
2309 (edi.vecs.radcon.neig != 0));
2315 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for flooding/
2316 * umbrella sampling simulations. */
2317 static void copyEvecReference(t_eigvec* floodvecs, real * initialReferenceProjection)
2322 if (nullptr == initialReferenceProjection)
2324 snew(initialReferenceProjection, floodvecs->neig);
2327 for (i = 0; i < floodvecs->neig; i++)
2329 initialReferenceProjection[i] = floodvecs->refproj[i];
2334 /* Call on MASTER only. Check whether the essential dynamics / flooding
2335 * groups of the checkpoint file are consistent with the provided .edi file. */
2336 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam &ed, edsamhistory_t *EDstate)
2338 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2340 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2341 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2342 "it must also continue with/without ED constraints when checkpointing.\n"
2343 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2344 "from without a checkpoint.\n");
2347 for (size_t edinum = 0; edinum < ed.edpar.size(); ++edinum)
2349 /* Check number of atoms in the reference and average structures */
2350 if (EDstate->nref[edinum] != ed.edpar[edinum].sref.nr)
2352 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2353 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2354 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], ed.edpar[edinum].sref.nr);
2356 if (EDstate->nav[edinum] != ed.edpar[edinum].sav.nr)
2358 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2359 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2360 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], ed.edpar[edinum].sav.nr);
2364 if (static_cast<int>(ed.edpar.size()) != EDstate->nED)
2366 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2367 "There are %d ED groups in the .cpt file, but %zu in the .edi file!\n"
2368 "Are you sure this is the correct .edi file?\n", EDstate->nED, ed.edpar.size());
2373 /* The edsamstate struct stores the information we need to make the ED group
2374 * whole again after restarts from a checkpoint file. Here we do the following:
2375 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2376 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2377 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2378 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2379 * all ED structures at the last time step. */
2380 static void init_edsamstate(const gmx_edsam &ed, edsamhistory_t *EDstate)
2382 snew(EDstate->old_sref_p, EDstate->nED);
2383 snew(EDstate->old_sav_p, EDstate->nED);
2385 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2386 if (!EDstate->bFromCpt)
2388 snew(EDstate->nref, EDstate->nED);
2389 snew(EDstate->nav, EDstate->nED);
2392 /* Loop over all ED/flooding data sets (usually only one, though) */
2393 for (int nr_edi = 0; nr_edi < EDstate->nED; nr_edi++)
2395 const auto &edi = ed.edpar[nr_edi];
2396 /* We always need the last reference and average positions such that
2397 * in the next time step we can make the ED group whole again
2398 * if the atoms do not have the correct PBC representation */
2399 if (EDstate->bFromCpt)
2401 /* Copy the last whole positions of reference and average group from .cpt */
2402 for (int i = 0; i < edi.sref.nr; i++)
2404 copy_rvec(EDstate->old_sref[nr_edi][i], edi.sref.x_old[i]);
2406 for (int i = 0; i < edi.sav.nr; i++)
2408 copy_rvec(EDstate->old_sav [nr_edi][i], edi.sav.x_old [i]);
2413 EDstate->nref[nr_edi] = edi.sref.nr;
2414 EDstate->nav [nr_edi] = edi.sav.nr;
2417 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2418 EDstate->old_sref_p[nr_edi] = edi.sref.x_old;
2419 EDstate->old_sav_p [nr_edi] = edi.sav.x_old;
2424 /* Adds 'buf' to 'str' */
2425 static void add_to_string(char **str, const char *buf)
2430 len = strlen(*str) + strlen(buf) + 1;
2436 static void add_to_string_aligned(char **str, const char *buf)
2438 char buf_aligned[STRLEN];
2440 sprintf(buf_aligned, EDcol_sfmt, buf);
2441 add_to_string(str, buf_aligned);
2445 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, const char *value, const char *unit, char EDgroupchar)
2447 auto tmp = gmx::formatString("%c %s", EDgroupchar, value);
2448 add_to_string_aligned(LegendStr, tmp.c_str());
2449 tmp += gmx::formatString(" (%s)", unit);
2450 (*setname)[*nsets] = gmx_strdup(tmp.c_str());
2455 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2461 for (i = 0; i < evec->neig; i++)
2463 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2464 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2469 /* Makes a legend for the xvg output file. Call on MASTER only! */
2470 static void write_edo_legend(gmx_edsam * ed, int nED, const gmx_output_env_t *oenv)
2473 int nr_edi, nsets, n_flood, n_edsam;
2474 const char **setname;
2476 char *LegendStr = nullptr;
2479 auto edi = ed->edpar.begin();
2481 fprintf(ed->edo, "# Output will be written every %d step%s\n", edi->outfrq, edi->outfrq != 1 ? "s" : "");
2483 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2485 fprintf(ed->edo, "#\n");
2486 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2487 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2488 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2489 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2490 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2491 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2492 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2493 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2494 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2496 if (edi->flood.vecs.neig)
2498 /* If in any of the groups we find a flooding vector, flooding is turned on */
2499 ed->eEDtype = EssentialDynamicsType::Flooding;
2501 /* Print what flavor of flooding we will do */
2502 if (0 == edi->flood.tau) /* constant flooding strength */
2504 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2505 if (edi->flood.bHarmonic)
2507 fprintf(ed->edo, ", harmonic");
2510 else /* adaptive flooding */
2512 fprintf(ed->edo, ", adaptive");
2515 fprintf(ed->edo, "\n");
2520 /* Print a nice legend */
2522 LegendStr[0] = '\0';
2523 sprintf(buf, "# %6s", "time");
2524 add_to_string(&LegendStr, buf);
2526 /* Calculate the maximum number of columns we could end up with */
2527 edi = ed->edpar.begin();
2529 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2531 nsets += 5 +edi->vecs.mon.neig
2532 +edi->vecs.linfix.neig
2533 +edi->vecs.linacc.neig
2534 +edi->vecs.radfix.neig
2535 +edi->vecs.radacc.neig
2536 +edi->vecs.radcon.neig
2537 + 6*edi->flood.vecs.neig;
2540 snew(setname, nsets);
2542 /* In the mdrun time step in a first function call (do_flood()) the flooding
2543 * forces are calculated and in a second function call (do_edsam()) the
2544 * ED constraints. To get a corresponding legend, we need to loop twice
2545 * over the edi groups and output first the flooding, then the ED part */
2547 /* The flooding-related legend entries, if flooding is done */
2549 if (EssentialDynamicsType::Flooding == ed->eEDtype)
2551 edi = ed->edpar.begin();
2552 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2554 /* Always write out the projection on the flooding EVs. Of course, this can also
2555 * be achieved with the monitoring option in do_edsam() (if switched on by the
2556 * user), but in that case the positions need to be communicated in do_edsam(),
2557 * which is not necessary when doing flooding only. */
2558 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2560 for (i = 0; i < edi->flood.vecs.neig; i++)
2562 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2563 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2565 /* Output the current reference projection if it changes with time;
2566 * this can happen when flooding is used as harmonic restraint */
2567 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2569 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2570 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2573 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2574 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2576 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2577 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2580 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2581 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2583 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2585 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2586 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2589 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2590 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2594 } /* End of flooding-related legend entries */
2598 /* Now the ED-related entries, if essential dynamics is done */
2599 edi = ed->edpar.begin();
2600 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2602 if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2604 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2606 /* Essential dynamics, projections on eigenvectors */
2607 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2608 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2609 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2610 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2611 if (edi->vecs.radfix.neig)
2613 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2615 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2616 if (edi->vecs.radacc.neig)
2618 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2620 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2621 if (edi->vecs.radcon.neig)
2623 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2627 } /* end of 'pure' essential dynamics legend entries */
2628 n_edsam = nsets - n_flood;
2630 xvgr_legend(ed->edo, nsets, setname, oenv);
2633 fprintf(ed->edo, "#\n"
2634 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2635 n_flood, 1 == n_flood ? "" : "s",
2636 n_edsam, 1 == n_edsam ? "" : "s");
2637 fprintf(ed->edo, "%s", LegendStr);
2644 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2645 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2646 std::unique_ptr<gmx::EssentialDynamics> init_edsam(
2647 const gmx::MDLogger &mdlog,
2648 const char *ediFileName,
2649 const char *edoFileName,
2650 const gmx_mtop_t *mtop,
2651 const t_inputrec *ir,
2652 const t_commrec *cr,
2653 gmx::Constraints *constr,
2654 const t_state *globalState,
2655 ObservablesHistory *oh,
2656 const gmx_output_env_t *oenv,
2660 rvec *x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2661 rvec *xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2662 rvec fit_transvec; /* translation ... */
2663 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2664 rvec *ref_x_old = nullptr; /* helper pointer */
2669 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2671 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName) )
2673 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or flooding\n"
2674 "simulation. Please also set the .edi file on the command line with -ei.\n");
2678 GMX_LOG(mdlog.info).asParagraph().
2679 appendText("Activating essential dynamics via files passed to "
2680 "gmx mdrun -edi will change in future to be files passed to "
2681 "gmx grompp and the related .mdp options may change also.");
2683 /* Open input and output files, allocate space for ED data structure */
2684 auto edHandle = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr);
2685 auto ed = edHandle->getLegacyED();
2686 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2687 constr->saveEdsamPointer(ed);
2689 /* Needed for initializing radacc radius in do_edsam */
2692 /* The input file is read by the master and the edi structures are
2693 * initialized here. Input is stored in ed->edpar. Then the edi
2694 * structures are transferred to the other nodes */
2697 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2698 * flooding vector, Essential dynamics can be applied to more than one structure
2699 * as well, but will be done in the order given in the edi file, so
2700 * expect different results for different order of edi file concatenation! */
2701 for (auto &edi : ed->edpar)
2703 init_edi(mtop, &edi);
2704 init_flood(&edi, ed, ir->delta_t);
2708 /* The master does the work here. The other nodes get the positions
2709 * not before dd_partition_system which is called after init_edsam */
2712 edsamhistory_t *EDstate = oh->edsamHistory.get();
2714 if (!EDstate->bFromCpt)
2716 /* Remove PBC, make molecule(s) subject to ED whole. */
2717 snew(x_pbc, mtop->natoms);
2718 copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop->natoms);
2719 do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
2721 /* Reset pointer to first ED data set which contains the actual ED data */
2722 auto edi = ed->edpar.begin();
2723 GMX_ASSERT(!ed->edpar.empty(), "Essential dynamics parameters is undefined");
2725 /* Loop over all ED/flooding data sets (usually only one, though) */
2726 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2728 /* For multiple ED groups we use the output frequency that was specified
2729 * in the first set */
2732 edi->outfrq = ed->edpar.begin()->outfrq;
2735 /* Extract the initial reference and average positions. When starting
2736 * from .cpt, these have already been read into sref.x_old
2737 * in init_edsamstate() */
2738 if (!EDstate->bFromCpt)
2740 /* If this is the first run (i.e. no checkpoint present) we assume
2741 * that the starting positions give us the correct PBC representation */
2742 for (i = 0; i < edi->sref.nr; i++)
2744 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2747 for (i = 0; i < edi->sav.nr; i++)
2749 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2753 /* Now we have the PBC-correct start positions of the reference and
2754 average structure. We copy that over to dummy arrays on which we
2755 can apply fitting to print out the RMSD. We srenew the memory since
2756 the size of the buffers is likely different for every ED group */
2757 srenew(xfit, edi->sref.nr );
2758 srenew(xstart, edi->sav.nr );
2761 /* Reference indices are the same as average indices */
2762 ref_x_old = edi->sav.x_old;
2766 ref_x_old = edi->sref.x_old;
2768 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2769 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2771 /* Make the fit to the REFERENCE structure, get translation and rotation */
2772 fit_to_reference(xfit, fit_transvec, fit_rotmat, &(*edi));
2774 /* Output how well we fit to the reference at the start */
2775 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2776 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2777 rmsd_from_structure(xfit, &edi->sref));
2778 if (EDstate->nED > 1)
2780 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2782 fprintf(stderr, "\n");
2784 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2785 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2787 /* calculate initial projections */
2788 project(xstart, &(*edi));
2790 /* For the target and origin structure both a reference (fit) and an
2791 * average structure can be provided in make_edi. If both structures
2792 * are the same, make_edi only stores one of them in the .edi file.
2793 * If they differ, first the fit and then the average structure is stored
2794 * in star (or sor), thus the number of entries in star/sor is
2795 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2796 * the size of the average group. */
2798 /* process target structure, if required */
2799 if (edi->star.nr > 0)
2801 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2803 /* get translation & rotation for fit of target structure to reference structure */
2804 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, &(*edi));
2806 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2807 if (edi->star.nr == edi->sav.nr)
2811 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2813 /* The last sav.nr indices of the target structure correspond to
2814 * the average structure, which must be projected */
2815 avindex = edi->star.nr - edi->sav.nr;
2817 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2821 rad_project(*edi, xstart, &edi->vecs.radcon);
2824 /* process structure that will serve as origin of expansion circle */
2825 if ( (EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce) )
2827 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2830 if (edi->sori.nr > 0)
2832 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2834 /* fit this structure to reference structure */
2835 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, &(*edi));
2837 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2838 if (edi->sori.nr == edi->sav.nr)
2842 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2844 /* For the projection, we need the last sav.nr indices of sori */
2845 avindex = edi->sori.nr - edi->sav.nr;
2848 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2849 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2850 if ( (EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce) )
2852 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2853 /* Set center of flooding potential to the ORIGIN structure */
2854 rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2855 /* We already know that no (moving) reference position was provided,
2856 * therefore we can overwrite refproj[0]*/
2857 copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
2860 else /* No origin structure given */
2862 rad_project(*edi, xstart, &edi->vecs.radacc);
2863 rad_project(*edi, xstart, &edi->vecs.radfix);
2864 if ( (EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce) )
2866 if (edi->flood.bHarmonic)
2868 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2869 for (i = 0; i < edi->flood.vecs.neig; i++)
2871 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
2876 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2877 /* Set center of flooding potential to the center of the covariance matrix,
2878 * i.e. the average structure, i.e. zero in the projected system */
2879 for (i = 0; i < edi->flood.vecs.neig; i++)
2881 edi->flood.vecs.refproj[i] = 0.0;
2886 /* For convenience, output the center of the flooding potential for the eigenvectors */
2887 if ( (EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce) )
2889 for (i = 0; i < edi->flood.vecs.neig; i++)
2891 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2892 if (edi->flood.bHarmonic)
2894 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.referenceProjectionSlope[i]);
2896 fprintf(stdout, "\n");
2900 /* set starting projections for linsam */
2901 rad_project(*edi, xstart, &edi->vecs.linacc);
2902 rad_project(*edi, xstart, &edi->vecs.linfix);
2904 /* Prepare for the next edi data set: */
2907 /* Cleaning up on the master node: */
2908 if (!EDstate->bFromCpt)
2915 } /* end of MASTER only section */
2919 /* Broadcast the essential dynamics / flooding data to all nodes */
2920 broadcast_ed_data(cr, ed);
2924 /* In the single-CPU case, point the local atom numbers pointers to the global
2925 * one, so that we can use the same notation in serial and parallel case: */
2926 /* Loop over all ED data sets (usually only one, though) */
2927 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
2929 edi->sref.anrs_loc = edi->sref.anrs;
2930 edi->sav.anrs_loc = edi->sav.anrs;
2931 edi->star.anrs_loc = edi->star.anrs;
2932 edi->sori.anrs_loc = edi->sori.anrs;
2933 /* For the same reason as above, make a dummy c_ind array: */
2934 snew(edi->sav.c_ind, edi->sav.nr);
2935 /* Initialize the array */
2936 for (i = 0; i < edi->sav.nr; i++)
2938 edi->sav.c_ind[i] = i;
2940 /* In the general case we will need a different-sized array for the reference indices: */
2943 snew(edi->sref.c_ind, edi->sref.nr);
2944 for (i = 0; i < edi->sref.nr; i++)
2946 edi->sref.c_ind[i] = i;
2949 /* Point to the very same array in case of other structures: */
2950 edi->star.c_ind = edi->sav.c_ind;
2951 edi->sori.c_ind = edi->sav.c_ind;
2952 /* In the serial case, the local number of atoms is the global one: */
2953 edi->sref.nr_loc = edi->sref.nr;
2954 edi->sav.nr_loc = edi->sav.nr;
2955 edi->star.nr_loc = edi->star.nr;
2956 edi->sori.nr_loc = edi->sori.nr;
2960 /* Allocate space for ED buffer variables */
2961 /* Again, loop over ED data sets */
2962 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
2964 /* Allocate space for ED buffer variables */
2965 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
2966 snew(edi->buf->do_edsam, 1);
2968 /* Space for collective ED buffer variables */
2970 /* Collective positions of atoms with the average indices */
2971 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2972 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
2973 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2974 /* Collective positions of atoms with the reference indices */
2977 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2978 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
2979 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2982 /* Get memory for flooding forces */
2983 snew(edi->flood.forces_cartesian, edi->sav.nr);
2986 /* Flush the edo file so that the user can check some things
2987 * when the simulation has started */
2997 void do_edsam(const t_inputrec *ir,
2999 const t_commrec *cr,
3005 int i, edinr, iupdate = 500;
3006 matrix rotmat; /* rotation matrix */
3007 rvec transvec; /* translation vector */
3008 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3009 real dt_1; /* 1/dt */
3010 struct t_do_edsam *buf;
3011 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3012 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3015 /* Check if ED sampling has to be performed */
3016 if (ed->eEDtype == EssentialDynamicsType::None)
3021 dt_1 = 1.0/ir->delta_t;
3023 /* Loop over all ED groups (usually one) */
3025 for (auto &edi : ed->edpar)
3028 if (bNeedDoEdsam(edi))
3031 buf = edi.buf->do_edsam;
3035 /* initialize radacc radius for slope criterion */
3036 buf->oldrad = calc_radius(edi.vecs.radacc);
3039 /* Copy the positions into buf->xc* arrays and after ED
3040 * feed back corrections to the official positions */
3042 /* Broadcast the ED positions such that every node has all of them
3043 * Every node contributes its local positions xs and stores it in
3044 * the collective buf->xcoll array. Note that for edinr > 1
3045 * xs could already have been modified by an earlier ED */
3047 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3048 edi.sav.nr, edi.sav.nr_loc, edi.sav.anrs_loc, edi.sav.c_ind, edi.sav.x_old, box);
3050 /* Only assembly reference positions if their indices differ from the average ones */
3053 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3054 edi.sref.nr, edi.sref.nr_loc, edi.sref.anrs_loc, edi.sref.c_ind, edi.sref.x_old, box);
3057 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3058 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3059 * set bUpdateShifts=TRUE in the parallel case. */
3060 buf->bUpdateShifts = FALSE;
3062 /* Now all nodes have all of the ED positions in edi.sav->xcoll,
3063 * as well as the indices in edi.sav.anrs */
3065 /* Fit the reference indices to the reference structure */
3068 fit_to_reference(buf->xcoll, transvec, rotmat, &edi);
3072 fit_to_reference(buf->xc_ref, transvec, rotmat, &edi);
3075 /* Now apply the translation and rotation to the ED structure */
3076 translate_and_rotate(buf->xcoll, edi.sav.nr, transvec, rotmat);
3078 /* Find out how well we fit to the reference (just for output steps) */
3079 if (do_per_step(step, edi.outfrq) && MASTER(cr))
3083 /* Indices of reference and average structures are identical,
3084 * thus we can calculate the rmsd to SREF using xcoll */
3085 rmsdev = rmsd_from_structure(buf->xcoll, &edi.sref);
3089 /* We have to translate & rotate the reference atoms first */
3090 translate_and_rotate(buf->xc_ref, edi.sref.nr, transvec, rotmat);
3091 rmsdev = rmsd_from_structure(buf->xc_ref, &edi.sref);
3095 /* update radsam references, when required */
3096 if (do_per_step(step, edi.maxedsteps) && step >= edi.presteps)
3098 project(buf->xcoll, &edi);
3099 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3100 rad_project(edi, buf->xcoll, &edi.vecs.radfix);
3101 buf->oldrad = -1.e5;
3104 /* update radacc references, when required */
3105 if (do_per_step(step, iupdate) && step >= edi.presteps)
3107 edi.vecs.radacc.radius = calc_radius(edi.vecs.radacc);
3108 if (edi.vecs.radacc.radius - buf->oldrad < edi.slope)
3110 project(buf->xcoll, &edi);
3111 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3116 buf->oldrad = edi.vecs.radacc.radius;
3120 /* apply the constraints */
3121 if (step >= edi.presteps && ed_constraints(ed->eEDtype, edi))
3123 /* ED constraints should be applied already in the first MD step
3124 * (which is step 0), therefore we pass step+1 to the routine */
3125 ed_apply_constraints(buf->xcoll, &edi, step+1 - ir->init_step);
3128 /* write to edo, when required */
3129 if (do_per_step(step, edi.outfrq))
3131 project(buf->xcoll, &edi);
3132 if (MASTER(cr) && !bSuppress)
3134 write_edo(edi, ed->edo, rmsdev);
3138 /* Copy back the positions unless monitoring only */
3139 if (ed_constraints(ed->eEDtype, edi))
3141 /* remove fitting */
3142 rmfit(edi.sav.nr, buf->xcoll, transvec, rotmat);
3144 /* Copy the ED corrected positions into the coordinate array */
3145 /* Each node copies its local part. In the serial case, nat_loc is the
3146 * total number of ED atoms */
3147 for (i = 0; i < edi.sav.nr_loc; i++)
3149 /* Unshift local ED coordinate and store in x_unsh */
3150 ed_unshift_single_coord(box, buf->xcoll[edi.sav.c_ind[i]],
3151 buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh);
3153 /* dx is the ED correction to the positions: */
3154 rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx);
3158 /* dv is the ED correction to the velocity: */
3159 svmul(dt_1, dx, dv);
3160 /* apply the velocity correction: */
3161 rvec_inc(v[edi.sav.anrs_loc[i]], dv);
3163 /* Finally apply the position correction due to ED: */
3164 copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]);
3167 } /* END of if ( bNeedDoEdsam(edi) ) */
3169 /* Prepare for the next ED group */
3171 } /* END of loop over ED groups */