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"
79 /* enum to identify the type of ED: none, normal ED, flooding */
81 eEDnone, eEDedsam, eEDflood, eEDnr
84 /* enum to identify operations on reference, average, origin, target structures */
86 eedREF, eedAV, eedORI, eedTAR, eedNR
92 int neig; /* nr of eigenvectors */
93 int *ieig; /* index nrs of eigenvectors */
94 real *stpsz; /* stepsizes (per eigenvector) */
95 rvec **vec; /* eigenvector components */
96 real *xproj; /* instantaneous x projections */
97 real *fproj; /* instantaneous f projections */
98 real radius; /* instantaneous radius */
99 real *refproj; /* starting or target projections */
105 t_eigvec mon; /* only monitored, no constraints */
106 t_eigvec linfix; /* fixed linear constraints */
107 t_eigvec linacc; /* acceptance linear constraints */
108 t_eigvec radfix; /* fixed radial constraints (exp) */
109 t_eigvec radacc; /* acceptance radial constraints (exp) */
110 t_eigvec radcon; /* acceptance rad. contraction constr. */
117 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
119 gmx_bool bConstForce; /* Do not calculate a flooding potential,
120 instead flood with a constant force */
129 rvec *forces_cartesian;
130 t_eigvec vecs; /* use flooding for these */
131 /* When using flooding as harmonic restraint: The current reference projection
132 * is at each step calculated from the initial initialReferenceProjection and the slope. */
133 real *initialReferenceProjection;
134 real *referenceProjectionSlope;
138 /* This type is for the average, reference, target, and origin structure */
139 typedef struct gmx_edx
141 int nr; /* number of atoms this structure contains */
142 int nr_loc; /* number of atoms on local node */
143 int *anrs; /* atom index numbers */
144 int *anrs_loc; /* local atom index numbers */
145 int nalloc_loc; /* allocation size of anrs_loc */
146 int *c_ind; /* at which position of the whole anrs
147 * array is a local atom?, i.e.
148 * c_ind[0...nr_loc-1] gives the atom index
149 * with respect to the collective
150 * anrs[0...nr-1] array */
151 rvec *x; /* positions for this structure */
152 rvec *x_old; /* Last positions which have the correct PBC
153 representation of the ED group. In
154 combination with keeping track of the
155 shift vectors, the ED group can always
157 real *m; /* masses */
158 real mtot; /* total mass (only used in sref) */
159 real *sqrtm; /* sqrt of the masses used for mass-
160 * weighting of analysis (only used in sav) */
166 int nini; /* total Nr of atoms */
167 gmx_bool fitmas; /* true if trans fit with cm */
168 gmx_bool pcamas; /* true if mass-weighted PCA */
169 int presteps; /* number of steps to run without any
170 * perturbations ... just monitoring */
171 int outfrq; /* freq (in steps) of writing to edo */
172 int maxedsteps; /* max nr of steps per cycle */
174 /* all gmx_edx datasets are copied to all nodes in the parallel case */
175 struct gmx_edx sref; /* reference positions, to these fitting
177 gmx_bool bRefEqAv; /* If true, reference & average indices
178 * are the same. Used for optimization */
179 struct gmx_edx sav; /* average positions */
180 struct gmx_edx star; /* target positions */
181 struct gmx_edx sori; /* origin positions */
183 t_edvecs vecs; /* eigenvectors */
184 real slope; /* minimal slope in acceptance radexp */
186 t_edflood flood; /* parameters especially for flooding */
187 struct t_ed_buffer *buf; /* handle to local buffers */
188 struct edpar *next_edi; /* Pointer to another ED group */
195 int eEDtype = eEDnone; /* Type of ED: see enums above */
196 FILE *edo = nullptr; /* output file pointer */
197 t_edpar *edpar = nullptr;
198 gmx_bool bFirst = false;
200 gmx_edsam::~gmx_edsam()
204 /* edo is opened sometimes with xvgropen, sometimes with
205 * gmx_fio_fopen, so we use the least common denominator for
215 rvec old_transvec, older_transvec, transvec_compact;
216 rvec *xcoll; /* Positions from all nodes, this is the
217 collective set we work on.
218 These are the positions of atoms with
219 average structure indices */
220 rvec *xc_ref; /* same but with reference structure indices */
221 ivec *shifts_xcoll; /* Shifts for xcoll */
222 ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
223 ivec *shifts_xc_ref; /* Shifts for xc_ref */
224 ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
225 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
226 ED shifts for this ED group need to
231 /* definition of ED buffer structure */
234 struct t_fit_to_ref * fit_to_ref;
235 struct t_do_edfit * do_edfit;
236 struct t_do_edsam * do_edsam;
237 struct t_do_radcon * do_radcon;
242 class EssentialDynamics::Impl
245 // TODO: move all fields from gmx_edsam here and remove gmx_edsam
246 gmx_edsam essentialDynamics_;
248 EssentialDynamics::EssentialDynamics() : impl_(new Impl)
251 EssentialDynamics::~EssentialDynamics() = default;
253 gmx_edsam *EssentialDynamics::getLegacyED()
255 return &impl_->essentialDynamics_;
259 /* Function declarations */
260 static void fit_to_reference(rvec *xcoll, rvec transvec, matrix rotmat, t_edpar *edi);
261 static void translate_and_rotate(rvec *x, int nat, rvec transvec, matrix rotmat);
262 static real rmsd_from_structure(rvec *x, struct gmx_edx *s);
263 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms);
264 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam &ed, edsamhistory_t *EDstate);
265 static void init_edsamstate(gmx_edsam * ed, edsamhistory_t *EDstate);
266 static void write_edo_legend(gmx_edsam * ed, int nED, const gmx_output_env_t *oenv);
267 /* End function declarations */
269 /* Define formats for the column width in the output file */
270 const char EDcol_sfmt[] = "%17s";
271 const char EDcol_efmt[] = "%17.5e";
272 const char EDcol_ffmt[] = "%17f";
274 /* Define a formats for reading, otherwise cppcheck complains for scanf without width limits */
275 const char max_ev_fmt_d[] = "%7d"; /* Number of eigenvectors. 9,999,999 should be enough */
276 const char max_ev_fmt_lf[] = "%12lf";
277 const char max_ev_fmt_dlf[] = "%7d%12lf";
278 const char max_ev_fmt_dlflflf[] = "%7d%12lf%12lf%12lf";
279 const char max_ev_fmt_lelele[] = "%12le%12le%12le";
283 /*!\brief Determines whether to perform essential dynamics constraints other than flooding.
284 * \param[in] edi the essential dynamics parameters
285 * \returns true if essential dyanmics constraints need to be performed
287 bool bNeedDoEdsam(const t_edpar &edi)
289 return edi.vecs.mon.neig
290 || edi.vecs.linfix.neig
291 || edi.vecs.linacc.neig
292 || edi.vecs.radfix.neig
293 || edi.vecs.radacc.neig
294 || edi.vecs.radcon.neig;
299 /* Multiple ED groups will be labeled with letters instead of numbers
300 * to avoid confusion with eigenvector indices */
301 static char get_EDgroupChar(int nr_edi, int nED)
309 * nr_edi = 2 -> B ...
311 return 'A' + nr_edi - 1;
316 /*! \brief The mass-weighted inner product of two coordinate vectors.
317 * Does not subtract average positions, projection on single eigenvector is returned
318 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
319 * Average position is subtracted in ed_apply_constraints prior to calling projectx
320 * \param[in] edi Essential dynamics parameters
321 * \param[in] xcoll vector of atom coordinates
322 * \param[in] vec vector of coordinates to project onto
323 * \return mass-weighted projection.
325 real projectx(const t_edpar &edi, rvec *xcoll, rvec *vec)
331 for (i = 0; i < edi.sav.nr; i++)
333 proj += edi.sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
338 /*!\brief Project coordinates onto vector after substracting average position.
339 * projection is stored in vec->refproj which is used for radacc, radfix,
340 * radcon and center of flooding potential.
341 * Average positions are first substracted from x, then added back again.
342 * \param[in] edi essential dynamics parameters with average position
343 * \param[in] x Coordinates to be projected
344 * \param[out] vec eigenvector, radius and refproj are overwritten here
346 void rad_project(const t_edpar &edi, rvec *x, t_eigvec *vec)
351 /* Subtract average positions */
352 for (i = 0; i < edi.sav.nr; i++)
354 rvec_dec(x[i], edi.sav.x[i]);
357 for (i = 0; i < vec->neig; i++)
359 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
360 rad += gmx::square((vec->refproj[i]-vec->xproj[i]));
362 vec->radius = sqrt(rad);
364 /* Add average positions */
365 for (i = 0; i < edi.sav.nr; i++)
367 rvec_inc(x[i], edi.sav.x[i]);
371 /*!\brief Projects coordinates onto eigenvectors and stores result in vec->xproj.
372 * Mass-weighting is applied. Subtracts average positions prior to projection and add
373 * them afterwards to retain the unchanged vector.
374 * \param[in] x The coordinates to project to an eigenvector
375 * \param[in,out] vec The eigenvectors
376 * \param[in] edi essential dynamics parameters holding average structure and masses
378 void project_to_eigvectors(rvec *x, t_eigvec *vec, const t_edpar &edi)
385 /* Subtract average positions */
386 for (int i = 0; i < edi.sav.nr; i++)
388 rvec_dec(x[i], edi.sav.x[i]);
391 for (int i = 0; i < vec->neig; i++)
393 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
396 /* Add average positions */
397 for (int i = 0; i < edi.sav.nr; i++)
399 rvec_inc(x[i], edi.sav.x[i]);
404 /* Project vector x onto all edi->vecs (mon, linfix,...) */
405 static void project(rvec *x, /* positions to project */
406 t_edpar *edi) /* edi data set */
408 /* It is not more work to subtract the average position in every
409 * subroutine again, because these routines are rarely used simultaneously */
410 project_to_eigvectors(x, &edi->vecs.mon, *edi);
411 project_to_eigvectors(x, &edi->vecs.linfix, *edi);
412 project_to_eigvectors(x, &edi->vecs.linacc, *edi);
413 project_to_eigvectors(x, &edi->vecs.radfix, *edi);
414 project_to_eigvectors(x, &edi->vecs.radacc, *edi);
415 project_to_eigvectors(x, &edi->vecs.radcon, *edi);
420 /*!\brief Evaluates the distance from reference to current eigenvector projection.
421 * \param[in] vec eigenvector
424 real calc_radius(const t_eigvec &vec)
428 for (int i = 0; i < vec.neig; i++)
430 rad += gmx::square((vec.refproj[i]-vec.xproj[i]));
433 return rad = sqrt(rad);
442 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
444 /* this is a copy of do_fit with some modifications */
445 int c, r, n, j, i, irot;
446 double d[6], xnr, xpc;
451 struct t_do_edfit *loc;
454 if (edi->buf->do_edfit != nullptr)
461 snew(edi->buf->do_edfit, 1);
463 loc = edi->buf->do_edfit;
467 snew(loc->omega, 2*DIM);
468 snew(loc->om, 2*DIM);
469 for (i = 0; i < 2*DIM; i++)
471 snew(loc->omega[i], 2*DIM);
472 snew(loc->om[i], 2*DIM);
476 for (i = 0; (i < 6); i++)
479 for (j = 0; (j < 6); j++)
481 loc->omega[i][j] = 0;
486 /* calculate the matrix U */
488 for (n = 0; (n < natoms); n++)
490 for (c = 0; (c < DIM); c++)
493 for (r = 0; (r < DIM); r++)
501 /* construct loc->omega */
502 /* loc->omega is symmetric -> loc->omega==loc->omega' */
503 for (r = 0; (r < 6); r++)
505 for (c = 0; (c <= r); c++)
507 if ((r >= 3) && (c < 3))
509 loc->omega[r][c] = u[r-3][c];
510 loc->omega[c][r] = u[r-3][c];
514 loc->omega[r][c] = 0;
515 loc->omega[c][r] = 0;
520 /* determine h and k */
521 jacobi(loc->omega, 6, d, loc->om, &irot);
525 fprintf(stderr, "IROT=0\n");
528 index = 0; /* For the compiler only */
530 for (j = 0; (j < 3); j++)
533 for (i = 0; (i < 6); i++)
542 for (i = 0; (i < 3); i++)
544 vh[j][i] = M_SQRT2*loc->om[i][index];
545 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
550 for (c = 0; (c < 3); c++)
552 for (r = 0; (r < 3); r++)
554 R[c][r] = vk[0][r]*vh[0][c]+
561 for (c = 0; (c < 3); c++)
563 for (r = 0; (r < 3); r++)
565 R[c][r] = vk[0][r]*vh[0][c]+
574 static void rmfit(int nat, rvec *xcoll, const rvec transvec, matrix rotmat)
581 * The inverse rotation is described by the transposed rotation matrix */
582 transpose(rotmat, tmat);
583 rotate_x(xcoll, nat, tmat);
585 /* Remove translation */
586 vec[XX] = -transvec[XX];
587 vec[YY] = -transvec[YY];
588 vec[ZZ] = -transvec[ZZ];
589 translate_x(xcoll, nat, vec);
593 /**********************************************************************************
594 ******************** FLOODING ****************************************************
595 **********************************************************************************
597 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
598 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
599 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
601 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
602 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
604 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
605 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
607 do_flood makes a copy of the positions,
608 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
609 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
610 fit. Then do_flood adds these forces to the forcefield-forces
611 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
613 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
614 structure is projected to the system of eigenvectors and then this position in the subspace is used as
615 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
616 i.e. the average structure as given in the make_edi file.
618 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
619 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
620 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
621 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
622 further adaption is applied, Efl will stay constant at zero.
624 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
625 used as spring constants for the harmonic potential.
626 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
627 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
629 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
630 the routine read_edi_file reads all of theses flooding files.
631 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
632 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
633 edi there is no interdependence whatsoever. The forces are added together.
635 To write energies into the .edr file, call the function
636 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
638 get_flood_energies(real Vfl[],int nnames);
641 - 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.
643 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
644 two edsam files from two peptide chains
647 // TODO split this into multiple files
650 /*!\brief Output flooding simulation settings and results to file.
651 * \param[in] edi Essential dynamics input parameters
652 * \param[in] fp output file
653 * \param[in] rmsd rmsd to reference structure
655 void write_edo_flood(const t_edpar &edi, FILE *fp, real rmsd)
657 /* Output how well we fit to the reference structure */
658 fprintf(fp, EDcol_ffmt, rmsd);
660 for (int i = 0; i < edi.flood.vecs.neig; i++)
662 fprintf(fp, EDcol_efmt, edi.flood.vecs.xproj[i]);
664 /* Check whether the reference projection changes with time (this can happen
665 * in case flooding is used as harmonic restraint). If so, output the
666 * current reference projection */
667 if (edi.flood.bHarmonic && edi.flood.referenceProjectionSlope[i] != 0.0)
669 fprintf(fp, EDcol_efmt, edi.flood.vecs.refproj[i]);
672 /* Output Efl if we are doing adaptive flooding */
673 if (0 != edi.flood.tau)
675 fprintf(fp, EDcol_efmt, edi.flood.Efl);
677 fprintf(fp, EDcol_efmt, edi.flood.Vfl);
679 /* Output deltaF if we are doing adaptive flooding */
680 if (0 != edi.flood.tau)
682 fprintf(fp, EDcol_efmt, edi.flood.deltaF);
684 fprintf(fp, EDcol_efmt, edi.flood.vecs.fproj[i]);
690 /* From flood.xproj compute the Vfl(x) at this point */
691 static real flood_energy(t_edpar *edi, int64_t step)
693 /* compute flooding energy Vfl
694 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
695 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
696 it is already computed by make_edi and stored in stpsz[i]
698 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
705 /* Each time this routine is called (i.e. each time step), we add a small
706 * value to the reference projection. This way a harmonic restraint towards
707 * a moving reference is realized. If no value for the additive constant
708 * is provided in the edi file, the reference will not change. */
709 if (edi->flood.bHarmonic)
711 for (i = 0; i < edi->flood.vecs.neig; i++)
713 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i] + step * edi->flood.referenceProjectionSlope[i];
718 /* Compute sum which will be the exponent of the exponential */
719 for (i = 0; i < edi->flood.vecs.neig; i++)
721 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
722 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]);
725 /* Compute the Gauss function*/
726 if (edi->flood.bHarmonic)
728 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
732 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*std::exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
739 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
740 static void flood_forces(t_edpar *edi)
742 /* compute the forces in the subspace of the flooding eigenvectors
743 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
746 real energy = edi->flood.Vfl;
749 if (edi->flood.bHarmonic)
751 for (i = 0; i < edi->flood.vecs.neig; i++)
753 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
758 for (i = 0; i < edi->flood.vecs.neig; i++)
760 /* if Efl is zero the forces are zero if not use the formula */
761 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;
768 /*!\brief Raise forces from subspace into cartesian space.
769 * This function lifts the forces from the subspace to the cartesian space
770 * all the values not contained in the subspace are assumed to be zero and then
771 * a coordinate transformation from eigenvector to cartesian vectors is performed
772 * The nonexistent values don't have to be set to zero explicitly, they would occur
773 * as zero valued summands, hence we just stop to compute this part of the sum.
774 * For every atom we add all the contributions to this atom from all the different eigenvectors.
775 * NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
776 * field forces_cart prior the computation, but we compute the forces separately
777 * to have them accessible for diagnostics
779 * \param[in] edi Essential dynamics input parameters
780 * \param[out] forces_cart The cartesian forces
783 void flood_blowup(const t_edpar &edi, rvec *forces_cart)
785 const real * forces_sub = edi.flood.vecs.fproj;
786 /* Calculate the cartesian forces for the local atoms */
788 /* Clear forces first */
789 for (int j = 0; j < edi.sav.nr_loc; j++)
791 clear_rvec(forces_cart[j]);
794 /* Now compute atomwise */
795 for (int j = 0; j < edi.sav.nr_loc; j++)
797 /* Compute forces_cart[edi.sav.anrs[j]] */
798 for (int eig = 0; eig < edi.flood.vecs.neig; eig++)
801 /* Force vector is force * eigenvector (compute only atom j) */
802 svmul(forces_sub[eig], edi.flood.vecs.vec[eig][edi.sav.c_ind[j]], addedForce);
803 /* Add this vector to the cartesian forces */
804 rvec_inc(forces_cart[j], addedForce);
812 /* Update the values of Efl, deltaF depending on tau and Vfl */
813 static void update_adaption(t_edpar *edi)
815 /* this function updates the parameter Efl and deltaF according to the rules given in
816 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
819 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
821 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
822 /* check if restrain (inverted flooding) -> don't let EFL become positive */
823 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
828 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
833 static void do_single_flood(
841 gmx_bool bNS) /* Are we in a neighbor searching step? */
844 matrix rotmat; /* rotation matrix */
845 matrix tmat; /* inverse rotation */
846 rvec transvec; /* translation vector */
848 struct t_do_edsam *buf;
851 buf = edi->buf->do_edsam;
854 /* Broadcast the positions of the AVERAGE structure such that they are known on
855 * every processor. Each node contributes its local positions x and stores them in
856 * the collective ED array buf->xcoll */
857 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
858 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
860 /* Only assembly REFERENCE positions if their indices differ from the average ones */
863 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
864 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
867 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
868 * We do not need to update the shifts until the next NS step */
869 buf->bUpdateShifts = FALSE;
871 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
872 * as well as the indices in edi->sav.anrs */
874 /* Fit the reference indices to the reference structure */
877 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
881 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
884 /* Now apply the translation and rotation to the ED structure */
885 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
887 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
888 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
890 if (FALSE == edi->flood.bConstForce)
892 /* Compute Vfl(x) from flood.xproj */
893 edi->flood.Vfl = flood_energy(edi, step);
895 update_adaption(edi);
897 /* Compute the flooding forces */
901 /* Translate them into cartesian positions */
902 flood_blowup(*edi, edi->flood.forces_cartesian);
904 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
905 /* Each node rotates back its local forces */
906 transpose(rotmat, tmat);
907 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
909 /* Finally add forces to the main force variable */
910 for (i = 0; i < edi->sav.nr_loc; i++)
912 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
915 /* Output is written by the master process */
916 if (do_per_step(step, edi->outfrq) && MASTER(cr))
918 /* Output how well we fit to the reference */
921 /* Indices of reference and average structures are identical,
922 * thus we can calculate the rmsd to SREF using xcoll */
923 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
927 /* We have to translate & rotate the reference atoms first */
928 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
929 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
932 write_edo_flood(*edi, edo, rmsdev);
937 /* Main flooding routine, called from do_force */
938 extern void do_flood(const t_commrec *cr,
939 const t_inputrec *ir,
952 /* Write time to edo, when required. Output the time anyhow since we need
953 * it in the output file for ED constraints. */
954 if (MASTER(cr) && do_per_step(step, edi->outfrq))
956 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
959 if (ed->eEDtype != eEDflood)
966 /* Call flooding for one matrix */
967 if (edi->flood.vecs.neig)
969 do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS);
976 /* Called by init_edi, configure some flooding related variables and structures,
977 * print headers to output files */
978 static void init_flood(t_edpar *edi, gmx_edsam * ed, real dt)
983 edi->flood.Efl = edi->flood.constEfl;
987 if (edi->flood.vecs.neig)
989 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
990 ed->eEDtype = eEDflood;
992 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
994 if (edi->flood.bConstForce)
996 /* We have used stpsz as a vehicle to carry the fproj values for constant
997 * force flooding. Now we copy that to flood.vecs.fproj. Note that
998 * in const force flooding, fproj is never changed. */
999 for (i = 0; i < edi->flood.vecs.neig; i++)
1001 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1003 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1004 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1012 /*********** Energy book keeping ******/
1013 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1022 srenew(names, count);
1023 sprintf(buf, "Vfl_%d", count);
1024 names[count-1] = gmx_strdup(buf);
1025 actual = actual->next_edi;
1032 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1034 /*fl has to be big enough to capture nnames-many entries*/
1043 Vfl[count-1] = actual->flood.Vfl;
1044 actual = actual->next_edi;
1047 if (nnames != count-1)
1049 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1052 /************* END of FLOODING IMPLEMENTATION ****************************/
1056 /* This function opens the ED input and output files, reads in all datasets it finds
1057 * in the input file, and cross-checks whether the .edi file information is consistent
1058 * with the essential dynamics data found in the checkpoint file (if present).
1059 * gmx make_edi can be used to create an .edi input file.
1061 static std::unique_ptr<gmx::EssentialDynamics> ed_open(
1063 ObservablesHistory *oh,
1064 const char *ediFileName,
1065 const char *edoFileName,
1067 const gmx_output_env_t *oenv,
1068 const t_commrec *cr)
1070 auto edHandle = gmx::compat::make_unique<gmx::EssentialDynamics>();
1071 auto ed = edHandle->getLegacyED();
1074 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1075 ed->eEDtype = eEDedsam;
1081 // If we start from a checkpoint file, we already have an edsamHistory struct
1082 if (oh->edsamHistory == nullptr)
1084 oh->edsamHistory = gmx::compat::make_unique<edsamhistory_t>(edsamhistory_t {});
1086 edsamhistory_t *EDstate = oh->edsamHistory.get();
1088 /* Read the edi input file: */
1089 nED = read_edi_file(ediFileName, ed->edpar, natoms);
1091 /* Make sure the checkpoint was produced in a run using this .edi file */
1092 if (EDstate->bFromCpt)
1094 crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1100 init_edsamstate(ed, EDstate);
1102 /* The master opens the ED output file */
1105 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1109 ed->edo = xvgropen(edoFileName,
1110 "Essential dynamics / flooding output",
1112 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1114 /* Make a descriptive legend */
1115 write_edo_legend(ed, EDstate->nED, oenv);
1122 /* Broadcasts the structure data */
1123 static void bc_ed_positions(const t_commrec *cr, struct gmx_edx *s, int stype)
1125 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1126 snew_bc(cr, s->x, s->nr ); /* Positions */
1127 nblock_bc(cr, s->nr, s->anrs );
1128 nblock_bc(cr, s->nr, s->x );
1130 /* For the average & reference structures we need an array for the collective indices,
1131 * and we need to broadcast the masses as well */
1132 if (stype == eedAV || stype == eedREF)
1134 /* We need these additional variables in the parallel case: */
1135 snew(s->c_ind, s->nr ); /* Collective indices */
1136 /* Local atom indices get assigned in dd_make_local_group_indices.
1137 * There, also memory is allocated */
1138 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1139 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1140 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1143 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1144 if (stype == eedREF)
1146 snew_bc(cr, s->m, s->nr);
1147 nblock_bc(cr, s->nr, s->m);
1150 /* For the average structure we might need the masses for mass-weighting */
1153 snew_bc(cr, s->sqrtm, s->nr);
1154 nblock_bc(cr, s->nr, s->sqrtm);
1155 snew_bc(cr, s->m, s->nr);
1156 nblock_bc(cr, s->nr, s->m);
1161 /* Broadcasts the eigenvector data */
1162 static void bc_ed_vecs(const t_commrec *cr, t_eigvec *ev, int length)
1166 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1167 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1168 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1169 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1170 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1172 nblock_bc(cr, ev->neig, ev->ieig );
1173 nblock_bc(cr, ev->neig, ev->stpsz );
1174 nblock_bc(cr, ev->neig, ev->xproj );
1175 nblock_bc(cr, ev->neig, ev->fproj );
1176 nblock_bc(cr, ev->neig, ev->refproj);
1178 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1179 for (i = 0; i < ev->neig; i++)
1181 snew_bc(cr, ev->vec[i], length);
1182 nblock_bc(cr, length, ev->vec[i]);
1188 /* Broadcasts the ED / flooding data to other nodes
1189 * and allocates memory where needed */
1190 static void broadcast_ed_data(const t_commrec *cr, gmx_edsam * ed, int numedis)
1196 /* Master lets the other nodes know if its ED only or also flooding */
1197 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1199 snew_bc(cr, ed->edpar, 1);
1200 /* Now transfer the ED data set(s) */
1202 for (nr = 0; nr < numedis; nr++)
1204 /* Broadcast a single ED data set */
1207 /* Broadcast positions */
1208 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1209 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1210 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1211 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1213 /* Broadcast eigenvectors */
1214 bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr);
1215 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr);
1216 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr);
1217 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr);
1218 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr);
1219 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr);
1220 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1221 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr);
1223 /* For harmonic restraints the reference projections can change with time */
1224 if (edi->flood.bHarmonic)
1226 snew_bc(cr, edi->flood.initialReferenceProjection, edi->flood.vecs.neig);
1227 snew_bc(cr, edi->flood.referenceProjectionSlope, edi->flood.vecs.neig);
1228 nblock_bc(cr, edi->flood.vecs.neig, edi->flood.initialReferenceProjection);
1229 nblock_bc(cr, edi->flood.vecs.neig, edi->flood.referenceProjectionSlope);
1233 /* Set the pointer to the next ED group */
1236 snew_bc(cr, edi->next_edi, 1);
1237 edi = edi->next_edi;
1243 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1244 static void init_edi(const gmx_mtop_t *mtop, t_edpar *edi)
1247 real totalmass = 0.0;
1250 /* NOTE Init_edi is executed on the master process only
1251 * The initialized data sets are then transmitted to the
1252 * other nodes in broadcast_ed_data */
1254 /* evaluate masses (reference structure) */
1255 snew(edi->sref.m, edi->sref.nr);
1257 for (i = 0; i < edi->sref.nr; i++)
1261 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1265 edi->sref.m[i] = 1.0;
1268 /* Check that every m > 0. Bad things will happen otherwise. */
1269 if (edi->sref.m[i] <= 0.0)
1271 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1272 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1273 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1274 "atoms from the reference structure by creating a proper index group.\n",
1275 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1278 totalmass += edi->sref.m[i];
1280 edi->sref.mtot = totalmass;
1282 /* Masses m and sqrt(m) for the average structure. Note that m
1283 * is needed if forces have to be evaluated in do_edsam */
1284 snew(edi->sav.sqrtm, edi->sav.nr );
1285 snew(edi->sav.m, edi->sav.nr );
1286 for (i = 0; i < edi->sav.nr; i++)
1288 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1291 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1295 edi->sav.sqrtm[i] = 1.0;
1298 /* Check that every m > 0. Bad things will happen otherwise. */
1299 if (edi->sav.sqrtm[i] <= 0.0)
1301 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1302 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1303 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1304 "atoms from the average structure by creating a proper index group.\n",
1305 i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1309 /* put reference structure in origin */
1310 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1314 translate_x(edi->sref.x, edi->sref.nr, com);
1316 /* Init ED buffer */
1321 static void check(const char *line, const char *label)
1323 if (!strstr(line, label))
1325 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1330 static int read_checked_edint(FILE *file, const char *label)
1332 char line[STRLEN+1];
1335 fgets2 (line, STRLEN, file);
1337 fgets2 (line, STRLEN, file);
1338 sscanf (line, max_ev_fmt_d, &idum);
1343 static int read_edint(FILE *file, bool *bEOF)
1345 char line[STRLEN+1];
1350 eof = fgets2 (line, STRLEN, file);
1356 eof = fgets2 (line, STRLEN, file);
1362 sscanf (line, max_ev_fmt_d, &idum);
1368 static real read_checked_edreal(FILE *file, const char *label)
1370 char line[STRLEN+1];
1374 fgets2 (line, STRLEN, file);
1376 fgets2 (line, STRLEN, file);
1377 sscanf (line, max_ev_fmt_lf, &rdum);
1378 return static_cast<real>(rdum); /* always read as double and convert to single */
1382 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1385 char line[STRLEN+1];
1389 for (i = 0; i < number; i++)
1391 fgets2 (line, STRLEN, file);
1392 sscanf (line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1393 anrs[i]--; /* we are reading FORTRAN indices */
1394 for (j = 0; j < 3; j++)
1396 x[i][j] = d[j]; /* always read as double and convert to single */
1403 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1404 * \param[in] in the file to read from
1405 * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1406 * \param[out] vec the eigen-vectors
1407 * \param[in] nEig the number of eigenvectors
1409 void scan_edvec(FILE *in, int numAtoms, rvec ***vec, int nEig)
1412 for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1414 snew((*vec)[iEigenvector], numAtoms);
1415 for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1417 char line[STRLEN+1];
1419 fgets2(line, STRLEN, in);
1420 sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1421 (*vec)[iEigenvector][iAtom][XX] = x;
1422 (*vec)[iEigenvector][iAtom][YY] = y;
1423 (*vec)[iEigenvector][iAtom][ZZ] = z;
1428 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1429 * \param[in] in input file
1430 * \param[in] nrAtoms number of atoms/coordinates
1431 * \param[out] tvec the eigenvector
1433 void read_edvec(FILE *in, int nrAtoms, t_eigvec *tvec)
1435 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1436 if (tvec->neig <= 0)
1441 snew(tvec->ieig, tvec->neig);
1442 snew(tvec->stpsz, tvec->neig);
1443 for (int i = 0; i < tvec->neig; i++)
1445 char line[STRLEN+1];
1446 fgets2(line, STRLEN, in);
1447 int numEigenVectors;
1449 const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1452 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1454 tvec->ieig[i] = numEigenVectors;
1455 tvec->stpsz[i] = stepSize;
1456 } /* end of loop over eigenvectors */
1458 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1460 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1462 * Eigenvector and its intial reference and slope are stored on the
1463 * same line in the .edi format. To avoid re-winding the .edi file,
1464 * the reference eigen-vector and reference data are read in one go.
1466 * \param[in] in input file
1467 * \param[in] nrAtoms number of atoms/coordinates
1468 * \param[out] tvec the eigenvector
1469 * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1470 * \param[out] referenceProjectionSlope The slope of the reference projections.
1472 bool readEdVecWithReferenceProjection(FILE *in, int nrAtoms, t_eigvec *tvec, real ** initialReferenceProjection, real ** referenceProjectionSlope)
1474 bool providesReference = false;
1475 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1476 if (tvec->neig <= 0)
1481 snew(tvec->ieig, tvec->neig);
1482 snew(tvec->stpsz, tvec->neig);
1483 snew(*initialReferenceProjection, tvec->neig);
1484 snew(*referenceProjectionSlope, tvec->neig);
1485 for (int i = 0; i < tvec->neig; i++)
1487 char line[STRLEN+1];
1488 fgets2 (line, STRLEN, in);
1489 int numEigenVectors;
1490 double stepSize = 0.;
1491 double referenceProjection = 0.;
1492 double referenceSlope = 0.;
1493 const int nscan = sscanf(line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize, &referenceProjection, &referenceSlope);
1494 /* Zero out values which were not scanned */
1498 /* Every 4 values read, including reference position */
1499 providesReference = true;
1502 /* A reference position is provided */
1503 providesReference = true;
1504 /* No value for slope, set to 0 */
1505 referenceSlope = 0.0;
1508 /* No values for reference projection and slope, set to 0 */
1509 referenceProjection = 0.0;
1510 referenceSlope = 0.0;
1513 gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1515 (*initialReferenceProjection)[i] = referenceProjection;
1516 (*referenceProjectionSlope)[i] = referenceSlope;
1518 tvec->ieig[i] = numEigenVectors;
1519 tvec->stpsz[i] = stepSize;
1520 } /* end of loop over eigenvectors */
1522 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1523 return providesReference;
1526 /*!\brief Allocate working memory for the eigenvectors.
1527 * \param[in,out] tvec eigenvector for which memory will be allocated
1529 void setup_edvec(t_eigvec *tvec)
1531 snew(tvec->xproj, tvec->neig);
1532 snew(tvec->fproj, tvec->neig);
1533 snew(tvec->refproj, tvec->neig);
1538 /* Check if the same atom indices are used for reference and average positions */
1539 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1544 /* If the number of atoms differs between the two structures,
1545 * they cannot be identical */
1546 if (sref.nr != sav.nr)
1551 /* Now that we know that both stuctures have the same number of atoms,
1552 * check if also the indices are identical */
1553 for (i = 0; i < sav.nr; i++)
1555 if (sref.anrs[i] != sav.anrs[i])
1560 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1568 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1569 constexpr int c_oldestUnsupportedMagicNumber = 666;
1570 //! Oldest (lowest) magic number indicating supported essential dynamics input
1571 constexpr int c_oldestSupportedMagicNumber = 669;
1572 //! Latest (highest) magic number indicating supported essential dynamics input
1573 constexpr int c_latestSupportedMagicNumber = 670;
1575 /*!\brief Set up essential dynamics work parameters.
1576 * \param[in] edi Essential dynamics working structure.
1578 void setup_edi(t_edpar *edi)
1580 snew(edi->sref.x_old, edi->sref.nr);
1581 edi->sref.sqrtm = nullptr;
1582 snew(edi->sav.x_old, edi->sav.nr);
1583 if (edi->star.nr > 0)
1585 edi->star.sqrtm = nullptr;
1587 if (edi->sori.nr > 0)
1589 edi->sori.sqrtm = nullptr;
1591 setup_edvec(&edi->vecs.linacc);
1592 setup_edvec(&edi->vecs.mon);
1593 setup_edvec(&edi->vecs.linfix);
1594 setup_edvec(&edi->vecs.linacc);
1595 setup_edvec(&edi->vecs.radfix);
1596 setup_edvec(&edi->vecs.radacc);
1597 setup_edvec(&edi->vecs.radcon);
1598 setup_edvec(&edi->flood.vecs);
1601 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1602 * \param[in] magicNumber indicates the essential dynamics input file version
1603 * \returns true if CONST_FORCE_FLOODING is supported
1605 bool ediFileHasConstForceFlooding(int magicNumber)
1607 return magicNumber > c_oldestSupportedMagicNumber;
1610 /*!\brief Reports reading success of the essential dynamics input file magic number.
1611 * \param[in] in pointer to input file
1612 * \param[out] magicNumber determines the edi file version
1613 * \returns true if setting file version from input file was successful.
1615 bool ediFileHasValidDataSet(FILE *in, int * magicNumber)
1617 /* the edi file is not free format, so expect problems if the input is corrupt. */
1618 bool reachedEndOfFile = true;
1619 /* check the magic number */
1620 *magicNumber = read_edint(in, &reachedEndOfFile);
1621 /* Check whether we have reached the end of the input file */
1622 return !reachedEndOfFile;
1625 /*!\brief Trigger fatal error if magic number is unsupported.
1626 * \param[in] magicNumber A magic number read from the edi file
1627 * \param[in] fn name of input file for error reporting
1629 void exitWhenMagicNumberIsInvalid(int magicNumber, const char * fn)
1632 if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1634 if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1636 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1640 gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1645 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1647 * \param[in] in input file
1648 * \param[in] edi essential dynamics parameters
1649 * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1650 * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1651 * \param[in] fn the file name of the input file for error reporting
1653 void read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, bool hasConstForceFlooding, const char *fn)
1656 /* check the number of atoms */
1657 edi->nini = read_edint(in, &bEOF);
1658 if (edi->nini != nr_mdatoms)
1660 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1663 /* Done checking. For the rest we blindly trust the input */
1664 edi->fitmas = read_checked_edint(in, "FITMAS");
1665 edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS");
1666 edi->outfrq = read_checked_edint(in, "OUTFRQ");
1667 edi->maxedsteps = read_checked_edint(in, "MAXLEN");
1668 edi->slope = read_checked_edreal(in, "SLOPECRIT");
1670 edi->presteps = read_checked_edint(in, "PRESTEPS");
1671 edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1672 edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1673 edi->flood.tau = read_checked_edreal(in, "TAU");
1674 edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1675 edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1676 edi->flood.kT = read_checked_edreal(in, "KT");
1677 edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1678 if (hasConstForceFlooding)
1680 edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1684 edi->flood.bConstForce = FALSE;
1686 edi->sref.nr = read_checked_edint(in, "NREF");
1688 /* allocate space for reference positions and read them */
1689 snew(edi->sref.anrs, edi->sref.nr);
1690 snew(edi->sref.x, edi->sref.nr);
1691 read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1693 /* average positions. they define which atoms will be used for ED sampling */
1694 edi->sav.nr = read_checked_edint(in, "NAV");
1695 snew(edi->sav.anrs, edi->sav.nr);
1696 snew(edi->sav.x, edi->sav.nr);
1697 read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1699 /* Check if the same atom indices are used for reference and average positions */
1700 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1704 read_edvec(in, edi->sav.nr, &edi->vecs.mon);
1705 read_edvec(in, edi->sav.nr, &edi->vecs.linfix);
1706 read_edvec(in, edi->sav.nr, &edi->vecs.linacc);
1707 read_edvec(in, edi->sav.nr, &edi->vecs.radfix);
1708 read_edvec(in, edi->sav.nr, &edi->vecs.radacc);
1709 read_edvec(in, edi->sav.nr, &edi->vecs.radcon);
1711 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1712 bool bHaveReference = false;
1713 if (edi->flood.bHarmonic)
1715 bHaveReference = readEdVecWithReferenceProjection(in, edi->sav.nr, &edi->flood.vecs, &edi->flood.initialReferenceProjection, &edi->flood.referenceProjectionSlope);
1719 read_edvec(in, edi->sav.nr, &edi->flood.vecs);
1722 /* target positions */
1723 edi->star.nr = read_edint(in, &bEOF);
1724 if (edi->star.nr > 0)
1726 snew(edi->star.anrs, edi->star.nr);
1727 snew(edi->star.x, edi->star.nr);
1728 read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1731 /* positions defining origin of expansion circle */
1732 edi->sori.nr = read_edint(in, &bEOF);
1733 if (edi->sori.nr > 0)
1737 /* Both an -ori structure and a at least one manual reference point have been
1738 * specified. That's ambiguous and probably not intentional. */
1739 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1740 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1742 snew(edi->sori.anrs, edi->sori.nr);
1743 snew(edi->sori.x, edi->sori.nr);
1744 read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1748 } // anonymous namespace
1750 /* Read in the edi input file. Note that it may contain several ED data sets which were
1751 * achieved by concatenating multiple edi files. The standard case would be a single ED
1752 * data set, though. */
1753 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1756 t_edpar *curr_edi, *last_edi;
1761 /* This routine is executed on the master only */
1763 /* Open the .edi parameter input file */
1764 in = gmx_fio_fopen(fn, "r");
1765 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1767 /* Now read a sequence of ED input parameter sets from the edi file */
1770 int ediFileMagicNumber;
1771 while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1773 exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1774 bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1775 read_edi(in, curr_edi, nr_mdatoms, hasConstForceFlooding, fn);
1776 setup_edi(curr_edi);
1779 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1780 /* We need to allocate space for the data: */
1782 /* Point the 'next_edi' entry to the next edi: */
1783 curr_edi->next_edi = edi_read;
1784 /* Keep the curr_edi pointer for the case that the next group is empty: */
1785 last_edi = curr_edi;
1786 /* Let's prepare to read in the next edi data set: */
1787 curr_edi = edi_read;
1791 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1794 /* Terminate the edi group list with a NULL pointer: */
1795 last_edi->next_edi = nullptr;
1797 fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1799 /* Close the .edi file again */
1806 struct t_fit_to_ref {
1807 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1810 /* Fit the current positions to the reference positions
1811 * Do not actually do the fit, just return rotation and translation.
1812 * Note that the COM of the reference structure was already put into
1813 * the origin by init_edi. */
1814 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1815 rvec transvec, /* The translation vector */
1816 matrix rotmat, /* The rotation matrix */
1817 t_edpar *edi) /* Just needed for do_edfit */
1819 rvec com; /* center of mass */
1821 struct t_fit_to_ref *loc;
1824 /* Allocate memory the first time this routine is called for each edi group */
1825 if (nullptr == edi->buf->fit_to_ref)
1827 snew(edi->buf->fit_to_ref, 1);
1828 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1830 loc = edi->buf->fit_to_ref;
1832 /* We do not touch the original positions but work on a copy. */
1833 for (i = 0; i < edi->sref.nr; i++)
1835 copy_rvec(xcoll[i], loc->xcopy[i]);
1838 /* Calculate the center of mass */
1839 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1841 transvec[XX] = -com[XX];
1842 transvec[YY] = -com[YY];
1843 transvec[ZZ] = -com[ZZ];
1845 /* Subtract the center of mass from the copy */
1846 translate_x(loc->xcopy, edi->sref.nr, transvec);
1848 /* Determine the rotation matrix */
1849 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1853 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1854 int nat, /* How many positions are there? */
1855 rvec transvec, /* The translation vector */
1856 matrix rotmat) /* The rotation matrix */
1859 translate_x(x, nat, transvec);
1862 rotate_x(x, nat, rotmat);
1866 /* Gets the rms deviation of the positions to the structure s */
1867 /* fit_to_structure has to be called before calling this routine! */
1868 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1869 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1875 for (i = 0; i < s->nr; i++)
1877 rmsd += distance2(s->x[i], x[i]);
1880 rmsd /= static_cast<real>(s->nr);
1887 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1892 if (ed->eEDtype != eEDnone)
1894 /* Loop over ED groups */
1898 /* Local atoms of the reference structure (for fitting), need only be assembled
1899 * if their indices differ from the average ones */
1902 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1903 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1906 /* Local atoms of the average structure (on these ED will be performed) */
1907 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1908 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1910 /* Indicate that the ED shift vectors for this structure need to be updated
1911 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1912 edi->buf->do_edsam->bUpdateShifts = TRUE;
1914 /* Set the pointer to the next ED group (if any) */
1915 edi = edi->next_edi;
1921 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1932 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1933 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1934 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1938 xu[XX] = x[XX]-tx*box[XX][XX];
1939 xu[YY] = x[YY]-ty*box[YY][YY];
1940 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1946 /*!\brief Apply fixed linear constraints to essential dynamics variable.
1947 * \param[in,out] xcoll The collected coordinates.
1948 * \param[in] edi the essential dynamics parameters
1949 * \param[in] step the current simulation step
1951 void do_linfix(rvec *xcoll, const t_edpar &edi, int64_t step)
1953 /* loop over linfix vectors */
1954 for (int i = 0; i < edi.vecs.linfix.neig; i++)
1956 /* calculate the projection */
1957 real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
1959 /* calculate the correction */
1960 real preFactor = edi.vecs.linfix.refproj[i] + step*edi.vecs.linfix.stpsz[i] - proj;
1962 /* apply the correction */
1963 preFactor /= edi.sav.sqrtm[i];
1964 for (int j = 0; j < edi.sav.nr; j++)
1966 rvec differenceVector;
1967 svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
1968 rvec_inc(xcoll[j], differenceVector);
1973 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
1974 * \param[in,out] xcoll The collected coordinates.
1975 * \param[in] edi the essential dynamics parameters
1977 void do_linacc(rvec *xcoll, t_edpar *edi)
1979 /* loop over linacc vectors */
1980 for (int i = 0; i < edi->vecs.linacc.neig; i++)
1982 /* calculate the projection */
1983 real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
1985 /* calculate the correction */
1986 real preFactor = 0.0;
1987 if (edi->vecs.linacc.stpsz[i] > 0.0)
1989 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1991 preFactor = edi->vecs.linacc.refproj[i] - proj;
1994 if (edi->vecs.linacc.stpsz[i] < 0.0)
1996 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1998 preFactor = edi->vecs.linacc.refproj[i] - proj;
2002 /* apply the correction */
2003 preFactor /= edi->sav.sqrtm[i];
2004 for (int j = 0; j < edi->sav.nr; j++)
2006 rvec differenceVector;
2007 svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2008 rvec_inc(xcoll[j], differenceVector);
2011 /* new positions will act as reference */
2012 edi->vecs.linacc.refproj[i] = proj + preFactor;
2018 static void do_radfix(rvec *xcoll, t_edpar *edi)
2021 real *proj, rad = 0.0, ratio;
2025 if (edi->vecs.radfix.neig == 0)
2030 snew(proj, edi->vecs.radfix.neig);
2032 /* loop over radfix vectors */
2033 for (i = 0; i < edi->vecs.radfix.neig; i++)
2035 /* calculate the projections, radius */
2036 proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2037 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2041 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2042 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2044 /* loop over radfix vectors */
2045 for (i = 0; i < edi->vecs.radfix.neig; i++)
2047 proj[i] -= edi->vecs.radfix.refproj[i];
2049 /* apply the correction */
2050 proj[i] /= edi->sav.sqrtm[i];
2052 for (j = 0; j < edi->sav.nr; j++)
2054 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2055 rvec_inc(xcoll[j], vec_dum);
2063 static void do_radacc(rvec *xcoll, t_edpar *edi)
2066 real *proj, rad = 0.0, ratio = 0.0;
2070 if (edi->vecs.radacc.neig == 0)
2075 snew(proj, edi->vecs.radacc.neig);
2077 /* loop over radacc vectors */
2078 for (i = 0; i < edi->vecs.radacc.neig; i++)
2080 /* calculate the projections, radius */
2081 proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2082 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2086 /* only correct when radius decreased */
2087 if (rad < edi->vecs.radacc.radius)
2089 ratio = edi->vecs.radacc.radius/rad - 1.0;
2093 edi->vecs.radacc.radius = rad;
2096 /* loop over radacc vectors */
2097 for (i = 0; i < edi->vecs.radacc.neig; i++)
2099 proj[i] -= edi->vecs.radacc.refproj[i];
2101 /* apply the correction */
2102 proj[i] /= edi->sav.sqrtm[i];
2104 for (j = 0; j < edi->sav.nr; j++)
2106 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2107 rvec_inc(xcoll[j], vec_dum);
2114 struct t_do_radcon {
2118 static void do_radcon(rvec *xcoll, t_edpar *edi)
2121 real rad = 0.0, ratio = 0.0;
2122 struct t_do_radcon *loc;
2127 if (edi->buf->do_radcon != nullptr)
2134 snew(edi->buf->do_radcon, 1);
2136 loc = edi->buf->do_radcon;
2138 if (edi->vecs.radcon.neig == 0)
2145 snew(loc->proj, edi->vecs.radcon.neig);
2148 /* loop over radcon vectors */
2149 for (i = 0; i < edi->vecs.radcon.neig; i++)
2151 /* calculate the projections, radius */
2152 loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2153 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2156 /* only correct when radius increased */
2157 if (rad > edi->vecs.radcon.radius)
2159 ratio = edi->vecs.radcon.radius/rad - 1.0;
2161 /* loop over radcon vectors */
2162 for (i = 0; i < edi->vecs.radcon.neig; i++)
2164 /* apply the correction */
2165 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2166 loc->proj[i] /= edi->sav.sqrtm[i];
2167 loc->proj[i] *= ratio;
2169 for (j = 0; j < edi->sav.nr; j++)
2171 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2172 rvec_inc(xcoll[j], vec_dum);
2179 edi->vecs.radcon.radius = rad;
2185 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, int64_t step)
2190 /* subtract the average positions */
2191 for (i = 0; i < edi->sav.nr; i++)
2193 rvec_dec(xcoll[i], edi->sav.x[i]);
2196 /* apply the constraints */
2199 do_linfix(xcoll, *edi, step);
2201 do_linacc(xcoll, edi);
2204 do_radfix(xcoll, edi);
2206 do_radacc(xcoll, edi);
2207 do_radcon(xcoll, edi);
2209 /* add back the average positions */
2210 for (i = 0; i < edi->sav.nr; i++)
2212 rvec_inc(xcoll[i], edi->sav.x[i]);
2219 /*!\brief Write out the projections onto the eigenvectors.
2220 * The order of output corresponds to ed_output_legend().
2221 * \param[in] edi The essential dyanmics parameters
2222 * \param[in] fp The output file
2223 * \param[in] rmsd the rmsd to the reference structure
2225 void write_edo(const t_edpar &edi, FILE *fp, real rmsd)
2227 /* Output how well we fit to the reference structure */
2228 fprintf(fp, EDcol_ffmt, rmsd);
2230 for (int i = 0; i < edi.vecs.mon.neig; i++)
2232 fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2235 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2237 fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2240 for (int i = 0; i < edi.vecs.linacc.neig; i++)
2242 fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2245 for (int i = 0; i < edi.vecs.radfix.neig; i++)
2247 fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2249 if (edi.vecs.radfix.neig)
2251 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2254 for (int i = 0; i < edi.vecs.radacc.neig; i++)
2256 fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2258 if (edi.vecs.radacc.neig)
2260 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2263 for (int i = 0; i < edi.vecs.radcon.neig; i++)
2265 fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2267 if (edi.vecs.radcon.neig)
2269 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2275 /* Returns if any constraints are switched on */
2276 static int ed_constraints(int edtype, t_edpar *edi)
2278 if (edtype == eEDedsam || edtype == eEDflood)
2280 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2281 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2282 edi->vecs.radcon.neig);
2288 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for flooding/
2289 * umbrella sampling simulations. */
2290 static void copyEvecReference(t_eigvec* floodvecs, real * initialReferenceProjection)
2295 if (nullptr == initialReferenceProjection)
2297 snew(initialReferenceProjection, floodvecs->neig);
2300 for (i = 0; i < floodvecs->neig; i++)
2302 initialReferenceProjection[i] = floodvecs->refproj[i];
2307 /* Call on MASTER only. Check whether the essential dynamics / flooding
2308 * groups of the checkpoint file are consistent with the provided .edi file. */
2309 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam &ed, edsamhistory_t *EDstate)
2311 t_edpar *edi = nullptr; /* points to a single edi data set */
2315 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2317 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2318 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2319 "it must also continue with/without ED constraints when checkpointing.\n"
2320 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2321 "from without a checkpoint.\n");
2326 while (edi != nullptr)
2328 /* Check number of atoms in the reference and average structures */
2329 if (EDstate->nref[edinum] != edi->sref.nr)
2331 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2332 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2333 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2335 if (EDstate->nav[edinum] != edi->sav.nr)
2337 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2338 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2339 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2341 edi = edi->next_edi;
2345 if (edinum != EDstate->nED)
2347 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2348 "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2349 "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2354 /* The edsamstate struct stores the information we need to make the ED group
2355 * whole again after restarts from a checkpoint file. Here we do the following:
2356 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2357 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2358 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2359 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2360 * all ED structures at the last time step. */
2361 static void init_edsamstate(gmx_edsam * ed, edsamhistory_t *EDstate)
2367 snew(EDstate->old_sref_p, EDstate->nED);
2368 snew(EDstate->old_sav_p, EDstate->nED);
2370 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2371 if (!EDstate->bFromCpt)
2373 snew(EDstate->nref, EDstate->nED);
2374 snew(EDstate->nav, EDstate->nED);
2377 /* Loop over all ED/flooding data sets (usually only one, though) */
2379 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2381 /* We always need the last reference and average positions such that
2382 * in the next time step we can make the ED group whole again
2383 * if the atoms do not have the correct PBC representation */
2384 if (EDstate->bFromCpt)
2386 /* Copy the last whole positions of reference and average group from .cpt */
2387 for (i = 0; i < edi->sref.nr; i++)
2389 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2391 for (i = 0; i < edi->sav.nr; i++)
2393 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2398 EDstate->nref[nr_edi-1] = edi->sref.nr;
2399 EDstate->nav [nr_edi-1] = edi->sav.nr;
2402 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2403 EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2404 EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2406 edi = edi->next_edi;
2411 /* Adds 'buf' to 'str' */
2412 static void add_to_string(char **str, char *buf)
2417 len = strlen(*str) + strlen(buf) + 1;
2423 static void add_to_string_aligned(char **str, char *buf)
2425 char buf_aligned[STRLEN];
2427 sprintf(buf_aligned, EDcol_sfmt, buf);
2428 add_to_string(str, buf_aligned);
2432 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, const char *value, const char *unit, char EDgroupchar)
2434 char tmp[STRLEN], tmp2[STRLEN];
2437 sprintf(tmp, "%c %s", EDgroupchar, value);
2438 add_to_string_aligned(LegendStr, tmp);
2439 sprintf(tmp2, "%s (%s)", tmp, unit);
2440 (*setname)[*nsets] = gmx_strdup(tmp2);
2445 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2451 for (i = 0; i < evec->neig; i++)
2453 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2454 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2459 /* Makes a legend for the xvg output file. Call on MASTER only! */
2460 static void write_edo_legend(gmx_edsam * ed, int nED, const gmx_output_env_t *oenv)
2462 t_edpar *edi = nullptr;
2464 int nr_edi, nsets, n_flood, n_edsam;
2465 const char **setname;
2467 char *LegendStr = nullptr;
2472 fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2474 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2476 fprintf(ed->edo, "#\n");
2477 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2478 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2479 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2480 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2481 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2482 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2483 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2484 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2485 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2487 if (edi->flood.vecs.neig)
2489 /* If in any of the groups we find a flooding vector, flooding is turned on */
2490 ed->eEDtype = eEDflood;
2492 /* Print what flavor of flooding we will do */
2493 if (0 == edi->flood.tau) /* constant flooding strength */
2495 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2496 if (edi->flood.bHarmonic)
2498 fprintf(ed->edo, ", harmonic");
2501 else /* adaptive flooding */
2503 fprintf(ed->edo, ", adaptive");
2506 fprintf(ed->edo, "\n");
2508 edi = edi->next_edi;
2511 /* Print a nice legend */
2513 LegendStr[0] = '\0';
2514 sprintf(buf, "# %6s", "time");
2515 add_to_string(&LegendStr, buf);
2517 /* Calculate the maximum number of columns we could end up with */
2520 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2522 nsets += 5 +edi->vecs.mon.neig
2523 +edi->vecs.linfix.neig
2524 +edi->vecs.linacc.neig
2525 +edi->vecs.radfix.neig
2526 +edi->vecs.radacc.neig
2527 +edi->vecs.radcon.neig
2528 + 6*edi->flood.vecs.neig;
2529 edi = edi->next_edi;
2531 snew(setname, nsets);
2533 /* In the mdrun time step in a first function call (do_flood()) the flooding
2534 * forces are calculated and in a second function call (do_edsam()) the
2535 * ED constraints. To get a corresponding legend, we need to loop twice
2536 * over the edi groups and output first the flooding, then the ED part */
2538 /* The flooding-related legend entries, if flooding is done */
2540 if (eEDflood == ed->eEDtype)
2543 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2545 /* Always write out the projection on the flooding EVs. Of course, this can also
2546 * be achieved with the monitoring option in do_edsam() (if switched on by the
2547 * user), but in that case the positions need to be communicated in do_edsam(),
2548 * which is not necessary when doing flooding only. */
2549 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2551 for (i = 0; i < edi->flood.vecs.neig; i++)
2553 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2554 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2556 /* Output the current reference projection if it changes with time;
2557 * this can happen when flooding is used as harmonic restraint */
2558 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2560 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2561 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2564 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2565 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2567 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2568 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2571 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2572 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2574 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2576 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2577 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2580 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2581 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2584 edi = edi->next_edi;
2585 } /* End of flooding-related legend entries */
2589 /* Now the ED-related entries, if essential dynamics is done */
2591 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2593 if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2595 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2597 /* Essential dynamics, projections on eigenvectors */
2598 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2599 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2600 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2601 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2602 if (edi->vecs.radfix.neig)
2604 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2606 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2607 if (edi->vecs.radacc.neig)
2609 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2611 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2612 if (edi->vecs.radcon.neig)
2614 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2617 edi = edi->next_edi;
2618 } /* end of 'pure' essential dynamics legend entries */
2619 n_edsam = nsets - n_flood;
2621 xvgr_legend(ed->edo, nsets, setname, oenv);
2624 fprintf(ed->edo, "#\n"
2625 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2626 n_flood, 1 == n_flood ? "" : "s",
2627 n_edsam, 1 == n_edsam ? "" : "s");
2628 fprintf(ed->edo, "%s", LegendStr);
2635 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2636 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2637 std::unique_ptr<gmx::EssentialDynamics> init_edsam(
2638 const char *ediFileName,
2639 const char *edoFileName,
2640 const gmx_mtop_t *mtop,
2641 const t_inputrec *ir,
2642 const t_commrec *cr,
2643 gmx::Constraints *constr,
2644 const t_state *globalState,
2645 ObservablesHistory *oh,
2646 const gmx_output_env_t *oenv,
2649 t_edpar *edi = nullptr; /* points to a single edi data set */
2651 int nED = 0; /* Number of ED data sets */
2652 rvec *x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2653 rvec *xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2654 rvec fit_transvec; /* translation ... */
2655 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2656 rvec *ref_x_old = nullptr; /* helper pointer */
2661 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2663 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName) )
2665 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or flooding\n"
2666 "simulation. Please also set the .edi file on the command line with -ei.\n");
2670 /* Open input and output files, allocate space for ED data structure */
2671 auto edHandle = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr);
2672 auto ed = edHandle->getLegacyED();
2673 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2674 constr->saveEdsamPointer(ed);
2676 /* Needed for initializing radacc radius in do_edsam */
2679 /* The input file is read by the master and the edi structures are
2680 * initialized here. Input is stored in ed->edpar. Then the edi
2681 * structures are transferred to the other nodes */
2684 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2685 * flooding vector, Essential dynamics can be applied to more than one structure
2686 * as well, but will be done in the order given in the edi file, so
2687 * expect different results for different order of edi file concatenation! */
2689 while (edi != nullptr)
2691 init_edi(mtop, edi);
2692 init_flood(edi, ed, ir->delta_t);
2693 edi = edi->next_edi;
2697 /* The master does the work here. The other nodes get the positions
2698 * not before dd_partition_system which is called after init_edsam */
2701 edsamhistory_t *EDstate = oh->edsamHistory.get();
2703 if (!EDstate->bFromCpt)
2705 /* Remove PBC, make molecule(s) subject to ED whole. */
2706 snew(x_pbc, mtop->natoms);
2707 copy_rvecn(as_rvec_array(globalState->x.data()), x_pbc, 0, mtop->natoms);
2708 do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
2710 /* Reset pointer to first ED data set which contains the actual ED data */
2712 GMX_ASSERT(nullptr != edi,
2713 "The pointer to the essential dynamics parameters is undefined");
2716 /* Loop over all ED/flooding data sets (usually only one, though) */
2717 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2719 /* For multiple ED groups we use the output frequency that was specified
2720 * in the first set */
2723 edi->outfrq = ed->edpar->outfrq;
2726 /* Extract the initial reference and average positions. When starting
2727 * from .cpt, these have already been read into sref.x_old
2728 * in init_edsamstate() */
2729 if (!EDstate->bFromCpt)
2731 /* If this is the first run (i.e. no checkpoint present) we assume
2732 * that the starting positions give us the correct PBC representation */
2733 for (i = 0; i < edi->sref.nr; i++)
2735 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2738 for (i = 0; i < edi->sav.nr; i++)
2740 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2744 /* Now we have the PBC-correct start positions of the reference and
2745 average structure. We copy that over to dummy arrays on which we
2746 can apply fitting to print out the RMSD. We srenew the memory since
2747 the size of the buffers is likely different for every ED group */
2748 srenew(xfit, edi->sref.nr );
2749 srenew(xstart, edi->sav.nr );
2752 /* Reference indices are the same as average indices */
2753 ref_x_old = edi->sav.x_old;
2757 ref_x_old = edi->sref.x_old;
2759 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2760 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2762 /* Make the fit to the REFERENCE structure, get translation and rotation */
2763 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2765 /* Output how well we fit to the reference at the start */
2766 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2767 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2768 rmsd_from_structure(xfit, &edi->sref));
2769 if (EDstate->nED > 1)
2771 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2773 fprintf(stderr, "\n");
2775 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2776 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2778 /* calculate initial projections */
2779 project(xstart, edi);
2781 /* For the target and origin structure both a reference (fit) and an
2782 * average structure can be provided in make_edi. If both structures
2783 * are the same, make_edi only stores one of them in the .edi file.
2784 * If they differ, first the fit and then the average structure is stored
2785 * in star (or sor), thus the number of entries in star/sor is
2786 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2787 * the size of the average group. */
2789 /* process target structure, if required */
2790 if (edi->star.nr > 0)
2792 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2794 /* get translation & rotation for fit of target structure to reference structure */
2795 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2797 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2798 if (edi->star.nr == edi->sav.nr)
2802 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2804 /* The last sav.nr indices of the target structure correspond to
2805 * the average structure, which must be projected */
2806 avindex = edi->star.nr - edi->sav.nr;
2808 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2812 rad_project(*edi, xstart, &edi->vecs.radcon);
2815 /* process structure that will serve as origin of expansion circle */
2816 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2818 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2821 if (edi->sori.nr > 0)
2823 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2825 /* fit this structure to reference structure */
2826 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2828 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2829 if (edi->sori.nr == edi->sav.nr)
2833 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2835 /* For the projection, we need the last sav.nr indices of sori */
2836 avindex = edi->sori.nr - edi->sav.nr;
2839 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2840 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2841 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2843 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2844 /* Set center of flooding potential to the ORIGIN structure */
2845 rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2846 /* We already know that no (moving) reference position was provided,
2847 * therefore we can overwrite refproj[0]*/
2848 copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
2851 else /* No origin structure given */
2853 rad_project(*edi, xstart, &edi->vecs.radacc);
2854 rad_project(*edi, xstart, &edi->vecs.radfix);
2855 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2857 if (edi->flood.bHarmonic)
2859 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2860 for (i = 0; i < edi->flood.vecs.neig; i++)
2862 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
2867 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2868 /* Set center of flooding potential to the center of the covariance matrix,
2869 * i.e. the average structure, i.e. zero in the projected system */
2870 for (i = 0; i < edi->flood.vecs.neig; i++)
2872 edi->flood.vecs.refproj[i] = 0.0;
2877 /* For convenience, output the center of the flooding potential for the eigenvectors */
2878 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2880 for (i = 0; i < edi->flood.vecs.neig; i++)
2882 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2883 if (edi->flood.bHarmonic)
2885 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.referenceProjectionSlope[i]);
2887 fprintf(stdout, "\n");
2891 /* set starting projections for linsam */
2892 rad_project(*edi, xstart, &edi->vecs.linacc);
2893 rad_project(*edi, xstart, &edi->vecs.linfix);
2895 /* Prepare for the next edi data set: */
2896 edi = edi->next_edi;
2898 /* Cleaning up on the master node: */
2899 if (!EDstate->bFromCpt)
2906 } /* end of MASTER only section */
2910 /* First let everybody know how many ED data sets to expect */
2911 gmx_bcast(sizeof(nED), &nED, cr);
2912 /* Broadcast the essential dynamics / flooding data to all nodes */
2913 broadcast_ed_data(cr, ed, nED);
2917 /* In the single-CPU case, point the local atom numbers pointers to the global
2918 * one, so that we can use the same notation in serial and parallel case: */
2919 /* Loop over all ED data sets (usually only one, though) */
2921 for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
2923 edi->sref.anrs_loc = edi->sref.anrs;
2924 edi->sav.anrs_loc = edi->sav.anrs;
2925 edi->star.anrs_loc = edi->star.anrs;
2926 edi->sori.anrs_loc = edi->sori.anrs;
2927 /* For the same reason as above, make a dummy c_ind array: */
2928 snew(edi->sav.c_ind, edi->sav.nr);
2929 /* Initialize the array */
2930 for (i = 0; i < edi->sav.nr; i++)
2932 edi->sav.c_ind[i] = i;
2934 /* In the general case we will need a different-sized array for the reference indices: */
2937 snew(edi->sref.c_ind, edi->sref.nr);
2938 for (i = 0; i < edi->sref.nr; i++)
2940 edi->sref.c_ind[i] = i;
2943 /* Point to the very same array in case of other structures: */
2944 edi->star.c_ind = edi->sav.c_ind;
2945 edi->sori.c_ind = edi->sav.c_ind;
2946 /* In the serial case, the local number of atoms is the global one: */
2947 edi->sref.nr_loc = edi->sref.nr;
2948 edi->sav.nr_loc = edi->sav.nr;
2949 edi->star.nr_loc = edi->star.nr;
2950 edi->sori.nr_loc = edi->sori.nr;
2952 /* An on we go to the next ED group */
2953 edi = edi->next_edi;
2957 /* Allocate space for ED buffer variables */
2958 /* Again, loop over ED data sets */
2960 for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
2962 /* Allocate space for ED buffer variables */
2963 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
2964 snew(edi->buf->do_edsam, 1);
2966 /* Space for collective ED buffer variables */
2968 /* Collective positions of atoms with the average indices */
2969 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2970 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
2971 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2972 /* Collective positions of atoms with the reference indices */
2975 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2976 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
2977 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2980 /* Get memory for flooding forces */
2981 snew(edi->flood.forces_cartesian, edi->sav.nr);
2984 edi = edi->next_edi;
2987 /* Flush the edo file so that the user can check some things
2988 * when the simulation has started */
2998 void do_edsam(const t_inputrec *ir,
3000 const t_commrec *cr,
3006 int i, edinr, iupdate = 500;
3007 matrix rotmat; /* rotation matrix */
3008 rvec transvec; /* translation vector */
3009 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3010 real dt_1; /* 1/dt */
3011 struct t_do_edsam *buf;
3013 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3014 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3017 /* Check if ED sampling has to be performed */
3018 if (ed->eEDtype == eEDnone)
3023 dt_1 = 1.0/ir->delta_t;
3025 /* Loop over all ED groups (usually one) */
3028 while (edi != nullptr)
3031 if (bNeedDoEdsam(*edi))
3034 buf = edi->buf->do_edsam;
3038 /* initialize radacc radius for slope criterion */
3039 buf->oldrad = calc_radius(edi->vecs.radacc);
3042 /* Copy the positions into buf->xc* arrays and after ED
3043 * feed back corrections to the official positions */
3045 /* Broadcast the ED positions such that every node has all of them
3046 * Every node contributes its local positions xs and stores it in
3047 * the collective buf->xcoll array. Note that for edinr > 1
3048 * xs could already have been modified by an earlier ED */
3050 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3051 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
3053 /* Only assembly reference positions if their indices differ from the average ones */
3056 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3057 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3060 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3061 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3062 * set bUpdateShifts=TRUE in the parallel case. */
3063 buf->bUpdateShifts = FALSE;
3065 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3066 * as well as the indices in edi->sav.anrs */
3068 /* Fit the reference indices to the reference structure */
3071 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3075 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3078 /* Now apply the translation and rotation to the ED structure */
3079 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3081 /* Find out how well we fit to the reference (just for output steps) */
3082 if (do_per_step(step, edi->outfrq) && MASTER(cr))
3086 /* Indices of reference and average structures are identical,
3087 * thus we can calculate the rmsd to SREF using xcoll */
3088 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3092 /* We have to translate & rotate the reference atoms first */
3093 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3094 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3098 /* update radsam references, when required */
3099 if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3101 project(buf->xcoll, edi);
3102 rad_project(*edi, buf->xcoll, &edi->vecs.radacc);
3103 rad_project(*edi, buf->xcoll, &edi->vecs.radfix);
3104 buf->oldrad = -1.e5;
3107 /* update radacc references, when required */
3108 if (do_per_step(step, iupdate) && step >= edi->presteps)
3110 edi->vecs.radacc.radius = calc_radius(edi->vecs.radacc);
3111 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3113 project(buf->xcoll, edi);
3114 rad_project(*edi, buf->xcoll, &edi->vecs.radacc);
3119 buf->oldrad = edi->vecs.radacc.radius;
3123 /* apply the constraints */
3124 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3126 /* ED constraints should be applied already in the first MD step
3127 * (which is step 0), therefore we pass step+1 to the routine */
3128 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3131 /* write to edo, when required */
3132 if (do_per_step(step, edi->outfrq))
3134 project(buf->xcoll, edi);
3135 if (MASTER(cr) && !bSuppress)
3137 write_edo(*edi, ed->edo, rmsdev);
3141 /* Copy back the positions unless monitoring only */
3142 if (ed_constraints(ed->eEDtype, edi))
3144 /* remove fitting */
3145 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3147 /* Copy the ED corrected positions into the coordinate array */
3148 /* Each node copies its local part. In the serial case, nat_loc is the
3149 * total number of ED atoms */
3150 for (i = 0; i < edi->sav.nr_loc; i++)
3152 /* Unshift local ED coordinate and store in x_unsh */
3153 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3154 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3156 /* dx is the ED correction to the positions: */
3157 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3161 /* dv is the ED correction to the velocity: */
3162 svmul(dt_1, dx, dv);
3163 /* apply the velocity correction: */
3164 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3166 /* Finally apply the position correction due to ED: */
3167 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3170 } /* END of if ( bNeedDoEdsam(edi) ) */
3172 /* Prepare for the next ED group */
3173 edi = edi->next_edi;
3175 } /* END of loop over ED groups */