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/smalloc.h"
77 #include "gromacs/utility/strconvert.h"
81 /*! \brief Identifies the type of ED: none, normal ED, flooding. */
82 enum class EssentialDynamicsType
84 //! No essential dynamics
86 //! Essential dynamics sampling
88 //! Essential dynamics flooding
92 /*! \brief Identify on which structure to operate. */
93 enum class EssentialDynamicsStructure
95 //! Reference structure
105 /*! \brief Essential dynamics vector.
106 * TODO: split into working data and input data
107 * NOTE: called eigvec, because traditionally eigenvectors from PCA analysis
108 * were used as essential dynamics vectors, however, vectors used for ED need
109 * not have any special properties
113 //! nr of eigenvectors
115 //! index nrs of eigenvectors
117 //! stepsizes (per eigenvector)
118 real *stpsz = nullptr;
119 //! eigenvector components
120 rvec **vec = nullptr;
121 //! instantaneous x projections
122 real *xproj = nullptr;
123 //! instantaneous f projections
124 real *fproj = nullptr;
125 //! instantaneous radius
127 //! starting or target projections
128 real *refproj = nullptr;
131 /*! \brief Essential dynamics vectors per method implementation.
135 //! only monitored, no constraints
137 //! fixed linear constraints
138 t_eigvec linfix = {};
139 //! acceptance linear constraints
140 t_eigvec linacc = {};
141 //! fixed radial constraints (exp)
142 t_eigvec radfix = {};
143 //! acceptance radial constraints (exp)
144 t_eigvec radacc = {};
145 //! acceptance rad. contraction constr.
146 t_eigvec radcon = {};
149 /*! \brief Essential dynamics flooding parameters and work data.
150 * TODO: split into working data and input parameters
151 * NOTE: The implementation here follows:
152 * O.E. Lange, L.V. Schafer, and H. Grubmuller,
153 * “Flooding in GROMACS: Accelerated barrier crossings in molecular dynamics,”
154 * J. Comp. Chem., 27 1693–1702 (2006)
158 /*! \brief Target destabilisation free energy.
159 * Controls the flooding potential strength.
160 * Linked to the expected speed-up of mean first passage time out of flooded minimum */
162 //! Do not calculate a flooding potential, instead flood with a constant force
163 bool bConstForce = false;
164 //! Relaxation time scale for slowest flooded degree of freedom
166 //! Current estimated destabilisation free energy
168 //! Flooding energy, acting as a proportionality factor for the flooding potential
170 //! Boltzmann constant times temperature, provided by user
172 //! The flooding potential
174 //! Integration time step
176 //! Inital flooding strenth
178 //! Empirical global scaling parameter of the essential dynamics vectors.
180 //! The forces from flooding in atom coordinate space (in contrast to projections onto essential dynamics vectors)
181 rvec *forces_cartesian = nullptr;
182 //! The vectors along which to flood
184 //! Use flooding for harmonic restraint on the eigenvector
185 bool bHarmonic = false;
186 //! The initial reference projection of the flooding vectors. Only with harmonic restraint.
187 real *initialReferenceProjection = nullptr;
188 //! The current reference projection is the initialReferenceProjection + step * slope. Only with harmonic restraint.
189 real *referenceProjectionSlope = nullptr;
194 /* This type is for the average, reference, target, and origin structure */
197 int nr = 0; /* number of atoms this structure contains */
198 int nr_loc = 0; /* number of atoms on local node */
199 int *anrs = nullptr; /* atom index numbers */
200 int *anrs_loc = nullptr; /* local atom index numbers */
201 int nalloc_loc = 0; /* allocation size of anrs_loc */
202 int *c_ind = nullptr; /* at which position of the whole anrs
203 * array is a local atom?, i.e.
204 * c_ind[0...nr_loc-1] gives the atom index
205 * with respect to the collective
206 * anrs[0...nr-1] array */
207 rvec *x = nullptr; /* positions for this structure */
208 rvec *x_old = nullptr; /* Last positions which have the correct PBC
209 representation of the ED group. In
210 combination with keeping track of the
211 shift vectors, the ED group can always
213 real *m = nullptr; /* masses */
214 real mtot = 0.; /* total mass (only used in sref) */
215 real *sqrtm = nullptr; /* sqrt of the masses used for mass-
216 * weighting of analysis (only used in sav) */
222 int nini = 0; /* total Nr of atoms */
223 gmx_bool fitmas = false; /* true if trans fit with cm */
224 gmx_bool pcamas = false; /* true if mass-weighted PCA */
225 int presteps = 0; /* number of steps to run without any
226 * perturbations ... just monitoring */
227 int outfrq = 0; /* freq (in steps) of writing to edo */
228 int maxedsteps = 0; /* max nr of steps per cycle */
230 /* all gmx_edx datasets are copied to all nodes in the parallel case */
231 gmx_edx sref = {}; /* reference positions, to these fitting
233 gmx_bool bRefEqAv = false; /* If true, reference & average indices
234 * are the same. Used for optimization */
235 gmx_edx sav = {}; /* average positions */
236 gmx_edx star = {}; /* target positions */
237 gmx_edx sori = {}; /* origin positions */
239 t_edvecs vecs = {}; /* eigenvectors */
240 real slope = 0; /* minimal slope in acceptance radexp */
242 t_edflood flood = {}; /* parameters especially for flooding */
243 struct t_ed_buffer *buf = nullptr; /* handle to local buffers */
251 EssentialDynamicsType eEDtype = EssentialDynamicsType::None;
252 //! output file pointer
254 std::vector<t_edpar> edpar;
255 gmx_bool bFirst = false;
257 gmx_edsam::~gmx_edsam()
261 /* edo is opened sometimes with xvgropen, sometimes with
262 * gmx_fio_fopen, so we use the least common denominator for
272 rvec old_transvec, older_transvec, transvec_compact;
273 rvec *xcoll; /* Positions from all nodes, this is the
274 collective set we work on.
275 These are the positions of atoms with
276 average structure indices */
277 rvec *xc_ref; /* same but with reference structure indices */
278 ivec *shifts_xcoll; /* Shifts for xcoll */
279 ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
280 ivec *shifts_xc_ref; /* Shifts for xc_ref */
281 ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
282 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
283 ED shifts for this ED group need to
288 /* definition of ED buffer structure */
291 struct t_fit_to_ref * fit_to_ref;
292 struct t_do_edfit * do_edfit;
293 struct t_do_edsam * do_edsam;
294 struct t_do_radcon * do_radcon;
299 class EssentialDynamics::Impl
302 // TODO: move all fields from gmx_edsam here and remove gmx_edsam
303 gmx_edsam essentialDynamics_;
305 EssentialDynamics::EssentialDynamics() : impl_(new Impl)
308 EssentialDynamics::~EssentialDynamics() = default;
310 gmx_edsam *EssentialDynamics::getLegacyED()
312 return &impl_->essentialDynamics_;
316 /* Function declarations */
317 static void fit_to_reference(rvec *xcoll, rvec transvec, matrix rotmat, t_edpar *edi);
318 static void translate_and_rotate(rvec *x, int nat, rvec transvec, matrix rotmat);
319 static real rmsd_from_structure(rvec *x, struct gmx_edx *s);
322 /*! \brief Read in the essential dynamics input file and return its contents.
323 * \param[in] fn the file name of the edi file to be read
324 * \param[in] nr_mdatoms the number of atoms in the simulation
325 * \returns A vector of containing the essentiail dyanmics parameters.
326 * NOTE: edi files may that it may contain several ED data sets from concatenated edi files.
327 * The standard case would be a single ED data set, though. */
328 std::vector<t_edpar> read_edi_file(const char *fn, int nr_mdatoms);
330 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam &ed, edsamhistory_t *EDstate);
331 static void init_edsamstate(const gmx_edsam &ed, edsamhistory_t *EDstate);
332 static void write_edo_legend(gmx_edsam * ed, int nED, const gmx_output_env_t *oenv);
333 /* End function declarations */
335 /* Define formats for the column width in the output file */
336 const char EDcol_sfmt[] = "%17s";
337 const char EDcol_efmt[] = "%17.5e";
338 const char EDcol_ffmt[] = "%17f";
340 /* Define a formats for reading, otherwise cppcheck complains for scanf without width limits */
341 const char max_ev_fmt_d[] = "%7d"; /* Number of eigenvectors. 9,999,999 should be enough */
342 const char max_ev_fmt_lf[] = "%12lf";
343 const char max_ev_fmt_dlf[] = "%7d%12lf";
344 const char max_ev_fmt_dlflflf[] = "%7d%12lf%12lf%12lf";
345 const char max_ev_fmt_lelele[] = "%12le%12le%12le";
349 /*!\brief Determines whether to perform essential dynamics constraints other than flooding.
350 * \param[in] edi the essential dynamics parameters
351 * \returns true if essential dyanmics constraints need to be performed
353 bool bNeedDoEdsam(const t_edpar &edi)
355 return (edi.vecs.mon.neig != 0)
356 || (edi.vecs.linfix.neig != 0)
357 || (edi.vecs.linacc.neig != 0)
358 || (edi.vecs.radfix.neig != 0)
359 || (edi.vecs.radacc.neig != 0)
360 || (edi.vecs.radcon.neig != 0);
365 /* Multiple ED groups will be labeled with letters instead of numbers
366 * to avoid confusion with eigenvector indices */
367 static char get_EDgroupChar(int nr_edi, int nED)
375 * nr_edi = 2 -> B ...
377 return 'A' + nr_edi - 1;
382 /*! \brief The mass-weighted inner product of two coordinate vectors.
383 * Does not subtract average positions, projection on single eigenvector is returned
384 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
385 * Average position is subtracted in ed_apply_constraints prior to calling projectx
386 * \param[in] edi Essential dynamics parameters
387 * \param[in] xcoll vector of atom coordinates
388 * \param[in] vec vector of coordinates to project onto
389 * \return mass-weighted projection.
391 real projectx(const t_edpar &edi, rvec *xcoll, rvec *vec)
397 for (i = 0; i < edi.sav.nr; i++)
399 proj += edi.sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
404 /*!\brief Project coordinates onto vector after substracting average position.
405 * projection is stored in vec->refproj which is used for radacc, radfix,
406 * radcon and center of flooding potential.
407 * Average positions are first substracted from x, then added back again.
408 * \param[in] edi essential dynamics parameters with average position
409 * \param[in] x Coordinates to be projected
410 * \param[out] vec eigenvector, radius and refproj are overwritten here
412 void rad_project(const t_edpar &edi, rvec *x, t_eigvec *vec)
417 /* Subtract average positions */
418 for (i = 0; i < edi.sav.nr; i++)
420 rvec_dec(x[i], edi.sav.x[i]);
423 for (i = 0; i < vec->neig; i++)
425 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
426 rad += gmx::square((vec->refproj[i]-vec->xproj[i]));
428 vec->radius = sqrt(rad);
430 /* Add average positions */
431 for (i = 0; i < edi.sav.nr; i++)
433 rvec_inc(x[i], edi.sav.x[i]);
437 /*!\brief Projects coordinates onto eigenvectors and stores result in vec->xproj.
438 * Mass-weighting is applied. Subtracts average positions prior to projection and add
439 * them afterwards to retain the unchanged vector.
440 * \param[in] x The coordinates to project to an eigenvector
441 * \param[in,out] vec The eigenvectors
442 * \param[in] edi essential dynamics parameters holding average structure and masses
444 void project_to_eigvectors(rvec *x, t_eigvec *vec, const t_edpar &edi)
451 /* Subtract average positions */
452 for (int i = 0; i < edi.sav.nr; i++)
454 rvec_dec(x[i], edi.sav.x[i]);
457 for (int i = 0; i < vec->neig; i++)
459 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
462 /* Add average positions */
463 for (int i = 0; i < edi.sav.nr; i++)
465 rvec_inc(x[i], edi.sav.x[i]);
470 /* Project vector x onto all edi->vecs (mon, linfix,...) */
471 static void project(rvec *x, /* positions to project */
472 t_edpar *edi) /* edi data set */
474 /* It is not more work to subtract the average position in every
475 * subroutine again, because these routines are rarely used simultaneously */
476 project_to_eigvectors(x, &edi->vecs.mon, *edi);
477 project_to_eigvectors(x, &edi->vecs.linfix, *edi);
478 project_to_eigvectors(x, &edi->vecs.linacc, *edi);
479 project_to_eigvectors(x, &edi->vecs.radfix, *edi);
480 project_to_eigvectors(x, &edi->vecs.radacc, *edi);
481 project_to_eigvectors(x, &edi->vecs.radcon, *edi);
486 /*!\brief Evaluates the distance from reference to current eigenvector projection.
487 * \param[in] vec eigenvector
490 real calc_radius(const t_eigvec &vec)
494 for (int i = 0; i < vec.neig; i++)
496 rad += gmx::square((vec.refproj[i]-vec.xproj[i]));
499 return rad = sqrt(rad);
508 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
510 /* this is a copy of do_fit with some modifications */
511 int c, r, n, j, i, irot;
512 double d[6], xnr, xpc;
517 struct t_do_edfit *loc;
520 if (edi->buf->do_edfit != nullptr)
527 snew(edi->buf->do_edfit, 1);
529 loc = edi->buf->do_edfit;
533 snew(loc->omega, 2*DIM);
534 snew(loc->om, 2*DIM);
535 for (i = 0; i < 2*DIM; i++)
537 snew(loc->omega[i], 2*DIM);
538 snew(loc->om[i], 2*DIM);
542 for (i = 0; (i < 6); i++)
545 for (j = 0; (j < 6); j++)
547 loc->omega[i][j] = 0;
552 /* calculate the matrix U */
554 for (n = 0; (n < natoms); n++)
556 for (c = 0; (c < DIM); c++)
559 for (r = 0; (r < DIM); r++)
567 /* construct loc->omega */
568 /* loc->omega is symmetric -> loc->omega==loc->omega' */
569 for (r = 0; (r < 6); r++)
571 for (c = 0; (c <= r); c++)
573 if ((r >= 3) && (c < 3))
575 loc->omega[r][c] = u[r-3][c];
576 loc->omega[c][r] = u[r-3][c];
580 loc->omega[r][c] = 0;
581 loc->omega[c][r] = 0;
586 /* determine h and k */
587 jacobi(loc->omega, 6, d, loc->om, &irot);
591 fprintf(stderr, "IROT=0\n");
594 index = 0; /* For the compiler only */
596 for (j = 0; (j < 3); j++)
599 for (i = 0; (i < 6); i++)
608 for (i = 0; (i < 3); i++)
610 vh[j][i] = M_SQRT2*loc->om[i][index];
611 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
616 for (c = 0; (c < 3); c++)
618 for (r = 0; (r < 3); r++)
620 R[c][r] = vk[0][r]*vh[0][c]+
627 for (c = 0; (c < 3); c++)
629 for (r = 0; (r < 3); r++)
631 R[c][r] = vk[0][r]*vh[0][c]+
640 static void rmfit(int nat, rvec *xcoll, const rvec transvec, matrix rotmat)
647 * The inverse rotation is described by the transposed rotation matrix */
648 transpose(rotmat, tmat);
649 rotate_x(xcoll, nat, tmat);
651 /* Remove translation */
652 vec[XX] = -transvec[XX];
653 vec[YY] = -transvec[YY];
654 vec[ZZ] = -transvec[ZZ];
655 translate_x(xcoll, nat, vec);
659 /**********************************************************************************
660 ******************** FLOODING ****************************************************
661 **********************************************************************************
663 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
664 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
665 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
667 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
668 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
670 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
671 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
673 do_flood makes a copy of the positions,
674 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
675 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
676 fit. Then do_flood adds these forces to the forcefield-forces
677 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
679 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
680 structure is projected to the system of eigenvectors and then this position in the subspace is used as
681 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
682 i.e. the average structure as given in the make_edi file.
684 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
685 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
686 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
687 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
688 further adaption is applied, Efl will stay constant at zero.
690 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
691 used as spring constants for the harmonic potential.
692 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
693 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
695 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
696 the routine read_edi_file reads all of theses flooding files.
697 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
698 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
699 edi there is no interdependence whatsoever. The forces are added together.
701 To write energies into the .edr file, call the function
702 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
704 get_flood_energies(real Vfl[],int nnames);
707 - 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.
709 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
710 two edsam files from two peptide chains
713 // TODO split this into multiple files
716 /*!\brief Output flooding simulation settings and results to file.
717 * \param[in] edi Essential dynamics input parameters
718 * \param[in] fp output file
719 * \param[in] rmsd rmsd to reference structure
721 void write_edo_flood(const t_edpar &edi, FILE *fp, real rmsd)
723 /* Output how well we fit to the reference structure */
724 fprintf(fp, EDcol_ffmt, rmsd);
726 for (int i = 0; i < edi.flood.vecs.neig; i++)
728 fprintf(fp, EDcol_efmt, edi.flood.vecs.xproj[i]);
730 /* Check whether the reference projection changes with time (this can happen
731 * in case flooding is used as harmonic restraint). If so, output the
732 * current reference projection */
733 if (edi.flood.bHarmonic && edi.flood.referenceProjectionSlope[i] != 0.0)
735 fprintf(fp, EDcol_efmt, edi.flood.vecs.refproj[i]);
738 /* Output Efl if we are doing adaptive flooding */
739 if (0 != edi.flood.tau)
741 fprintf(fp, EDcol_efmt, edi.flood.Efl);
743 fprintf(fp, EDcol_efmt, edi.flood.Vfl);
745 /* Output deltaF if we are doing adaptive flooding */
746 if (0 != edi.flood.tau)
748 fprintf(fp, EDcol_efmt, edi.flood.deltaF);
750 fprintf(fp, EDcol_efmt, edi.flood.vecs.fproj[i]);
756 /* From flood.xproj compute the Vfl(x) at this point */
757 static real flood_energy(t_edpar *edi, int64_t step)
759 /* compute flooding energy Vfl
760 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
761 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
762 it is already computed by make_edi and stored in stpsz[i]
764 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
771 /* Each time this routine is called (i.e. each time step), we add a small
772 * value to the reference projection. This way a harmonic restraint towards
773 * a moving reference is realized. If no value for the additive constant
774 * is provided in the edi file, the reference will not change. */
775 if (edi->flood.bHarmonic)
777 for (i = 0; i < edi->flood.vecs.neig; i++)
779 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i] + step * edi->flood.referenceProjectionSlope[i];
784 /* Compute sum which will be the exponent of the exponential */
785 for (i = 0; i < edi->flood.vecs.neig; i++)
787 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
788 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]);
791 /* Compute the Gauss function*/
792 if (edi->flood.bHarmonic)
794 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
798 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*std::exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
805 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
806 static void flood_forces(t_edpar *edi)
808 /* compute the forces in the subspace of the flooding eigenvectors
809 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
812 real energy = edi->flood.Vfl;
815 if (edi->flood.bHarmonic)
817 for (i = 0; i < edi->flood.vecs.neig; i++)
819 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
824 for (i = 0; i < edi->flood.vecs.neig; i++)
826 /* if Efl is zero the forces are zero if not use the formula */
827 edi->flood.vecs.fproj[i] = 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;
834 /*!\brief Raise forces from subspace into cartesian space.
835 * This function lifts the forces from the subspace to the cartesian space
836 * all the values not contained in the subspace are assumed to be zero and then
837 * a coordinate transformation from eigenvector to cartesian vectors is performed
838 * The nonexistent values don't have to be set to zero explicitly, they would occur
839 * as zero valued summands, hence we just stop to compute this part of the sum.
840 * For every atom we add all the contributions to this atom from all the different eigenvectors.
841 * NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
842 * field forces_cart prior the computation, but we compute the forces separately
843 * to have them accessible for diagnostics
845 * \param[in] edi Essential dynamics input parameters
846 * \param[out] forces_cart The cartesian forces
849 void flood_blowup(const t_edpar &edi, rvec *forces_cart)
851 const real * forces_sub = edi.flood.vecs.fproj;
852 /* Calculate the cartesian forces for the local atoms */
854 /* Clear forces first */
855 for (int j = 0; j < edi.sav.nr_loc; j++)
857 clear_rvec(forces_cart[j]);
860 /* Now compute atomwise */
861 for (int j = 0; j < edi.sav.nr_loc; j++)
863 /* Compute forces_cart[edi.sav.anrs[j]] */
864 for (int eig = 0; eig < edi.flood.vecs.neig; eig++)
867 /* Force vector is force * eigenvector (compute only atom j) */
868 svmul(forces_sub[eig], edi.flood.vecs.vec[eig][edi.sav.c_ind[j]], addedForce);
869 /* Add this vector to the cartesian forces */
870 rvec_inc(forces_cart[j], addedForce);
878 /* Update the values of Efl, deltaF depending on tau and Vfl */
879 static void update_adaption(t_edpar *edi)
881 /* this function updates the parameter Efl and deltaF according to the rules given in
882 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
885 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
887 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
888 /* check if restrain (inverted flooding) -> don't let EFL become positive */
889 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
894 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
899 static void do_single_flood(
907 gmx_bool bNS) /* Are we in a neighbor searching step? */
910 matrix rotmat; /* rotation matrix */
911 matrix tmat; /* inverse rotation */
912 rvec transvec; /* translation vector */
914 struct t_do_edsam *buf;
917 buf = edi->buf->do_edsam;
920 /* Broadcast the positions of the AVERAGE structure such that they are known on
921 * every processor. Each node contributes its local positions x and stores them in
922 * the collective ED array buf->xcoll */
923 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
924 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
926 /* Only assembly REFERENCE positions if their indices differ from the average ones */
929 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
930 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
933 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
934 * We do not need to update the shifts until the next NS step */
935 buf->bUpdateShifts = FALSE;
937 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
938 * as well as the indices in edi->sav.anrs */
940 /* Fit the reference indices to the reference structure */
943 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
947 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
950 /* Now apply the translation and rotation to the ED structure */
951 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
953 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
954 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
956 if (!edi->flood.bConstForce)
958 /* Compute Vfl(x) from flood.xproj */
959 edi->flood.Vfl = flood_energy(edi, step);
961 update_adaption(edi);
963 /* Compute the flooding forces */
967 /* Translate them into cartesian positions */
968 flood_blowup(*edi, edi->flood.forces_cartesian);
970 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
971 /* Each node rotates back its local forces */
972 transpose(rotmat, tmat);
973 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
975 /* Finally add forces to the main force variable */
976 for (i = 0; i < edi->sav.nr_loc; i++)
978 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
981 /* Output is written by the master process */
982 if (do_per_step(step, edi->outfrq) && MASTER(cr))
984 /* Output how well we fit to the reference */
987 /* Indices of reference and average structures are identical,
988 * thus we can calculate the rmsd to SREF using xcoll */
989 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
993 /* We have to translate & rotate the reference atoms first */
994 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
995 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
998 write_edo_flood(*edi, edo, rmsdev);
1003 /* Main flooding routine, called from do_force */
1004 extern void do_flood(const t_commrec *cr,
1005 const t_inputrec *ir,
1013 /* Write time to edo, when required. Output the time anyhow since we need
1014 * it in the output file for ED constraints. */
1015 if (MASTER(cr) && do_per_step(step, ed->edpar.begin()->outfrq))
1017 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
1020 if (ed->eEDtype != EssentialDynamicsType::Flooding)
1025 for (auto &edi : ed->edpar)
1027 /* Call flooding for one matrix */
1028 if (edi.flood.vecs.neig)
1030 do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS);
1036 /* Called by init_edi, configure some flooding related variables and structures,
1037 * print headers to output files */
1038 static void init_flood(t_edpar *edi, gmx_edsam * ed, real dt)
1043 edi->flood.Efl = edi->flood.constEfl;
1047 if (edi->flood.vecs.neig)
1049 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1050 ed->eEDtype = EssentialDynamicsType::Flooding;
1052 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
1054 if (edi->flood.bConstForce)
1056 /* We have used stpsz as a vehicle to carry the fproj values for constant
1057 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1058 * in const force flooding, fproj is never changed. */
1059 for (i = 0; i < edi->flood.vecs.neig; i++)
1061 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1063 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1064 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1072 /*********** Energy book keeping ******/
1073 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1082 srenew(names, count);
1083 sprintf(buf, "Vfl_%d", count);
1084 names[count-1] = gmx_strdup(buf);
1085 actual = actual->next_edi;
1092 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1094 /*fl has to be big enough to capture nnames-many entries*/
1103 Vfl[count-1] = actual->flood.Vfl;
1104 actual = actual->next_edi;
1107 if (nnames != count-1)
1109 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1112 /************* END of FLOODING IMPLEMENTATION ****************************/
1116 /* This function opens the ED input and output files, reads in all datasets it finds
1117 * in the input file, and cross-checks whether the .edi file information is consistent
1118 * with the essential dynamics data found in the checkpoint file (if present).
1119 * gmx make_edi can be used to create an .edi input file.
1121 static std::unique_ptr<gmx::EssentialDynamics> ed_open(
1123 ObservablesHistory *oh,
1124 const char *ediFileName,
1125 const char *edoFileName,
1127 const gmx_output_env_t *oenv,
1128 const t_commrec *cr)
1130 auto edHandle = gmx::compat::make_unique<gmx::EssentialDynamics>();
1131 auto ed = edHandle->getLegacyED();
1132 /* We want to perform ED (this switch might later be upgraded to EssentialDynamicsType::Flooding) */
1133 ed->eEDtype = EssentialDynamicsType::EDSampling;
1137 // If we start from a checkpoint file, we already have an edsamHistory struct
1138 if (oh->edsamHistory == nullptr)
1140 oh->edsamHistory = gmx::compat::make_unique<edsamhistory_t>(edsamhistory_t {});
1142 edsamhistory_t *EDstate = oh->edsamHistory.get();
1144 /* Read the edi input file: */
1145 ed->edpar = read_edi_file(ediFileName, natoms);
1147 /* Make sure the checkpoint was produced in a run using this .edi file */
1148 if (EDstate->bFromCpt)
1150 crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1154 EDstate->nED = ed->edpar.size();
1156 init_edsamstate(*ed, EDstate);
1158 /* The master opens the ED output file */
1161 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1165 ed->edo = xvgropen(edoFileName,
1166 "Essential dynamics / flooding output",
1168 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1170 /* Make a descriptive legend */
1171 write_edo_legend(ed, EDstate->nED, oenv);
1178 /* Broadcasts the structure data */
1179 static void bc_ed_positions(const t_commrec *cr, struct gmx_edx *s, EssentialDynamicsStructure stype)
1181 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1182 snew_bc(cr, s->x, s->nr ); /* Positions */
1183 nblock_bc(cr, s->nr, s->anrs );
1184 nblock_bc(cr, s->nr, s->x );
1186 /* For the average & reference structures we need an array for the collective indices,
1187 * and we need to broadcast the masses as well */
1188 if (stype == EssentialDynamicsStructure::Average || stype == EssentialDynamicsStructure::Reference)
1190 /* We need these additional variables in the parallel case: */
1191 snew(s->c_ind, s->nr ); /* Collective indices */
1192 /* Local atom indices get assigned in dd_make_local_group_indices.
1193 * There, also memory is allocated */
1194 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1195 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1196 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1199 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1200 if (stype == EssentialDynamicsStructure::Reference)
1202 snew_bc(cr, s->m, s->nr);
1203 nblock_bc(cr, s->nr, s->m);
1206 /* For the average structure we might need the masses for mass-weighting */
1207 if (stype == EssentialDynamicsStructure::Average)
1209 snew_bc(cr, s->sqrtm, s->nr);
1210 nblock_bc(cr, s->nr, s->sqrtm);
1211 snew_bc(cr, s->m, s->nr);
1212 nblock_bc(cr, s->nr, s->m);
1217 /* Broadcasts the eigenvector data */
1218 static void bc_ed_vecs(const t_commrec *cr, t_eigvec *ev, int length)
1222 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1223 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1224 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1225 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1226 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1228 nblock_bc(cr, ev->neig, ev->ieig );
1229 nblock_bc(cr, ev->neig, ev->stpsz );
1230 nblock_bc(cr, ev->neig, ev->xproj );
1231 nblock_bc(cr, ev->neig, ev->fproj );
1232 nblock_bc(cr, ev->neig, ev->refproj);
1234 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1235 for (i = 0; i < ev->neig; i++)
1237 snew_bc(cr, ev->vec[i], length);
1238 nblock_bc(cr, length, ev->vec[i]);
1244 /* Broadcasts the ED / flooding data to other nodes
1245 * and allocates memory where needed */
1246 static void broadcast_ed_data(const t_commrec *cr, gmx_edsam * ed)
1248 /* Master lets the other nodes know if its ED only or also flooding */
1249 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1251 int numedis = ed->edpar.size();
1252 /* First let everybody know how many ED data sets to expect */
1253 gmx_bcast(sizeof(numedis), &numedis, cr);
1254 nblock_abc(cr, numedis, &(ed->edpar));
1255 /* Now transfer the ED data set(s) */
1256 for (auto &edi : ed->edpar)
1258 /* Broadcast a single ED data set */
1261 /* Broadcast positions */
1262 bc_ed_positions(cr, &(edi.sref), EssentialDynamicsStructure::Reference); /* reference positions (don't broadcast masses) */
1263 bc_ed_positions(cr, &(edi.sav ), EssentialDynamicsStructure::Average ); /* average positions (do broadcast masses as well) */
1264 bc_ed_positions(cr, &(edi.star), EssentialDynamicsStructure::Target); /* target positions */
1265 bc_ed_positions(cr, &(edi.sori), EssentialDynamicsStructure::Origin); /* origin positions */
1267 /* Broadcast eigenvectors */
1268 bc_ed_vecs(cr, &edi.vecs.mon, edi.sav.nr);
1269 bc_ed_vecs(cr, &edi.vecs.linfix, edi.sav.nr);
1270 bc_ed_vecs(cr, &edi.vecs.linacc, edi.sav.nr);
1271 bc_ed_vecs(cr, &edi.vecs.radfix, edi.sav.nr);
1272 bc_ed_vecs(cr, &edi.vecs.radacc, edi.sav.nr);
1273 bc_ed_vecs(cr, &edi.vecs.radcon, edi.sav.nr);
1274 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1275 bc_ed_vecs(cr, &edi.flood.vecs, edi.sav.nr);
1277 /* For harmonic restraints the reference projections can change with time */
1278 if (edi.flood.bHarmonic)
1280 snew_bc(cr, edi.flood.initialReferenceProjection, edi.flood.vecs.neig);
1281 snew_bc(cr, edi.flood.referenceProjectionSlope, edi.flood.vecs.neig);
1282 nblock_bc(cr, edi.flood.vecs.neig, edi.flood.initialReferenceProjection);
1283 nblock_bc(cr, edi.flood.vecs.neig, edi.flood.referenceProjectionSlope);
1289 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1290 static void init_edi(const gmx_mtop_t *mtop, t_edpar *edi)
1293 real totalmass = 0.0;
1296 /* NOTE Init_edi is executed on the master process only
1297 * The initialized data sets are then transmitted to the
1298 * other nodes in broadcast_ed_data */
1300 /* evaluate masses (reference structure) */
1301 snew(edi->sref.m, edi->sref.nr);
1303 for (i = 0; i < edi->sref.nr; i++)
1307 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1311 edi->sref.m[i] = 1.0;
1314 /* Check that every m > 0. Bad things will happen otherwise. */
1315 if (edi->sref.m[i] <= 0.0)
1317 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1318 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1319 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1320 "atoms from the reference structure by creating a proper index group.\n",
1321 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1324 totalmass += edi->sref.m[i];
1326 edi->sref.mtot = totalmass;
1328 /* Masses m and sqrt(m) for the average structure. Note that m
1329 * is needed if forces have to be evaluated in do_edsam */
1330 snew(edi->sav.sqrtm, edi->sav.nr );
1331 snew(edi->sav.m, edi->sav.nr );
1332 for (i = 0; i < edi->sav.nr; i++)
1334 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1337 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1341 edi->sav.sqrtm[i] = 1.0;
1344 /* Check that every m > 0. Bad things will happen otherwise. */
1345 if (edi->sav.sqrtm[i] <= 0.0)
1347 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1348 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1349 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1350 "atoms from the average structure by creating a proper index group.\n",
1351 i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1355 /* put reference structure in origin */
1356 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1360 translate_x(edi->sref.x, edi->sref.nr, com);
1362 /* Init ED buffer */
1367 static void check(const char *line, const char *label)
1369 if (!strstr(line, label))
1371 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1376 static int read_checked_edint(FILE *file, const char *label)
1378 char line[STRLEN+1];
1381 fgets2 (line, STRLEN, file);
1383 fgets2 (line, STRLEN, file);
1384 sscanf (line, max_ev_fmt_d, &idum);
1388 static bool read_checked_edbool(FILE *file, const char *label)
1390 return read_checked_edint(file, label) != 0;
1393 static int read_edint(FILE *file, bool *bEOF)
1395 char line[STRLEN+1];
1400 eof = fgets2 (line, STRLEN, file);
1406 eof = fgets2 (line, STRLEN, file);
1412 sscanf (line, max_ev_fmt_d, &idum);
1418 static real read_checked_edreal(FILE *file, const char *label)
1420 char line[STRLEN+1];
1424 fgets2 (line, STRLEN, file);
1426 fgets2 (line, STRLEN, file);
1427 sscanf (line, max_ev_fmt_lf, &rdum);
1428 return static_cast<real>(rdum); /* always read as double and convert to single */
1432 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1435 char line[STRLEN+1];
1439 for (i = 0; i < number; i++)
1441 fgets2 (line, STRLEN, file);
1442 sscanf (line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1443 anrs[i]--; /* we are reading FORTRAN indices */
1444 for (j = 0; j < 3; j++)
1446 x[i][j] = d[j]; /* always read as double and convert to single */
1453 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1454 * \param[in] in the file to read from
1455 * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1456 * \param[out] vec the eigen-vectors
1457 * \param[in] nEig the number of eigenvectors
1459 void scan_edvec(FILE *in, int numAtoms, rvec ***vec, int nEig)
1462 for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1464 snew((*vec)[iEigenvector], numAtoms);
1465 for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1467 char line[STRLEN+1];
1469 fgets2(line, STRLEN, in);
1470 sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1471 (*vec)[iEigenvector][iAtom][XX] = x;
1472 (*vec)[iEigenvector][iAtom][YY] = y;
1473 (*vec)[iEigenvector][iAtom][ZZ] = z;
1478 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1479 * \param[in] in input file
1480 * \param[in] nrAtoms number of atoms/coordinates
1481 * \param[out] tvec the eigenvector
1483 void read_edvec(FILE *in, int nrAtoms, t_eigvec *tvec)
1485 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1486 if (tvec->neig <= 0)
1491 snew(tvec->ieig, tvec->neig);
1492 snew(tvec->stpsz, tvec->neig);
1493 for (int i = 0; i < tvec->neig; i++)
1495 char line[STRLEN+1];
1496 fgets2(line, STRLEN, in);
1497 int numEigenVectors;
1499 const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1502 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1504 tvec->ieig[i] = numEigenVectors;
1505 tvec->stpsz[i] = stepSize;
1506 } /* end of loop over eigenvectors */
1508 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1510 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1512 * Eigenvector and its intial reference and slope are stored on the
1513 * same line in the .edi format. To avoid re-winding the .edi file,
1514 * the reference eigen-vector and reference data are read in one go.
1516 * \param[in] in input file
1517 * \param[in] nrAtoms number of atoms/coordinates
1518 * \param[out] tvec the eigenvector
1519 * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1520 * \param[out] referenceProjectionSlope The slope of the reference projections.
1522 bool readEdVecWithReferenceProjection(FILE *in, int nrAtoms, t_eigvec *tvec, real ** initialReferenceProjection, real ** referenceProjectionSlope)
1524 bool providesReference = false;
1525 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1526 if (tvec->neig <= 0)
1531 snew(tvec->ieig, tvec->neig);
1532 snew(tvec->stpsz, tvec->neig);
1533 snew(*initialReferenceProjection, tvec->neig);
1534 snew(*referenceProjectionSlope, tvec->neig);
1535 for (int i = 0; i < tvec->neig; i++)
1537 char line[STRLEN+1];
1538 fgets2 (line, STRLEN, in);
1539 int numEigenVectors;
1540 double stepSize = 0.;
1541 double referenceProjection = 0.;
1542 double referenceSlope = 0.;
1543 const int nscan = sscanf(line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize, &referenceProjection, &referenceSlope);
1544 /* Zero out values which were not scanned */
1548 /* Every 4 values read, including reference position */
1549 providesReference = true;
1552 /* A reference position is provided */
1553 providesReference = true;
1554 /* No value for slope, set to 0 */
1555 referenceSlope = 0.0;
1558 /* No values for reference projection and slope, set to 0 */
1559 referenceProjection = 0.0;
1560 referenceSlope = 0.0;
1563 gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1565 (*initialReferenceProjection)[i] = referenceProjection;
1566 (*referenceProjectionSlope)[i] = referenceSlope;
1568 tvec->ieig[i] = numEigenVectors;
1569 tvec->stpsz[i] = stepSize;
1570 } /* end of loop over eigenvectors */
1572 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1573 return providesReference;
1576 /*!\brief Allocate working memory for the eigenvectors.
1577 * \param[in,out] tvec eigenvector for which memory will be allocated
1579 void setup_edvec(t_eigvec *tvec)
1581 snew(tvec->xproj, tvec->neig);
1582 snew(tvec->fproj, tvec->neig);
1583 snew(tvec->refproj, tvec->neig);
1588 /* Check if the same atom indices are used for reference and average positions */
1589 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1594 /* If the number of atoms differs between the two structures,
1595 * they cannot be identical */
1596 if (sref.nr != sav.nr)
1601 /* Now that we know that both stuctures have the same number of atoms,
1602 * check if also the indices are identical */
1603 for (i = 0; i < sav.nr; i++)
1605 if (sref.anrs[i] != sav.anrs[i])
1610 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1618 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1619 constexpr int c_oldestUnsupportedMagicNumber = 666;
1620 //! Oldest (lowest) magic number indicating supported essential dynamics input
1621 constexpr int c_oldestSupportedMagicNumber = 669;
1622 //! Latest (highest) magic number indicating supported essential dynamics input
1623 constexpr int c_latestSupportedMagicNumber = 670;
1625 /*!\brief Set up essential dynamics work parameters.
1626 * \param[in] edi Essential dynamics working structure.
1628 void setup_edi(t_edpar *edi)
1630 snew(edi->sref.x_old, edi->sref.nr);
1631 edi->sref.sqrtm = nullptr;
1632 snew(edi->sav.x_old, edi->sav.nr);
1633 if (edi->star.nr > 0)
1635 edi->star.sqrtm = nullptr;
1637 if (edi->sori.nr > 0)
1639 edi->sori.sqrtm = nullptr;
1641 setup_edvec(&edi->vecs.linacc);
1642 setup_edvec(&edi->vecs.mon);
1643 setup_edvec(&edi->vecs.linfix);
1644 setup_edvec(&edi->vecs.linacc);
1645 setup_edvec(&edi->vecs.radfix);
1646 setup_edvec(&edi->vecs.radacc);
1647 setup_edvec(&edi->vecs.radcon);
1648 setup_edvec(&edi->flood.vecs);
1651 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1652 * \param[in] magicNumber indicates the essential dynamics input file version
1653 * \returns true if CONST_FORCE_FLOODING is supported
1655 bool ediFileHasConstForceFlooding(int magicNumber)
1657 return magicNumber > c_oldestSupportedMagicNumber;
1660 /*!\brief Reports reading success of the essential dynamics input file magic number.
1661 * \param[in] in pointer to input file
1662 * \param[out] magicNumber determines the edi file version
1663 * \returns true if setting file version from input file was successful.
1665 bool ediFileHasValidDataSet(FILE *in, int * magicNumber)
1667 /* the edi file is not free format, so expect problems if the input is corrupt. */
1668 bool reachedEndOfFile = true;
1669 /* check the magic number */
1670 *magicNumber = read_edint(in, &reachedEndOfFile);
1671 /* Check whether we have reached the end of the input file */
1672 return !reachedEndOfFile;
1675 /*!\brief Trigger fatal error if magic number is unsupported.
1676 * \param[in] magicNumber A magic number read from the edi file
1677 * \param[in] fn name of input file for error reporting
1679 void exitWhenMagicNumberIsInvalid(int magicNumber, const char * fn)
1682 if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1684 if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1686 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1690 gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1695 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1697 * \param[in] in input file
1698 * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1699 * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1700 * \param[in] fn the file name of the input file for error reporting
1701 * \returns edi essential dynamics parameters
1703 t_edpar read_edi(FILE* in, int nr_mdatoms, bool hasConstForceFlooding, const char *fn)
1707 /* check the number of atoms */
1708 edi.nini = read_edint(in, &bEOF);
1709 if (edi.nini != nr_mdatoms)
1711 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi.nini, nr_mdatoms);
1714 /* Done checking. For the rest we blindly trust the input */
1715 edi.fitmas = read_checked_edbool(in, "FITMAS");
1716 edi.pcamas = read_checked_edbool(in, "ANALYSIS_MAS");
1717 edi.outfrq = read_checked_edint(in, "OUTFRQ");
1718 edi.maxedsteps = read_checked_edint(in, "MAXLEN");
1719 edi.slope = read_checked_edreal(in, "SLOPECRIT");
1721 edi.presteps = read_checked_edint(in, "PRESTEPS");
1722 edi.flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1723 edi.flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1724 edi.flood.tau = read_checked_edreal(in, "TAU");
1725 edi.flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1726 edi.flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1727 edi.flood.kT = read_checked_edreal(in, "KT");
1728 edi.flood.bHarmonic = read_checked_edbool(in, "HARMONIC");
1729 if (hasConstForceFlooding)
1731 edi.flood.bConstForce = read_checked_edbool(in, "CONST_FORCE_FLOODING");
1735 edi.flood.bConstForce = FALSE;
1737 edi.sref.nr = read_checked_edint(in, "NREF");
1739 /* allocate space for reference positions and read them */
1740 snew(edi.sref.anrs, edi.sref.nr);
1741 snew(edi.sref.x, edi.sref.nr);
1742 read_edx(in, edi.sref.nr, edi.sref.anrs, edi.sref.x);
1744 /* average positions. they define which atoms will be used for ED sampling */
1745 edi.sav.nr = read_checked_edint(in, "NAV");
1746 snew(edi.sav.anrs, edi.sav.nr);
1747 snew(edi.sav.x, edi.sav.nr);
1748 read_edx(in, edi.sav.nr, edi.sav.anrs, edi.sav.x);
1750 /* Check if the same atom indices are used for reference and average positions */
1751 edi.bRefEqAv = check_if_same(edi.sref, edi.sav);
1755 read_edvec(in, edi.sav.nr, &edi.vecs.mon);
1756 read_edvec(in, edi.sav.nr, &edi.vecs.linfix);
1757 read_edvec(in, edi.sav.nr, &edi.vecs.linacc);
1758 read_edvec(in, edi.sav.nr, &edi.vecs.radfix);
1759 read_edvec(in, edi.sav.nr, &edi.vecs.radacc);
1760 read_edvec(in, edi.sav.nr, &edi.vecs.radcon);
1762 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1763 bool bHaveReference = false;
1764 if (edi.flood.bHarmonic)
1766 bHaveReference = readEdVecWithReferenceProjection(in, edi.sav.nr, &edi.flood.vecs, &edi.flood.initialReferenceProjection, &edi.flood.referenceProjectionSlope);
1770 read_edvec(in, edi.sav.nr, &edi.flood.vecs);
1773 /* target positions */
1774 edi.star.nr = read_edint(in, &bEOF);
1775 if (edi.star.nr > 0)
1777 snew(edi.star.anrs, edi.star.nr);
1778 snew(edi.star.x, edi.star.nr);
1779 read_edx(in, edi.star.nr, edi.star.anrs, edi.star.x);
1782 /* positions defining origin of expansion circle */
1783 edi.sori.nr = read_edint(in, &bEOF);
1784 if (edi.sori.nr > 0)
1788 /* Both an -ori structure and a at least one manual reference point have been
1789 * specified. That's ambiguous and probably not intentional. */
1790 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1791 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1793 snew(edi.sori.anrs, edi.sori.nr);
1794 snew(edi.sori.x, edi.sori.nr);
1795 read_edx(in, edi.sori.nr, edi.sori.anrs, edi.sori.x);
1800 std::vector<t_edpar> read_edi_file(const char *fn, int nr_mdatoms)
1802 std::vector<t_edpar> essentialDynamicsParameters;
1804 /* This routine is executed on the master only */
1806 /* Open the .edi parameter input file */
1807 in = gmx_fio_fopen(fn, "r");
1808 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1810 /* Now read a sequence of ED input parameter sets from the edi file */
1811 int ediFileMagicNumber;
1812 while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1814 exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1815 bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1816 auto edi = read_edi(in, nr_mdatoms, hasConstForceFlooding, fn);
1818 essentialDynamicsParameters.emplace_back(edi);
1820 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1821 /* We need to allocate space for the data: */
1823 if (essentialDynamicsParameters.empty())
1825 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1828 fprintf(stderr, "ED: Found %zu ED group%s.\n", essentialDynamicsParameters.size(), essentialDynamicsParameters.size() > 1 ? "s" : "");
1830 /* Close the .edi file again */
1833 return essentialDynamicsParameters;
1835 } // anonymous namespace
1838 struct t_fit_to_ref {
1839 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1842 /* Fit the current positions to the reference positions
1843 * Do not actually do the fit, just return rotation and translation.
1844 * Note that the COM of the reference structure was already put into
1845 * the origin by init_edi. */
1846 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1847 rvec transvec, /* The translation vector */
1848 matrix rotmat, /* The rotation matrix */
1849 t_edpar *edi) /* Just needed for do_edfit */
1851 rvec com; /* center of mass */
1853 struct t_fit_to_ref *loc;
1856 /* Allocate memory the first time this routine is called for each edi group */
1857 if (nullptr == edi->buf->fit_to_ref)
1859 snew(edi->buf->fit_to_ref, 1);
1860 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1862 loc = edi->buf->fit_to_ref;
1864 /* We do not touch the original positions but work on a copy. */
1865 for (i = 0; i < edi->sref.nr; i++)
1867 copy_rvec(xcoll[i], loc->xcopy[i]);
1870 /* Calculate the center of mass */
1871 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1873 transvec[XX] = -com[XX];
1874 transvec[YY] = -com[YY];
1875 transvec[ZZ] = -com[ZZ];
1877 /* Subtract the center of mass from the copy */
1878 translate_x(loc->xcopy, edi->sref.nr, transvec);
1880 /* Determine the rotation matrix */
1881 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1885 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1886 int nat, /* How many positions are there? */
1887 rvec transvec, /* The translation vector */
1888 matrix rotmat) /* The rotation matrix */
1891 translate_x(x, nat, transvec);
1894 rotate_x(x, nat, rotmat);
1898 /* Gets the rms deviation of the positions to the structure s */
1899 /* fit_to_structure has to be called before calling this routine! */
1900 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1901 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1907 for (i = 0; i < s->nr; i++)
1909 rmsd += distance2(s->x[i], x[i]);
1912 rmsd /= static_cast<real>(s->nr);
1919 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1921 if (ed->eEDtype != EssentialDynamicsType::None)
1923 /* Loop over ED groups */
1925 for (auto &edi : ed->edpar)
1927 /* Local atoms of the reference structure (for fitting), need only be assembled
1928 * if their indices differ from the average ones */
1931 dd_make_local_group_indices(dd->ga2la, edi.sref.nr, edi.sref.anrs,
1932 &edi.sref.nr_loc, &edi.sref.anrs_loc, &edi.sref.nalloc_loc, edi.sref.c_ind);
1935 /* Local atoms of the average structure (on these ED will be performed) */
1936 dd_make_local_group_indices(dd->ga2la, edi.sav.nr, edi.sav.anrs,
1937 &edi.sav.nr_loc, &edi.sav.anrs_loc, &edi.sav.nalloc_loc, edi.sav.c_ind);
1939 /* Indicate that the ED shift vectors for this structure need to be updated
1940 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1941 edi.buf->do_edsam->bUpdateShifts = TRUE;
1947 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1958 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1959 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1960 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1964 xu[XX] = x[XX]-tx*box[XX][XX];
1965 xu[YY] = x[YY]-ty*box[YY][YY];
1966 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1972 /*!\brief Apply fixed linear constraints to essential dynamics variable.
1973 * \param[in,out] xcoll The collected coordinates.
1974 * \param[in] edi the essential dynamics parameters
1975 * \param[in] step the current simulation step
1977 void do_linfix(rvec *xcoll, const t_edpar &edi, int64_t step)
1979 /* loop over linfix vectors */
1980 for (int i = 0; i < edi.vecs.linfix.neig; i++)
1982 /* calculate the projection */
1983 real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
1985 /* calculate the correction */
1986 real preFactor = edi.vecs.linfix.refproj[i] + step*edi.vecs.linfix.stpsz[i] - proj;
1988 /* apply the correction */
1989 preFactor /= edi.sav.sqrtm[i];
1990 for (int j = 0; j < edi.sav.nr; j++)
1992 rvec differenceVector;
1993 svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
1994 rvec_inc(xcoll[j], differenceVector);
1999 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
2000 * \param[in,out] xcoll The collected coordinates.
2001 * \param[in] edi the essential dynamics parameters
2003 void do_linacc(rvec *xcoll, t_edpar *edi)
2005 /* loop over linacc vectors */
2006 for (int i = 0; i < edi->vecs.linacc.neig; i++)
2008 /* calculate the projection */
2009 real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
2011 /* calculate the correction */
2012 real preFactor = 0.0;
2013 if (edi->vecs.linacc.stpsz[i] > 0.0)
2015 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
2017 preFactor = edi->vecs.linacc.refproj[i] - proj;
2020 if (edi->vecs.linacc.stpsz[i] < 0.0)
2022 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
2024 preFactor = edi->vecs.linacc.refproj[i] - proj;
2028 /* apply the correction */
2029 preFactor /= edi->sav.sqrtm[i];
2030 for (int j = 0; j < edi->sav.nr; j++)
2032 rvec differenceVector;
2033 svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2034 rvec_inc(xcoll[j], differenceVector);
2037 /* new positions will act as reference */
2038 edi->vecs.linacc.refproj[i] = proj + preFactor;
2044 static void do_radfix(rvec *xcoll, t_edpar *edi)
2047 real *proj, rad = 0.0, ratio;
2051 if (edi->vecs.radfix.neig == 0)
2056 snew(proj, edi->vecs.radfix.neig);
2058 /* loop over radfix vectors */
2059 for (i = 0; i < edi->vecs.radfix.neig; i++)
2061 /* calculate the projections, radius */
2062 proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2063 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2067 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2068 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2070 /* loop over radfix vectors */
2071 for (i = 0; i < edi->vecs.radfix.neig; i++)
2073 proj[i] -= edi->vecs.radfix.refproj[i];
2075 /* apply the correction */
2076 proj[i] /= edi->sav.sqrtm[i];
2078 for (j = 0; j < edi->sav.nr; j++)
2080 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2081 rvec_inc(xcoll[j], vec_dum);
2089 static void do_radacc(rvec *xcoll, t_edpar *edi)
2092 real *proj, rad = 0.0, ratio = 0.0;
2096 if (edi->vecs.radacc.neig == 0)
2101 snew(proj, edi->vecs.radacc.neig);
2103 /* loop over radacc vectors */
2104 for (i = 0; i < edi->vecs.radacc.neig; i++)
2106 /* calculate the projections, radius */
2107 proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2108 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2112 /* only correct when radius decreased */
2113 if (rad < edi->vecs.radacc.radius)
2115 ratio = edi->vecs.radacc.radius/rad - 1.0;
2119 edi->vecs.radacc.radius = rad;
2122 /* loop over radacc vectors */
2123 for (i = 0; i < edi->vecs.radacc.neig; i++)
2125 proj[i] -= edi->vecs.radacc.refproj[i];
2127 /* apply the correction */
2128 proj[i] /= edi->sav.sqrtm[i];
2130 for (j = 0; j < edi->sav.nr; j++)
2132 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2133 rvec_inc(xcoll[j], vec_dum);
2140 struct t_do_radcon {
2144 static void do_radcon(rvec *xcoll, t_edpar *edi)
2147 real rad = 0.0, ratio = 0.0;
2148 struct t_do_radcon *loc;
2153 if (edi->buf->do_radcon != nullptr)
2160 snew(edi->buf->do_radcon, 1);
2162 loc = edi->buf->do_radcon;
2164 if (edi->vecs.radcon.neig == 0)
2171 snew(loc->proj, edi->vecs.radcon.neig);
2174 /* loop over radcon vectors */
2175 for (i = 0; i < edi->vecs.radcon.neig; i++)
2177 /* calculate the projections, radius */
2178 loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2179 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2182 /* only correct when radius increased */
2183 if (rad > edi->vecs.radcon.radius)
2185 ratio = edi->vecs.radcon.radius/rad - 1.0;
2187 /* loop over radcon vectors */
2188 for (i = 0; i < edi->vecs.radcon.neig; i++)
2190 /* apply the correction */
2191 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2192 loc->proj[i] /= edi->sav.sqrtm[i];
2193 loc->proj[i] *= ratio;
2195 for (j = 0; j < edi->sav.nr; j++)
2197 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2198 rvec_inc(xcoll[j], vec_dum);
2205 edi->vecs.radcon.radius = rad;
2211 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, int64_t step)
2216 /* subtract the average positions */
2217 for (i = 0; i < edi->sav.nr; i++)
2219 rvec_dec(xcoll[i], edi->sav.x[i]);
2222 /* apply the constraints */
2225 do_linfix(xcoll, *edi, step);
2227 do_linacc(xcoll, edi);
2230 do_radfix(xcoll, edi);
2232 do_radacc(xcoll, edi);
2233 do_radcon(xcoll, edi);
2235 /* add back the average positions */
2236 for (i = 0; i < edi->sav.nr; i++)
2238 rvec_inc(xcoll[i], edi->sav.x[i]);
2245 /*!\brief Write out the projections onto the eigenvectors.
2246 * The order of output corresponds to ed_output_legend().
2247 * \param[in] edi The essential dyanmics parameters
2248 * \param[in] fp The output file
2249 * \param[in] rmsd the rmsd to the reference structure
2251 void write_edo(const t_edpar &edi, FILE *fp, real rmsd)
2253 /* Output how well we fit to the reference structure */
2254 fprintf(fp, EDcol_ffmt, rmsd);
2256 for (int i = 0; i < edi.vecs.mon.neig; i++)
2258 fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2261 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2263 fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2266 for (int i = 0; i < edi.vecs.linacc.neig; i++)
2268 fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2271 for (int i = 0; i < edi.vecs.radfix.neig; i++)
2273 fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2275 if (edi.vecs.radfix.neig)
2277 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2280 for (int i = 0; i < edi.vecs.radacc.neig; i++)
2282 fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2284 if (edi.vecs.radacc.neig)
2286 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2289 for (int i = 0; i < edi.vecs.radcon.neig; i++)
2291 fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2293 if (edi.vecs.radcon.neig)
2295 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2301 /* Returns if any constraints are switched on */
2302 static bool ed_constraints(EssentialDynamicsType edtype, const t_edpar &edi)
2304 if (edtype == EssentialDynamicsType::EDSampling || edtype == EssentialDynamicsType::Flooding)
2306 return ((edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0) ||
2307 (edi.vecs.radfix.neig != 0) || (edi.vecs.radacc.neig != 0) ||
2308 (edi.vecs.radcon.neig != 0));
2314 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for flooding/
2315 * umbrella sampling simulations. */
2316 static void copyEvecReference(t_eigvec* floodvecs, real * initialReferenceProjection)
2321 if (nullptr == initialReferenceProjection)
2323 snew(initialReferenceProjection, floodvecs->neig);
2326 for (i = 0; i < floodvecs->neig; i++)
2328 initialReferenceProjection[i] = floodvecs->refproj[i];
2333 /* Call on MASTER only. Check whether the essential dynamics / flooding
2334 * groups of the checkpoint file are consistent with the provided .edi file. */
2335 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam &ed, edsamhistory_t *EDstate)
2337 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2339 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2340 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2341 "it must also continue with/without ED constraints when checkpointing.\n"
2342 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2343 "from without a checkpoint.\n");
2346 for (size_t edinum = 0; edinum < ed.edpar.size(); ++edinum)
2348 /* Check number of atoms in the reference and average structures */
2349 if (EDstate->nref[edinum] != ed.edpar[edinum].sref.nr)
2351 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2352 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2353 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], ed.edpar[edinum].sref.nr);
2355 if (EDstate->nav[edinum] != ed.edpar[edinum].sav.nr)
2357 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2358 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2359 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], ed.edpar[edinum].sav.nr);
2363 if (static_cast<int>(ed.edpar.size()) != EDstate->nED)
2365 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2366 "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2367 "Are you sure this is the correct .edi file?\n", EDstate->nED, ed.edpar.size());
2372 /* The edsamstate struct stores the information we need to make the ED group
2373 * whole again after restarts from a checkpoint file. Here we do the following:
2374 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2375 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2376 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2377 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2378 * all ED structures at the last time step. */
2379 static void init_edsamstate(const gmx_edsam &ed, edsamhistory_t *EDstate)
2381 snew(EDstate->old_sref_p, EDstate->nED);
2382 snew(EDstate->old_sav_p, EDstate->nED);
2384 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2385 if (!EDstate->bFromCpt)
2387 snew(EDstate->nref, EDstate->nED);
2388 snew(EDstate->nav, EDstate->nED);
2391 /* Loop over all ED/flooding data sets (usually only one, though) */
2392 for (int nr_edi = 0; nr_edi < EDstate->nED; nr_edi++)
2394 const auto &edi = ed.edpar[nr_edi];
2395 /* We always need the last reference and average positions such that
2396 * in the next time step we can make the ED group whole again
2397 * if the atoms do not have the correct PBC representation */
2398 if (EDstate->bFromCpt)
2400 /* Copy the last whole positions of reference and average group from .cpt */
2401 for (int i = 0; i < edi.sref.nr; i++)
2403 copy_rvec(EDstate->old_sref[nr_edi][i], edi.sref.x_old[i]);
2405 for (int i = 0; i < edi.sav.nr; i++)
2407 copy_rvec(EDstate->old_sav [nr_edi][i], edi.sav.x_old [i]);
2412 EDstate->nref[nr_edi] = edi.sref.nr;
2413 EDstate->nav [nr_edi] = edi.sav.nr;
2416 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2417 EDstate->old_sref_p[nr_edi] = edi.sref.x_old;
2418 EDstate->old_sav_p [nr_edi] = edi.sav.x_old;
2423 /* Adds 'buf' to 'str' */
2424 static void add_to_string(char **str, char *buf)
2429 len = strlen(*str) + strlen(buf) + 1;
2435 static void add_to_string_aligned(char **str, char *buf)
2437 char buf_aligned[STRLEN];
2439 sprintf(buf_aligned, EDcol_sfmt, buf);
2440 add_to_string(str, buf_aligned);
2444 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, const char *value, const char *unit, char EDgroupchar)
2446 char tmp[STRLEN], tmp2[STRLEN];
2449 sprintf(tmp, "%c %s", EDgroupchar, value);
2450 add_to_string_aligned(LegendStr, tmp);
2451 sprintf(tmp2, "%s (%s)", tmp, unit);
2452 (*setname)[*nsets] = gmx_strdup(tmp2);
2457 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2463 for (i = 0; i < evec->neig; i++)
2465 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2466 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2471 /* Makes a legend for the xvg output file. Call on MASTER only! */
2472 static void write_edo_legend(gmx_edsam * ed, int nED, const gmx_output_env_t *oenv)
2475 int nr_edi, nsets, n_flood, n_edsam;
2476 const char **setname;
2478 char *LegendStr = nullptr;
2481 auto edi = ed->edpar.begin();
2483 fprintf(ed->edo, "# Output will be written every %d step%s\n", edi->outfrq, edi->outfrq != 1 ? "s" : "");
2485 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2487 fprintf(ed->edo, "#\n");
2488 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2489 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2490 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2491 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2492 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2493 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2494 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2495 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2496 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2498 if (edi->flood.vecs.neig)
2500 /* If in any of the groups we find a flooding vector, flooding is turned on */
2501 ed->eEDtype = EssentialDynamicsType::Flooding;
2503 /* Print what flavor of flooding we will do */
2504 if (0 == edi->flood.tau) /* constant flooding strength */
2506 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2507 if (edi->flood.bHarmonic)
2509 fprintf(ed->edo, ", harmonic");
2512 else /* adaptive flooding */
2514 fprintf(ed->edo, ", adaptive");
2517 fprintf(ed->edo, "\n");
2522 /* Print a nice legend */
2524 LegendStr[0] = '\0';
2525 sprintf(buf, "# %6s", "time");
2526 add_to_string(&LegendStr, buf);
2528 /* Calculate the maximum number of columns we could end up with */
2529 edi = ed->edpar.begin();
2531 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2533 nsets += 5 +edi->vecs.mon.neig
2534 +edi->vecs.linfix.neig
2535 +edi->vecs.linacc.neig
2536 +edi->vecs.radfix.neig
2537 +edi->vecs.radacc.neig
2538 +edi->vecs.radcon.neig
2539 + 6*edi->flood.vecs.neig;
2542 snew(setname, nsets);
2544 /* In the mdrun time step in a first function call (do_flood()) the flooding
2545 * forces are calculated and in a second function call (do_edsam()) the
2546 * ED constraints. To get a corresponding legend, we need to loop twice
2547 * over the edi groups and output first the flooding, then the ED part */
2549 /* The flooding-related legend entries, if flooding is done */
2551 if (EssentialDynamicsType::Flooding == ed->eEDtype)
2553 edi = ed->edpar.begin();
2554 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2556 /* Always write out the projection on the flooding EVs. Of course, this can also
2557 * be achieved with the monitoring option in do_edsam() (if switched on by the
2558 * user), but in that case the positions need to be communicated in do_edsam(),
2559 * which is not necessary when doing flooding only. */
2560 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2562 for (i = 0; i < edi->flood.vecs.neig; i++)
2564 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2565 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2567 /* Output the current reference projection if it changes with time;
2568 * this can happen when flooding is used as harmonic restraint */
2569 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2571 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2572 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2575 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2576 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2578 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2579 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2582 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2583 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2585 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2587 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2588 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2591 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2592 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2596 } /* End of flooding-related legend entries */
2600 /* Now the ED-related entries, if essential dynamics is done */
2601 edi = ed->edpar.begin();
2602 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2604 if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2606 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2608 /* Essential dynamics, projections on eigenvectors */
2609 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2610 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2611 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2612 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2613 if (edi->vecs.radfix.neig)
2615 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2617 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2618 if (edi->vecs.radacc.neig)
2620 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2622 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2623 if (edi->vecs.radcon.neig)
2625 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2629 } /* end of 'pure' essential dynamics legend entries */
2630 n_edsam = nsets - n_flood;
2632 xvgr_legend(ed->edo, nsets, setname, oenv);
2635 fprintf(ed->edo, "#\n"
2636 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2637 n_flood, 1 == n_flood ? "" : "s",
2638 n_edsam, 1 == n_edsam ? "" : "s");
2639 fprintf(ed->edo, "%s", LegendStr);
2646 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2647 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2648 std::unique_ptr<gmx::EssentialDynamics> init_edsam(
2649 const char *ediFileName,
2650 const char *edoFileName,
2651 const gmx_mtop_t *mtop,
2652 const t_inputrec *ir,
2653 const t_commrec *cr,
2654 gmx::Constraints *constr,
2655 const t_state *globalState,
2656 ObservablesHistory *oh,
2657 const gmx_output_env_t *oenv,
2661 rvec *x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2662 rvec *xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2663 rvec fit_transvec; /* translation ... */
2664 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2665 rvec *ref_x_old = nullptr; /* helper pointer */
2670 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2672 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName) )
2674 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or flooding\n"
2675 "simulation. Please also set the .edi file on the command line with -ei.\n");
2679 /* Open input and output files, allocate space for ED data structure */
2680 auto edHandle = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr);
2681 auto ed = edHandle->getLegacyED();
2682 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2683 constr->saveEdsamPointer(ed);
2685 /* Needed for initializing radacc radius in do_edsam */
2688 /* The input file is read by the master and the edi structures are
2689 * initialized here. Input is stored in ed->edpar. Then the edi
2690 * structures are transferred to the other nodes */
2693 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2694 * flooding vector, Essential dynamics can be applied to more than one structure
2695 * as well, but will be done in the order given in the edi file, so
2696 * expect different results for different order of edi file concatenation! */
2697 for (auto &edi : ed->edpar)
2699 init_edi(mtop, &edi);
2700 init_flood(&edi, ed, ir->delta_t);
2704 /* The master does the work here. The other nodes get the positions
2705 * not before dd_partition_system which is called after init_edsam */
2708 edsamhistory_t *EDstate = oh->edsamHistory.get();
2710 if (!EDstate->bFromCpt)
2712 /* Remove PBC, make molecule(s) subject to ED whole. */
2713 snew(x_pbc, mtop->natoms);
2714 copy_rvecn(as_rvec_array(globalState->x.data()), x_pbc, 0, mtop->natoms);
2715 do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
2717 /* Reset pointer to first ED data set which contains the actual ED data */
2718 auto edi = ed->edpar.begin();
2719 GMX_ASSERT(!ed->edpar.empty(), "Essential dynamics parameters is undefined");
2721 /* Loop over all ED/flooding data sets (usually only one, though) */
2722 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2724 /* For multiple ED groups we use the output frequency that was specified
2725 * in the first set */
2728 edi->outfrq = ed->edpar.begin()->outfrq;
2731 /* Extract the initial reference and average positions. When starting
2732 * from .cpt, these have already been read into sref.x_old
2733 * in init_edsamstate() */
2734 if (!EDstate->bFromCpt)
2736 /* If this is the first run (i.e. no checkpoint present) we assume
2737 * that the starting positions give us the correct PBC representation */
2738 for (i = 0; i < edi->sref.nr; i++)
2740 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2743 for (i = 0; i < edi->sav.nr; i++)
2745 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2749 /* Now we have the PBC-correct start positions of the reference and
2750 average structure. We copy that over to dummy arrays on which we
2751 can apply fitting to print out the RMSD. We srenew the memory since
2752 the size of the buffers is likely different for every ED group */
2753 srenew(xfit, edi->sref.nr );
2754 srenew(xstart, edi->sav.nr );
2757 /* Reference indices are the same as average indices */
2758 ref_x_old = edi->sav.x_old;
2762 ref_x_old = edi->sref.x_old;
2764 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2765 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2767 /* Make the fit to the REFERENCE structure, get translation and rotation */
2768 fit_to_reference(xfit, fit_transvec, fit_rotmat, &(*edi));
2770 /* Output how well we fit to the reference at the start */
2771 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2772 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2773 rmsd_from_structure(xfit, &edi->sref));
2774 if (EDstate->nED > 1)
2776 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2778 fprintf(stderr, "\n");
2780 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2781 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2783 /* calculate initial projections */
2784 project(xstart, &(*edi));
2786 /* For the target and origin structure both a reference (fit) and an
2787 * average structure can be provided in make_edi. If both structures
2788 * are the same, make_edi only stores one of them in the .edi file.
2789 * If they differ, first the fit and then the average structure is stored
2790 * in star (or sor), thus the number of entries in star/sor is
2791 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2792 * the size of the average group. */
2794 /* process target structure, if required */
2795 if (edi->star.nr > 0)
2797 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2799 /* get translation & rotation for fit of target structure to reference structure */
2800 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, &(*edi));
2802 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2803 if (edi->star.nr == edi->sav.nr)
2807 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2809 /* The last sav.nr indices of the target structure correspond to
2810 * the average structure, which must be projected */
2811 avindex = edi->star.nr - edi->sav.nr;
2813 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2817 rad_project(*edi, xstart, &edi->vecs.radcon);
2820 /* process structure that will serve as origin of expansion circle */
2821 if ( (EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce) )
2823 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2826 if (edi->sori.nr > 0)
2828 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2830 /* fit this structure to reference structure */
2831 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, &(*edi));
2833 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2834 if (edi->sori.nr == edi->sav.nr)
2838 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2840 /* For the projection, we need the last sav.nr indices of sori */
2841 avindex = edi->sori.nr - edi->sav.nr;
2844 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2845 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2846 if ( (EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce) )
2848 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2849 /* Set center of flooding potential to the ORIGIN structure */
2850 rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2851 /* We already know that no (moving) reference position was provided,
2852 * therefore we can overwrite refproj[0]*/
2853 copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
2856 else /* No origin structure given */
2858 rad_project(*edi, xstart, &edi->vecs.radacc);
2859 rad_project(*edi, xstart, &edi->vecs.radfix);
2860 if ( (EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce) )
2862 if (edi->flood.bHarmonic)
2864 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2865 for (i = 0; i < edi->flood.vecs.neig; i++)
2867 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
2872 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2873 /* Set center of flooding potential to the center of the covariance matrix,
2874 * i.e. the average structure, i.e. zero in the projected system */
2875 for (i = 0; i < edi->flood.vecs.neig; i++)
2877 edi->flood.vecs.refproj[i] = 0.0;
2882 /* For convenience, output the center of the flooding potential for the eigenvectors */
2883 if ( (EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce) )
2885 for (i = 0; i < edi->flood.vecs.neig; i++)
2887 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2888 if (edi->flood.bHarmonic)
2890 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.referenceProjectionSlope[i]);
2892 fprintf(stdout, "\n");
2896 /* set starting projections for linsam */
2897 rad_project(*edi, xstart, &edi->vecs.linacc);
2898 rad_project(*edi, xstart, &edi->vecs.linfix);
2900 /* Prepare for the next edi data set: */
2903 /* Cleaning up on the master node: */
2904 if (!EDstate->bFromCpt)
2911 } /* end of MASTER only section */
2915 /* Broadcast the essential dynamics / flooding data to all nodes */
2916 broadcast_ed_data(cr, ed);
2920 /* In the single-CPU case, point the local atom numbers pointers to the global
2921 * one, so that we can use the same notation in serial and parallel case: */
2922 /* Loop over all ED data sets (usually only one, though) */
2923 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
2925 edi->sref.anrs_loc = edi->sref.anrs;
2926 edi->sav.anrs_loc = edi->sav.anrs;
2927 edi->star.anrs_loc = edi->star.anrs;
2928 edi->sori.anrs_loc = edi->sori.anrs;
2929 /* For the same reason as above, make a dummy c_ind array: */
2930 snew(edi->sav.c_ind, edi->sav.nr);
2931 /* Initialize the array */
2932 for (i = 0; i < edi->sav.nr; i++)
2934 edi->sav.c_ind[i] = i;
2936 /* In the general case we will need a different-sized array for the reference indices: */
2939 snew(edi->sref.c_ind, edi->sref.nr);
2940 for (i = 0; i < edi->sref.nr; i++)
2942 edi->sref.c_ind[i] = i;
2945 /* Point to the very same array in case of other structures: */
2946 edi->star.c_ind = edi->sav.c_ind;
2947 edi->sori.c_ind = edi->sav.c_ind;
2948 /* In the serial case, the local number of atoms is the global one: */
2949 edi->sref.nr_loc = edi->sref.nr;
2950 edi->sav.nr_loc = edi->sav.nr;
2951 edi->star.nr_loc = edi->star.nr;
2952 edi->sori.nr_loc = edi->sori.nr;
2956 /* Allocate space for ED buffer variables */
2957 /* Again, loop over ED data sets */
2958 for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
2960 /* Allocate space for ED buffer variables */
2961 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
2962 snew(edi->buf->do_edsam, 1);
2964 /* Space for collective ED buffer variables */
2966 /* Collective positions of atoms with the average indices */
2967 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2968 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
2969 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2970 /* Collective positions of atoms with the reference indices */
2973 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2974 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
2975 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2978 /* Get memory for flooding forces */
2979 snew(edi->flood.forces_cartesian, edi->sav.nr);
2982 /* Flush the edo file so that the user can check some things
2983 * when the simulation has started */
2993 void do_edsam(const t_inputrec *ir,
2995 const t_commrec *cr,
3001 int i, edinr, iupdate = 500;
3002 matrix rotmat; /* rotation matrix */
3003 rvec transvec; /* translation vector */
3004 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3005 real dt_1; /* 1/dt */
3006 struct t_do_edsam *buf;
3007 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3008 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3011 /* Check if ED sampling has to be performed */
3012 if (ed->eEDtype == EssentialDynamicsType::None)
3017 dt_1 = 1.0/ir->delta_t;
3019 /* Loop over all ED groups (usually one) */
3021 for (auto &edi : ed->edpar)
3024 if (bNeedDoEdsam(edi))
3027 buf = edi.buf->do_edsam;
3031 /* initialize radacc radius for slope criterion */
3032 buf->oldrad = calc_radius(edi.vecs.radacc);
3035 /* Copy the positions into buf->xc* arrays and after ED
3036 * feed back corrections to the official positions */
3038 /* Broadcast the ED positions such that every node has all of them
3039 * Every node contributes its local positions xs and stores it in
3040 * the collective buf->xcoll array. Note that for edinr > 1
3041 * xs could already have been modified by an earlier ED */
3043 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3044 edi.sav.nr, edi.sav.nr_loc, edi.sav.anrs_loc, edi.sav.c_ind, edi.sav.x_old, box);
3046 /* Only assembly reference positions if their indices differ from the average ones */
3049 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3050 edi.sref.nr, edi.sref.nr_loc, edi.sref.anrs_loc, edi.sref.c_ind, edi.sref.x_old, box);
3053 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3054 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3055 * set bUpdateShifts=TRUE in the parallel case. */
3056 buf->bUpdateShifts = FALSE;
3058 /* Now all nodes have all of the ED positions in edi.sav->xcoll,
3059 * as well as the indices in edi.sav.anrs */
3061 /* Fit the reference indices to the reference structure */
3064 fit_to_reference(buf->xcoll, transvec, rotmat, &edi);
3068 fit_to_reference(buf->xc_ref, transvec, rotmat, &edi);
3071 /* Now apply the translation and rotation to the ED structure */
3072 translate_and_rotate(buf->xcoll, edi.sav.nr, transvec, rotmat);
3074 /* Find out how well we fit to the reference (just for output steps) */
3075 if (do_per_step(step, edi.outfrq) && MASTER(cr))
3079 /* Indices of reference and average structures are identical,
3080 * thus we can calculate the rmsd to SREF using xcoll */
3081 rmsdev = rmsd_from_structure(buf->xcoll, &edi.sref);
3085 /* We have to translate & rotate the reference atoms first */
3086 translate_and_rotate(buf->xc_ref, edi.sref.nr, transvec, rotmat);
3087 rmsdev = rmsd_from_structure(buf->xc_ref, &edi.sref);
3091 /* update radsam references, when required */
3092 if (do_per_step(step, edi.maxedsteps) && step >= edi.presteps)
3094 project(buf->xcoll, &edi);
3095 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3096 rad_project(edi, buf->xcoll, &edi.vecs.radfix);
3097 buf->oldrad = -1.e5;
3100 /* update radacc references, when required */
3101 if (do_per_step(step, iupdate) && step >= edi.presteps)
3103 edi.vecs.radacc.radius = calc_radius(edi.vecs.radacc);
3104 if (edi.vecs.radacc.radius - buf->oldrad < edi.slope)
3106 project(buf->xcoll, &edi);
3107 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3112 buf->oldrad = edi.vecs.radacc.radius;
3116 /* apply the constraints */
3117 if (step >= edi.presteps && ed_constraints(ed->eEDtype, edi))
3119 /* ED constraints should be applied already in the first MD step
3120 * (which is step 0), therefore we pass step+1 to the routine */
3121 ed_apply_constraints(buf->xcoll, &edi, step+1 - ir->init_step);
3124 /* write to edo, when required */
3125 if (do_per_step(step, edi.outfrq))
3127 project(buf->xcoll, &edi);
3128 if (MASTER(cr) && !bSuppress)
3130 write_edo(edi, ed->edo, rmsdev);
3134 /* Copy back the positions unless monitoring only */
3135 if (ed_constraints(ed->eEDtype, edi))
3137 /* remove fitting */
3138 rmfit(edi.sav.nr, buf->xcoll, transvec, rotmat);
3140 /* Copy the ED corrected positions into the coordinate array */
3141 /* Each node copies its local part. In the serial case, nat_loc is the
3142 * total number of ED atoms */
3143 for (i = 0; i < edi.sav.nr_loc; i++)
3145 /* Unshift local ED coordinate and store in x_unsh */
3146 ed_unshift_single_coord(box, buf->xcoll[edi.sav.c_ind[i]],
3147 buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh);
3149 /* dx is the ED correction to the positions: */
3150 rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx);
3154 /* dv is the ED correction to the velocity: */
3155 svmul(dt_1, dx, dv);
3156 /* apply the velocity correction: */
3157 rvec_inc(v[edi.sav.anrs_loc[i]], dv);
3159 /* Finally apply the position correction due to ED: */
3160 copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]);
3163 } /* END of if ( bNeedDoEdsam(edi) ) */
3165 /* Prepare for the next ED group */
3167 } /* END of loop over ED groups */