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 */
101 /* When using flooding as harmonic restraint: The current reference projection
102 * is at each step calculated from the initial refproj0 and the slope. */
103 real *refproj0, *refprojslope;
109 t_eigvec mon; /* only monitored, no constraints */
110 t_eigvec linfix; /* fixed linear constraints */
111 t_eigvec linacc; /* acceptance linear constraints */
112 t_eigvec radfix; /* fixed radial constraints (exp) */
113 t_eigvec radacc; /* acceptance radial constraints (exp) */
114 t_eigvec radcon; /* acceptance rad. contraction constr. */
121 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
123 gmx_bool bConstForce; /* Do not calculate a flooding potential,
124 instead flood with a constant force */
133 rvec *forces_cartesian;
134 t_eigvec vecs; /* use flooding for these */
138 /* This type is for the average, reference, target, and origin structure */
139 typedef struct gmx_edx
141 int nr; /* number of atoms this structure contains */
142 int nr_loc; /* number of atoms on local node */
143 int *anrs; /* atom index numbers */
144 int *anrs_loc; /* local atom index numbers */
145 int nalloc_loc; /* allocation size of anrs_loc */
146 int *c_ind; /* at which position of the whole anrs
147 * array is a local atom?, i.e.
148 * c_ind[0...nr_loc-1] gives the atom index
149 * with respect to the collective
150 * anrs[0...nr-1] array */
151 rvec *x; /* positions for this structure */
152 rvec *x_old; /* Last positions which have the correct PBC
153 representation of the ED group. In
154 combination with keeping track of the
155 shift vectors, the ED group can always
157 real *m; /* masses */
158 real mtot; /* total mass (only used in sref) */
159 real *sqrtm; /* sqrt of the masses used for mass-
160 * weighting of analysis (only used in sav) */
166 int nini; /* total Nr of atoms */
167 gmx_bool fitmas; /* true if trans fit with cm */
168 gmx_bool pcamas; /* true if mass-weighted PCA */
169 int presteps; /* number of steps to run without any
170 * perturbations ... just monitoring */
171 int outfrq; /* freq (in steps) of writing to edo */
172 int maxedsteps; /* max nr of steps per cycle */
174 /* all gmx_edx datasets are copied to all nodes in the parallel case */
175 struct gmx_edx sref; /* reference positions, to these fitting
177 gmx_bool bRefEqAv; /* If true, reference & average indices
178 * are the same. Used for optimization */
179 struct gmx_edx sav; /* average positions */
180 struct gmx_edx star; /* target positions */
181 struct gmx_edx sori; /* origin positions */
183 t_edvecs vecs; /* eigenvectors */
184 real slope; /* minimal slope in acceptance radexp */
186 t_edflood flood; /* parameters especially for flooding */
187 struct t_ed_buffer *buf; /* handle to local buffers */
188 struct edpar *next_edi; /* Pointer to another ED group */
192 typedef struct gmx_edsam
194 int eEDtype; /* Type of ED: see enums above */
195 FILE *edo; /* output file pointer */
205 rvec old_transvec, older_transvec, transvec_compact;
206 rvec *xcoll; /* Positions from all nodes, this is the
207 collective set we work on.
208 These are the positions of atoms with
209 average structure indices */
210 rvec *xc_ref; /* same but with reference structure indices */
211 ivec *shifts_xcoll; /* Shifts for xcoll */
212 ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
213 ivec *shifts_xc_ref; /* Shifts for xc_ref */
214 ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
215 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
216 ED shifts for this ED group need to
221 /* definition of ED buffer structure */
224 struct t_fit_to_ref * fit_to_ref;
225 struct t_do_edfit * do_edfit;
226 struct t_do_edsam * do_edsam;
227 struct t_do_radcon * do_radcon;
231 /* Function declarations */
232 static void fit_to_reference(rvec *xcoll, rvec transvec, matrix rotmat, t_edpar *edi);
233 static void translate_and_rotate(rvec *x, int nat, rvec transvec, matrix rotmat);
234 static real rmsd_from_structure(rvec *x, struct gmx_edx *s);
235 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms);
236 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamhistory_t *EDstate);
237 static void init_edsamstate(gmx_edsam_t ed, edsamhistory_t *EDstate);
238 static void write_edo_legend(gmx_edsam_t ed, int nED, const gmx_output_env_t *oenv);
239 /* End function declarations */
241 /* Define formats for the column width in the output file */
242 const char EDcol_sfmt[] = "%17s";
243 const char EDcol_efmt[] = "%17.5e";
244 const char EDcol_ffmt[] = "%17f";
246 /* Define a formats for reading, otherwise cppcheck complains for scanf without width limits */
247 const char max_ev_fmt_d[] = "%7d"; /* Number of eigenvectors. 9,999,999 should be enough */
248 const char max_ev_fmt_lf[] = "%12lf";
249 const char max_ev_fmt_dlf[] = "%7d%12lf";
250 const char max_ev_fmt_dlflflf[] = "%7d%12lf%12lf%12lf";
251 const char max_ev_fmt_lelele[] = "%12le%12le%12le";
253 /* Do we have to perform essential dynamics constraints or possibly only flooding
254 * for any of the ED groups? */
255 static gmx_bool bNeedDoEdsam(t_edpar *edi)
257 return edi->vecs.mon.neig
258 || edi->vecs.linfix.neig
259 || edi->vecs.linacc.neig
260 || edi->vecs.radfix.neig
261 || edi->vecs.radacc.neig
262 || edi->vecs.radcon.neig;
266 /* Multiple ED groups will be labeled with letters instead of numbers
267 * to avoid confusion with eigenvector indices */
268 static char get_EDgroupChar(int nr_edi, int nED)
276 * nr_edi = 2 -> B ...
278 return 'A' + nr_edi - 1;
282 /* Does not subtract average positions, projection on single eigenvector is returned
283 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
284 * Average position is subtracted in ed_apply_constraints prior to calling projectx
286 static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec)
292 for (i = 0; i < edi->sav.nr; i++)
294 proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
301 /* Specialized: projection is stored in vec->refproj
302 * -> used for radacc, radfix, radcon and center of flooding potential
303 * subtracts average positions, projects vector x */
304 static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec)
309 /* Subtract average positions */
310 for (i = 0; i < edi->sav.nr; i++)
312 rvec_dec(x[i], edi->sav.x[i]);
315 for (i = 0; i < vec->neig; i++)
317 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
318 rad += gmx::square((vec->refproj[i]-vec->xproj[i]));
320 vec->radius = sqrt(rad);
322 /* Add average positions */
323 for (i = 0; i < edi->sav.nr; i++)
325 rvec_inc(x[i], edi->sav.x[i]);
330 /* Project vector x, subtract average positions prior to projection and add
331 * them afterwards to retain the unchanged vector. Store in xproj. Mass-weighting
333 static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
334 t_eigvec *vec, /* The eigenvectors */
345 /* Subtract average positions */
346 for (i = 0; i < edi->sav.nr; i++)
348 rvec_dec(x[i], edi->sav.x[i]);
351 for (i = 0; i < vec->neig; i++)
353 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
356 /* Add average positions */
357 for (i = 0; i < edi->sav.nr; i++)
359 rvec_inc(x[i], edi->sav.x[i]);
364 /* Project vector x onto all edi->vecs (mon, linfix,...) */
365 static void project(rvec *x, /* positions to project */
366 t_edpar *edi) /* edi data set */
368 /* It is not more work to subtract the average position in every
369 * subroutine again, because these routines are rarely used simultaneously */
370 project_to_eigvectors(x, &edi->vecs.mon, edi);
371 project_to_eigvectors(x, &edi->vecs.linfix, edi);
372 project_to_eigvectors(x, &edi->vecs.linacc, edi);
373 project_to_eigvectors(x, &edi->vecs.radfix, edi);
374 project_to_eigvectors(x, &edi->vecs.radacc, edi);
375 project_to_eigvectors(x, &edi->vecs.radcon, edi);
379 static real calc_radius(t_eigvec *vec)
385 for (i = 0; i < vec->neig; i++)
387 rad += gmx::square((vec->refproj[i]-vec->xproj[i]));
390 return rad = sqrt(rad);
395 static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
400 fprintf(out, "#%s positions:\n%d\n", name, s->nr);
406 fprintf(out, "#index, x, y, z");
409 fprintf(out, ", sqrt(m)");
411 for (i = 0; i < s->nr; i++)
413 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]);
416 fprintf(out, "%9.3f", s->sqrtm[i]);
424 static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
425 const char name[], int length)
430 fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig);
431 /* Dump the data for every eigenvector: */
432 for (i = 0; i < ev->neig; i++)
434 fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
435 ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
436 for (j = 0; j < length; j++)
438 fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
445 static void dump_edi(t_edpar *edpars, const t_commrec *cr, int nr_edi)
451 sprintf(fn, "EDdump_rank%d_edi%d", cr->nodeid, nr_edi);
452 out = gmx_ffopen(fn, "w");
454 fprintf(out, "#NINI\n %d\n#FITMAS\n %s\n#ANALYSIS_MAS\n %s\n",
455 edpars->nini, gmx::boolToString(edpars->fitmas), gmx::boolToString(edpars->pcamas));
456 fprintf(out, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
457 edpars->outfrq, edpars->maxedsteps, edpars->slope);
458 fprintf(out, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n",
459 edpars->presteps, edpars->flood.deltaF0, edpars->flood.tau,
460 edpars->flood.constEfl, edpars->flood.alpha2);
462 /* Dump reference, average, target, origin positions */
463 dump_edi_positions(out, &edpars->sref, "REFERENCE");
464 dump_edi_positions(out, &edpars->sav, "AVERAGE" );
465 dump_edi_positions(out, &edpars->star, "TARGET" );
466 dump_edi_positions(out, &edpars->sori, "ORIGIN" );
468 /* Dump eigenvectors */
469 dump_edi_eigenvecs(out, &edpars->vecs.mon, "MONITORED", edpars->sav.nr);
470 dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX", edpars->sav.nr);
471 dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC", edpars->sav.nr);
472 dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX", edpars->sav.nr);
473 dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC", edpars->sav.nr);
474 dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON", edpars->sav.nr);
476 /* Dump flooding eigenvectors */
477 dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING", edpars->sav.nr);
479 /* Dump ed local buffer */
480 fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit );
481 fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam );
482 fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon );
493 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
495 /* this is a copy of do_fit with some modifications */
496 int c, r, n, j, i, irot;
497 double d[6], xnr, xpc;
502 struct t_do_edfit *loc;
505 if (edi->buf->do_edfit != nullptr)
512 snew(edi->buf->do_edfit, 1);
514 loc = edi->buf->do_edfit;
518 snew(loc->omega, 2*DIM);
519 snew(loc->om, 2*DIM);
520 for (i = 0; i < 2*DIM; i++)
522 snew(loc->omega[i], 2*DIM);
523 snew(loc->om[i], 2*DIM);
527 for (i = 0; (i < 6); i++)
530 for (j = 0; (j < 6); j++)
532 loc->omega[i][j] = 0;
537 /* calculate the matrix U */
539 for (n = 0; (n < natoms); n++)
541 for (c = 0; (c < DIM); c++)
544 for (r = 0; (r < DIM); r++)
552 /* construct loc->omega */
553 /* loc->omega is symmetric -> loc->omega==loc->omega' */
554 for (r = 0; (r < 6); r++)
556 for (c = 0; (c <= r); c++)
558 if ((r >= 3) && (c < 3))
560 loc->omega[r][c] = u[r-3][c];
561 loc->omega[c][r] = u[r-3][c];
565 loc->omega[r][c] = 0;
566 loc->omega[c][r] = 0;
571 /* determine h and k */
572 jacobi(loc->omega, 6, d, loc->om, &irot);
576 fprintf(stderr, "IROT=0\n");
579 index = 0; /* For the compiler only */
581 for (j = 0; (j < 3); j++)
584 for (i = 0; (i < 6); i++)
593 for (i = 0; (i < 3); i++)
595 vh[j][i] = M_SQRT2*loc->om[i][index];
596 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
601 for (c = 0; (c < 3); c++)
603 for (r = 0; (r < 3); r++)
605 R[c][r] = vk[0][r]*vh[0][c]+
612 for (c = 0; (c < 3); c++)
614 for (r = 0; (r < 3); r++)
616 R[c][r] = vk[0][r]*vh[0][c]+
625 static void rmfit(int nat, rvec *xcoll, const rvec transvec, matrix rotmat)
632 * The inverse rotation is described by the transposed rotation matrix */
633 transpose(rotmat, tmat);
634 rotate_x(xcoll, nat, tmat);
636 /* Remove translation */
637 vec[XX] = -transvec[XX];
638 vec[YY] = -transvec[YY];
639 vec[ZZ] = -transvec[ZZ];
640 translate_x(xcoll, nat, vec);
644 /**********************************************************************************
645 ******************** FLOODING ****************************************************
646 **********************************************************************************
648 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
649 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
650 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
652 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
653 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
655 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
656 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
658 do_flood makes a copy of the positions,
659 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
660 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
661 fit. Then do_flood adds these forces to the forcefield-forces
662 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
664 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
665 structure is projected to the system of eigenvectors and then this position in the subspace is used as
666 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
667 i.e. the average structure as given in the make_edi file.
669 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
670 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
671 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
672 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
673 further adaption is applied, Efl will stay constant at zero.
675 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
676 used as spring constants for the harmonic potential.
677 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
678 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
680 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
681 the routine read_edi_file reads all of theses flooding files.
682 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
683 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
684 edi there is no interdependence whatsoever. The forces are added together.
686 To write energies into the .edr file, call the function
687 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
689 get_flood_energies(real Vfl[],int nnames);
692 - 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.
694 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
695 two edsam files from two peptide chains
698 static void write_edo_flood(t_edpar *edi, FILE *fp, real rmsd)
703 /* Output how well we fit to the reference structure */
704 fprintf(fp, EDcol_ffmt, rmsd);
706 for (i = 0; i < edi->flood.vecs.neig; i++)
708 fprintf(fp, EDcol_efmt, edi->flood.vecs.xproj[i]);
710 /* Check whether the reference projection changes with time (this can happen
711 * in case flooding is used as harmonic restraint). If so, output the
712 * current reference projection */
713 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
715 fprintf(fp, EDcol_efmt, edi->flood.vecs.refproj[i]);
718 /* Output Efl if we are doing adaptive flooding */
719 if (0 != edi->flood.tau)
721 fprintf(fp, EDcol_efmt, edi->flood.Efl);
723 fprintf(fp, EDcol_efmt, edi->flood.Vfl);
725 /* Output deltaF if we are doing adaptive flooding */
726 if (0 != edi->flood.tau)
728 fprintf(fp, EDcol_efmt, edi->flood.deltaF);
730 fprintf(fp, EDcol_efmt, edi->flood.vecs.fproj[i]);
735 /* From flood.xproj compute the Vfl(x) at this point */
736 static real flood_energy(t_edpar *edi, int64_t step)
738 /* compute flooding energy Vfl
739 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
740 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
741 it is already computed by make_edi and stored in stpsz[i]
743 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
750 /* Each time this routine is called (i.e. each time step), we add a small
751 * value to the reference projection. This way a harmonic restraint towards
752 * a moving reference is realized. If no value for the additive constant
753 * is provided in the edi file, the reference will not change. */
754 if (edi->flood.bHarmonic)
756 for (i = 0; i < edi->flood.vecs.neig; i++)
758 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
763 /* Compute sum which will be the exponent of the exponential */
764 for (i = 0; i < edi->flood.vecs.neig; i++)
766 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
767 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]);
770 /* Compute the Gauss function*/
771 if (edi->flood.bHarmonic)
773 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
777 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*std::exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
784 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
785 static void flood_forces(t_edpar *edi)
787 /* compute the forces in the subspace of the flooding eigenvectors
788 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
791 real energy = edi->flood.Vfl;
794 if (edi->flood.bHarmonic)
796 for (i = 0; i < edi->flood.vecs.neig; i++)
798 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
803 for (i = 0; i < edi->flood.vecs.neig; i++)
805 /* if Efl is zero the forces are zero if not use the formula */
806 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;
812 /* Raise forces from subspace into cartesian space */
813 static void flood_blowup(t_edpar *edi, rvec *forces_cart)
815 /* this function lifts the forces from the subspace to the cartesian space
816 all the values not contained in the subspace are assumed to be zero and then
817 a coordinate transformation from eigenvector to cartesian vectors is performed
818 The nonexistent values don't have to be set to zero explicitly, they would occur
819 as zero valued summands, hence we just stop to compute this part of the sum.
821 for every atom we add all the contributions to this atom from all the different eigenvectors.
823 NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
824 field forces_cart prior the computation, but we compute the forces separately
825 to have them accessible for diagnostics
832 forces_sub = edi->flood.vecs.fproj;
835 /* Calculate the cartesian forces for the local atoms */
837 /* Clear forces first */
838 for (j = 0; j < edi->sav.nr_loc; j++)
840 clear_rvec(forces_cart[j]);
843 /* Now compute atomwise */
844 for (j = 0; j < edi->sav.nr_loc; j++)
846 /* Compute forces_cart[edi->sav.anrs[j]] */
847 for (eig = 0; eig < edi->flood.vecs.neig; eig++)
849 /* Force vector is force * eigenvector (compute only atom j) */
850 svmul(forces_sub[eig], edi->flood.vecs.vec[eig][edi->sav.c_ind[j]], dum);
851 /* Add this vector to the cartesian forces */
852 rvec_inc(forces_cart[j], dum);
858 /* Update the values of Efl, deltaF depending on tau and Vfl */
859 static void update_adaption(t_edpar *edi)
861 /* this function updates the parameter Efl and deltaF according to the rules given in
862 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
865 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
867 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
868 /* check if restrain (inverted flooding) -> don't let EFL become positive */
869 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
874 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
879 static void do_single_flood(
887 gmx_bool bNS) /* Are we in a neighbor searching step? */
890 matrix rotmat; /* rotation matrix */
891 matrix tmat; /* inverse rotation */
892 rvec transvec; /* translation vector */
894 struct t_do_edsam *buf;
897 buf = edi->buf->do_edsam;
900 /* Broadcast the positions of the AVERAGE structure such that they are known on
901 * every processor. Each node contributes its local positions x and stores them in
902 * the collective ED array buf->xcoll */
903 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
904 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
906 /* Only assembly REFERENCE positions if their indices differ from the average ones */
909 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
910 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
913 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
914 * We do not need to update the shifts until the next NS step */
915 buf->bUpdateShifts = FALSE;
917 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
918 * as well as the indices in edi->sav.anrs */
920 /* Fit the reference indices to the reference structure */
923 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
927 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
930 /* Now apply the translation and rotation to the ED structure */
931 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
933 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
934 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, edi);
936 if (FALSE == edi->flood.bConstForce)
938 /* Compute Vfl(x) from flood.xproj */
939 edi->flood.Vfl = flood_energy(edi, step);
941 update_adaption(edi);
943 /* Compute the flooding forces */
947 /* Translate them into cartesian positions */
948 flood_blowup(edi, edi->flood.forces_cartesian);
950 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
951 /* Each node rotates back its local forces */
952 transpose(rotmat, tmat);
953 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
955 /* Finally add forces to the main force variable */
956 for (i = 0; i < edi->sav.nr_loc; i++)
958 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
961 /* Output is written by the master process */
962 if (do_per_step(step, edi->outfrq) && MASTER(cr))
964 /* Output how well we fit to the reference */
967 /* Indices of reference and average structures are identical,
968 * thus we can calculate the rmsd to SREF using xcoll */
969 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
973 /* We have to translate & rotate the reference atoms first */
974 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
975 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
978 write_edo_flood(edi, edo, rmsdev);
983 /* Main flooding routine, called from do_force */
984 extern void do_flood(const t_commrec *cr,
985 const t_inputrec *ir,
998 /* Write time to edo, when required. Output the time anyhow since we need
999 * it in the output file for ED constraints. */
1000 if (MASTER(cr) && do_per_step(step, edi->outfrq))
1002 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
1005 if (ed->eEDtype != eEDflood)
1012 /* Call flooding for one matrix */
1013 if (edi->flood.vecs.neig)
1015 do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS);
1017 edi = edi->next_edi;
1022 /* Called by init_edi, configure some flooding related variables and structures,
1023 * print headers to output files */
1024 static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt)
1029 edi->flood.Efl = edi->flood.constEfl;
1033 if (edi->flood.vecs.neig)
1035 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1036 ed->eEDtype = eEDflood;
1038 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
1040 if (edi->flood.bConstForce)
1042 /* We have used stpsz as a vehicle to carry the fproj values for constant
1043 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1044 * in const force flooding, fproj is never changed. */
1045 for (i = 0; i < edi->flood.vecs.neig; i++)
1047 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1049 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1050 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1058 /*********** Energy book keeping ******/
1059 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1068 srenew(names, count);
1069 sprintf(buf, "Vfl_%d", count);
1070 names[count-1] = gmx_strdup(buf);
1071 actual = actual->next_edi;
1078 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1080 /*fl has to be big enough to capture nnames-many entries*/
1089 Vfl[count-1] = actual->flood.Vfl;
1090 actual = actual->next_edi;
1093 if (nnames != count-1)
1095 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1098 /************* END of FLOODING IMPLEMENTATION ****************************/
1102 /* This function opens the ED input and output files, reads in all datasets it finds
1103 * in the input file, and cross-checks whether the .edi file information is consistent
1104 * with the essential dynamics data found in the checkpoint file (if present).
1105 * gmx make_edi can be used to create an .edi input file.
1107 static gmx_edsam_t ed_open(
1109 ObservablesHistory *oh,
1110 const char *ediFileName,
1111 const char *edoFileName,
1113 const gmx_output_env_t *oenv,
1114 const t_commrec *cr)
1119 /* Allocate space for the ED data structure */
1122 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1123 ed->eEDtype = eEDedsam;
1129 // If we start from a checkpoint file, we already have an edsamHistory struct
1130 if (oh->edsamHistory == nullptr)
1132 oh->edsamHistory = gmx::compat::make_unique<edsamhistory_t>(edsamhistory_t {});
1134 edsamhistory_t *EDstate = oh->edsamHistory.get();
1136 /* Read the edi input file: */
1137 nED = read_edi_file(ediFileName, ed->edpar, natoms);
1139 /* Make sure the checkpoint was produced in a run using this .edi file */
1140 if (EDstate->bFromCpt)
1142 crosscheck_edi_file_vs_checkpoint(ed, EDstate);
1148 init_edsamstate(ed, EDstate);
1150 /* The master opens the ED output file */
1153 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1157 ed->edo = xvgropen(edoFileName,
1158 "Essential dynamics / flooding output",
1160 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1162 /* Make a descriptive legend */
1163 write_edo_legend(ed, EDstate->nED, oenv);
1170 /* Broadcasts the structure data */
1171 static void bc_ed_positions(const t_commrec *cr, struct gmx_edx *s, int stype)
1173 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1174 snew_bc(cr, s->x, s->nr ); /* Positions */
1175 nblock_bc(cr, s->nr, s->anrs );
1176 nblock_bc(cr, s->nr, s->x );
1178 /* For the average & reference structures we need an array for the collective indices,
1179 * and we need to broadcast the masses as well */
1180 if (stype == eedAV || stype == eedREF)
1182 /* We need these additional variables in the parallel case: */
1183 snew(s->c_ind, s->nr ); /* Collective indices */
1184 /* Local atom indices get assigned in dd_make_local_group_indices.
1185 * There, also memory is allocated */
1186 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1187 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1188 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1191 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1192 if (stype == eedREF)
1194 snew_bc(cr, s->m, s->nr);
1195 nblock_bc(cr, s->nr, s->m);
1198 /* For the average structure we might need the masses for mass-weighting */
1201 snew_bc(cr, s->sqrtm, s->nr);
1202 nblock_bc(cr, s->nr, s->sqrtm);
1203 snew_bc(cr, s->m, s->nr);
1204 nblock_bc(cr, s->nr, s->m);
1209 /* Broadcasts the eigenvector data */
1210 static void bc_ed_vecs(const t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
1214 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1215 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1216 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1217 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1218 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1220 nblock_bc(cr, ev->neig, ev->ieig );
1221 nblock_bc(cr, ev->neig, ev->stpsz );
1222 nblock_bc(cr, ev->neig, ev->xproj );
1223 nblock_bc(cr, ev->neig, ev->fproj );
1224 nblock_bc(cr, ev->neig, ev->refproj);
1226 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1227 for (i = 0; i < ev->neig; i++)
1229 snew_bc(cr, ev->vec[i], length);
1230 nblock_bc(cr, length, ev->vec[i]);
1233 /* For harmonic restraints the reference projections can change with time */
1236 snew_bc(cr, ev->refproj0, ev->neig);
1237 snew_bc(cr, ev->refprojslope, ev->neig);
1238 nblock_bc(cr, ev->neig, ev->refproj0 );
1239 nblock_bc(cr, ev->neig, ev->refprojslope);
1244 /* Broadcasts the ED / flooding data to other nodes
1245 * and allocates memory where needed */
1246 static void broadcast_ed_data(const t_commrec *cr, gmx_edsam_t ed, int numedis)
1252 /* Master lets the other nodes know if its ED only or also flooding */
1253 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1255 snew_bc(cr, ed->edpar, 1);
1256 /* Now transfer the ED data set(s) */
1258 for (nr = 0; nr < numedis; nr++)
1260 /* Broadcast a single ED data set */
1263 /* Broadcast positions */
1264 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1265 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1266 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1267 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1269 /* Broadcast eigenvectors */
1270 bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr, FALSE);
1271 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
1272 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
1273 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
1274 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
1275 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
1276 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1277 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
1279 /* Set the pointer to the next ED group */
1282 snew_bc(cr, edi->next_edi, 1);
1283 edi = edi->next_edi;
1289 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1290 static void init_edi(const gmx_mtop_t *mtop, t_edpar *edi)
1293 real totalmass = 0.0;
1296 /* NOTE Init_edi is executed on the master process only
1297 * The initialized data sets are then transmitted to the
1298 * other nodes in broadcast_ed_data */
1300 /* evaluate masses (reference structure) */
1301 snew(edi->sref.m, edi->sref.nr);
1303 for (i = 0; i < edi->sref.nr; i++)
1307 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1311 edi->sref.m[i] = 1.0;
1314 /* Check that every m > 0. Bad things will happen otherwise. */
1315 if (edi->sref.m[i] <= 0.0)
1317 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1318 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1319 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1320 "atoms from the reference structure by creating a proper index group.\n",
1321 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1324 totalmass += edi->sref.m[i];
1326 edi->sref.mtot = totalmass;
1328 /* Masses m and sqrt(m) for the average structure. Note that m
1329 * is needed if forces have to be evaluated in do_edsam */
1330 snew(edi->sav.sqrtm, edi->sav.nr );
1331 snew(edi->sav.m, edi->sav.nr );
1332 for (i = 0; i < edi->sav.nr; i++)
1334 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1337 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1341 edi->sav.sqrtm[i] = 1.0;
1344 /* Check that every m > 0. Bad things will happen otherwise. */
1345 if (edi->sav.sqrtm[i] <= 0.0)
1347 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1348 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1349 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1350 "atoms from the average structure by creating a proper index group.\n",
1351 i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1355 /* put reference structure in origin */
1356 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1360 translate_x(edi->sref.x, edi->sref.nr, com);
1362 /* Init ED buffer */
1367 static void check(const char *line, const char *label)
1369 if (!strstr(line, label))
1371 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1376 static int read_checked_edint(FILE *file, const char *label)
1378 char line[STRLEN+1];
1381 fgets2 (line, STRLEN, file);
1383 fgets2 (line, STRLEN, file);
1384 sscanf (line, max_ev_fmt_d, &idum);
1389 static int read_edint(FILE *file, bool *bEOF)
1391 char line[STRLEN+1];
1396 eof = fgets2 (line, STRLEN, file);
1402 eof = fgets2 (line, STRLEN, file);
1408 sscanf (line, max_ev_fmt_d, &idum);
1414 static real read_checked_edreal(FILE *file, const char *label)
1416 char line[STRLEN+1];
1420 fgets2 (line, STRLEN, file);
1422 fgets2 (line, STRLEN, file);
1423 sscanf (line, max_ev_fmt_lf, &rdum);
1424 return static_cast<real>(rdum); /* always read as double and convert to single */
1428 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1431 char line[STRLEN+1];
1435 for (i = 0; i < number; i++)
1437 fgets2 (line, STRLEN, file);
1438 sscanf (line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1439 anrs[i]--; /* we are reading FORTRAN indices */
1440 for (j = 0; j < 3; j++)
1442 x[i][j] = d[j]; /* always read as double and convert to single */
1448 static void scan_edvec(FILE *in, int nr, rvec *vec)
1450 char line[STRLEN+1];
1455 for (i = 0; (i < nr); i++)
1457 fgets2 (line, STRLEN, in);
1458 sscanf (line, max_ev_fmt_lelele, &x, &y, &z);
1466 static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1469 double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0;
1470 char line[STRLEN+1];
1473 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1476 snew(tvec->ieig, tvec->neig);
1477 snew(tvec->stpsz, tvec->neig);
1478 snew(tvec->vec, tvec->neig);
1479 snew(tvec->xproj, tvec->neig);
1480 snew(tvec->fproj, tvec->neig);
1481 snew(tvec->refproj, tvec->neig);
1484 snew(tvec->refproj0, tvec->neig);
1485 snew(tvec->refprojslope, tvec->neig);
1488 for (i = 0; (i < tvec->neig); i++)
1490 fgets2 (line, STRLEN, in);
1491 if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1493 nscan = sscanf(line, max_ev_fmt_dlflflf, &idum, &rdum, &refproj_dum, &refprojslope_dum);
1494 /* Zero out values which were not scanned */
1498 /* Every 4 values read, including reference position */
1499 *bHaveReference = TRUE;
1502 /* A reference position is provided */
1503 *bHaveReference = TRUE;
1504 /* No value for slope, set to 0 */
1505 refprojslope_dum = 0.0;
1508 /* No values for reference projection and slope, set to 0 */
1510 refprojslope_dum = 0.0;
1513 gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1515 tvec->refproj[i] = refproj_dum;
1516 tvec->refproj0[i] = refproj_dum;
1517 tvec->refprojslope[i] = refprojslope_dum;
1519 else /* Normal flooding */
1521 nscan = sscanf(line, max_ev_fmt_dlf, &idum, &rdum);
1524 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1527 tvec->ieig[i] = idum;
1528 tvec->stpsz[i] = rdum;
1529 } /* end of loop over eigenvectors */
1531 for (i = 0; (i < tvec->neig); i++)
1533 snew(tvec->vec[i], nr);
1534 scan_edvec(in, nr, tvec->vec[i]);
1540 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1541 static void read_edvecs(FILE *in, int nr, t_edvecs *vecs)
1543 gmx_bool bHaveReference = FALSE;
1546 read_edvec(in, nr, &vecs->mon, FALSE, &bHaveReference);
1547 read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1548 read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1549 read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1550 read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1551 read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1555 /* Check if the same atom indices are used for reference and average positions */
1556 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1561 /* If the number of atoms differs between the two structures,
1562 * they cannot be identical */
1563 if (sref.nr != sav.nr)
1568 /* Now that we know that both stuctures have the same number of atoms,
1569 * check if also the indices are identical */
1570 for (i = 0; i < sav.nr; i++)
1572 if (sref.anrs[i] != sav.anrs[i])
1577 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1585 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1586 constexpr int c_oldestUnsupportedMagicNumber = 666;
1587 //! Oldest (lowest) magic number indicating supported essential dynamics input
1588 constexpr int c_oldestSupportedMagicNumber = 669;
1589 //! Latest (highest) magic number indicating supported essential dynamics input
1590 constexpr int c_latestSupportedMagicNumber = 670;
1592 /*!\brief Set up essential dynamics work parameters.
1593 * \param[in] edi Essential dynamics working structure.
1595 void setup_edi(t_edpar *edi)
1597 snew(edi->sref.x_old, edi->sref.nr);
1598 edi->sref.sqrtm = nullptr;
1599 snew(edi->sav.x_old, edi->sav.nr);
1600 if (edi->star.nr > 0)
1602 edi->star.sqrtm = nullptr;
1604 if (edi->sori.nr > 0)
1606 edi->sori.sqrtm = nullptr;
1610 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1611 * \param[in] magicNumber indicates the essential dynamics input file version
1612 * \returns true if CONST_FORCE_FLOODING is supported
1614 bool ediFileHasConstForceFlooding(int magicNumber)
1616 return magicNumber > c_oldestSupportedMagicNumber;
1619 /*!\brief Reports reading success of the essential dynamics input file magic number.
1620 * \param[in] in pointer to input file
1621 * \param[out] magicNumber determines the edi file version
1622 * \returns true if setting file version from input file was successful.
1624 bool ediFileHasValidDataSet(FILE *in, int * magicNumber)
1626 /* the edi file is not free format, so expect problems if the input is corrupt. */
1627 bool reachedEndOfFile = true;
1628 /* check the magic number */
1629 *magicNumber = read_edint(in, &reachedEndOfFile);
1630 /* Check whether we have reached the end of the input file */
1631 return !reachedEndOfFile;
1634 /*!\brief Trigger fatal error if magic number is unsupported.
1635 * \param[in] magicNumber A magic number read from the edi file
1636 * \param[in] fn name of input file for error reporting
1638 void exitWhenMagicNumberIsInvalid(int magicNumber, const char * fn)
1641 if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1643 if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1645 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1649 gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1654 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1656 * \param[in] in input file
1657 * \param[in] edi essential dynamics parameters
1658 * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1659 * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1660 * \param[in] fn the file name of the input file for error reporting
1662 void read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, bool hasConstForceFlooding, const char *fn)
1664 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1665 gmx_bool bHaveReference = FALSE;
1668 /* check the number of atoms */
1669 edi->nini = read_edint(in, &bEOF);
1670 if (edi->nini != nr_mdatoms)
1672 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1675 /* Done checking. For the rest we blindly trust the input */
1676 edi->fitmas = read_checked_edint(in, "FITMAS");
1677 edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS");
1678 edi->outfrq = read_checked_edint(in, "OUTFRQ");
1679 edi->maxedsteps = read_checked_edint(in, "MAXLEN");
1680 edi->slope = read_checked_edreal(in, "SLOPECRIT");
1682 edi->presteps = read_checked_edint(in, "PRESTEPS");
1683 edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1684 edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1685 edi->flood.tau = read_checked_edreal(in, "TAU");
1686 edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1687 edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1688 edi->flood.kT = read_checked_edreal(in, "KT");
1689 edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1690 if (hasConstForceFlooding)
1692 edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1696 edi->flood.bConstForce = FALSE;
1698 edi->sref.nr = read_checked_edint(in, "NREF");
1700 /* allocate space for reference positions and read them */
1701 snew(edi->sref.anrs, edi->sref.nr);
1702 snew(edi->sref.x, edi->sref.nr);
1703 read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1705 /* average positions. they define which atoms will be used for ED sampling */
1706 edi->sav.nr = read_checked_edint(in, "NAV");
1707 snew(edi->sav.anrs, edi->sav.nr);
1708 snew(edi->sav.x, edi->sav.nr);
1709 read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1711 /* Check if the same atom indices are used for reference and average positions */
1712 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1715 read_edvecs(in, edi->sav.nr, &edi->vecs);
1716 read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
1718 /* target positions */
1719 edi->star.nr = read_edint(in, &bEOF);
1720 if (edi->star.nr > 0)
1722 snew(edi->star.anrs, edi->star.nr);
1723 snew(edi->star.x, edi->star.nr);
1724 read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1727 /* positions defining origin of expansion circle */
1728 edi->sori.nr = read_edint(in, &bEOF);
1729 if (edi->sori.nr > 0)
1733 /* Both an -ori structure and a at least one manual reference point have been
1734 * specified. That's ambiguous and probably not intentional. */
1735 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1736 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1738 snew(edi->sori.anrs, edi->sori.nr);
1739 snew(edi->sori.x, edi->sori.nr);
1740 read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1744 } // anonymous namespace
1746 /* Read in the edi input file. Note that it may contain several ED data sets which were
1747 * achieved by concatenating multiple edi files. The standard case would be a single ED
1748 * data set, though. */
1749 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1752 t_edpar *curr_edi, *last_edi;
1757 /* This routine is executed on the master only */
1759 /* Open the .edi parameter input file */
1760 in = gmx_fio_fopen(fn, "r");
1761 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1763 /* Now read a sequence of ED input parameter sets from the edi file */
1766 int ediFileMagicNumber;
1767 while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1769 exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1770 bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1771 read_edi(in, curr_edi, nr_mdatoms, hasConstForceFlooding, fn);
1772 setup_edi(curr_edi);
1775 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1776 /* We need to allocate space for the data: */
1778 /* Point the 'next_edi' entry to the next edi: */
1779 curr_edi->next_edi = edi_read;
1780 /* Keep the curr_edi pointer for the case that the next group is empty: */
1781 last_edi = curr_edi;
1782 /* Let's prepare to read in the next edi data set: */
1783 curr_edi = edi_read;
1787 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1790 /* Terminate the edi group list with a NULL pointer: */
1791 last_edi->next_edi = nullptr;
1793 fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1795 /* Close the .edi file again */
1802 struct t_fit_to_ref {
1803 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1806 /* Fit the current positions to the reference positions
1807 * Do not actually do the fit, just return rotation and translation.
1808 * Note that the COM of the reference structure was already put into
1809 * the origin by init_edi. */
1810 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1811 rvec transvec, /* The translation vector */
1812 matrix rotmat, /* The rotation matrix */
1813 t_edpar *edi) /* Just needed for do_edfit */
1815 rvec com; /* center of mass */
1817 struct t_fit_to_ref *loc;
1820 /* Allocate memory the first time this routine is called for each edi group */
1821 if (nullptr == edi->buf->fit_to_ref)
1823 snew(edi->buf->fit_to_ref, 1);
1824 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1826 loc = edi->buf->fit_to_ref;
1828 /* We do not touch the original positions but work on a copy. */
1829 for (i = 0; i < edi->sref.nr; i++)
1831 copy_rvec(xcoll[i], loc->xcopy[i]);
1834 /* Calculate the center of mass */
1835 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1837 transvec[XX] = -com[XX];
1838 transvec[YY] = -com[YY];
1839 transvec[ZZ] = -com[ZZ];
1841 /* Subtract the center of mass from the copy */
1842 translate_x(loc->xcopy, edi->sref.nr, transvec);
1844 /* Determine the rotation matrix */
1845 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1849 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1850 int nat, /* How many positions are there? */
1851 rvec transvec, /* The translation vector */
1852 matrix rotmat) /* The rotation matrix */
1855 translate_x(x, nat, transvec);
1858 rotate_x(x, nat, rotmat);
1862 /* Gets the rms deviation of the positions to the structure s */
1863 /* fit_to_structure has to be called before calling this routine! */
1864 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1865 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1871 for (i = 0; i < s->nr; i++)
1873 rmsd += distance2(s->x[i], x[i]);
1876 rmsd /= static_cast<real>(s->nr);
1883 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1888 if (ed->eEDtype != eEDnone)
1890 /* Loop over ED groups */
1894 /* Local atoms of the reference structure (for fitting), need only be assembled
1895 * if their indices differ from the average ones */
1898 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1899 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1902 /* Local atoms of the average structure (on these ED will be performed) */
1903 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1904 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1906 /* Indicate that the ED shift vectors for this structure need to be updated
1907 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1908 edi->buf->do_edsam->bUpdateShifts = TRUE;
1910 /* Set the pointer to the next ED group (if any) */
1911 edi = edi->next_edi;
1917 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1928 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1929 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1930 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1934 xu[XX] = x[XX]-tx*box[XX][XX];
1935 xu[YY] = x[YY]-ty*box[YY][YY];
1936 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1941 static void do_linfix(rvec *xcoll, t_edpar *edi, int64_t step)
1948 /* loop over linfix vectors */
1949 for (i = 0; i < edi->vecs.linfix.neig; i++)
1951 /* calculate the projection */
1952 proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1954 /* calculate the correction */
1955 add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1957 /* apply the correction */
1958 add /= edi->sav.sqrtm[i];
1959 for (j = 0; j < edi->sav.nr; j++)
1961 svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1962 rvec_inc(xcoll[j], vec_dum);
1968 static void do_linacc(rvec *xcoll, t_edpar *edi)
1975 /* loop over linacc vectors */
1976 for (i = 0; i < edi->vecs.linacc.neig; i++)
1978 /* calculate the projection */
1979 proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1981 /* calculate the correction */
1983 if (edi->vecs.linacc.stpsz[i] > 0.0)
1985 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1987 add = edi->vecs.linacc.refproj[i] - proj;
1990 if (edi->vecs.linacc.stpsz[i] < 0.0)
1992 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1994 add = edi->vecs.linacc.refproj[i] - proj;
1998 /* apply the correction */
1999 add /= edi->sav.sqrtm[i];
2000 for (j = 0; j < edi->sav.nr; j++)
2002 svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
2003 rvec_inc(xcoll[j], vec_dum);
2006 /* new positions will act as reference */
2007 edi->vecs.linacc.refproj[i] = proj + add;
2012 static void do_radfix(rvec *xcoll, t_edpar *edi)
2015 real *proj, rad = 0.0, ratio;
2019 if (edi->vecs.radfix.neig == 0)
2024 snew(proj, edi->vecs.radfix.neig);
2026 /* loop over radfix vectors */
2027 for (i = 0; i < edi->vecs.radfix.neig; i++)
2029 /* calculate the projections, radius */
2030 proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
2031 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2035 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2036 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2038 /* loop over radfix vectors */
2039 for (i = 0; i < edi->vecs.radfix.neig; i++)
2041 proj[i] -= edi->vecs.radfix.refproj[i];
2043 /* apply the correction */
2044 proj[i] /= edi->sav.sqrtm[i];
2046 for (j = 0; j < edi->sav.nr; j++)
2048 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2049 rvec_inc(xcoll[j], vec_dum);
2057 static void do_radacc(rvec *xcoll, t_edpar *edi)
2060 real *proj, rad = 0.0, ratio = 0.0;
2064 if (edi->vecs.radacc.neig == 0)
2069 snew(proj, edi->vecs.radacc.neig);
2071 /* loop over radacc vectors */
2072 for (i = 0; i < edi->vecs.radacc.neig; i++)
2074 /* calculate the projections, radius */
2075 proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2076 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2080 /* only correct when radius decreased */
2081 if (rad < edi->vecs.radacc.radius)
2083 ratio = edi->vecs.radacc.radius/rad - 1.0;
2087 edi->vecs.radacc.radius = rad;
2090 /* loop over radacc vectors */
2091 for (i = 0; i < edi->vecs.radacc.neig; i++)
2093 proj[i] -= edi->vecs.radacc.refproj[i];
2095 /* apply the correction */
2096 proj[i] /= edi->sav.sqrtm[i];
2098 for (j = 0; j < edi->sav.nr; j++)
2100 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2101 rvec_inc(xcoll[j], vec_dum);
2108 struct t_do_radcon {
2112 static void do_radcon(rvec *xcoll, t_edpar *edi)
2115 real rad = 0.0, ratio = 0.0;
2116 struct t_do_radcon *loc;
2121 if (edi->buf->do_radcon != nullptr)
2128 snew(edi->buf->do_radcon, 1);
2130 loc = edi->buf->do_radcon;
2132 if (edi->vecs.radcon.neig == 0)
2139 snew(loc->proj, edi->vecs.radcon.neig);
2142 /* loop over radcon vectors */
2143 for (i = 0; i < edi->vecs.radcon.neig; i++)
2145 /* calculate the projections, radius */
2146 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2147 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2150 /* only correct when radius increased */
2151 if (rad > edi->vecs.radcon.radius)
2153 ratio = edi->vecs.radcon.radius/rad - 1.0;
2155 /* loop over radcon vectors */
2156 for (i = 0; i < edi->vecs.radcon.neig; i++)
2158 /* apply the correction */
2159 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2160 loc->proj[i] /= edi->sav.sqrtm[i];
2161 loc->proj[i] *= ratio;
2163 for (j = 0; j < edi->sav.nr; j++)
2165 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2166 rvec_inc(xcoll[j], vec_dum);
2173 edi->vecs.radcon.radius = rad;
2179 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, int64_t step)
2184 /* subtract the average positions */
2185 for (i = 0; i < edi->sav.nr; i++)
2187 rvec_dec(xcoll[i], edi->sav.x[i]);
2190 /* apply the constraints */
2193 do_linfix(xcoll, edi, step);
2195 do_linacc(xcoll, edi);
2198 do_radfix(xcoll, edi);
2200 do_radacc(xcoll, edi);
2201 do_radcon(xcoll, edi);
2203 /* add back the average positions */
2204 for (i = 0; i < edi->sav.nr; i++)
2206 rvec_inc(xcoll[i], edi->sav.x[i]);
2211 /* Write out the projections onto the eigenvectors. The order of output
2212 * corresponds to ed_output_legend() */
2213 static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
2218 /* Output how well we fit to the reference structure */
2219 fprintf(fp, EDcol_ffmt, rmsd);
2221 for (i = 0; i < edi->vecs.mon.neig; i++)
2223 fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
2226 for (i = 0; i < edi->vecs.linfix.neig; i++)
2228 fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
2231 for (i = 0; i < edi->vecs.linacc.neig; i++)
2233 fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
2236 for (i = 0; i < edi->vecs.radfix.neig; i++)
2238 fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
2240 if (edi->vecs.radfix.neig)
2242 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
2245 for (i = 0; i < edi->vecs.radacc.neig; i++)
2247 fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
2249 if (edi->vecs.radacc.neig)
2251 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
2254 for (i = 0; i < edi->vecs.radcon.neig; i++)
2256 fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
2258 if (edi->vecs.radcon.neig)
2260 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
2264 /* Returns if any constraints are switched on */
2265 static int ed_constraints(int edtype, t_edpar *edi)
2267 if (edtype == eEDedsam || edtype == eEDflood)
2269 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2270 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2271 edi->vecs.radcon.neig);
2277 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2278 * umbrella sampling simulations. */
2279 static void copyEvecReference(t_eigvec* floodvecs)
2284 if (nullptr == floodvecs->refproj0)
2286 snew(floodvecs->refproj0, floodvecs->neig);
2289 for (i = 0; i < floodvecs->neig; i++)
2291 floodvecs->refproj0[i] = floodvecs->refproj[i];
2296 /* Call on MASTER only. Check whether the essential dynamics / flooding
2297 * groups of the checkpoint file are consistent with the provided .edi file. */
2298 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamhistory_t *EDstate)
2300 t_edpar *edi = nullptr; /* points to a single edi data set */
2304 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2306 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2307 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2308 "it must also continue with/without ED constraints when checkpointing.\n"
2309 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2310 "from without a checkpoint.\n");
2315 while (edi != nullptr)
2317 /* Check number of atoms in the reference and average structures */
2318 if (EDstate->nref[edinum] != edi->sref.nr)
2320 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2321 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2322 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2324 if (EDstate->nav[edinum] != edi->sav.nr)
2326 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2327 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2328 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2330 edi = edi->next_edi;
2334 if (edinum != EDstate->nED)
2336 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2337 "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2338 "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2343 /* The edsamstate struct stores the information we need to make the ED group
2344 * whole again after restarts from a checkpoint file. Here we do the following:
2345 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2346 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2347 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2348 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2349 * all ED structures at the last time step. */
2350 static void init_edsamstate(gmx_edsam_t ed, edsamhistory_t *EDstate)
2356 snew(EDstate->old_sref_p, EDstate->nED);
2357 snew(EDstate->old_sav_p, EDstate->nED);
2359 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2360 if (!EDstate->bFromCpt)
2362 snew(EDstate->nref, EDstate->nED);
2363 snew(EDstate->nav, EDstate->nED);
2366 /* Loop over all ED/flooding data sets (usually only one, though) */
2368 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2370 /* We always need the last reference and average positions such that
2371 * in the next time step we can make the ED group whole again
2372 * if the atoms do not have the correct PBC representation */
2373 if (EDstate->bFromCpt)
2375 /* Copy the last whole positions of reference and average group from .cpt */
2376 for (i = 0; i < edi->sref.nr; i++)
2378 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2380 for (i = 0; i < edi->sav.nr; i++)
2382 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2387 EDstate->nref[nr_edi-1] = edi->sref.nr;
2388 EDstate->nav [nr_edi-1] = edi->sav.nr;
2391 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2392 EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2393 EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2395 edi = edi->next_edi;
2400 /* Adds 'buf' to 'str' */
2401 static void add_to_string(char **str, char *buf)
2406 len = strlen(*str) + strlen(buf) + 1;
2412 static void add_to_string_aligned(char **str, char *buf)
2414 char buf_aligned[STRLEN];
2416 sprintf(buf_aligned, EDcol_sfmt, buf);
2417 add_to_string(str, buf_aligned);
2421 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, const char *value, const char *unit, char EDgroupchar)
2423 char tmp[STRLEN], tmp2[STRLEN];
2426 sprintf(tmp, "%c %s", EDgroupchar, value);
2427 add_to_string_aligned(LegendStr, tmp);
2428 sprintf(tmp2, "%s (%s)", tmp, unit);
2429 (*setname)[*nsets] = gmx_strdup(tmp2);
2434 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2440 for (i = 0; i < evec->neig; i++)
2442 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2443 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2448 /* Makes a legend for the xvg output file. Call on MASTER only! */
2449 static void write_edo_legend(gmx_edsam_t ed, int nED, const gmx_output_env_t *oenv)
2451 t_edpar *edi = nullptr;
2453 int nr_edi, nsets, n_flood, n_edsam;
2454 const char **setname;
2456 char *LegendStr = nullptr;
2461 fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2463 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2465 fprintf(ed->edo, "#\n");
2466 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2467 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2468 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2469 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2470 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2471 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2472 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2473 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2474 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2476 if (edi->flood.vecs.neig)
2478 /* If in any of the groups we find a flooding vector, flooding is turned on */
2479 ed->eEDtype = eEDflood;
2481 /* Print what flavor of flooding we will do */
2482 if (0 == edi->flood.tau) /* constant flooding strength */
2484 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2485 if (edi->flood.bHarmonic)
2487 fprintf(ed->edo, ", harmonic");
2490 else /* adaptive flooding */
2492 fprintf(ed->edo, ", adaptive");
2495 fprintf(ed->edo, "\n");
2497 edi = edi->next_edi;
2500 /* Print a nice legend */
2502 LegendStr[0] = '\0';
2503 sprintf(buf, "# %6s", "time");
2504 add_to_string(&LegendStr, buf);
2506 /* Calculate the maximum number of columns we could end up with */
2509 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2511 nsets += 5 +edi->vecs.mon.neig
2512 +edi->vecs.linfix.neig
2513 +edi->vecs.linacc.neig
2514 +edi->vecs.radfix.neig
2515 +edi->vecs.radacc.neig
2516 +edi->vecs.radcon.neig
2517 + 6*edi->flood.vecs.neig;
2518 edi = edi->next_edi;
2520 snew(setname, nsets);
2522 /* In the mdrun time step in a first function call (do_flood()) the flooding
2523 * forces are calculated and in a second function call (do_edsam()) the
2524 * ED constraints. To get a corresponding legend, we need to loop twice
2525 * over the edi groups and output first the flooding, then the ED part */
2527 /* The flooding-related legend entries, if flooding is done */
2529 if (eEDflood == ed->eEDtype)
2532 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2534 /* Always write out the projection on the flooding EVs. Of course, this can also
2535 * be achieved with the monitoring option in do_edsam() (if switched on by the
2536 * user), but in that case the positions need to be communicated in do_edsam(),
2537 * which is not necessary when doing flooding only. */
2538 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2540 for (i = 0; i < edi->flood.vecs.neig; i++)
2542 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2543 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2545 /* Output the current reference projection if it changes with time;
2546 * this can happen when flooding is used as harmonic restraint */
2547 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
2549 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2550 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2553 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2554 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2556 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2557 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2560 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2561 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2563 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2565 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2566 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2569 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2570 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2573 edi = edi->next_edi;
2574 } /* End of flooding-related legend entries */
2578 /* Now the ED-related entries, if essential dynamics is done */
2580 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2582 if (bNeedDoEdsam(edi)) /* Only print ED legend if at least one ED option is on */
2584 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2586 /* Essential dynamics, projections on eigenvectors */
2587 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2588 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2589 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2590 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2591 if (edi->vecs.radfix.neig)
2593 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2595 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2596 if (edi->vecs.radacc.neig)
2598 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2600 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2601 if (edi->vecs.radcon.neig)
2603 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2606 edi = edi->next_edi;
2607 } /* end of 'pure' essential dynamics legend entries */
2608 n_edsam = nsets - n_flood;
2610 xvgr_legend(ed->edo, nsets, setname, oenv);
2613 fprintf(ed->edo, "#\n"
2614 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2615 n_flood, 1 == n_flood ? "" : "s",
2616 n_edsam, 1 == n_edsam ? "" : "s");
2617 fprintf(ed->edo, "%s", LegendStr);
2624 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2625 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2626 gmx_edsam_t init_edsam(
2627 const char *ediFileName,
2628 const char *edoFileName,
2629 const gmx_mtop_t *mtop,
2630 const t_inputrec *ir,
2631 const t_commrec *cr,
2632 gmx::Constraints *constr,
2633 const t_state *globalState,
2634 ObservablesHistory *oh,
2635 const gmx_output_env_t *oenv,
2638 t_edpar *edi = nullptr; /* points to a single edi data set */
2640 int nED = 0; /* Number of ED data sets */
2641 rvec *x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2642 rvec *xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2643 rvec fit_transvec; /* translation ... */
2644 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2645 rvec *ref_x_old = nullptr; /* helper pointer */
2650 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2652 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName) )
2654 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or flooding\n"
2655 "simulation. Please also set the .edi file on the command line with -ei.\n");
2659 /* Open input and output files, allocate space for ED data structure */
2660 gmx_edsam_t ed = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr);
2661 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2662 constr->saveEdsamPointer(ed);
2664 /* Needed for initializing radacc radius in do_edsam */
2667 /* The input file is read by the master and the edi structures are
2668 * initialized here. Input is stored in ed->edpar. Then the edi
2669 * structures are transferred to the other nodes */
2672 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2673 * flooding vector, Essential dynamics can be applied to more than one structure
2674 * as well, but will be done in the order given in the edi file, so
2675 * expect different results for different order of edi file concatenation! */
2677 while (edi != nullptr)
2679 init_edi(mtop, edi);
2680 init_flood(edi, ed, ir->delta_t);
2681 edi = edi->next_edi;
2685 /* The master does the work here. The other nodes get the positions
2686 * not before dd_partition_system which is called after init_edsam */
2689 edsamhistory_t *EDstate = oh->edsamHistory.get();
2691 if (!EDstate->bFromCpt)
2693 /* Remove PBC, make molecule(s) subject to ED whole. */
2694 snew(x_pbc, mtop->natoms);
2695 copy_rvecn(as_rvec_array(globalState->x.data()), x_pbc, 0, mtop->natoms);
2696 do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
2698 /* Reset pointer to first ED data set which contains the actual ED data */
2700 GMX_ASSERT(nullptr != edi,
2701 "The pointer to the essential dynamics parameters is undefined");
2704 /* Loop over all ED/flooding data sets (usually only one, though) */
2705 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2707 /* For multiple ED groups we use the output frequency that was specified
2708 * in the first set */
2711 edi->outfrq = ed->edpar->outfrq;
2714 /* Extract the initial reference and average positions. When starting
2715 * from .cpt, these have already been read into sref.x_old
2716 * in init_edsamstate() */
2717 if (!EDstate->bFromCpt)
2719 /* If this is the first run (i.e. no checkpoint present) we assume
2720 * that the starting positions give us the correct PBC representation */
2721 for (i = 0; i < edi->sref.nr; i++)
2723 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2726 for (i = 0; i < edi->sav.nr; i++)
2728 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2732 /* Now we have the PBC-correct start positions of the reference and
2733 average structure. We copy that over to dummy arrays on which we
2734 can apply fitting to print out the RMSD. We srenew the memory since
2735 the size of the buffers is likely different for every ED group */
2736 srenew(xfit, edi->sref.nr );
2737 srenew(xstart, edi->sav.nr );
2740 /* Reference indices are the same as average indices */
2741 ref_x_old = edi->sav.x_old;
2745 ref_x_old = edi->sref.x_old;
2747 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2748 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2750 /* Make the fit to the REFERENCE structure, get translation and rotation */
2751 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2753 /* Output how well we fit to the reference at the start */
2754 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2755 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2756 rmsd_from_structure(xfit, &edi->sref));
2757 if (EDstate->nED > 1)
2759 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2761 fprintf(stderr, "\n");
2763 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2764 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2766 /* calculate initial projections */
2767 project(xstart, edi);
2769 /* For the target and origin structure both a reference (fit) and an
2770 * average structure can be provided in make_edi. If both structures
2771 * are the same, make_edi only stores one of them in the .edi file.
2772 * If they differ, first the fit and then the average structure is stored
2773 * in star (or sor), thus the number of entries in star/sor is
2774 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2775 * the size of the average group. */
2777 /* process target structure, if required */
2778 if (edi->star.nr > 0)
2780 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2782 /* get translation & rotation for fit of target structure to reference structure */
2783 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2785 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2786 if (edi->star.nr == edi->sav.nr)
2790 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2792 /* The last sav.nr indices of the target structure correspond to
2793 * the average structure, which must be projected */
2794 avindex = edi->star.nr - edi->sav.nr;
2796 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
2800 rad_project(edi, xstart, &edi->vecs.radcon);
2803 /* process structure that will serve as origin of expansion circle */
2804 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2806 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2809 if (edi->sori.nr > 0)
2811 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2813 /* fit this structure to reference structure */
2814 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2816 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2817 if (edi->sori.nr == edi->sav.nr)
2821 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2823 /* For the projection, we need the last sav.nr indices of sori */
2824 avindex = edi->sori.nr - edi->sav.nr;
2827 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2828 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2829 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2831 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2832 /* Set center of flooding potential to the ORIGIN structure */
2833 rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
2834 /* We already know that no (moving) reference position was provided,
2835 * therefore we can overwrite refproj[0]*/
2836 copyEvecReference(&edi->flood.vecs);
2839 else /* No origin structure given */
2841 rad_project(edi, xstart, &edi->vecs.radacc);
2842 rad_project(edi, xstart, &edi->vecs.radfix);
2843 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2845 if (edi->flood.bHarmonic)
2847 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2848 for (i = 0; i < edi->flood.vecs.neig; i++)
2850 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2855 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2856 /* Set center of flooding potential to the center of the covariance matrix,
2857 * i.e. the average structure, i.e. zero in the projected system */
2858 for (i = 0; i < edi->flood.vecs.neig; i++)
2860 edi->flood.vecs.refproj[i] = 0.0;
2865 /* For convenience, output the center of the flooding potential for the eigenvectors */
2866 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2868 for (i = 0; i < edi->flood.vecs.neig; i++)
2870 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2871 if (edi->flood.bHarmonic)
2873 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2875 fprintf(stdout, "\n");
2879 /* set starting projections for linsam */
2880 rad_project(edi, xstart, &edi->vecs.linacc);
2881 rad_project(edi, xstart, &edi->vecs.linfix);
2883 /* Prepare for the next edi data set: */
2884 edi = edi->next_edi;
2886 /* Cleaning up on the master node: */
2887 if (!EDstate->bFromCpt)
2894 } /* end of MASTER only section */
2898 /* First let everybody know how many ED data sets to expect */
2899 gmx_bcast(sizeof(nED), &nED, cr);
2900 /* Broadcast the essential dynamics / flooding data to all nodes */
2901 broadcast_ed_data(cr, ed, nED);
2905 /* In the single-CPU case, point the local atom numbers pointers to the global
2906 * one, so that we can use the same notation in serial and parallel case: */
2907 /* Loop over all ED data sets (usually only one, though) */
2909 for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
2911 edi->sref.anrs_loc = edi->sref.anrs;
2912 edi->sav.anrs_loc = edi->sav.anrs;
2913 edi->star.anrs_loc = edi->star.anrs;
2914 edi->sori.anrs_loc = edi->sori.anrs;
2915 /* For the same reason as above, make a dummy c_ind array: */
2916 snew(edi->sav.c_ind, edi->sav.nr);
2917 /* Initialize the array */
2918 for (i = 0; i < edi->sav.nr; i++)
2920 edi->sav.c_ind[i] = i;
2922 /* In the general case we will need a different-sized array for the reference indices: */
2925 snew(edi->sref.c_ind, edi->sref.nr);
2926 for (i = 0; i < edi->sref.nr; i++)
2928 edi->sref.c_ind[i] = i;
2931 /* Point to the very same array in case of other structures: */
2932 edi->star.c_ind = edi->sav.c_ind;
2933 edi->sori.c_ind = edi->sav.c_ind;
2934 /* In the serial case, the local number of atoms is the global one: */
2935 edi->sref.nr_loc = edi->sref.nr;
2936 edi->sav.nr_loc = edi->sav.nr;
2937 edi->star.nr_loc = edi->star.nr;
2938 edi->sori.nr_loc = edi->sori.nr;
2940 /* An on we go to the next ED group */
2941 edi = edi->next_edi;
2945 /* Allocate space for ED buffer variables */
2946 /* Again, loop over ED data sets */
2948 for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
2950 /* Allocate space for ED buffer variables */
2951 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
2952 snew(edi->buf->do_edsam, 1);
2954 /* Space for collective ED buffer variables */
2956 /* Collective positions of atoms with the average indices */
2957 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2958 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
2959 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2960 /* Collective positions of atoms with the reference indices */
2963 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2964 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
2965 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2968 /* Get memory for flooding forces */
2969 snew(edi->flood.forces_cartesian, edi->sav.nr);
2973 /* Dump it all into one file per process */
2974 dump_edi(edi, cr, nr_edi);
2978 edi = edi->next_edi;
2981 /* Flush the edo file so that the user can check some things
2982 * when the simulation has started */
2992 void do_edsam(const t_inputrec *ir,
2994 const t_commrec *cr,
3000 int i, edinr, iupdate = 500;
3001 matrix rotmat; /* rotation matrix */
3002 rvec transvec; /* translation vector */
3003 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3004 real dt_1; /* 1/dt */
3005 struct t_do_edsam *buf;
3007 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3008 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3011 /* Check if ED sampling has to be performed */
3012 if (ed->eEDtype == eEDnone)
3017 dt_1 = 1.0/ir->delta_t;
3019 /* Loop over all ED groups (usually one) */
3022 while (edi != nullptr)
3025 if (bNeedDoEdsam(edi))
3028 buf = edi->buf->do_edsam;
3032 /* initialize radacc radius for slope criterion */
3033 buf->oldrad = calc_radius(&edi->vecs.radacc);
3036 /* Copy the positions into buf->xc* arrays and after ED
3037 * feed back corrections to the official positions */
3039 /* Broadcast the ED positions such that every node has all of them
3040 * Every node contributes its local positions xs and stores it in
3041 * the collective buf->xcoll array. Note that for edinr > 1
3042 * xs could already have been modified by an earlier ED */
3044 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3045 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
3047 /* Only assembly reference positions if their indices differ from the average ones */
3050 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3051 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3054 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3055 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3056 * set bUpdateShifts=TRUE in the parallel case. */
3057 buf->bUpdateShifts = FALSE;
3059 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3060 * as well as the indices in edi->sav.anrs */
3062 /* Fit the reference indices to the reference structure */
3065 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3069 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3072 /* Now apply the translation and rotation to the ED structure */
3073 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3075 /* Find out how well we fit to the reference (just for output steps) */
3076 if (do_per_step(step, edi->outfrq) && MASTER(cr))
3080 /* Indices of reference and average structures are identical,
3081 * thus we can calculate the rmsd to SREF using xcoll */
3082 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3086 /* We have to translate & rotate the reference atoms first */
3087 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3088 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3092 /* update radsam references, when required */
3093 if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3095 project(buf->xcoll, edi);
3096 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3097 rad_project(edi, buf->xcoll, &edi->vecs.radfix);
3098 buf->oldrad = -1.e5;
3101 /* update radacc references, when required */
3102 if (do_per_step(step, iupdate) && step >= edi->presteps)
3104 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
3105 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3107 project(buf->xcoll, edi);
3108 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3113 buf->oldrad = edi->vecs.radacc.radius;
3117 /* apply the constraints */
3118 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3120 /* ED constraints should be applied already in the first MD step
3121 * (which is step 0), therefore we pass step+1 to the routine */
3122 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3125 /* write to edo, when required */
3126 if (do_per_step(step, edi->outfrq))
3128 project(buf->xcoll, edi);
3129 if (MASTER(cr) && !bSuppress)
3131 write_edo(edi, ed->edo, rmsdev);
3135 /* Copy back the positions unless monitoring only */
3136 if (ed_constraints(ed->eEDtype, edi))
3138 /* remove fitting */
3139 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3141 /* Copy the ED corrected positions into the coordinate array */
3142 /* Each node copies its local part. In the serial case, nat_loc is the
3143 * total number of ED atoms */
3144 for (i = 0; i < edi->sav.nr_loc; i++)
3146 /* Unshift local ED coordinate and store in x_unsh */
3147 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3148 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3150 /* dx is the ED correction to the positions: */
3151 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3155 /* dv is the ED correction to the velocity: */
3156 svmul(dt_1, dx, dv);
3157 /* apply the velocity correction: */
3158 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3160 /* Finally apply the position correction due to ED: */
3161 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3164 } /* END of if ( bNeedDoEdsam(edi) ) */
3166 /* Prepare for the next ED group */
3167 edi = edi->next_edi;
3169 } /* END of loop over ED groups */
3174 void done_ed(gmx_edsam_t *ed)
3178 /* ed->edo is opened sometimes with xvgropen, sometimes with
3179 * gmx_fio_fopen, so we use the least common denominator for
3181 gmx_fio_fclose((*ed)->edo);
3184 /* TODO deallocate ed and set pointer to NULL */