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.
47 #include "gromacs/commandline/filenm.h"
48 #include "gromacs/compat/make_unique.h"
49 #include "gromacs/domdec/domdec_struct.h"
50 #include "gromacs/fileio/gmxfio.h"
51 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/gmxlib/network.h"
53 #include "gromacs/gmxlib/nrnb.h"
54 #include "gromacs/linearalgebra/nrjac.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/utilities.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/math/vectypes.h"
59 #include "gromacs/mdlib/broadcaststructs.h"
60 #include "gromacs/mdlib/constr.h"
61 #include "gromacs/mdlib/groupcoord.h"
62 #include "gromacs/mdlib/mdrun.h"
63 #include "gromacs/mdlib/sim_util.h"
64 #include "gromacs/mdlib/update.h"
65 #include "gromacs/mdtypes/commrec.h"
66 #include "gromacs/mdtypes/edsamhistory.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/md_enums.h"
69 #include "gromacs/mdtypes/observableshistory.h"
70 #include "gromacs/mdtypes/state.h"
71 #include "gromacs/pbcutil/pbc.h"
72 #include "gromacs/topology/mtop_lookup.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/gmxassert.h"
77 #include "gromacs/utility/smalloc.h"
78 #include "gromacs/utility/strconvert.h"
80 /* enum to identify the type of ED: none, normal ED, flooding */
82 eEDnone, eEDedsam, eEDflood, eEDnr
85 /* enum to identify operations on reference, average, origin, target structures */
87 eedREF, eedAV, eedORI, eedTAR, eedNR
93 int neig; /* nr of eigenvectors */
94 int *ieig; /* index nrs of eigenvectors */
95 real *stpsz; /* stepsizes (per eigenvector) */
96 rvec **vec; /* eigenvector components */
97 real *xproj; /* instantaneous x projections */
98 real *fproj; /* instantaneous f projections */
99 real radius; /* instantaneous radius */
100 real *refproj; /* starting or target projections */
106 t_eigvec mon; /* only monitored, no constraints */
107 t_eigvec linfix; /* fixed linear constraints */
108 t_eigvec linacc; /* acceptance linear constraints */
109 t_eigvec radfix; /* fixed radial constraints (exp) */
110 t_eigvec radacc; /* acceptance radial constraints (exp) */
111 t_eigvec radcon; /* acceptance rad. contraction constr. */
118 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
120 gmx_bool bConstForce; /* Do not calculate a flooding potential,
121 instead flood with a constant force */
130 rvec *forces_cartesian;
131 t_eigvec vecs; /* use flooding for these */
132 /* When using flooding as harmonic restraint: The current reference projection
133 * is at each step calculated from the initial initialReferenceProjection and the slope. */
134 real *initialReferenceProjection;
135 real *referenceProjectionSlope;
139 /* This type is for the average, reference, target, and origin structure */
140 typedef struct gmx_edx
142 int nr; /* number of atoms this structure contains */
143 int nr_loc; /* number of atoms on local node */
144 int *anrs; /* atom index numbers */
145 int *anrs_loc; /* local atom index numbers */
146 int nalloc_loc; /* allocation size of anrs_loc */
147 int *c_ind; /* at which position of the whole anrs
148 * array is a local atom?, i.e.
149 * c_ind[0...nr_loc-1] gives the atom index
150 * with respect to the collective
151 * anrs[0...nr-1] array */
152 rvec *x; /* positions for this structure */
153 rvec *x_old; /* Last positions which have the correct PBC
154 representation of the ED group. In
155 combination with keeping track of the
156 shift vectors, the ED group can always
158 real *m; /* masses */
159 real mtot; /* total mass (only used in sref) */
160 real *sqrtm; /* sqrt of the masses used for mass-
161 * weighting of analysis (only used in sav) */
167 int nini; /* total Nr of atoms */
168 gmx_bool fitmas; /* true if trans fit with cm */
169 gmx_bool pcamas; /* true if mass-weighted PCA */
170 int presteps; /* number of steps to run without any
171 * perturbations ... just monitoring */
172 int outfrq; /* freq (in steps) of writing to edo */
173 int maxedsteps; /* max nr of steps per cycle */
175 /* all gmx_edx datasets are copied to all nodes in the parallel case */
176 struct gmx_edx sref; /* reference positions, to these fitting
178 gmx_bool bRefEqAv; /* If true, reference & average indices
179 * are the same. Used for optimization */
180 struct gmx_edx sav; /* average positions */
181 struct gmx_edx star; /* target positions */
182 struct gmx_edx sori; /* origin positions */
184 t_edvecs vecs; /* eigenvectors */
185 real slope; /* minimal slope in acceptance radexp */
187 t_edflood flood; /* parameters especially for flooding */
188 struct t_ed_buffer *buf; /* handle to local buffers */
189 struct edpar *next_edi; /* Pointer to another ED group */
196 int eEDtype = eEDnone; /* Type of ED: see enums above */
197 FILE *edo = nullptr; /* output file pointer */
198 t_edpar *edpar = nullptr;
199 gmx_bool bFirst = false;
201 gmx_edsam::~gmx_edsam()
205 /* edo is opened sometimes with xvgropen, sometimes with
206 * gmx_fio_fopen, so we use the least common denominator for
216 rvec old_transvec, older_transvec, transvec_compact;
217 rvec *xcoll; /* Positions from all nodes, this is the
218 collective set we work on.
219 These are the positions of atoms with
220 average structure indices */
221 rvec *xc_ref; /* same but with reference structure indices */
222 ivec *shifts_xcoll; /* Shifts for xcoll */
223 ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
224 ivec *shifts_xc_ref; /* Shifts for xc_ref */
225 ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
226 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
227 ED shifts for this ED group need to
232 /* definition of ED buffer structure */
235 struct t_fit_to_ref * fit_to_ref;
236 struct t_do_edfit * do_edfit;
237 struct t_do_edsam * do_edsam;
238 struct t_do_radcon * do_radcon;
243 class EssentialDynamics::Impl
246 // TODO: move all fields from gmx_edsam here and remove gmx_edsam
247 gmx_edsam essentialDynamics_;
249 EssentialDynamics::EssentialDynamics() : impl_(new Impl)
252 EssentialDynamics::~EssentialDynamics() = default;
254 gmx_edsam *EssentialDynamics::getLegacyED()
256 return &impl_->essentialDynamics_;
260 /* Function declarations */
261 static void fit_to_reference(rvec *xcoll, rvec transvec, matrix rotmat, t_edpar *edi);
262 static void translate_and_rotate(rvec *x, int nat, rvec transvec, matrix rotmat);
263 static real rmsd_from_structure(rvec *x, struct gmx_edx *s);
264 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms);
265 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam &ed, edsamhistory_t *EDstate);
266 static void init_edsamstate(gmx_edsam * ed, edsamhistory_t *EDstate);
267 static void write_edo_legend(gmx_edsam * ed, int nED, const gmx_output_env_t *oenv);
268 /* End function declarations */
270 /* Define formats for the column width in the output file */
271 const char EDcol_sfmt[] = "%17s";
272 const char EDcol_efmt[] = "%17.5e";
273 const char EDcol_ffmt[] = "%17f";
275 /* Define a formats for reading, otherwise cppcheck complains for scanf without width limits */
276 const char max_ev_fmt_d[] = "%7d"; /* Number of eigenvectors. 9,999,999 should be enough */
277 const char max_ev_fmt_lf[] = "%12lf";
278 const char max_ev_fmt_dlf[] = "%7d%12lf";
279 const char max_ev_fmt_dlflflf[] = "%7d%12lf%12lf%12lf";
280 const char max_ev_fmt_lelele[] = "%12le%12le%12le";
284 /*!\brief Determines whether to perform essential dynamics constraints other than flooding.
285 * \param[in] edi the essential dynamics parameters
286 * \returns true if essential dyanmics constraints need to be performed
288 bool bNeedDoEdsam(const t_edpar &edi)
290 return edi.vecs.mon.neig
291 || edi.vecs.linfix.neig
292 || edi.vecs.linacc.neig
293 || edi.vecs.radfix.neig
294 || edi.vecs.radacc.neig
295 || edi.vecs.radcon.neig;
300 /* Multiple ED groups will be labeled with letters instead of numbers
301 * to avoid confusion with eigenvector indices */
302 static char get_EDgroupChar(int nr_edi, int nED)
310 * nr_edi = 2 -> B ...
312 return 'A' + nr_edi - 1;
317 /*! \brief The mass-weighted inner product of two coordinate vectors.
318 * Does not subtract average positions, projection on single eigenvector is returned
319 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
320 * Average position is subtracted in ed_apply_constraints prior to calling projectx
321 * \param[in] edi Essential dynamics parameters
322 * \param[in] xcoll vector of atom coordinates
323 * \param[in] vec vector of coordinates to project onto
324 * \return mass-weighted projection.
326 real projectx(const t_edpar &edi, rvec *xcoll, rvec *vec)
332 for (i = 0; i < edi.sav.nr; i++)
334 proj += edi.sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
339 /*!\brief Project coordinates onto vector after substracting average position.
340 * projection is stored in vec->refproj which is used for radacc, radfix,
341 * radcon and center of flooding potential.
342 * Average positions are first substracted from x, then added back again.
343 * \param[in] edi essential dynamics parameters with average position
344 * \param[in] x Coordinates to be projected
345 * \param[out] vec eigenvector, radius and refproj are overwritten here
347 void rad_project(const t_edpar &edi, rvec *x, t_eigvec *vec)
352 /* Subtract average positions */
353 for (i = 0; i < edi.sav.nr; i++)
355 rvec_dec(x[i], edi.sav.x[i]);
358 for (i = 0; i < vec->neig; i++)
360 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
361 rad += gmx::square((vec->refproj[i]-vec->xproj[i]));
363 vec->radius = sqrt(rad);
365 /* Add average positions */
366 for (i = 0; i < edi.sav.nr; i++)
368 rvec_inc(x[i], edi.sav.x[i]);
372 /*!\brief Projects coordinates onto eigenvectors and stores result in vec->xproj.
373 * Mass-weighting is applied. Subtracts average positions prior to projection and add
374 * them afterwards to retain the unchanged vector.
375 * \param[in] x The coordinates to project to an eigenvector
376 * \param[in,out] vec The eigenvectors
377 * \param[in] edi essential dynamics parameters holding average structure and masses
379 void project_to_eigvectors(rvec *x, t_eigvec *vec, const t_edpar &edi)
386 /* Subtract average positions */
387 for (int i = 0; i < edi.sav.nr; i++)
389 rvec_dec(x[i], edi.sav.x[i]);
392 for (int i = 0; i < vec->neig; i++)
394 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
397 /* Add average positions */
398 for (int i = 0; i < edi.sav.nr; i++)
400 rvec_inc(x[i], edi.sav.x[i]);
405 /* Project vector x onto all edi->vecs (mon, linfix,...) */
406 static void project(rvec *x, /* positions to project */
407 t_edpar *edi) /* edi data set */
409 /* It is not more work to subtract the average position in every
410 * subroutine again, because these routines are rarely used simultaneously */
411 project_to_eigvectors(x, &edi->vecs.mon, *edi);
412 project_to_eigvectors(x, &edi->vecs.linfix, *edi);
413 project_to_eigvectors(x, &edi->vecs.linacc, *edi);
414 project_to_eigvectors(x, &edi->vecs.radfix, *edi);
415 project_to_eigvectors(x, &edi->vecs.radacc, *edi);
416 project_to_eigvectors(x, &edi->vecs.radcon, *edi);
421 /*!\brief Evaluates the distance from reference to current eigenvector projection.
422 * \param[in] vec eigenvector
425 real calc_radius(const t_eigvec &vec)
429 for (int i = 0; i < vec.neig; i++)
431 rad += gmx::square((vec.refproj[i]-vec.xproj[i]));
434 return rad = sqrt(rad);
440 static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
445 fprintf(out, "#%s positions:\n%d\n", name, s->nr);
451 fprintf(out, "#index, x, y, z");
454 fprintf(out, ", sqrt(m)");
456 for (i = 0; i < s->nr; i++)
458 fprintf(out, "\n%6d %11.6f %11.6f %11.6f", s->anrs[i], s->x[i][XX], s->x[i][YY], s->x[i][ZZ]);
461 fprintf(out, "%9.3f", s->sqrtm[i]);
469 static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
470 const char name[], int length)
475 fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig);
476 /* Dump the data for every eigenvector: */
477 for (i = 0; i < ev->neig; i++)
479 fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
480 ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
481 for (j = 0; j < length; j++)
483 fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
490 static void dump_edi(t_edpar *edpars, const t_commrec *cr, int nr_edi)
496 sprintf(fn, "EDdump_rank%d_edi%d", cr->nodeid, nr_edi);
497 out = gmx_ffopen(fn, "w");
499 fprintf(out, "#NINI\n %d\n#FITMAS\n %s\n#ANALYSIS_MAS\n %s\n",
500 edpars->nini, gmx::boolToString(edpars->fitmas), gmx::boolToString(edpars->pcamas));
501 fprintf(out, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
502 edpars->outfrq, edpars->maxedsteps, edpars->slope);
503 fprintf(out, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n",
504 edpars->presteps, edpars->flood.deltaF0, edpars->flood.tau,
505 edpars->flood.constEfl, edpars->flood.alpha2);
507 /* Dump reference, average, target, origin positions */
508 dump_edi_positions(out, &edpars->sref, "REFERENCE");
509 dump_edi_positions(out, &edpars->sav, "AVERAGE" );
510 dump_edi_positions(out, &edpars->star, "TARGET" );
511 dump_edi_positions(out, &edpars->sori, "ORIGIN" );
513 /* Dump eigenvectors */
514 dump_edi_eigenvecs(out, &edpars->vecs.mon, "MONITORED", edpars->sav.nr);
515 dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX", edpars->sav.nr);
516 dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC", edpars->sav.nr);
517 dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX", edpars->sav.nr);
518 dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC", edpars->sav.nr);
519 dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON", edpars->sav.nr);
521 /* Dump flooding eigenvectors */
522 dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING", edpars->sav.nr);
524 /* Dump ed local buffer */
525 fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit );
526 fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam );
527 fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon );
538 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
540 /* this is a copy of do_fit with some modifications */
541 int c, r, n, j, i, irot;
542 double d[6], xnr, xpc;
547 struct t_do_edfit *loc;
550 if (edi->buf->do_edfit != nullptr)
557 snew(edi->buf->do_edfit, 1);
559 loc = edi->buf->do_edfit;
563 snew(loc->omega, 2*DIM);
564 snew(loc->om, 2*DIM);
565 for (i = 0; i < 2*DIM; i++)
567 snew(loc->omega[i], 2*DIM);
568 snew(loc->om[i], 2*DIM);
572 for (i = 0; (i < 6); i++)
575 for (j = 0; (j < 6); j++)
577 loc->omega[i][j] = 0;
582 /* calculate the matrix U */
584 for (n = 0; (n < natoms); n++)
586 for (c = 0; (c < DIM); c++)
589 for (r = 0; (r < DIM); r++)
597 /* construct loc->omega */
598 /* loc->omega is symmetric -> loc->omega==loc->omega' */
599 for (r = 0; (r < 6); r++)
601 for (c = 0; (c <= r); c++)
603 if ((r >= 3) && (c < 3))
605 loc->omega[r][c] = u[r-3][c];
606 loc->omega[c][r] = u[r-3][c];
610 loc->omega[r][c] = 0;
611 loc->omega[c][r] = 0;
616 /* determine h and k */
617 jacobi(loc->omega, 6, d, loc->om, &irot);
621 fprintf(stderr, "IROT=0\n");
624 index = 0; /* For the compiler only */
626 for (j = 0; (j < 3); j++)
629 for (i = 0; (i < 6); i++)
638 for (i = 0; (i < 3); i++)
640 vh[j][i] = M_SQRT2*loc->om[i][index];
641 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
646 for (c = 0; (c < 3); c++)
648 for (r = 0; (r < 3); r++)
650 R[c][r] = vk[0][r]*vh[0][c]+
657 for (c = 0; (c < 3); c++)
659 for (r = 0; (r < 3); r++)
661 R[c][r] = vk[0][r]*vh[0][c]+
670 static void rmfit(int nat, rvec *xcoll, const rvec transvec, matrix rotmat)
677 * The inverse rotation is described by the transposed rotation matrix */
678 transpose(rotmat, tmat);
679 rotate_x(xcoll, nat, tmat);
681 /* Remove translation */
682 vec[XX] = -transvec[XX];
683 vec[YY] = -transvec[YY];
684 vec[ZZ] = -transvec[ZZ];
685 translate_x(xcoll, nat, vec);
689 /**********************************************************************************
690 ******************** FLOODING ****************************************************
691 **********************************************************************************
693 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
694 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
695 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
697 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
698 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
700 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
701 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
703 do_flood makes a copy of the positions,
704 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
705 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
706 fit. Then do_flood adds these forces to the forcefield-forces
707 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
709 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
710 structure is projected to the system of eigenvectors and then this position in the subspace is used as
711 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
712 i.e. the average structure as given in the make_edi file.
714 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
715 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
716 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
717 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
718 further adaption is applied, Efl will stay constant at zero.
720 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
721 used as spring constants for the harmonic potential.
722 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
723 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
725 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
726 the routine read_edi_file reads all of theses flooding files.
727 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
728 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
729 edi there is no interdependence whatsoever. The forces are added together.
731 To write energies into the .edr file, call the function
732 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
734 get_flood_energies(real Vfl[],int nnames);
737 - 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.
739 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
740 two edsam files from two peptide chains
743 // TODO split this into multiple files
746 /*!\brief Output flooding simulation settings and results to file.
747 * \param[in] edi Essential dynamics input parameters
748 * \param[in] fp output file
749 * \param[in] rmsd rmsd to reference structure
751 void write_edo_flood(const t_edpar &edi, FILE *fp, real rmsd)
753 /* Output how well we fit to the reference structure */
754 fprintf(fp, EDcol_ffmt, rmsd);
756 for (int i = 0; i < edi.flood.vecs.neig; i++)
758 fprintf(fp, EDcol_efmt, edi.flood.vecs.xproj[i]);
760 /* Check whether the reference projection changes with time (this can happen
761 * in case flooding is used as harmonic restraint). If so, output the
762 * current reference projection */
763 if (edi.flood.bHarmonic && edi.flood.referenceProjectionSlope[i] != 0.0)
765 fprintf(fp, EDcol_efmt, edi.flood.vecs.refproj[i]);
768 /* Output Efl if we are doing adaptive flooding */
769 if (0 != edi.flood.tau)
771 fprintf(fp, EDcol_efmt, edi.flood.Efl);
773 fprintf(fp, EDcol_efmt, edi.flood.Vfl);
775 /* Output deltaF if we are doing adaptive flooding */
776 if (0 != edi.flood.tau)
778 fprintf(fp, EDcol_efmt, edi.flood.deltaF);
780 fprintf(fp, EDcol_efmt, edi.flood.vecs.fproj[i]);
786 /* From flood.xproj compute the Vfl(x) at this point */
787 static real flood_energy(t_edpar *edi, int64_t step)
789 /* compute flooding energy Vfl
790 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
791 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
792 it is already computed by make_edi and stored in stpsz[i]
794 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
801 /* Each time this routine is called (i.e. each time step), we add a small
802 * value to the reference projection. This way a harmonic restraint towards
803 * a moving reference is realized. If no value for the additive constant
804 * is provided in the edi file, the reference will not change. */
805 if (edi->flood.bHarmonic)
807 for (i = 0; i < edi->flood.vecs.neig; i++)
809 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i] + step * edi->flood.referenceProjectionSlope[i];
814 /* Compute sum which will be the exponent of the exponential */
815 for (i = 0; i < edi->flood.vecs.neig; i++)
817 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
818 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]);
821 /* Compute the Gauss function*/
822 if (edi->flood.bHarmonic)
824 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
828 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*std::exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
835 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
836 static void flood_forces(t_edpar *edi)
838 /* compute the forces in the subspace of the flooding eigenvectors
839 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
842 real energy = edi->flood.Vfl;
845 if (edi->flood.bHarmonic)
847 for (i = 0; i < edi->flood.vecs.neig; i++)
849 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
854 for (i = 0; i < edi->flood.vecs.neig; i++)
856 /* if Efl is zero the forces are zero if not use the formula */
857 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;
864 /*!\brief Raise forces from subspace into cartesian space.
865 * This function lifts the forces from the subspace to the cartesian space
866 * all the values not contained in the subspace are assumed to be zero and then
867 * a coordinate transformation from eigenvector to cartesian vectors is performed
868 * The nonexistent values don't have to be set to zero explicitly, they would occur
869 * as zero valued summands, hence we just stop to compute this part of the sum.
870 * For every atom we add all the contributions to this atom from all the different eigenvectors.
871 * NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
872 * field forces_cart prior the computation, but we compute the forces separately
873 * to have them accessible for diagnostics
875 * \param[in] edi Essential dynamics input parameters
876 * \param[out] forces_cart The cartesian forces
879 void flood_blowup(const t_edpar &edi, rvec *forces_cart)
881 const real * forces_sub = edi.flood.vecs.fproj;
882 /* Calculate the cartesian forces for the local atoms */
884 /* Clear forces first */
885 for (int j = 0; j < edi.sav.nr_loc; j++)
887 clear_rvec(forces_cart[j]);
890 /* Now compute atomwise */
891 for (int j = 0; j < edi.sav.nr_loc; j++)
893 /* Compute forces_cart[edi.sav.anrs[j]] */
894 for (int eig = 0; eig < edi.flood.vecs.neig; eig++)
897 /* Force vector is force * eigenvector (compute only atom j) */
898 svmul(forces_sub[eig], edi.flood.vecs.vec[eig][edi.sav.c_ind[j]], addedForce);
899 /* Add this vector to the cartesian forces */
900 rvec_inc(forces_cart[j], addedForce);
908 /* Update the values of Efl, deltaF depending on tau and Vfl */
909 static void update_adaption(t_edpar *edi)
911 /* this function updates the parameter Efl and deltaF according to the rules given in
912 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
915 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
917 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
918 /* check if restrain (inverted flooding) -> don't let EFL become positive */
919 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
924 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
929 static void do_single_flood(
937 gmx_bool bNS) /* Are we in a neighbor searching step? */
940 matrix rotmat; /* rotation matrix */
941 matrix tmat; /* inverse rotation */
942 rvec transvec; /* translation vector */
944 struct t_do_edsam *buf;
947 buf = edi->buf->do_edsam;
950 /* Broadcast the positions of the AVERAGE structure such that they are known on
951 * every processor. Each node contributes its local positions x and stores them in
952 * the collective ED array buf->xcoll */
953 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
954 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
956 /* Only assembly REFERENCE positions if their indices differ from the average ones */
959 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
960 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
963 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
964 * We do not need to update the shifts until the next NS step */
965 buf->bUpdateShifts = FALSE;
967 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
968 * as well as the indices in edi->sav.anrs */
970 /* Fit the reference indices to the reference structure */
973 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
977 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
980 /* Now apply the translation and rotation to the ED structure */
981 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
983 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
984 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
986 if (FALSE == edi->flood.bConstForce)
988 /* Compute Vfl(x) from flood.xproj */
989 edi->flood.Vfl = flood_energy(edi, step);
991 update_adaption(edi);
993 /* Compute the flooding forces */
997 /* Translate them into cartesian positions */
998 flood_blowup(*edi, edi->flood.forces_cartesian);
1000 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
1001 /* Each node rotates back its local forces */
1002 transpose(rotmat, tmat);
1003 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
1005 /* Finally add forces to the main force variable */
1006 for (i = 0; i < edi->sav.nr_loc; i++)
1008 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
1011 /* Output is written by the master process */
1012 if (do_per_step(step, edi->outfrq) && MASTER(cr))
1014 /* Output how well we fit to the reference */
1017 /* Indices of reference and average structures are identical,
1018 * thus we can calculate the rmsd to SREF using xcoll */
1019 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
1023 /* We have to translate & rotate the reference atoms first */
1024 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1025 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1028 write_edo_flood(*edi, edo, rmsdev);
1033 /* Main flooding routine, called from do_force */
1034 extern void do_flood(const t_commrec *cr,
1035 const t_inputrec *ir,
1038 const gmx_edsam *ed,
1048 /* Write time to edo, when required. Output the time anyhow since we need
1049 * it in the output file for ED constraints. */
1050 if (MASTER(cr) && do_per_step(step, edi->outfrq))
1052 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
1055 if (ed->eEDtype != eEDflood)
1062 /* Call flooding for one matrix */
1063 if (edi->flood.vecs.neig)
1065 do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS);
1067 edi = edi->next_edi;
1072 /* Called by init_edi, configure some flooding related variables and structures,
1073 * print headers to output files */
1074 static void init_flood(t_edpar *edi, gmx_edsam * ed, real dt)
1079 edi->flood.Efl = edi->flood.constEfl;
1083 if (edi->flood.vecs.neig)
1085 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1086 ed->eEDtype = eEDflood;
1088 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
1090 if (edi->flood.bConstForce)
1092 /* We have used stpsz as a vehicle to carry the fproj values for constant
1093 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1094 * in const force flooding, fproj is never changed. */
1095 for (i = 0; i < edi->flood.vecs.neig; i++)
1097 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1099 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1100 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1108 /*********** Energy book keeping ******/
1109 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1118 srenew(names, count);
1119 sprintf(buf, "Vfl_%d", count);
1120 names[count-1] = gmx_strdup(buf);
1121 actual = actual->next_edi;
1128 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1130 /*fl has to be big enough to capture nnames-many entries*/
1139 Vfl[count-1] = actual->flood.Vfl;
1140 actual = actual->next_edi;
1143 if (nnames != count-1)
1145 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1148 /************* END of FLOODING IMPLEMENTATION ****************************/
1152 /* This function opens the ED input and output files, reads in all datasets it finds
1153 * in the input file, and cross-checks whether the .edi file information is consistent
1154 * with the essential dynamics data found in the checkpoint file (if present).
1155 * gmx make_edi can be used to create an .edi input file.
1157 static std::unique_ptr<gmx::EssentialDynamics> ed_open(
1159 ObservablesHistory *oh,
1160 const char *ediFileName,
1161 const char *edoFileName,
1163 const gmx_output_env_t *oenv,
1164 const t_commrec *cr)
1166 auto edHandle = gmx::compat::make_unique<gmx::EssentialDynamics>();
1167 auto ed = edHandle->getLegacyED();
1170 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1171 ed->eEDtype = eEDedsam;
1177 // If we start from a checkpoint file, we already have an edsamHistory struct
1178 if (oh->edsamHistory == nullptr)
1180 oh->edsamHistory = gmx::compat::make_unique<edsamhistory_t>(edsamhistory_t {});
1182 edsamhistory_t *EDstate = oh->edsamHistory.get();
1184 /* Read the edi input file: */
1185 nED = read_edi_file(ediFileName, ed->edpar, natoms);
1187 /* Make sure the checkpoint was produced in a run using this .edi file */
1188 if (EDstate->bFromCpt)
1190 crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1196 init_edsamstate(ed, EDstate);
1198 /* The master opens the ED output file */
1201 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1205 ed->edo = xvgropen(edoFileName,
1206 "Essential dynamics / flooding output",
1208 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1210 /* Make a descriptive legend */
1211 write_edo_legend(ed, EDstate->nED, oenv);
1218 /* Broadcasts the structure data */
1219 static void bc_ed_positions(const t_commrec *cr, struct gmx_edx *s, int stype)
1221 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1222 snew_bc(cr, s->x, s->nr ); /* Positions */
1223 nblock_bc(cr, s->nr, s->anrs );
1224 nblock_bc(cr, s->nr, s->x );
1226 /* For the average & reference structures we need an array for the collective indices,
1227 * and we need to broadcast the masses as well */
1228 if (stype == eedAV || stype == eedREF)
1230 /* We need these additional variables in the parallel case: */
1231 snew(s->c_ind, s->nr ); /* Collective indices */
1232 /* Local atom indices get assigned in dd_make_local_group_indices.
1233 * There, also memory is allocated */
1234 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1235 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1236 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1239 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1240 if (stype == eedREF)
1242 snew_bc(cr, s->m, s->nr);
1243 nblock_bc(cr, s->nr, s->m);
1246 /* For the average structure we might need the masses for mass-weighting */
1249 snew_bc(cr, s->sqrtm, s->nr);
1250 nblock_bc(cr, s->nr, s->sqrtm);
1251 snew_bc(cr, s->m, s->nr);
1252 nblock_bc(cr, s->nr, s->m);
1257 /* Broadcasts the eigenvector data */
1258 static void bc_ed_vecs(const t_commrec *cr, t_eigvec *ev, int length)
1262 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1263 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1264 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1265 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1266 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1268 nblock_bc(cr, ev->neig, ev->ieig );
1269 nblock_bc(cr, ev->neig, ev->stpsz );
1270 nblock_bc(cr, ev->neig, ev->xproj );
1271 nblock_bc(cr, ev->neig, ev->fproj );
1272 nblock_bc(cr, ev->neig, ev->refproj);
1274 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1275 for (i = 0; i < ev->neig; i++)
1277 snew_bc(cr, ev->vec[i], length);
1278 nblock_bc(cr, length, ev->vec[i]);
1284 /* Broadcasts the ED / flooding data to other nodes
1285 * and allocates memory where needed */
1286 static void broadcast_ed_data(const t_commrec *cr, gmx_edsam * ed, int numedis)
1292 /* Master lets the other nodes know if its ED only or also flooding */
1293 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1295 snew_bc(cr, ed->edpar, 1);
1296 /* Now transfer the ED data set(s) */
1298 for (nr = 0; nr < numedis; nr++)
1300 /* Broadcast a single ED data set */
1303 /* Broadcast positions */
1304 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1305 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1306 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1307 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1309 /* Broadcast eigenvectors */
1310 bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr);
1311 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr);
1312 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr);
1313 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr);
1314 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr);
1315 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr);
1316 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1317 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr);
1319 /* For harmonic restraints the reference projections can change with time */
1320 if (edi->flood.bHarmonic)
1322 snew_bc(cr, edi->flood.initialReferenceProjection, edi->flood.vecs.neig);
1323 snew_bc(cr, edi->flood.referenceProjectionSlope, edi->flood.vecs.neig);
1324 nblock_bc(cr, edi->flood.vecs.neig, edi->flood.initialReferenceProjection);
1325 nblock_bc(cr, edi->flood.vecs.neig, edi->flood.referenceProjectionSlope);
1329 /* Set the pointer to the next ED group */
1332 snew_bc(cr, edi->next_edi, 1);
1333 edi = edi->next_edi;
1339 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1340 static void init_edi(const gmx_mtop_t *mtop, t_edpar *edi)
1343 real totalmass = 0.0;
1346 /* NOTE Init_edi is executed on the master process only
1347 * The initialized data sets are then transmitted to the
1348 * other nodes in broadcast_ed_data */
1350 /* evaluate masses (reference structure) */
1351 snew(edi->sref.m, edi->sref.nr);
1353 for (i = 0; i < edi->sref.nr; i++)
1357 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1361 edi->sref.m[i] = 1.0;
1364 /* Check that every m > 0. Bad things will happen otherwise. */
1365 if (edi->sref.m[i] <= 0.0)
1367 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1368 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1369 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1370 "atoms from the reference structure by creating a proper index group.\n",
1371 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1374 totalmass += edi->sref.m[i];
1376 edi->sref.mtot = totalmass;
1378 /* Masses m and sqrt(m) for the average structure. Note that m
1379 * is needed if forces have to be evaluated in do_edsam */
1380 snew(edi->sav.sqrtm, edi->sav.nr );
1381 snew(edi->sav.m, edi->sav.nr );
1382 for (i = 0; i < edi->sav.nr; i++)
1384 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1387 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1391 edi->sav.sqrtm[i] = 1.0;
1394 /* Check that every m > 0. Bad things will happen otherwise. */
1395 if (edi->sav.sqrtm[i] <= 0.0)
1397 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1398 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1399 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1400 "atoms from the average structure by creating a proper index group.\n",
1401 i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1405 /* put reference structure in origin */
1406 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1410 translate_x(edi->sref.x, edi->sref.nr, com);
1412 /* Init ED buffer */
1417 static void check(const char *line, const char *label)
1419 if (!strstr(line, label))
1421 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1426 static int read_checked_edint(FILE *file, const char *label)
1428 char line[STRLEN+1];
1431 fgets2 (line, STRLEN, file);
1433 fgets2 (line, STRLEN, file);
1434 sscanf (line, max_ev_fmt_d, &idum);
1439 static int read_edint(FILE *file, bool *bEOF)
1441 char line[STRLEN+1];
1446 eof = fgets2 (line, STRLEN, file);
1452 eof = fgets2 (line, STRLEN, file);
1458 sscanf (line, max_ev_fmt_d, &idum);
1464 static real read_checked_edreal(FILE *file, const char *label)
1466 char line[STRLEN+1];
1470 fgets2 (line, STRLEN, file);
1472 fgets2 (line, STRLEN, file);
1473 sscanf (line, max_ev_fmt_lf, &rdum);
1474 return static_cast<real>(rdum); /* always read as double and convert to single */
1478 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1481 char line[STRLEN+1];
1485 for (i = 0; i < number; i++)
1487 fgets2 (line, STRLEN, file);
1488 sscanf (line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1489 anrs[i]--; /* we are reading FORTRAN indices */
1490 for (j = 0; j < 3; j++)
1492 x[i][j] = d[j]; /* always read as double and convert to single */
1499 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1500 * \param[in] in the file to read from
1501 * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1502 * \param[out] vec the eigen-vectors
1503 * \param[in] nEig the number of eigenvectors
1505 void scan_edvec(FILE *in, int numAtoms, rvec ***vec, int nEig)
1508 for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1510 snew((*vec)[iEigenvector], numAtoms);
1511 for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1513 char line[STRLEN+1];
1515 fgets2(line, STRLEN, in);
1516 sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1517 (*vec)[iEigenvector][iAtom][XX] = x;
1518 (*vec)[iEigenvector][iAtom][YY] = y;
1519 (*vec)[iEigenvector][iAtom][ZZ] = z;
1524 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1525 * \param[in] in input file
1526 * \param[in] nrAtoms number of atoms/coordinates
1527 * \param[out] tvec the eigenvector
1529 void read_edvec(FILE *in, int nrAtoms, t_eigvec *tvec)
1531 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1532 if (tvec->neig <= 0)
1537 snew(tvec->ieig, tvec->neig);
1538 snew(tvec->stpsz, tvec->neig);
1539 for (int i = 0; i < tvec->neig; i++)
1541 char line[STRLEN+1];
1542 fgets2(line, STRLEN, in);
1543 int numEigenVectors;
1545 const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1548 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1550 tvec->ieig[i] = numEigenVectors;
1551 tvec->stpsz[i] = stepSize;
1552 } /* end of loop over eigenvectors */
1554 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1556 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1558 * Eigenvector and its intial reference and slope are stored on the
1559 * same line in the .edi format. To avoid re-winding the .edi file,
1560 * the reference eigen-vector and reference data are read in one go.
1562 * \param[in] in input file
1563 * \param[in] nrAtoms number of atoms/coordinates
1564 * \param[out] tvec the eigenvector
1565 * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1566 * \param[out] referenceProjectionSlope The slope of the reference projections.
1568 bool readEdVecWithReferenceProjection(FILE *in, int nrAtoms, t_eigvec *tvec, real ** initialReferenceProjection, real ** referenceProjectionSlope)
1570 bool providesReference = false;
1571 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1572 if (tvec->neig <= 0)
1577 snew(tvec->ieig, tvec->neig);
1578 snew(tvec->stpsz, tvec->neig);
1579 snew(*initialReferenceProjection, tvec->neig);
1580 snew(*referenceProjectionSlope, tvec->neig);
1581 for (int i = 0; i < tvec->neig; i++)
1583 char line[STRLEN+1];
1584 fgets2 (line, STRLEN, in);
1585 int numEigenVectors;
1586 double stepSize = 0.;
1587 double referenceProjection = 0.;
1588 double referenceSlope = 0.;
1589 const int nscan = sscanf(line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize, &referenceProjection, &referenceSlope);
1590 /* Zero out values which were not scanned */
1594 /* Every 4 values read, including reference position */
1595 providesReference = true;
1598 /* A reference position is provided */
1599 providesReference = true;
1600 /* No value for slope, set to 0 */
1601 referenceSlope = 0.0;
1604 /* No values for reference projection and slope, set to 0 */
1605 referenceProjection = 0.0;
1606 referenceSlope = 0.0;
1609 gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1611 (*initialReferenceProjection)[i] = referenceProjection;
1612 (*referenceProjectionSlope)[i] = referenceSlope;
1614 tvec->ieig[i] = numEigenVectors;
1615 tvec->stpsz[i] = stepSize;
1616 } /* end of loop over eigenvectors */
1618 scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1619 return providesReference;
1622 /*!\brief Allocate working memory for the eigenvectors.
1623 * \param[in,out] tvec eigenvector for which memory will be allocated
1625 void setup_edvec(t_eigvec *tvec)
1627 snew(tvec->xproj, tvec->neig);
1628 snew(tvec->fproj, tvec->neig);
1629 snew(tvec->refproj, tvec->neig);
1634 /* Check if the same atom indices are used for reference and average positions */
1635 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1640 /* If the number of atoms differs between the two structures,
1641 * they cannot be identical */
1642 if (sref.nr != sav.nr)
1647 /* Now that we know that both stuctures have the same number of atoms,
1648 * check if also the indices are identical */
1649 for (i = 0; i < sav.nr; i++)
1651 if (sref.anrs[i] != sav.anrs[i])
1656 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1664 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1665 constexpr int c_oldestUnsupportedMagicNumber = 666;
1666 //! Oldest (lowest) magic number indicating supported essential dynamics input
1667 constexpr int c_oldestSupportedMagicNumber = 669;
1668 //! Latest (highest) magic number indicating supported essential dynamics input
1669 constexpr int c_latestSupportedMagicNumber = 670;
1671 /*!\brief Set up essential dynamics work parameters.
1672 * \param[in] edi Essential dynamics working structure.
1674 void setup_edi(t_edpar *edi)
1676 snew(edi->sref.x_old, edi->sref.nr);
1677 edi->sref.sqrtm = nullptr;
1678 snew(edi->sav.x_old, edi->sav.nr);
1679 if (edi->star.nr > 0)
1681 edi->star.sqrtm = nullptr;
1683 if (edi->sori.nr > 0)
1685 edi->sori.sqrtm = nullptr;
1687 setup_edvec(&edi->vecs.linacc);
1688 setup_edvec(&edi->vecs.mon);
1689 setup_edvec(&edi->vecs.linfix);
1690 setup_edvec(&edi->vecs.linacc);
1691 setup_edvec(&edi->vecs.radfix);
1692 setup_edvec(&edi->vecs.radacc);
1693 setup_edvec(&edi->vecs.radcon);
1694 setup_edvec(&edi->flood.vecs);
1697 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1698 * \param[in] magicNumber indicates the essential dynamics input file version
1699 * \returns true if CONST_FORCE_FLOODING is supported
1701 bool ediFileHasConstForceFlooding(int magicNumber)
1703 return magicNumber > c_oldestSupportedMagicNumber;
1706 /*!\brief Reports reading success of the essential dynamics input file magic number.
1707 * \param[in] in pointer to input file
1708 * \param[out] magicNumber determines the edi file version
1709 * \returns true if setting file version from input file was successful.
1711 bool ediFileHasValidDataSet(FILE *in, int * magicNumber)
1713 /* the edi file is not free format, so expect problems if the input is corrupt. */
1714 bool reachedEndOfFile = true;
1715 /* check the magic number */
1716 *magicNumber = read_edint(in, &reachedEndOfFile);
1717 /* Check whether we have reached the end of the input file */
1718 return !reachedEndOfFile;
1721 /*!\brief Trigger fatal error if magic number is unsupported.
1722 * \param[in] magicNumber A magic number read from the edi file
1723 * \param[in] fn name of input file for error reporting
1725 void exitWhenMagicNumberIsInvalid(int magicNumber, const char * fn)
1728 if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1730 if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1732 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1736 gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1741 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1743 * \param[in] in input file
1744 * \param[in] edi essential dynamics parameters
1745 * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1746 * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1747 * \param[in] fn the file name of the input file for error reporting
1749 void read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, bool hasConstForceFlooding, const char *fn)
1752 /* check the number of atoms */
1753 edi->nini = read_edint(in, &bEOF);
1754 if (edi->nini != nr_mdatoms)
1756 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1759 /* Done checking. For the rest we blindly trust the input */
1760 edi->fitmas = read_checked_edint(in, "FITMAS");
1761 edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS");
1762 edi->outfrq = read_checked_edint(in, "OUTFRQ");
1763 edi->maxedsteps = read_checked_edint(in, "MAXLEN");
1764 edi->slope = read_checked_edreal(in, "SLOPECRIT");
1766 edi->presteps = read_checked_edint(in, "PRESTEPS");
1767 edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1768 edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1769 edi->flood.tau = read_checked_edreal(in, "TAU");
1770 edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1771 edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1772 edi->flood.kT = read_checked_edreal(in, "KT");
1773 edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1774 if (hasConstForceFlooding)
1776 edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1780 edi->flood.bConstForce = FALSE;
1782 edi->sref.nr = read_checked_edint(in, "NREF");
1784 /* allocate space for reference positions and read them */
1785 snew(edi->sref.anrs, edi->sref.nr);
1786 snew(edi->sref.x, edi->sref.nr);
1787 read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1789 /* average positions. they define which atoms will be used for ED sampling */
1790 edi->sav.nr = read_checked_edint(in, "NAV");
1791 snew(edi->sav.anrs, edi->sav.nr);
1792 snew(edi->sav.x, edi->sav.nr);
1793 read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1795 /* Check if the same atom indices are used for reference and average positions */
1796 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1800 read_edvec(in, edi->sav.nr, &edi->vecs.mon);
1801 read_edvec(in, edi->sav.nr, &edi->vecs.linfix);
1802 read_edvec(in, edi->sav.nr, &edi->vecs.linacc);
1803 read_edvec(in, edi->sav.nr, &edi->vecs.radfix);
1804 read_edvec(in, edi->sav.nr, &edi->vecs.radacc);
1805 read_edvec(in, edi->sav.nr, &edi->vecs.radcon);
1807 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1808 bool bHaveReference = false;
1809 if (edi->flood.bHarmonic)
1811 bHaveReference = readEdVecWithReferenceProjection(in, edi->sav.nr, &edi->flood.vecs, &edi->flood.initialReferenceProjection, &edi->flood.referenceProjectionSlope);
1815 read_edvec(in, edi->sav.nr, &edi->flood.vecs);
1818 /* target positions */
1819 edi->star.nr = read_edint(in, &bEOF);
1820 if (edi->star.nr > 0)
1822 snew(edi->star.anrs, edi->star.nr);
1823 snew(edi->star.x, edi->star.nr);
1824 read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1827 /* positions defining origin of expansion circle */
1828 edi->sori.nr = read_edint(in, &bEOF);
1829 if (edi->sori.nr > 0)
1833 /* Both an -ori structure and a at least one manual reference point have been
1834 * specified. That's ambiguous and probably not intentional. */
1835 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1836 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1838 snew(edi->sori.anrs, edi->sori.nr);
1839 snew(edi->sori.x, edi->sori.nr);
1840 read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1844 } // anonymous namespace
1846 /* Read in the edi input file. Note that it may contain several ED data sets which were
1847 * achieved by concatenating multiple edi files. The standard case would be a single ED
1848 * data set, though. */
1849 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1852 t_edpar *curr_edi, *last_edi;
1857 /* This routine is executed on the master only */
1859 /* Open the .edi parameter input file */
1860 in = gmx_fio_fopen(fn, "r");
1861 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1863 /* Now read a sequence of ED input parameter sets from the edi file */
1866 int ediFileMagicNumber;
1867 while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1869 exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1870 bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1871 read_edi(in, curr_edi, nr_mdatoms, hasConstForceFlooding, fn);
1872 setup_edi(curr_edi);
1875 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1876 /* We need to allocate space for the data: */
1878 /* Point the 'next_edi' entry to the next edi: */
1879 curr_edi->next_edi = edi_read;
1880 /* Keep the curr_edi pointer for the case that the next group is empty: */
1881 last_edi = curr_edi;
1882 /* Let's prepare to read in the next edi data set: */
1883 curr_edi = edi_read;
1887 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1890 /* Terminate the edi group list with a NULL pointer: */
1891 last_edi->next_edi = nullptr;
1893 fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1895 /* Close the .edi file again */
1902 struct t_fit_to_ref {
1903 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1906 /* Fit the current positions to the reference positions
1907 * Do not actually do the fit, just return rotation and translation.
1908 * Note that the COM of the reference structure was already put into
1909 * the origin by init_edi. */
1910 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1911 rvec transvec, /* The translation vector */
1912 matrix rotmat, /* The rotation matrix */
1913 t_edpar *edi) /* Just needed for do_edfit */
1915 rvec com; /* center of mass */
1917 struct t_fit_to_ref *loc;
1920 /* Allocate memory the first time this routine is called for each edi group */
1921 if (nullptr == edi->buf->fit_to_ref)
1923 snew(edi->buf->fit_to_ref, 1);
1924 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1926 loc = edi->buf->fit_to_ref;
1928 /* We do not touch the original positions but work on a copy. */
1929 for (i = 0; i < edi->sref.nr; i++)
1931 copy_rvec(xcoll[i], loc->xcopy[i]);
1934 /* Calculate the center of mass */
1935 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1937 transvec[XX] = -com[XX];
1938 transvec[YY] = -com[YY];
1939 transvec[ZZ] = -com[ZZ];
1941 /* Subtract the center of mass from the copy */
1942 translate_x(loc->xcopy, edi->sref.nr, transvec);
1944 /* Determine the rotation matrix */
1945 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1949 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1950 int nat, /* How many positions are there? */
1951 rvec transvec, /* The translation vector */
1952 matrix rotmat) /* The rotation matrix */
1955 translate_x(x, nat, transvec);
1958 rotate_x(x, nat, rotmat);
1962 /* Gets the rms deviation of the positions to the structure s */
1963 /* fit_to_structure has to be called before calling this routine! */
1964 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1965 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1971 for (i = 0; i < s->nr; i++)
1973 rmsd += distance2(s->x[i], x[i]);
1976 rmsd /= static_cast<real>(s->nr);
1983 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1988 if (ed->eEDtype != eEDnone)
1990 /* Loop over ED groups */
1994 /* Local atoms of the reference structure (for fitting), need only be assembled
1995 * if their indices differ from the average ones */
1998 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1999 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
2002 /* Local atoms of the average structure (on these ED will be performed) */
2003 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
2004 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
2006 /* Indicate that the ED shift vectors for this structure need to be updated
2007 * at the next call to communicate_group_positions, since obviously we are in a NS step */
2008 edi->buf->do_edsam->bUpdateShifts = TRUE;
2010 /* Set the pointer to the next ED group (if any) */
2011 edi = edi->next_edi;
2017 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
2028 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
2029 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
2030 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
2034 xu[XX] = x[XX]-tx*box[XX][XX];
2035 xu[YY] = x[YY]-ty*box[YY][YY];
2036 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
2042 /*!\brief Apply fixed linear constraints to essential dynamics variable.
2043 * \param[in,out] xcoll The collected coordinates.
2044 * \param[in] edi the essential dynamics parameters
2045 * \param[in] step the current simulation step
2047 void do_linfix(rvec *xcoll, const t_edpar &edi, int64_t step)
2049 /* loop over linfix vectors */
2050 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2052 /* calculate the projection */
2053 real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
2055 /* calculate the correction */
2056 real preFactor = edi.vecs.linfix.refproj[i] + step*edi.vecs.linfix.stpsz[i] - proj;
2058 /* apply the correction */
2059 preFactor /= edi.sav.sqrtm[i];
2060 for (int j = 0; j < edi.sav.nr; j++)
2062 rvec differenceVector;
2063 svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
2064 rvec_inc(xcoll[j], differenceVector);
2069 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
2070 * \param[in,out] xcoll The collected coordinates.
2071 * \param[in] edi the essential dynamics parameters
2073 void do_linacc(rvec *xcoll, t_edpar *edi)
2075 /* loop over linacc vectors */
2076 for (int i = 0; i < edi->vecs.linacc.neig; i++)
2078 /* calculate the projection */
2079 real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
2081 /* calculate the correction */
2082 real preFactor = 0.0;
2083 if (edi->vecs.linacc.stpsz[i] > 0.0)
2085 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
2087 preFactor = edi->vecs.linacc.refproj[i] - proj;
2090 if (edi->vecs.linacc.stpsz[i] < 0.0)
2092 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
2094 preFactor = edi->vecs.linacc.refproj[i] - proj;
2098 /* apply the correction */
2099 preFactor /= edi->sav.sqrtm[i];
2100 for (int j = 0; j < edi->sav.nr; j++)
2102 rvec differenceVector;
2103 svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2104 rvec_inc(xcoll[j], differenceVector);
2107 /* new positions will act as reference */
2108 edi->vecs.linacc.refproj[i] = proj + preFactor;
2114 static void do_radfix(rvec *xcoll, t_edpar *edi)
2117 real *proj, rad = 0.0, ratio;
2121 if (edi->vecs.radfix.neig == 0)
2126 snew(proj, edi->vecs.radfix.neig);
2128 /* loop over radfix vectors */
2129 for (i = 0; i < edi->vecs.radfix.neig; i++)
2131 /* calculate the projections, radius */
2132 proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2133 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2137 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2138 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2140 /* loop over radfix vectors */
2141 for (i = 0; i < edi->vecs.radfix.neig; i++)
2143 proj[i] -= edi->vecs.radfix.refproj[i];
2145 /* apply the correction */
2146 proj[i] /= edi->sav.sqrtm[i];
2148 for (j = 0; j < edi->sav.nr; j++)
2150 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2151 rvec_inc(xcoll[j], vec_dum);
2159 static void do_radacc(rvec *xcoll, t_edpar *edi)
2162 real *proj, rad = 0.0, ratio = 0.0;
2166 if (edi->vecs.radacc.neig == 0)
2171 snew(proj, edi->vecs.radacc.neig);
2173 /* loop over radacc vectors */
2174 for (i = 0; i < edi->vecs.radacc.neig; i++)
2176 /* calculate the projections, radius */
2177 proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2178 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2182 /* only correct when radius decreased */
2183 if (rad < edi->vecs.radacc.radius)
2185 ratio = edi->vecs.radacc.radius/rad - 1.0;
2189 edi->vecs.radacc.radius = rad;
2192 /* loop over radacc vectors */
2193 for (i = 0; i < edi->vecs.radacc.neig; i++)
2195 proj[i] -= edi->vecs.radacc.refproj[i];
2197 /* apply the correction */
2198 proj[i] /= edi->sav.sqrtm[i];
2200 for (j = 0; j < edi->sav.nr; j++)
2202 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2203 rvec_inc(xcoll[j], vec_dum);
2210 struct t_do_radcon {
2214 static void do_radcon(rvec *xcoll, t_edpar *edi)
2217 real rad = 0.0, ratio = 0.0;
2218 struct t_do_radcon *loc;
2223 if (edi->buf->do_radcon != nullptr)
2230 snew(edi->buf->do_radcon, 1);
2232 loc = edi->buf->do_radcon;
2234 if (edi->vecs.radcon.neig == 0)
2241 snew(loc->proj, edi->vecs.radcon.neig);
2244 /* loop over radcon vectors */
2245 for (i = 0; i < edi->vecs.radcon.neig; i++)
2247 /* calculate the projections, radius */
2248 loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2249 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2252 /* only correct when radius increased */
2253 if (rad > edi->vecs.radcon.radius)
2255 ratio = edi->vecs.radcon.radius/rad - 1.0;
2257 /* loop over radcon vectors */
2258 for (i = 0; i < edi->vecs.radcon.neig; i++)
2260 /* apply the correction */
2261 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2262 loc->proj[i] /= edi->sav.sqrtm[i];
2263 loc->proj[i] *= ratio;
2265 for (j = 0; j < edi->sav.nr; j++)
2267 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2268 rvec_inc(xcoll[j], vec_dum);
2275 edi->vecs.radcon.radius = rad;
2281 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, int64_t step)
2286 /* subtract the average positions */
2287 for (i = 0; i < edi->sav.nr; i++)
2289 rvec_dec(xcoll[i], edi->sav.x[i]);
2292 /* apply the constraints */
2295 do_linfix(xcoll, *edi, step);
2297 do_linacc(xcoll, edi);
2300 do_radfix(xcoll, edi);
2302 do_radacc(xcoll, edi);
2303 do_radcon(xcoll, edi);
2305 /* add back the average positions */
2306 for (i = 0; i < edi->sav.nr; i++)
2308 rvec_inc(xcoll[i], edi->sav.x[i]);
2315 /*!\brief Write out the projections onto the eigenvectors.
2316 * The order of output corresponds to ed_output_legend().
2317 * \param[in] edi The essential dyanmics parameters
2318 * \param[in] fp The output file
2319 * \param[in] rmsd the rmsd to the reference structure
2321 void write_edo(const t_edpar &edi, FILE *fp, real rmsd)
2323 /* Output how well we fit to the reference structure */
2324 fprintf(fp, EDcol_ffmt, rmsd);
2326 for (int i = 0; i < edi.vecs.mon.neig; i++)
2328 fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2331 for (int i = 0; i < edi.vecs.linfix.neig; i++)
2333 fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2336 for (int i = 0; i < edi.vecs.linacc.neig; i++)
2338 fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2341 for (int i = 0; i < edi.vecs.radfix.neig; i++)
2343 fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2345 if (edi.vecs.radfix.neig)
2347 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2350 for (int i = 0; i < edi.vecs.radacc.neig; i++)
2352 fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2354 if (edi.vecs.radacc.neig)
2356 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2359 for (int i = 0; i < edi.vecs.radcon.neig; i++)
2361 fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2363 if (edi.vecs.radcon.neig)
2365 fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2371 /* Returns if any constraints are switched on */
2372 static int ed_constraints(int edtype, t_edpar *edi)
2374 if (edtype == eEDedsam || edtype == eEDflood)
2376 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2377 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2378 edi->vecs.radcon.neig);
2384 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for flooding/
2385 * umbrella sampling simulations. */
2386 static void copyEvecReference(t_eigvec* floodvecs, real * initialReferenceProjection)
2391 if (nullptr == initialReferenceProjection)
2393 snew(initialReferenceProjection, floodvecs->neig);
2396 for (i = 0; i < floodvecs->neig; i++)
2398 initialReferenceProjection[i] = floodvecs->refproj[i];
2403 /* Call on MASTER only. Check whether the essential dynamics / flooding
2404 * groups of the checkpoint file are consistent with the provided .edi file. */
2405 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam &ed, edsamhistory_t *EDstate)
2407 t_edpar *edi = nullptr; /* points to a single edi data set */
2411 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2413 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2414 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2415 "it must also continue with/without ED constraints when checkpointing.\n"
2416 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2417 "from without a checkpoint.\n");
2422 while (edi != nullptr)
2424 /* Check number of atoms in the reference and average structures */
2425 if (EDstate->nref[edinum] != edi->sref.nr)
2427 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2428 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2429 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2431 if (EDstate->nav[edinum] != edi->sav.nr)
2433 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2434 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2435 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2437 edi = edi->next_edi;
2441 if (edinum != EDstate->nED)
2443 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2444 "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2445 "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2450 /* The edsamstate struct stores the information we need to make the ED group
2451 * whole again after restarts from a checkpoint file. Here we do the following:
2452 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2453 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2454 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2455 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2456 * all ED structures at the last time step. */
2457 static void init_edsamstate(gmx_edsam * ed, edsamhistory_t *EDstate)
2463 snew(EDstate->old_sref_p, EDstate->nED);
2464 snew(EDstate->old_sav_p, EDstate->nED);
2466 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2467 if (!EDstate->bFromCpt)
2469 snew(EDstate->nref, EDstate->nED);
2470 snew(EDstate->nav, EDstate->nED);
2473 /* Loop over all ED/flooding data sets (usually only one, though) */
2475 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2477 /* We always need the last reference and average positions such that
2478 * in the next time step we can make the ED group whole again
2479 * if the atoms do not have the correct PBC representation */
2480 if (EDstate->bFromCpt)
2482 /* Copy the last whole positions of reference and average group from .cpt */
2483 for (i = 0; i < edi->sref.nr; i++)
2485 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2487 for (i = 0; i < edi->sav.nr; i++)
2489 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2494 EDstate->nref[nr_edi-1] = edi->sref.nr;
2495 EDstate->nav [nr_edi-1] = edi->sav.nr;
2498 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2499 EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2500 EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2502 edi = edi->next_edi;
2507 /* Adds 'buf' to 'str' */
2508 static void add_to_string(char **str, char *buf)
2513 len = strlen(*str) + strlen(buf) + 1;
2519 static void add_to_string_aligned(char **str, char *buf)
2521 char buf_aligned[STRLEN];
2523 sprintf(buf_aligned, EDcol_sfmt, buf);
2524 add_to_string(str, buf_aligned);
2528 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, const char *value, const char *unit, char EDgroupchar)
2530 char tmp[STRLEN], tmp2[STRLEN];
2533 sprintf(tmp, "%c %s", EDgroupchar, value);
2534 add_to_string_aligned(LegendStr, tmp);
2535 sprintf(tmp2, "%s (%s)", tmp, unit);
2536 (*setname)[*nsets] = gmx_strdup(tmp2);
2541 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2547 for (i = 0; i < evec->neig; i++)
2549 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2550 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2555 /* Makes a legend for the xvg output file. Call on MASTER only! */
2556 static void write_edo_legend(gmx_edsam * ed, int nED, const gmx_output_env_t *oenv)
2558 t_edpar *edi = nullptr;
2560 int nr_edi, nsets, n_flood, n_edsam;
2561 const char **setname;
2563 char *LegendStr = nullptr;
2568 fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2570 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2572 fprintf(ed->edo, "#\n");
2573 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2574 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2575 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2576 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2577 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2578 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2579 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2580 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2581 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2583 if (edi->flood.vecs.neig)
2585 /* If in any of the groups we find a flooding vector, flooding is turned on */
2586 ed->eEDtype = eEDflood;
2588 /* Print what flavor of flooding we will do */
2589 if (0 == edi->flood.tau) /* constant flooding strength */
2591 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2592 if (edi->flood.bHarmonic)
2594 fprintf(ed->edo, ", harmonic");
2597 else /* adaptive flooding */
2599 fprintf(ed->edo, ", adaptive");
2602 fprintf(ed->edo, "\n");
2604 edi = edi->next_edi;
2607 /* Print a nice legend */
2609 LegendStr[0] = '\0';
2610 sprintf(buf, "# %6s", "time");
2611 add_to_string(&LegendStr, buf);
2613 /* Calculate the maximum number of columns we could end up with */
2616 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2618 nsets += 5 +edi->vecs.mon.neig
2619 +edi->vecs.linfix.neig
2620 +edi->vecs.linacc.neig
2621 +edi->vecs.radfix.neig
2622 +edi->vecs.radacc.neig
2623 +edi->vecs.radcon.neig
2624 + 6*edi->flood.vecs.neig;
2625 edi = edi->next_edi;
2627 snew(setname, nsets);
2629 /* In the mdrun time step in a first function call (do_flood()) the flooding
2630 * forces are calculated and in a second function call (do_edsam()) the
2631 * ED constraints. To get a corresponding legend, we need to loop twice
2632 * over the edi groups and output first the flooding, then the ED part */
2634 /* The flooding-related legend entries, if flooding is done */
2636 if (eEDflood == ed->eEDtype)
2639 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2641 /* Always write out the projection on the flooding EVs. Of course, this can also
2642 * be achieved with the monitoring option in do_edsam() (if switched on by the
2643 * user), but in that case the positions need to be communicated in do_edsam(),
2644 * which is not necessary when doing flooding only. */
2645 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2647 for (i = 0; i < edi->flood.vecs.neig; i++)
2649 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2650 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2652 /* Output the current reference projection if it changes with time;
2653 * this can happen when flooding is used as harmonic restraint */
2654 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2656 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2657 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2660 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2661 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2663 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2664 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2667 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2668 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2670 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2672 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2673 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2676 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2677 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2680 edi = edi->next_edi;
2681 } /* End of flooding-related legend entries */
2685 /* Now the ED-related entries, if essential dynamics is done */
2687 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2689 if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2691 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2693 /* Essential dynamics, projections on eigenvectors */
2694 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2695 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2696 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2697 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2698 if (edi->vecs.radfix.neig)
2700 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2702 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2703 if (edi->vecs.radacc.neig)
2705 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2707 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2708 if (edi->vecs.radcon.neig)
2710 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2713 edi = edi->next_edi;
2714 } /* end of 'pure' essential dynamics legend entries */
2715 n_edsam = nsets - n_flood;
2717 xvgr_legend(ed->edo, nsets, setname, oenv);
2720 fprintf(ed->edo, "#\n"
2721 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2722 n_flood, 1 == n_flood ? "" : "s",
2723 n_edsam, 1 == n_edsam ? "" : "s");
2724 fprintf(ed->edo, "%s", LegendStr);
2731 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2732 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2733 std::unique_ptr<gmx::EssentialDynamics> init_edsam(
2734 const char *ediFileName,
2735 const char *edoFileName,
2736 const gmx_mtop_t *mtop,
2737 const t_inputrec *ir,
2738 const t_commrec *cr,
2739 gmx::Constraints *constr,
2740 const t_state *globalState,
2741 ObservablesHistory *oh,
2742 const gmx_output_env_t *oenv,
2745 t_edpar *edi = nullptr; /* points to a single edi data set */
2747 int nED = 0; /* Number of ED data sets */
2748 rvec *x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2749 rvec *xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2750 rvec fit_transvec; /* translation ... */
2751 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2752 rvec *ref_x_old = nullptr; /* helper pointer */
2757 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2759 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName) )
2761 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or flooding\n"
2762 "simulation. Please also set the .edi file on the command line with -ei.\n");
2766 /* Open input and output files, allocate space for ED data structure */
2767 auto edHandle = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr);
2768 auto ed = edHandle->getLegacyED();
2769 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2770 constr->saveEdsamPointer(ed);
2772 /* Needed for initializing radacc radius in do_edsam */
2775 /* The input file is read by the master and the edi structures are
2776 * initialized here. Input is stored in ed->edpar. Then the edi
2777 * structures are transferred to the other nodes */
2780 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2781 * flooding vector, Essential dynamics can be applied to more than one structure
2782 * as well, but will be done in the order given in the edi file, so
2783 * expect different results for different order of edi file concatenation! */
2785 while (edi != nullptr)
2787 init_edi(mtop, edi);
2788 init_flood(edi, ed, ir->delta_t);
2789 edi = edi->next_edi;
2793 /* The master does the work here. The other nodes get the positions
2794 * not before dd_partition_system which is called after init_edsam */
2797 edsamhistory_t *EDstate = oh->edsamHistory.get();
2799 if (!EDstate->bFromCpt)
2801 /* Remove PBC, make molecule(s) subject to ED whole. */
2802 snew(x_pbc, mtop->natoms);
2803 copy_rvecn(as_rvec_array(globalState->x.data()), x_pbc, 0, mtop->natoms);
2804 do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
2806 /* Reset pointer to first ED data set which contains the actual ED data */
2808 GMX_ASSERT(nullptr != edi,
2809 "The pointer to the essential dynamics parameters is undefined");
2812 /* Loop over all ED/flooding data sets (usually only one, though) */
2813 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2815 /* For multiple ED groups we use the output frequency that was specified
2816 * in the first set */
2819 edi->outfrq = ed->edpar->outfrq;
2822 /* Extract the initial reference and average positions. When starting
2823 * from .cpt, these have already been read into sref.x_old
2824 * in init_edsamstate() */
2825 if (!EDstate->bFromCpt)
2827 /* If this is the first run (i.e. no checkpoint present) we assume
2828 * that the starting positions give us the correct PBC representation */
2829 for (i = 0; i < edi->sref.nr; i++)
2831 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2834 for (i = 0; i < edi->sav.nr; i++)
2836 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2840 /* Now we have the PBC-correct start positions of the reference and
2841 average structure. We copy that over to dummy arrays on which we
2842 can apply fitting to print out the RMSD. We srenew the memory since
2843 the size of the buffers is likely different for every ED group */
2844 srenew(xfit, edi->sref.nr );
2845 srenew(xstart, edi->sav.nr );
2848 /* Reference indices are the same as average indices */
2849 ref_x_old = edi->sav.x_old;
2853 ref_x_old = edi->sref.x_old;
2855 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2856 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2858 /* Make the fit to the REFERENCE structure, get translation and rotation */
2859 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2861 /* Output how well we fit to the reference at the start */
2862 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2863 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2864 rmsd_from_structure(xfit, &edi->sref));
2865 if (EDstate->nED > 1)
2867 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2869 fprintf(stderr, "\n");
2871 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2872 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2874 /* calculate initial projections */
2875 project(xstart, edi);
2877 /* For the target and origin structure both a reference (fit) and an
2878 * average structure can be provided in make_edi. If both structures
2879 * are the same, make_edi only stores one of them in the .edi file.
2880 * If they differ, first the fit and then the average structure is stored
2881 * in star (or sor), thus the number of entries in star/sor is
2882 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2883 * the size of the average group. */
2885 /* process target structure, if required */
2886 if (edi->star.nr > 0)
2888 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2890 /* get translation & rotation for fit of target structure to reference structure */
2891 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2893 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2894 if (edi->star.nr == edi->sav.nr)
2898 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2900 /* The last sav.nr indices of the target structure correspond to
2901 * the average structure, which must be projected */
2902 avindex = edi->star.nr - edi->sav.nr;
2904 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2908 rad_project(*edi, xstart, &edi->vecs.radcon);
2911 /* process structure that will serve as origin of expansion circle */
2912 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2914 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2917 if (edi->sori.nr > 0)
2919 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2921 /* fit this structure to reference structure */
2922 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2924 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2925 if (edi->sori.nr == edi->sav.nr)
2929 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2931 /* For the projection, we need the last sav.nr indices of sori */
2932 avindex = edi->sori.nr - edi->sav.nr;
2935 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2936 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2937 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2939 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2940 /* Set center of flooding potential to the ORIGIN structure */
2941 rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2942 /* We already know that no (moving) reference position was provided,
2943 * therefore we can overwrite refproj[0]*/
2944 copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
2947 else /* No origin structure given */
2949 rad_project(*edi, xstart, &edi->vecs.radacc);
2950 rad_project(*edi, xstart, &edi->vecs.radfix);
2951 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2953 if (edi->flood.bHarmonic)
2955 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2956 for (i = 0; i < edi->flood.vecs.neig; i++)
2958 edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
2963 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2964 /* Set center of flooding potential to the center of the covariance matrix,
2965 * i.e. the average structure, i.e. zero in the projected system */
2966 for (i = 0; i < edi->flood.vecs.neig; i++)
2968 edi->flood.vecs.refproj[i] = 0.0;
2973 /* For convenience, output the center of the flooding potential for the eigenvectors */
2974 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2976 for (i = 0; i < edi->flood.vecs.neig; i++)
2978 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2979 if (edi->flood.bHarmonic)
2981 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.referenceProjectionSlope[i]);
2983 fprintf(stdout, "\n");
2987 /* set starting projections for linsam */
2988 rad_project(*edi, xstart, &edi->vecs.linacc);
2989 rad_project(*edi, xstart, &edi->vecs.linfix);
2991 /* Prepare for the next edi data set: */
2992 edi = edi->next_edi;
2994 /* Cleaning up on the master node: */
2995 if (!EDstate->bFromCpt)
3002 } /* end of MASTER only section */
3006 /* First let everybody know how many ED data sets to expect */
3007 gmx_bcast(sizeof(nED), &nED, cr);
3008 /* Broadcast the essential dynamics / flooding data to all nodes */
3009 broadcast_ed_data(cr, ed, nED);
3013 /* In the single-CPU case, point the local atom numbers pointers to the global
3014 * one, so that we can use the same notation in serial and parallel case: */
3015 /* Loop over all ED data sets (usually only one, though) */
3017 for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
3019 edi->sref.anrs_loc = edi->sref.anrs;
3020 edi->sav.anrs_loc = edi->sav.anrs;
3021 edi->star.anrs_loc = edi->star.anrs;
3022 edi->sori.anrs_loc = edi->sori.anrs;
3023 /* For the same reason as above, make a dummy c_ind array: */
3024 snew(edi->sav.c_ind, edi->sav.nr);
3025 /* Initialize the array */
3026 for (i = 0; i < edi->sav.nr; i++)
3028 edi->sav.c_ind[i] = i;
3030 /* In the general case we will need a different-sized array for the reference indices: */
3033 snew(edi->sref.c_ind, edi->sref.nr);
3034 for (i = 0; i < edi->sref.nr; i++)
3036 edi->sref.c_ind[i] = i;
3039 /* Point to the very same array in case of other structures: */
3040 edi->star.c_ind = edi->sav.c_ind;
3041 edi->sori.c_ind = edi->sav.c_ind;
3042 /* In the serial case, the local number of atoms is the global one: */
3043 edi->sref.nr_loc = edi->sref.nr;
3044 edi->sav.nr_loc = edi->sav.nr;
3045 edi->star.nr_loc = edi->star.nr;
3046 edi->sori.nr_loc = edi->sori.nr;
3048 /* An on we go to the next ED group */
3049 edi = edi->next_edi;
3053 /* Allocate space for ED buffer variables */
3054 /* Again, loop over ED data sets */
3056 for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
3058 /* Allocate space for ED buffer variables */
3059 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
3060 snew(edi->buf->do_edsam, 1);
3062 /* Space for collective ED buffer variables */
3064 /* Collective positions of atoms with the average indices */
3065 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
3066 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
3067 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
3068 /* Collective positions of atoms with the reference indices */
3071 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
3072 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
3073 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
3076 /* Get memory for flooding forces */
3077 snew(edi->flood.forces_cartesian, edi->sav.nr);
3081 /* Dump it all into one file per process */
3082 dump_edi(edi, cr, nr_edi);
3086 edi = edi->next_edi;
3089 /* Flush the edo file so that the user can check some things
3090 * when the simulation has started */
3100 void do_edsam(const t_inputrec *ir,
3102 const t_commrec *cr,
3108 int i, edinr, iupdate = 500;
3109 matrix rotmat; /* rotation matrix */
3110 rvec transvec; /* translation vector */
3111 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3112 real dt_1; /* 1/dt */
3113 struct t_do_edsam *buf;
3115 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3116 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3119 /* Check if ED sampling has to be performed */
3120 if (ed->eEDtype == eEDnone)
3125 dt_1 = 1.0/ir->delta_t;
3127 /* Loop over all ED groups (usually one) */
3130 while (edi != nullptr)
3133 if (bNeedDoEdsam(*edi))
3136 buf = edi->buf->do_edsam;
3140 /* initialize radacc radius for slope criterion */
3141 buf->oldrad = calc_radius(edi->vecs.radacc);
3144 /* Copy the positions into buf->xc* arrays and after ED
3145 * feed back corrections to the official positions */
3147 /* Broadcast the ED positions such that every node has all of them
3148 * Every node contributes its local positions xs and stores it in
3149 * the collective buf->xcoll array. Note that for edinr > 1
3150 * xs could already have been modified by an earlier ED */
3152 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3153 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
3155 /* Only assembly reference positions if their indices differ from the average ones */
3158 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3159 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3162 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3163 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3164 * set bUpdateShifts=TRUE in the parallel case. */
3165 buf->bUpdateShifts = FALSE;
3167 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3168 * as well as the indices in edi->sav.anrs */
3170 /* Fit the reference indices to the reference structure */
3173 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3177 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3180 /* Now apply the translation and rotation to the ED structure */
3181 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3183 /* Find out how well we fit to the reference (just for output steps) */
3184 if (do_per_step(step, edi->outfrq) && MASTER(cr))
3188 /* Indices of reference and average structures are identical,
3189 * thus we can calculate the rmsd to SREF using xcoll */
3190 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3194 /* We have to translate & rotate the reference atoms first */
3195 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3196 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3200 /* update radsam references, when required */
3201 if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3203 project(buf->xcoll, edi);
3204 rad_project(*edi, buf->xcoll, &edi->vecs.radacc);
3205 rad_project(*edi, buf->xcoll, &edi->vecs.radfix);
3206 buf->oldrad = -1.e5;
3209 /* update radacc references, when required */
3210 if (do_per_step(step, iupdate) && step >= edi->presteps)
3212 edi->vecs.radacc.radius = calc_radius(edi->vecs.radacc);
3213 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3215 project(buf->xcoll, edi);
3216 rad_project(*edi, buf->xcoll, &edi->vecs.radacc);
3221 buf->oldrad = edi->vecs.radacc.radius;
3225 /* apply the constraints */
3226 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3228 /* ED constraints should be applied already in the first MD step
3229 * (which is step 0), therefore we pass step+1 to the routine */
3230 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3233 /* write to edo, when required */
3234 if (do_per_step(step, edi->outfrq))
3236 project(buf->xcoll, edi);
3237 if (MASTER(cr) && !bSuppress)
3239 write_edo(*edi, ed->edo, rmsdev);
3243 /* Copy back the positions unless monitoring only */
3244 if (ed_constraints(ed->eEDtype, edi))
3246 /* remove fitting */
3247 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3249 /* Copy the ED corrected positions into the coordinate array */
3250 /* Each node copies its local part. In the serial case, nat_loc is the
3251 * total number of ED atoms */
3252 for (i = 0; i < edi->sav.nr_loc; i++)
3254 /* Unshift local ED coordinate and store in x_unsh */
3255 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3256 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3258 /* dx is the ED correction to the positions: */
3259 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3263 /* dv is the ED correction to the velocity: */
3264 svmul(dt_1, dx, dv);
3265 /* apply the velocity correction: */
3266 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3268 /* Finally apply the position correction due to ED: */
3269 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3272 } /* END of if ( bNeedDoEdsam(edi) ) */
3274 /* Prepare for the next ED group */
3275 edi = edi->next_edi;
3277 } /* END of loop over ED groups */