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/domdec/domdec_struct.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxlib/network.h"
52 #include "gromacs/gmxlib/nrnb.h"
53 #include "gromacs/linearalgebra/nrjac.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/math/vectypes.h"
58 #include "gromacs/mdlib/broadcaststructs.h"
59 #include "gromacs/mdlib/constr.h"
60 #include "gromacs/mdlib/groupcoord.h"
61 #include "gromacs/mdlib/mdrun.h"
62 #include "gromacs/mdlib/sim_util.h"
63 #include "gromacs/mdlib/update.h"
64 #include "gromacs/mdtypes/commrec.h"
65 #include "gromacs/mdtypes/edsamhistory.h"
66 #include "gromacs/mdtypes/inputrec.h"
67 #include "gromacs/mdtypes/md_enums.h"
68 #include "gromacs/mdtypes/observableshistory.h"
69 #include "gromacs/mdtypes/state.h"
70 #include "gromacs/pbcutil/pbc.h"
71 #include "gromacs/topology/mtop_lookup.h"
72 #include "gromacs/topology/topology.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/gmxassert.h"
76 #include "gromacs/utility/smalloc.h"
79 /* enum to identify the type of ED: none, normal ED, flooding */
81 eEDnone, eEDedsam, eEDflood, eEDnr
84 /* enum to identify operations on reference, average, origin, target structures */
86 eedREF, eedAV, eedORI, eedTAR, eedNR
92 int neig; /* nr of eigenvectors */
93 int *ieig; /* index nrs of eigenvectors */
94 real *stpsz; /* stepsizes (per eigenvector) */
95 rvec **vec; /* eigenvector components */
96 real *xproj; /* instantaneous x projections */
97 real *fproj; /* instantaneous f projections */
98 real radius; /* instantaneous radius */
99 real *refproj; /* starting or target projections */
100 /* When using flooding as harmonic restraint: The current reference projection
101 * is at each step calculated from the initial refproj0 and the slope. */
102 real *refproj0, *refprojslope;
108 t_eigvec mon; /* only monitored, no constraints */
109 t_eigvec linfix; /* fixed linear constraints */
110 t_eigvec linacc; /* acceptance linear constraints */
111 t_eigvec radfix; /* fixed radial constraints (exp) */
112 t_eigvec radacc; /* acceptance radial constraints (exp) */
113 t_eigvec radcon; /* acceptance rad. contraction constr. */
120 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
122 gmx_bool bConstForce; /* Do not calculate a flooding potential,
123 instead flood with a constant force */
132 rvec *forces_cartesian;
133 t_eigvec vecs; /* use flooding for these */
137 /* This type is for the average, reference, target, and origin structure */
138 typedef struct gmx_edx
140 int nr; /* number of atoms this structure contains */
141 int nr_loc; /* number of atoms on local node */
142 int *anrs; /* atom index numbers */
143 int *anrs_loc; /* local atom index numbers */
144 int nalloc_loc; /* allocation size of anrs_loc */
145 int *c_ind; /* at which position of the whole anrs
146 * array is a local atom?, i.e.
147 * c_ind[0...nr_loc-1] gives the atom index
148 * with respect to the collective
149 * anrs[0...nr-1] array */
150 rvec *x; /* positions for this structure */
151 rvec *x_old; /* Last positions which have the correct PBC
152 representation of the ED group. In
153 combination with keeping track of the
154 shift vectors, the ED group can always
156 real *m; /* masses */
157 real mtot; /* total mass (only used in sref) */
158 real *sqrtm; /* sqrt of the masses used for mass-
159 * weighting of analysis (only used in sav) */
165 int nini; /* total Nr of atoms */
166 gmx_bool fitmas; /* true if trans fit with cm */
167 gmx_bool pcamas; /* true if mass-weighted PCA */
168 int presteps; /* number of steps to run without any
169 * perturbations ... just monitoring */
170 int outfrq; /* freq (in steps) of writing to edo */
171 int maxedsteps; /* max nr of steps per cycle */
173 /* all gmx_edx datasets are copied to all nodes in the parallel case */
174 struct gmx_edx sref; /* reference positions, to these fitting
176 gmx_bool bRefEqAv; /* If true, reference & average indices
177 * are the same. Used for optimization */
178 struct gmx_edx sav; /* average positions */
179 struct gmx_edx star; /* target positions */
180 struct gmx_edx sori; /* origin positions */
182 t_edvecs vecs; /* eigenvectors */
183 real slope; /* minimal slope in acceptance radexp */
185 t_edflood flood; /* parameters especially for flooding */
186 struct t_ed_buffer *buf; /* handle to local buffers */
187 struct edpar *next_edi; /* Pointer to another ED group */
191 typedef struct gmx_edsam
193 int eEDtype; /* Type of ED: see enums above */
194 FILE *edo; /* output file pointer */
204 rvec old_transvec, older_transvec, transvec_compact;
205 rvec *xcoll; /* Positions from all nodes, this is the
206 collective set we work on.
207 These are the positions of atoms with
208 average structure indices */
209 rvec *xc_ref; /* same but with reference structure indices */
210 ivec *shifts_xcoll; /* Shifts for xcoll */
211 ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
212 ivec *shifts_xc_ref; /* Shifts for xc_ref */
213 ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
214 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
215 ED shifts for this ED group need to
220 /* definition of ED buffer structure */
223 struct t_fit_to_ref * fit_to_ref;
224 struct t_do_edfit * do_edfit;
225 struct t_do_edsam * do_edsam;
226 struct t_do_radcon * do_radcon;
230 /* Function declarations */
231 static void fit_to_reference(rvec *xcoll, rvec transvec, matrix rotmat, t_edpar *edi);
232 static void translate_and_rotate(rvec *x, int nat, rvec transvec, matrix rotmat);
233 static real rmsd_from_structure(rvec *x, struct gmx_edx *s);
234 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms);
235 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamhistory_t *EDstate);
236 static void init_edsamstate(gmx_edsam_t ed, edsamhistory_t *EDstate);
237 static void write_edo_legend(gmx_edsam_t ed, int nED, const gmx_output_env_t *oenv);
238 /* End function declarations */
240 /* Define formats for the column width in the output file */
241 const char EDcol_sfmt[] = "%17s";
242 const char EDcol_efmt[] = "%17.5e";
243 const char EDcol_ffmt[] = "%17f";
245 /* Define a formats for reading, otherwise cppcheck complains for scanf without width limits */
246 const char max_ev_fmt_d[] = "%7d"; /* Number of eigenvectors. 9,999,999 should be enough */
247 const char max_ev_fmt_lf[] = "%12lf";
248 const char max_ev_fmt_dlf[] = "%7d%12lf";
249 const char max_ev_fmt_dlflflf[] = "%7d%12lf%12lf%12lf";
250 const char max_ev_fmt_lelele[] = "%12le%12le%12le";
252 /* Do we have to perform essential dynamics constraints or possibly only flooding
253 * for any of the ED groups? */
254 static gmx_bool bNeedDoEdsam(t_edpar *edi)
256 return edi->vecs.mon.neig
257 || edi->vecs.linfix.neig
258 || edi->vecs.linacc.neig
259 || edi->vecs.radfix.neig
260 || edi->vecs.radacc.neig
261 || edi->vecs.radcon.neig;
265 /* Multiple ED groups will be labeled with letters instead of numbers
266 * to avoid confusion with eigenvector indices */
267 static char get_EDgroupChar(int nr_edi, int nED)
275 * nr_edi = 2 -> B ...
277 return 'A' + nr_edi - 1;
281 /* Does not subtract average positions, projection on single eigenvector is returned
282 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
283 * Average position is subtracted in ed_apply_constraints prior to calling projectx
285 static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec)
291 for (i = 0; i < edi->sav.nr; i++)
293 proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
300 /* Specialized: projection is stored in vec->refproj
301 * -> used for radacc, radfix, radcon and center of flooding potential
302 * subtracts average positions, projects vector x */
303 static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec)
308 /* Subtract average positions */
309 for (i = 0; i < edi->sav.nr; i++)
311 rvec_dec(x[i], edi->sav.x[i]);
314 for (i = 0; i < vec->neig; i++)
316 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
317 rad += gmx::square((vec->refproj[i]-vec->xproj[i]));
319 vec->radius = sqrt(rad);
321 /* Add average positions */
322 for (i = 0; i < edi->sav.nr; i++)
324 rvec_inc(x[i], edi->sav.x[i]);
329 /* Project vector x, subtract average positions prior to projection and add
330 * them afterwards to retain the unchanged vector. Store in xproj. Mass-weighting
332 static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
333 t_eigvec *vec, /* The eigenvectors */
344 /* Subtract average positions */
345 for (i = 0; i < edi->sav.nr; i++)
347 rvec_dec(x[i], edi->sav.x[i]);
350 for (i = 0; i < vec->neig; i++)
352 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
355 /* Add average positions */
356 for (i = 0; i < edi->sav.nr; i++)
358 rvec_inc(x[i], edi->sav.x[i]);
363 /* Project vector x onto all edi->vecs (mon, linfix,...) */
364 static void project(rvec *x, /* positions to project */
365 t_edpar *edi) /* edi data set */
367 /* It is not more work to subtract the average position in every
368 * subroutine again, because these routines are rarely used simultaneously */
369 project_to_eigvectors(x, &edi->vecs.mon, edi);
370 project_to_eigvectors(x, &edi->vecs.linfix, edi);
371 project_to_eigvectors(x, &edi->vecs.linacc, edi);
372 project_to_eigvectors(x, &edi->vecs.radfix, edi);
373 project_to_eigvectors(x, &edi->vecs.radacc, edi);
374 project_to_eigvectors(x, &edi->vecs.radcon, edi);
378 static real calc_radius(t_eigvec *vec)
384 for (i = 0; i < vec->neig; i++)
386 rad += gmx::square((vec->refproj[i]-vec->xproj[i]));
389 return rad = sqrt(rad);
394 static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
399 fprintf(out, "#%s positions:\n%d\n", name, s->nr);
405 fprintf(out, "#index, x, y, z");
408 fprintf(out, ", sqrt(m)");
410 for (i = 0; i < s->nr; i++)
412 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]);
415 fprintf(out, "%9.3f", s->sqrtm[i]);
423 static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
424 const char name[], int length)
429 fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig);
430 /* Dump the data for every eigenvector: */
431 for (i = 0; i < ev->neig; i++)
433 fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
434 ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
435 for (j = 0; j < length; j++)
437 fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
444 static void dump_edi(t_edpar *edpars, const t_commrec *cr, int nr_edi)
450 sprintf(fn, "EDdump_rank%d_edi%d", cr->nodeid, nr_edi);
451 out = gmx_ffopen(fn, "w");
453 fprintf(out, "#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
454 edpars->nini, edpars->fitmas, edpars->pcamas);
455 fprintf(out, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
456 edpars->outfrq, edpars->maxedsteps, edpars->slope);
457 fprintf(out, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n",
458 edpars->presteps, edpars->flood.deltaF0, edpars->flood.tau,
459 edpars->flood.constEfl, edpars->flood.alpha2);
461 /* Dump reference, average, target, origin positions */
462 dump_edi_positions(out, &edpars->sref, "REFERENCE");
463 dump_edi_positions(out, &edpars->sav, "AVERAGE" );
464 dump_edi_positions(out, &edpars->star, "TARGET" );
465 dump_edi_positions(out, &edpars->sori, "ORIGIN" );
467 /* Dump eigenvectors */
468 dump_edi_eigenvecs(out, &edpars->vecs.mon, "MONITORED", edpars->sav.nr);
469 dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX", edpars->sav.nr);
470 dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC", edpars->sav.nr);
471 dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX", edpars->sav.nr);
472 dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC", edpars->sav.nr);
473 dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON", edpars->sav.nr);
475 /* Dump flooding eigenvectors */
476 dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING", edpars->sav.nr);
478 /* Dump ed local buffer */
479 fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit );
480 fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam );
481 fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon );
492 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
494 /* this is a copy of do_fit with some modifications */
495 int c, r, n, j, i, irot;
496 double d[6], xnr, xpc;
501 struct t_do_edfit *loc;
504 if (edi->buf->do_edfit != nullptr)
511 snew(edi->buf->do_edfit, 1);
513 loc = edi->buf->do_edfit;
517 snew(loc->omega, 2*DIM);
518 snew(loc->om, 2*DIM);
519 for (i = 0; i < 2*DIM; i++)
521 snew(loc->omega[i], 2*DIM);
522 snew(loc->om[i], 2*DIM);
526 for (i = 0; (i < 6); i++)
529 for (j = 0; (j < 6); j++)
531 loc->omega[i][j] = 0;
536 /* calculate the matrix U */
538 for (n = 0; (n < natoms); n++)
540 for (c = 0; (c < DIM); c++)
543 for (r = 0; (r < DIM); r++)
551 /* construct loc->omega */
552 /* loc->omega is symmetric -> loc->omega==loc->omega' */
553 for (r = 0; (r < 6); r++)
555 for (c = 0; (c <= r); c++)
557 if ((r >= 3) && (c < 3))
559 loc->omega[r][c] = u[r-3][c];
560 loc->omega[c][r] = u[r-3][c];
564 loc->omega[r][c] = 0;
565 loc->omega[c][r] = 0;
570 /* determine h and k */
571 jacobi(loc->omega, 6, d, loc->om, &irot);
575 fprintf(stderr, "IROT=0\n");
578 index = 0; /* For the compiler only */
580 for (j = 0; (j < 3); j++)
583 for (i = 0; (i < 6); i++)
592 for (i = 0; (i < 3); i++)
594 vh[j][i] = M_SQRT2*loc->om[i][index];
595 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
600 for (c = 0; (c < 3); c++)
602 for (r = 0; (r < 3); r++)
604 R[c][r] = vk[0][r]*vh[0][c]+
611 for (c = 0; (c < 3); c++)
613 for (r = 0; (r < 3); r++)
615 R[c][r] = vk[0][r]*vh[0][c]+
624 static void rmfit(int nat, rvec *xcoll, const rvec transvec, matrix rotmat)
631 * The inverse rotation is described by the transposed rotation matrix */
632 transpose(rotmat, tmat);
633 rotate_x(xcoll, nat, tmat);
635 /* Remove translation */
636 vec[XX] = -transvec[XX];
637 vec[YY] = -transvec[YY];
638 vec[ZZ] = -transvec[ZZ];
639 translate_x(xcoll, nat, vec);
643 /**********************************************************************************
644 ******************** FLOODING ****************************************************
645 **********************************************************************************
647 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
648 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
649 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
651 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
652 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
654 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
655 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
657 do_flood makes a copy of the positions,
658 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
659 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
660 fit. Then do_flood adds these forces to the forcefield-forces
661 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
663 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
664 structure is projected to the system of eigenvectors and then this position in the subspace is used as
665 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
666 i.e. the average structure as given in the make_edi file.
668 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
669 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
670 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
671 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
672 further adaption is applied, Efl will stay constant at zero.
674 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
675 used as spring constants for the harmonic potential.
676 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
677 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
679 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
680 the routine read_edi_file reads all of theses flooding files.
681 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
682 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
683 edi there is no interdependence whatsoever. The forces are added together.
685 To write energies into the .edr file, call the function
686 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
688 get_flood_energies(real Vfl[],int nnames);
691 - 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.
693 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
694 two edsam files from two peptide chains
697 static void write_edo_flood(t_edpar *edi, FILE *fp, real rmsd)
702 /* Output how well we fit to the reference structure */
703 fprintf(fp, EDcol_ffmt, rmsd);
705 for (i = 0; i < edi->flood.vecs.neig; i++)
707 fprintf(fp, EDcol_efmt, edi->flood.vecs.xproj[i]);
709 /* Check whether the reference projection changes with time (this can happen
710 * in case flooding is used as harmonic restraint). If so, output the
711 * current reference projection */
712 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
714 fprintf(fp, EDcol_efmt, edi->flood.vecs.refproj[i]);
717 /* Output Efl if we are doing adaptive flooding */
718 if (0 != edi->flood.tau)
720 fprintf(fp, EDcol_efmt, edi->flood.Efl);
722 fprintf(fp, EDcol_efmt, edi->flood.Vfl);
724 /* Output deltaF if we are doing adaptive flooding */
725 if (0 != edi->flood.tau)
727 fprintf(fp, EDcol_efmt, edi->flood.deltaF);
729 fprintf(fp, EDcol_efmt, edi->flood.vecs.fproj[i]);
734 /* From flood.xproj compute the Vfl(x) at this point */
735 static real flood_energy(t_edpar *edi, gmx_int64_t step)
737 /* compute flooding energy Vfl
738 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
739 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
740 it is already computed by make_edi and stored in stpsz[i]
742 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
749 /* Each time this routine is called (i.e. each time step), we add a small
750 * value to the reference projection. This way a harmonic restraint towards
751 * a moving reference is realized. If no value for the additive constant
752 * is provided in the edi file, the reference will not change. */
753 if (edi->flood.bHarmonic)
755 for (i = 0; i < edi->flood.vecs.neig; i++)
757 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
762 /* Compute sum which will be the exponent of the exponential */
763 for (i = 0; i < edi->flood.vecs.neig; i++)
765 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
766 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]);
769 /* Compute the Gauss function*/
770 if (edi->flood.bHarmonic)
772 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
776 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*std::exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
783 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
784 static void flood_forces(t_edpar *edi)
786 /* compute the forces in the subspace of the flooding eigenvectors
787 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
790 real energy = edi->flood.Vfl;
793 if (edi->flood.bHarmonic)
795 for (i = 0; i < edi->flood.vecs.neig; i++)
797 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
802 for (i = 0; i < edi->flood.vecs.neig; i++)
804 /* if Efl is zero the forces are zero if not use the formula */
805 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;
811 /* Raise forces from subspace into cartesian space */
812 static void flood_blowup(t_edpar *edi, rvec *forces_cart)
814 /* this function lifts the forces from the subspace to the cartesian space
815 all the values not contained in the subspace are assumed to be zero and then
816 a coordinate transformation from eigenvector to cartesian vectors is performed
817 The nonexistent values don't have to be set to zero explicitly, they would occur
818 as zero valued summands, hence we just stop to compute this part of the sum.
820 for every atom we add all the contributions to this atom from all the different eigenvectors.
822 NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
823 field forces_cart prior the computation, but we compute the forces separately
824 to have them accessible for diagnostics
831 forces_sub = edi->flood.vecs.fproj;
834 /* Calculate the cartesian forces for the local atoms */
836 /* Clear forces first */
837 for (j = 0; j < edi->sav.nr_loc; j++)
839 clear_rvec(forces_cart[j]);
842 /* Now compute atomwise */
843 for (j = 0; j < edi->sav.nr_loc; j++)
845 /* Compute forces_cart[edi->sav.anrs[j]] */
846 for (eig = 0; eig < edi->flood.vecs.neig; eig++)
848 /* Force vector is force * eigenvector (compute only atom j) */
849 svmul(forces_sub[eig], edi->flood.vecs.vec[eig][edi->sav.c_ind[j]], dum);
850 /* Add this vector to the cartesian forces */
851 rvec_inc(forces_cart[j], dum);
857 /* Update the values of Efl, deltaF depending on tau and Vfl */
858 static void update_adaption(t_edpar *edi)
860 /* this function updates the parameter Efl and deltaF according to the rules given in
861 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
864 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
866 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
867 /* check if restrain (inverted flooding) -> don't let EFL become positive */
868 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
873 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
878 static void do_single_flood(
886 gmx_bool bNS) /* Are we in a neighbor searching step? */
889 matrix rotmat; /* rotation matrix */
890 matrix tmat; /* inverse rotation */
891 rvec transvec; /* translation vector */
893 struct t_do_edsam *buf;
896 buf = edi->buf->do_edsam;
899 /* Broadcast the positions of the AVERAGE structure such that they are known on
900 * every processor. Each node contributes its local positions x and stores them in
901 * the collective ED array buf->xcoll */
902 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
903 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
905 /* Only assembly REFERENCE positions if their indices differ from the average ones */
908 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
909 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
912 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
913 * We do not need to update the shifts until the next NS step */
914 buf->bUpdateShifts = FALSE;
916 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
917 * as well as the indices in edi->sav.anrs */
919 /* Fit the reference indices to the reference structure */
922 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
926 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
929 /* Now apply the translation and rotation to the ED structure */
930 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
932 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
933 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, edi);
935 if (FALSE == edi->flood.bConstForce)
937 /* Compute Vfl(x) from flood.xproj */
938 edi->flood.Vfl = flood_energy(edi, step);
940 update_adaption(edi);
942 /* Compute the flooding forces */
946 /* Translate them into cartesian positions */
947 flood_blowup(edi, edi->flood.forces_cartesian);
949 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
950 /* Each node rotates back its local forces */
951 transpose(rotmat, tmat);
952 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
954 /* Finally add forces to the main force variable */
955 for (i = 0; i < edi->sav.nr_loc; i++)
957 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
960 /* Output is written by the master process */
961 if (do_per_step(step, edi->outfrq) && MASTER(cr))
963 /* Output how well we fit to the reference */
966 /* Indices of reference and average structures are identical,
967 * thus we can calculate the rmsd to SREF using xcoll */
968 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
972 /* We have to translate & rotate the reference atoms first */
973 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
974 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
977 write_edo_flood(edi, edo, rmsdev);
982 /* Main flooding routine, called from do_force */
983 extern void do_flood(const t_commrec *cr,
984 const t_inputrec *ir,
997 /* Write time to edo, when required. Output the time anyhow since we need
998 * it in the output file for ED constraints. */
999 if (MASTER(cr) && do_per_step(step, edi->outfrq))
1001 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
1004 if (ed->eEDtype != eEDflood)
1011 /* Call flooding for one matrix */
1012 if (edi->flood.vecs.neig)
1014 do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS);
1016 edi = edi->next_edi;
1021 /* Called by init_edi, configure some flooding related variables and structures,
1022 * print headers to output files */
1023 static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt)
1028 edi->flood.Efl = edi->flood.constEfl;
1032 if (edi->flood.vecs.neig)
1034 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1035 ed->eEDtype = eEDflood;
1037 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
1039 if (edi->flood.bConstForce)
1041 /* We have used stpsz as a vehicle to carry the fproj values for constant
1042 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1043 * in const force flooding, fproj is never changed. */
1044 for (i = 0; i < edi->flood.vecs.neig; i++)
1046 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1048 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1049 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1057 /*********** Energy book keeping ******/
1058 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1067 srenew(names, count);
1068 sprintf(buf, "Vfl_%d", count);
1069 names[count-1] = gmx_strdup(buf);
1070 actual = actual->next_edi;
1077 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1079 /*fl has to be big enough to capture nnames-many entries*/
1088 Vfl[count-1] = actual->flood.Vfl;
1089 actual = actual->next_edi;
1092 if (nnames != count-1)
1094 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1097 /************* END of FLOODING IMPLEMENTATION ****************************/
1101 /* This function opens the ED input and output files, reads in all datasets it finds
1102 * in the input file, and cross-checks whether the .edi file information is consistent
1103 * with the essential dynamics data found in the checkpoint file (if present).
1104 * gmx make_edi can be used to create an .edi input file.
1106 static gmx_edsam_t ed_open(
1108 ObservablesHistory *oh,
1109 const char *ediFileName,
1110 const char *edoFileName,
1112 const gmx_output_env_t *oenv,
1113 const t_commrec *cr)
1118 /* Allocate space for the ED data structure */
1121 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1122 ed->eEDtype = eEDedsam;
1128 // If we start from a checkpoint file, we already have an edsamHistory struct
1129 if (oh->edsamHistory == nullptr)
1131 oh->edsamHistory = std::unique_ptr<edsamhistory_t>(new edsamhistory_t {});
1133 edsamhistory_t *EDstate = oh->edsamHistory.get();
1135 /* Read the edi input file: */
1136 nED = read_edi_file(ediFileName, ed->edpar, natoms);
1138 /* Make sure the checkpoint was produced in a run using this .edi file */
1139 if (EDstate->bFromCpt)
1141 crosscheck_edi_file_vs_checkpoint(ed, EDstate);
1147 init_edsamstate(ed, EDstate);
1149 /* The master opens the ED output file */
1152 ed->edo = gmx_fio_fopen(edoFileName, "a+");
1156 ed->edo = xvgropen(edoFileName,
1157 "Essential dynamics / flooding output",
1159 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1161 /* Make a descriptive legend */
1162 write_edo_legend(ed, EDstate->nED, oenv);
1169 /* Broadcasts the structure data */
1170 static void bc_ed_positions(const t_commrec *cr, struct gmx_edx *s, int stype)
1172 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1173 snew_bc(cr, s->x, s->nr ); /* Positions */
1174 nblock_bc(cr, s->nr, s->anrs );
1175 nblock_bc(cr, s->nr, s->x );
1177 /* For the average & reference structures we need an array for the collective indices,
1178 * and we need to broadcast the masses as well */
1179 if (stype == eedAV || stype == eedREF)
1181 /* We need these additional variables in the parallel case: */
1182 snew(s->c_ind, s->nr ); /* Collective indices */
1183 /* Local atom indices get assigned in dd_make_local_group_indices.
1184 * There, also memory is allocated */
1185 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1186 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1187 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1190 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1191 if (stype == eedREF)
1193 snew_bc(cr, s->m, s->nr);
1194 nblock_bc(cr, s->nr, s->m);
1197 /* For the average structure we might need the masses for mass-weighting */
1200 snew_bc(cr, s->sqrtm, s->nr);
1201 nblock_bc(cr, s->nr, s->sqrtm);
1202 snew_bc(cr, s->m, s->nr);
1203 nblock_bc(cr, s->nr, s->m);
1208 /* Broadcasts the eigenvector data */
1209 static void bc_ed_vecs(const t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
1213 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1214 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1215 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1216 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1217 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1219 nblock_bc(cr, ev->neig, ev->ieig );
1220 nblock_bc(cr, ev->neig, ev->stpsz );
1221 nblock_bc(cr, ev->neig, ev->xproj );
1222 nblock_bc(cr, ev->neig, ev->fproj );
1223 nblock_bc(cr, ev->neig, ev->refproj);
1225 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1226 for (i = 0; i < ev->neig; i++)
1228 snew_bc(cr, ev->vec[i], length);
1229 nblock_bc(cr, length, ev->vec[i]);
1232 /* For harmonic restraints the reference projections can change with time */
1235 snew_bc(cr, ev->refproj0, ev->neig);
1236 snew_bc(cr, ev->refprojslope, ev->neig);
1237 nblock_bc(cr, ev->neig, ev->refproj0 );
1238 nblock_bc(cr, ev->neig, ev->refprojslope);
1243 /* Broadcasts the ED / flooding data to other nodes
1244 * and allocates memory where needed */
1245 static void broadcast_ed_data(const t_commrec *cr, gmx_edsam_t ed, int numedis)
1251 /* Master lets the other nodes know if its ED only or also flooding */
1252 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1254 snew_bc(cr, ed->edpar, 1);
1255 /* Now transfer the ED data set(s) */
1257 for (nr = 0; nr < numedis; nr++)
1259 /* Broadcast a single ED data set */
1262 /* Broadcast positions */
1263 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1264 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1265 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1266 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1268 /* Broadcast eigenvectors */
1269 bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr, FALSE);
1270 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
1271 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
1272 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
1273 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
1274 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
1275 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1276 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
1278 /* Set the pointer to the next ED group */
1281 snew_bc(cr, edi->next_edi, 1);
1282 edi = edi->next_edi;
1288 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1289 static void init_edi(const gmx_mtop_t *mtop, t_edpar *edi)
1292 real totalmass = 0.0;
1295 /* NOTE Init_edi is executed on the master process only
1296 * The initialized data sets are then transmitted to the
1297 * other nodes in broadcast_ed_data */
1299 /* evaluate masses (reference structure) */
1300 snew(edi->sref.m, edi->sref.nr);
1302 for (i = 0; i < edi->sref.nr; i++)
1306 edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1310 edi->sref.m[i] = 1.0;
1313 /* Check that every m > 0. Bad things will happen otherwise. */
1314 if (edi->sref.m[i] <= 0.0)
1316 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1317 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1318 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1319 "atoms from the reference structure by creating a proper index group.\n",
1320 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1323 totalmass += edi->sref.m[i];
1325 edi->sref.mtot = totalmass;
1327 /* Masses m and sqrt(m) for the average structure. Note that m
1328 * is needed if forces have to be evaluated in do_edsam */
1329 snew(edi->sav.sqrtm, edi->sav.nr );
1330 snew(edi->sav.m, edi->sav.nr );
1331 for (i = 0; i < edi->sav.nr; i++)
1333 edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1336 edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1340 edi->sav.sqrtm[i] = 1.0;
1343 /* Check that every m > 0. Bad things will happen otherwise. */
1344 if (edi->sav.sqrtm[i] <= 0.0)
1346 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1347 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1348 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1349 "atoms from the average structure by creating a proper index group.\n",
1350 i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1354 /* put reference structure in origin */
1355 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1359 translate_x(edi->sref.x, edi->sref.nr, com);
1361 /* Init ED buffer */
1366 static void check(const char *line, const char *label)
1368 if (!strstr(line, label))
1370 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1375 static int read_checked_edint(FILE *file, const char *label)
1377 char line[STRLEN+1];
1380 fgets2 (line, STRLEN, file);
1382 fgets2 (line, STRLEN, file);
1383 sscanf (line, max_ev_fmt_d, &idum);
1388 static int read_edint(FILE *file, gmx_bool *bEOF)
1390 char line[STRLEN+1];
1395 eof = fgets2 (line, STRLEN, file);
1401 eof = fgets2 (line, STRLEN, file);
1407 sscanf (line, max_ev_fmt_d, &idum);
1413 static real read_checked_edreal(FILE *file, const char *label)
1415 char line[STRLEN+1];
1419 fgets2 (line, STRLEN, file);
1421 fgets2 (line, STRLEN, file);
1422 sscanf (line, max_ev_fmt_lf, &rdum);
1423 return (real) rdum; /* always read as double and convert to single */
1427 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1430 char line[STRLEN+1];
1434 for (i = 0; i < number; i++)
1436 fgets2 (line, STRLEN, file);
1437 sscanf (line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1438 anrs[i]--; /* we are reading FORTRAN indices */
1439 for (j = 0; j < 3; j++)
1441 x[i][j] = d[j]; /* always read as double and convert to single */
1447 static void scan_edvec(FILE *in, int nr, rvec *vec)
1449 char line[STRLEN+1];
1454 for (i = 0; (i < nr); i++)
1456 fgets2 (line, STRLEN, in);
1457 sscanf (line, max_ev_fmt_lelele, &x, &y, &z);
1465 static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1468 double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0;
1469 char line[STRLEN+1];
1472 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1475 snew(tvec->ieig, tvec->neig);
1476 snew(tvec->stpsz, tvec->neig);
1477 snew(tvec->vec, tvec->neig);
1478 snew(tvec->xproj, tvec->neig);
1479 snew(tvec->fproj, tvec->neig);
1480 snew(tvec->refproj, tvec->neig);
1483 snew(tvec->refproj0, tvec->neig);
1484 snew(tvec->refprojslope, tvec->neig);
1487 for (i = 0; (i < tvec->neig); i++)
1489 fgets2 (line, STRLEN, in);
1490 if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1492 nscan = sscanf(line, max_ev_fmt_dlflflf, &idum, &rdum, &refproj_dum, &refprojslope_dum);
1493 /* Zero out values which were not scanned */
1497 /* Every 4 values read, including reference position */
1498 *bHaveReference = TRUE;
1501 /* A reference position is provided */
1502 *bHaveReference = TRUE;
1503 /* No value for slope, set to 0 */
1504 refprojslope_dum = 0.0;
1507 /* No values for reference projection and slope, set to 0 */
1509 refprojslope_dum = 0.0;
1512 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");
1583 static int read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, const char *fn)
1586 const int magic = 670;
1589 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1590 gmx_bool bHaveReference = FALSE;
1593 /* the edi file is not free format, so expect problems if the input is corrupt. */
1595 /* check the magic number */
1596 readmagic = read_edint(in, &bEOF);
1597 /* Check whether we have reached the end of the input file */
1603 if (readmagic != magic)
1605 if (readmagic == 666 || readmagic == 667 || readmagic == 668)
1607 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1609 else if (readmagic != 669)
1611 gmx_fatal(FARGS, "Wrong magic number %d in %s", readmagic, fn);
1615 /* check the number of atoms */
1616 edi->nini = read_edint(in, &bEOF);
1617 if (edi->nini != nr_mdatoms)
1619 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1622 /* Done checking. For the rest we blindly trust the input */
1623 edi->fitmas = read_checked_edint(in, "FITMAS");
1624 edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS");
1625 edi->outfrq = read_checked_edint(in, "OUTFRQ");
1626 edi->maxedsteps = read_checked_edint(in, "MAXLEN");
1627 edi->slope = read_checked_edreal(in, "SLOPECRIT");
1629 edi->presteps = read_checked_edint(in, "PRESTEPS");
1630 edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1631 edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1632 edi->flood.tau = read_checked_edreal(in, "TAU");
1633 edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1634 edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1635 edi->flood.kT = read_checked_edreal(in, "KT");
1636 edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1637 if (readmagic > 669)
1639 edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1643 edi->flood.bConstForce = FALSE;
1645 edi->sref.nr = read_checked_edint(in, "NREF");
1647 /* allocate space for reference positions and read them */
1648 snew(edi->sref.anrs, edi->sref.nr);
1649 snew(edi->sref.x, edi->sref.nr);
1650 snew(edi->sref.x_old, edi->sref.nr);
1651 edi->sref.sqrtm = nullptr;
1652 read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1654 /* average positions. they define which atoms will be used for ED sampling */
1655 edi->sav.nr = read_checked_edint(in, "NAV");
1656 snew(edi->sav.anrs, edi->sav.nr);
1657 snew(edi->sav.x, edi->sav.nr);
1658 snew(edi->sav.x_old, edi->sav.nr);
1659 read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1661 /* Check if the same atom indices are used for reference and average positions */
1662 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1665 read_edvecs(in, edi->sav.nr, &edi->vecs);
1666 read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
1668 /* target positions */
1669 edi->star.nr = read_edint(in, &bEOF);
1670 if (edi->star.nr > 0)
1672 snew(edi->star.anrs, edi->star.nr);
1673 snew(edi->star.x, edi->star.nr);
1674 edi->star.sqrtm = nullptr;
1675 read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1678 /* positions defining origin of expansion circle */
1679 edi->sori.nr = read_edint(in, &bEOF);
1680 if (edi->sori.nr > 0)
1684 /* Both an -ori structure and a at least one manual reference point have been
1685 * specified. That's ambiguous and probably not intentional. */
1686 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1687 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1689 snew(edi->sori.anrs, edi->sori.nr);
1690 snew(edi->sori.x, edi->sori.nr);
1691 edi->sori.sqrtm = nullptr;
1692 read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1701 /* Read in the edi input file. Note that it may contain several ED data sets which were
1702 * achieved by concatenating multiple edi files. The standard case would be a single ED
1703 * data set, though. */
1704 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1707 t_edpar *curr_edi, *last_edi;
1712 /* This routine is executed on the master only */
1714 /* Open the .edi parameter input file */
1715 in = gmx_fio_fopen(fn, "r");
1716 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1718 /* Now read a sequence of ED input parameter sets from the edi file */
1721 while (read_edi(in, curr_edi, nr_mdatoms, fn) )
1725 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1726 /* We need to allocate space for the data: */
1728 /* Point the 'next_edi' entry to the next edi: */
1729 curr_edi->next_edi = edi_read;
1730 /* Keep the curr_edi pointer for the case that the next group is empty: */
1731 last_edi = curr_edi;
1732 /* Let's prepare to read in the next edi data set: */
1733 curr_edi = edi_read;
1737 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1740 /* Terminate the edi group list with a NULL pointer: */
1741 last_edi->next_edi = nullptr;
1743 fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1745 /* Close the .edi file again */
1752 struct t_fit_to_ref {
1753 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1756 /* Fit the current positions to the reference positions
1757 * Do not actually do the fit, just return rotation and translation.
1758 * Note that the COM of the reference structure was already put into
1759 * the origin by init_edi. */
1760 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1761 rvec transvec, /* The translation vector */
1762 matrix rotmat, /* The rotation matrix */
1763 t_edpar *edi) /* Just needed for do_edfit */
1765 rvec com; /* center of mass */
1767 struct t_fit_to_ref *loc;
1770 /* Allocate memory the first time this routine is called for each edi group */
1771 if (nullptr == edi->buf->fit_to_ref)
1773 snew(edi->buf->fit_to_ref, 1);
1774 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1776 loc = edi->buf->fit_to_ref;
1778 /* We do not touch the original positions but work on a copy. */
1779 for (i = 0; i < edi->sref.nr; i++)
1781 copy_rvec(xcoll[i], loc->xcopy[i]);
1784 /* Calculate the center of mass */
1785 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1787 transvec[XX] = -com[XX];
1788 transvec[YY] = -com[YY];
1789 transvec[ZZ] = -com[ZZ];
1791 /* Subtract the center of mass from the copy */
1792 translate_x(loc->xcopy, edi->sref.nr, transvec);
1794 /* Determine the rotation matrix */
1795 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1799 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1800 int nat, /* How many positions are there? */
1801 rvec transvec, /* The translation vector */
1802 matrix rotmat) /* The rotation matrix */
1805 translate_x(x, nat, transvec);
1808 rotate_x(x, nat, rotmat);
1812 /* Gets the rms deviation of the positions to the structure s */
1813 /* fit_to_structure has to be called before calling this routine! */
1814 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1815 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1821 for (i = 0; i < s->nr; i++)
1823 rmsd += distance2(s->x[i], x[i]);
1826 rmsd /= (real) s->nr;
1833 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1838 if (ed->eEDtype != eEDnone)
1840 /* Loop over ED groups */
1844 /* Local atoms of the reference structure (for fitting), need only be assembled
1845 * if their indices differ from the average ones */
1848 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1849 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1852 /* Local atoms of the average structure (on these ED will be performed) */
1853 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1854 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1856 /* Indicate that the ED shift vectors for this structure need to be updated
1857 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1858 edi->buf->do_edsam->bUpdateShifts = TRUE;
1860 /* Set the pointer to the next ED group (if any) */
1861 edi = edi->next_edi;
1867 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1878 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1879 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1880 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1884 xu[XX] = x[XX]-tx*box[XX][XX];
1885 xu[YY] = x[YY]-ty*box[YY][YY];
1886 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1891 static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
1898 /* loop over linfix vectors */
1899 for (i = 0; i < edi->vecs.linfix.neig; i++)
1901 /* calculate the projection */
1902 proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1904 /* calculate the correction */
1905 add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1907 /* apply the correction */
1908 add /= edi->sav.sqrtm[i];
1909 for (j = 0; j < edi->sav.nr; j++)
1911 svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1912 rvec_inc(xcoll[j], vec_dum);
1918 static void do_linacc(rvec *xcoll, t_edpar *edi)
1925 /* loop over linacc vectors */
1926 for (i = 0; i < edi->vecs.linacc.neig; i++)
1928 /* calculate the projection */
1929 proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1931 /* calculate the correction */
1933 if (edi->vecs.linacc.stpsz[i] > 0.0)
1935 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1937 add = edi->vecs.linacc.refproj[i] - proj;
1940 if (edi->vecs.linacc.stpsz[i] < 0.0)
1942 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1944 add = edi->vecs.linacc.refproj[i] - proj;
1948 /* apply the correction */
1949 add /= edi->sav.sqrtm[i];
1950 for (j = 0; j < edi->sav.nr; j++)
1952 svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
1953 rvec_inc(xcoll[j], vec_dum);
1956 /* new positions will act as reference */
1957 edi->vecs.linacc.refproj[i] = proj + add;
1962 static void do_radfix(rvec *xcoll, t_edpar *edi)
1965 real *proj, rad = 0.0, ratio;
1969 if (edi->vecs.radfix.neig == 0)
1974 snew(proj, edi->vecs.radfix.neig);
1976 /* loop over radfix vectors */
1977 for (i = 0; i < edi->vecs.radfix.neig; i++)
1979 /* calculate the projections, radius */
1980 proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
1981 rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
1985 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
1986 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
1988 /* loop over radfix vectors */
1989 for (i = 0; i < edi->vecs.radfix.neig; i++)
1991 proj[i] -= edi->vecs.radfix.refproj[i];
1993 /* apply the correction */
1994 proj[i] /= edi->sav.sqrtm[i];
1996 for (j = 0; j < edi->sav.nr; j++)
1998 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
1999 rvec_inc(xcoll[j], vec_dum);
2007 static void do_radacc(rvec *xcoll, t_edpar *edi)
2010 real *proj, rad = 0.0, ratio = 0.0;
2014 if (edi->vecs.radacc.neig == 0)
2019 snew(proj, edi->vecs.radacc.neig);
2021 /* loop over radacc vectors */
2022 for (i = 0; i < edi->vecs.radacc.neig; i++)
2024 /* calculate the projections, radius */
2025 proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2026 rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2030 /* only correct when radius decreased */
2031 if (rad < edi->vecs.radacc.radius)
2033 ratio = edi->vecs.radacc.radius/rad - 1.0;
2037 edi->vecs.radacc.radius = rad;
2040 /* loop over radacc vectors */
2041 for (i = 0; i < edi->vecs.radacc.neig; i++)
2043 proj[i] -= edi->vecs.radacc.refproj[i];
2045 /* apply the correction */
2046 proj[i] /= edi->sav.sqrtm[i];
2048 for (j = 0; j < edi->sav.nr; j++)
2050 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2051 rvec_inc(xcoll[j], vec_dum);
2058 struct t_do_radcon {
2062 static void do_radcon(rvec *xcoll, t_edpar *edi)
2065 real rad = 0.0, ratio = 0.0;
2066 struct t_do_radcon *loc;
2071 if (edi->buf->do_radcon != nullptr)
2078 snew(edi->buf->do_radcon, 1);
2080 loc = edi->buf->do_radcon;
2082 if (edi->vecs.radcon.neig == 0)
2089 snew(loc->proj, edi->vecs.radcon.neig);
2092 /* loop over radcon vectors */
2093 for (i = 0; i < edi->vecs.radcon.neig; i++)
2095 /* calculate the projections, radius */
2096 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2097 rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2100 /* only correct when radius increased */
2101 if (rad > edi->vecs.radcon.radius)
2103 ratio = edi->vecs.radcon.radius/rad - 1.0;
2105 /* loop over radcon vectors */
2106 for (i = 0; i < edi->vecs.radcon.neig; i++)
2108 /* apply the correction */
2109 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2110 loc->proj[i] /= edi->sav.sqrtm[i];
2111 loc->proj[i] *= ratio;
2113 for (j = 0; j < edi->sav.nr; j++)
2115 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2116 rvec_inc(xcoll[j], vec_dum);
2123 edi->vecs.radcon.radius = rad;
2129 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
2134 /* subtract the average positions */
2135 for (i = 0; i < edi->sav.nr; i++)
2137 rvec_dec(xcoll[i], edi->sav.x[i]);
2140 /* apply the constraints */
2143 do_linfix(xcoll, edi, step);
2145 do_linacc(xcoll, edi);
2148 do_radfix(xcoll, edi);
2150 do_radacc(xcoll, edi);
2151 do_radcon(xcoll, edi);
2153 /* add back the average positions */
2154 for (i = 0; i < edi->sav.nr; i++)
2156 rvec_inc(xcoll[i], edi->sav.x[i]);
2161 /* Write out the projections onto the eigenvectors. The order of output
2162 * corresponds to ed_output_legend() */
2163 static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
2168 /* Output how well we fit to the reference structure */
2169 fprintf(fp, EDcol_ffmt, rmsd);
2171 for (i = 0; i < edi->vecs.mon.neig; i++)
2173 fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
2176 for (i = 0; i < edi->vecs.linfix.neig; i++)
2178 fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
2181 for (i = 0; i < edi->vecs.linacc.neig; i++)
2183 fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
2186 for (i = 0; i < edi->vecs.radfix.neig; i++)
2188 fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
2190 if (edi->vecs.radfix.neig)
2192 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
2195 for (i = 0; i < edi->vecs.radacc.neig; i++)
2197 fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
2199 if (edi->vecs.radacc.neig)
2201 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
2204 for (i = 0; i < edi->vecs.radcon.neig; i++)
2206 fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
2208 if (edi->vecs.radcon.neig)
2210 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
2214 /* Returns if any constraints are switched on */
2215 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2217 if (edtype == eEDedsam || edtype == eEDflood)
2219 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2220 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2221 edi->vecs.radcon.neig);
2227 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2228 * umbrella sampling simulations. */
2229 static void copyEvecReference(t_eigvec* floodvecs)
2234 if (nullptr == floodvecs->refproj0)
2236 snew(floodvecs->refproj0, floodvecs->neig);
2239 for (i = 0; i < floodvecs->neig; i++)
2241 floodvecs->refproj0[i] = floodvecs->refproj[i];
2246 /* Call on MASTER only. Check whether the essential dynamics / flooding
2247 * groups of the checkpoint file are consistent with the provided .edi file. */
2248 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamhistory_t *EDstate)
2250 t_edpar *edi = nullptr; /* points to a single edi data set */
2254 if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2256 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2257 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2258 "it must also continue with/without ED constraints when checkpointing.\n"
2259 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2260 "from without a checkpoint.\n");
2265 while (edi != nullptr)
2267 /* Check number of atoms in the reference and average structures */
2268 if (EDstate->nref[edinum] != edi->sref.nr)
2270 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2271 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2272 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2274 if (EDstate->nav[edinum] != edi->sav.nr)
2276 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2277 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2278 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2280 edi = edi->next_edi;
2284 if (edinum != EDstate->nED)
2286 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2287 "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2288 "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2293 /* The edsamstate struct stores the information we need to make the ED group
2294 * whole again after restarts from a checkpoint file. Here we do the following:
2295 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2296 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2297 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2298 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2299 * all ED structures at the last time step. */
2300 static void init_edsamstate(gmx_edsam_t ed, edsamhistory_t *EDstate)
2306 snew(EDstate->old_sref_p, EDstate->nED);
2307 snew(EDstate->old_sav_p, EDstate->nED);
2309 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2310 if (!EDstate->bFromCpt)
2312 snew(EDstate->nref, EDstate->nED);
2313 snew(EDstate->nav, EDstate->nED);
2316 /* Loop over all ED/flooding data sets (usually only one, though) */
2318 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2320 /* We always need the last reference and average positions such that
2321 * in the next time step we can make the ED group whole again
2322 * if the atoms do not have the correct PBC representation */
2323 if (EDstate->bFromCpt)
2325 /* Copy the last whole positions of reference and average group from .cpt */
2326 for (i = 0; i < edi->sref.nr; i++)
2328 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2330 for (i = 0; i < edi->sav.nr; i++)
2332 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2337 EDstate->nref[nr_edi-1] = edi->sref.nr;
2338 EDstate->nav [nr_edi-1] = edi->sav.nr;
2341 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2342 EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2343 EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2345 edi = edi->next_edi;
2350 /* Adds 'buf' to 'str' */
2351 static void add_to_string(char **str, char *buf)
2356 len = strlen(*str) + strlen(buf) + 1;
2362 static void add_to_string_aligned(char **str, char *buf)
2364 char buf_aligned[STRLEN];
2366 sprintf(buf_aligned, EDcol_sfmt, buf);
2367 add_to_string(str, buf_aligned);
2371 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, const char *value, const char *unit, char EDgroupchar)
2373 char tmp[STRLEN], tmp2[STRLEN];
2376 sprintf(tmp, "%c %s", EDgroupchar, value);
2377 add_to_string_aligned(LegendStr, tmp);
2378 sprintf(tmp2, "%s (%s)", tmp, unit);
2379 (*setname)[*nsets] = gmx_strdup(tmp2);
2384 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2390 for (i = 0; i < evec->neig; i++)
2392 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2393 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2398 /* Makes a legend for the xvg output file. Call on MASTER only! */
2399 static void write_edo_legend(gmx_edsam_t ed, int nED, const gmx_output_env_t *oenv)
2401 t_edpar *edi = nullptr;
2403 int nr_edi, nsets, n_flood, n_edsam;
2404 const char **setname;
2406 char *LegendStr = nullptr;
2411 fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2413 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2415 fprintf(ed->edo, "#\n");
2416 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2417 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2418 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2419 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2420 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2421 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2422 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2423 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2424 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2426 if (edi->flood.vecs.neig)
2428 /* If in any of the groups we find a flooding vector, flooding is turned on */
2429 ed->eEDtype = eEDflood;
2431 /* Print what flavor of flooding we will do */
2432 if (0 == edi->flood.tau) /* constant flooding strength */
2434 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2435 if (edi->flood.bHarmonic)
2437 fprintf(ed->edo, ", harmonic");
2440 else /* adaptive flooding */
2442 fprintf(ed->edo, ", adaptive");
2445 fprintf(ed->edo, "\n");
2447 edi = edi->next_edi;
2450 /* Print a nice legend */
2452 LegendStr[0] = '\0';
2453 sprintf(buf, "# %6s", "time");
2454 add_to_string(&LegendStr, buf);
2456 /* Calculate the maximum number of columns we could end up with */
2459 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2461 nsets += 5 +edi->vecs.mon.neig
2462 +edi->vecs.linfix.neig
2463 +edi->vecs.linacc.neig
2464 +edi->vecs.radfix.neig
2465 +edi->vecs.radacc.neig
2466 +edi->vecs.radcon.neig
2467 + 6*edi->flood.vecs.neig;
2468 edi = edi->next_edi;
2470 snew(setname, nsets);
2472 /* In the mdrun time step in a first function call (do_flood()) the flooding
2473 * forces are calculated and in a second function call (do_edsam()) the
2474 * ED constraints. To get a corresponding legend, we need to loop twice
2475 * over the edi groups and output first the flooding, then the ED part */
2477 /* The flooding-related legend entries, if flooding is done */
2479 if (eEDflood == ed->eEDtype)
2482 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2484 /* Always write out the projection on the flooding EVs. Of course, this can also
2485 * be achieved with the monitoring option in do_edsam() (if switched on by the
2486 * user), but in that case the positions need to be communicated in do_edsam(),
2487 * which is not necessary when doing flooding only. */
2488 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2490 for (i = 0; i < edi->flood.vecs.neig; i++)
2492 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2493 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2495 /* Output the current reference projection if it changes with time;
2496 * this can happen when flooding is used as harmonic restraint */
2497 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
2499 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2500 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2503 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2504 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2506 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2507 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2510 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2511 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2513 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2515 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2516 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2519 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2520 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2523 edi = edi->next_edi;
2524 } /* End of flooding-related legend entries */
2528 /* Now the ED-related entries, if essential dynamics is done */
2530 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2532 if (bNeedDoEdsam(edi)) /* Only print ED legend if at least one ED option is on */
2534 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2536 /* Essential dynamics, projections on eigenvectors */
2537 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2538 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2539 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2540 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2541 if (edi->vecs.radfix.neig)
2543 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2545 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2546 if (edi->vecs.radacc.neig)
2548 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2550 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2551 if (edi->vecs.radcon.neig)
2553 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2556 edi = edi->next_edi;
2557 } /* end of 'pure' essential dynamics legend entries */
2558 n_edsam = nsets - n_flood;
2560 xvgr_legend(ed->edo, nsets, setname, oenv);
2563 fprintf(ed->edo, "#\n"
2564 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2565 n_flood, 1 == n_flood ? "" : "s",
2566 n_edsam, 1 == n_edsam ? "" : "s");
2567 fprintf(ed->edo, "%s", LegendStr);
2574 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2575 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2576 gmx_edsam_t init_edsam(
2577 const char *ediFileName,
2578 const char *edoFileName,
2579 const gmx_mtop_t *mtop,
2580 const t_inputrec *ir,
2581 const t_commrec *cr,
2582 gmx::Constraints *constr,
2583 const t_state *globalState,
2584 ObservablesHistory *oh,
2585 const gmx_output_env_t *oenv,
2588 t_edpar *edi = nullptr; /* points to a single edi data set */
2590 int nED = 0; /* Number of ED data sets */
2591 rvec *x_pbc = nullptr; /* positions of the whole MD system with pbc removed */
2592 rvec *xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs */
2593 rvec fit_transvec; /* translation ... */
2594 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2595 rvec *ref_x_old = nullptr; /* helper pointer */
2600 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2602 if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName) )
2604 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or flooding\n"
2605 "simulation. Please also set the .edi file on the command line with -ei.\n");
2609 /* Open input and output files, allocate space for ED data structure */
2610 gmx_edsam_t ed = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr);
2611 GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2612 constr->saveEdsamPointer(ed);
2614 /* Needed for initializing radacc radius in do_edsam */
2617 /* The input file is read by the master and the edi structures are
2618 * initialized here. Input is stored in ed->edpar. Then the edi
2619 * structures are transferred to the other nodes */
2622 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2623 * flooding vector, Essential dynamics can be applied to more than one structure
2624 * as well, but will be done in the order given in the edi file, so
2625 * expect different results for different order of edi file concatenation! */
2627 while (edi != nullptr)
2629 init_edi(mtop, edi);
2630 init_flood(edi, ed, ir->delta_t);
2631 edi = edi->next_edi;
2635 /* The master does the work here. The other nodes get the positions
2636 * not before dd_partition_system which is called after init_edsam */
2639 edsamhistory_t *EDstate = oh->edsamHistory.get();
2641 if (!EDstate->bFromCpt)
2643 /* Remove PBC, make molecule(s) subject to ED whole. */
2644 snew(x_pbc, mtop->natoms);
2645 copy_rvecn(as_rvec_array(globalState->x.data()), x_pbc, 0, mtop->natoms);
2646 do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
2648 /* Reset pointer to first ED data set which contains the actual ED data */
2650 GMX_ASSERT(nullptr != edi,
2651 "The pointer to the essential dynamics parameters is undefined");
2654 /* Loop over all ED/flooding data sets (usually only one, though) */
2655 for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2657 /* For multiple ED groups we use the output frequency that was specified
2658 * in the first set */
2661 edi->outfrq = ed->edpar->outfrq;
2664 /* Extract the initial reference and average positions. When starting
2665 * from .cpt, these have already been read into sref.x_old
2666 * in init_edsamstate() */
2667 if (!EDstate->bFromCpt)
2669 /* If this is the first run (i.e. no checkpoint present) we assume
2670 * that the starting positions give us the correct PBC representation */
2671 for (i = 0; i < edi->sref.nr; i++)
2673 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2676 for (i = 0; i < edi->sav.nr; i++)
2678 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2682 /* Now we have the PBC-correct start positions of the reference and
2683 average structure. We copy that over to dummy arrays on which we
2684 can apply fitting to print out the RMSD. We srenew the memory since
2685 the size of the buffers is likely different for every ED group */
2686 srenew(xfit, edi->sref.nr );
2687 srenew(xstart, edi->sav.nr );
2690 /* Reference indices are the same as average indices */
2691 ref_x_old = edi->sav.x_old;
2695 ref_x_old = edi->sref.x_old;
2697 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2698 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2700 /* Make the fit to the REFERENCE structure, get translation and rotation */
2701 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2703 /* Output how well we fit to the reference at the start */
2704 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2705 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2706 rmsd_from_structure(xfit, &edi->sref));
2707 if (EDstate->nED > 1)
2709 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2711 fprintf(stderr, "\n");
2713 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2714 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2716 /* calculate initial projections */
2717 project(xstart, edi);
2719 /* For the target and origin structure both a reference (fit) and an
2720 * average structure can be provided in make_edi. If both structures
2721 * are the same, make_edi only stores one of them in the .edi file.
2722 * If they differ, first the fit and then the average structure is stored
2723 * in star (or sor), thus the number of entries in star/sor is
2724 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2725 * the size of the average group. */
2727 /* process target structure, if required */
2728 if (edi->star.nr > 0)
2730 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2732 /* get translation & rotation for fit of target structure to reference structure */
2733 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2735 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2736 if (edi->star.nr == edi->sav.nr)
2740 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2742 /* The last sav.nr indices of the target structure correspond to
2743 * the average structure, which must be projected */
2744 avindex = edi->star.nr - edi->sav.nr;
2746 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
2750 rad_project(edi, xstart, &edi->vecs.radcon);
2753 /* process structure that will serve as origin of expansion circle */
2754 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2756 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2759 if (edi->sori.nr > 0)
2761 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2763 /* fit this structure to reference structure */
2764 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2766 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2767 if (edi->sori.nr == edi->sav.nr)
2771 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2773 /* For the projection, we need the last sav.nr indices of sori */
2774 avindex = edi->sori.nr - edi->sav.nr;
2777 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2778 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2779 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2781 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2782 /* Set center of flooding potential to the ORIGIN structure */
2783 rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
2784 /* We already know that no (moving) reference position was provided,
2785 * therefore we can overwrite refproj[0]*/
2786 copyEvecReference(&edi->flood.vecs);
2789 else /* No origin structure given */
2791 rad_project(edi, xstart, &edi->vecs.radacc);
2792 rad_project(edi, xstart, &edi->vecs.radfix);
2793 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2795 if (edi->flood.bHarmonic)
2797 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2798 for (i = 0; i < edi->flood.vecs.neig; i++)
2800 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2805 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2806 /* Set center of flooding potential to the center of the covariance matrix,
2807 * i.e. the average structure, i.e. zero in the projected system */
2808 for (i = 0; i < edi->flood.vecs.neig; i++)
2810 edi->flood.vecs.refproj[i] = 0.0;
2815 /* For convenience, output the center of the flooding potential for the eigenvectors */
2816 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2818 for (i = 0; i < edi->flood.vecs.neig; i++)
2820 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2821 if (edi->flood.bHarmonic)
2823 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2825 fprintf(stdout, "\n");
2829 /* set starting projections for linsam */
2830 rad_project(edi, xstart, &edi->vecs.linacc);
2831 rad_project(edi, xstart, &edi->vecs.linfix);
2833 /* Prepare for the next edi data set: */
2834 edi = edi->next_edi;
2836 /* Cleaning up on the master node: */
2837 if (!EDstate->bFromCpt)
2844 } /* end of MASTER only section */
2848 /* First let everybody know how many ED data sets to expect */
2849 gmx_bcast(sizeof(nED), &nED, cr);
2850 /* Broadcast the essential dynamics / flooding data to all nodes */
2851 broadcast_ed_data(cr, ed, nED);
2855 /* In the single-CPU case, point the local atom numbers pointers to the global
2856 * one, so that we can use the same notation in serial and parallel case: */
2857 /* Loop over all ED data sets (usually only one, though) */
2859 for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
2861 edi->sref.anrs_loc = edi->sref.anrs;
2862 edi->sav.anrs_loc = edi->sav.anrs;
2863 edi->star.anrs_loc = edi->star.anrs;
2864 edi->sori.anrs_loc = edi->sori.anrs;
2865 /* For the same reason as above, make a dummy c_ind array: */
2866 snew(edi->sav.c_ind, edi->sav.nr);
2867 /* Initialize the array */
2868 for (i = 0; i < edi->sav.nr; i++)
2870 edi->sav.c_ind[i] = i;
2872 /* In the general case we will need a different-sized array for the reference indices: */
2875 snew(edi->sref.c_ind, edi->sref.nr);
2876 for (i = 0; i < edi->sref.nr; i++)
2878 edi->sref.c_ind[i] = i;
2881 /* Point to the very same array in case of other structures: */
2882 edi->star.c_ind = edi->sav.c_ind;
2883 edi->sori.c_ind = edi->sav.c_ind;
2884 /* In the serial case, the local number of atoms is the global one: */
2885 edi->sref.nr_loc = edi->sref.nr;
2886 edi->sav.nr_loc = edi->sav.nr;
2887 edi->star.nr_loc = edi->star.nr;
2888 edi->sori.nr_loc = edi->sori.nr;
2890 /* An on we go to the next ED group */
2891 edi = edi->next_edi;
2895 /* Allocate space for ED buffer variables */
2896 /* Again, loop over ED data sets */
2898 for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
2900 /* Allocate space for ED buffer variables */
2901 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
2902 snew(edi->buf->do_edsam, 1);
2904 /* Space for collective ED buffer variables */
2906 /* Collective positions of atoms with the average indices */
2907 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2908 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
2909 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2910 /* Collective positions of atoms with the reference indices */
2913 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2914 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
2915 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2918 /* Get memory for flooding forces */
2919 snew(edi->flood.forces_cartesian, edi->sav.nr);
2923 /* Dump it all into one file per process */
2924 dump_edi(edi, cr, nr_edi);
2928 edi = edi->next_edi;
2931 /* Flush the edo file so that the user can check some things
2932 * when the simulation has started */
2942 void do_edsam(const t_inputrec *ir,
2944 const t_commrec *cr,
2950 int i, edinr, iupdate = 500;
2951 matrix rotmat; /* rotation matrix */
2952 rvec transvec; /* translation vector */
2953 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
2954 real dt_1; /* 1/dt */
2955 struct t_do_edsam *buf;
2957 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
2958 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
2961 /* Check if ED sampling has to be performed */
2962 if (ed->eEDtype == eEDnone)
2967 dt_1 = 1.0/ir->delta_t;
2969 /* Loop over all ED groups (usually one) */
2972 while (edi != nullptr)
2975 if (bNeedDoEdsam(edi))
2978 buf = edi->buf->do_edsam;
2982 /* initialize radacc radius for slope criterion */
2983 buf->oldrad = calc_radius(&edi->vecs.radacc);
2986 /* Copy the positions into buf->xc* arrays and after ED
2987 * feed back corrections to the official positions */
2989 /* Broadcast the ED positions such that every node has all of them
2990 * Every node contributes its local positions xs and stores it in
2991 * the collective buf->xcoll array. Note that for edinr > 1
2992 * xs could already have been modified by an earlier ED */
2994 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
2995 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
2997 /* Only assembly reference positions if their indices differ from the average ones */
3000 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3001 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3004 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3005 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3006 * set bUpdateShifts=TRUE in the parallel case. */
3007 buf->bUpdateShifts = FALSE;
3009 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3010 * as well as the indices in edi->sav.anrs */
3012 /* Fit the reference indices to the reference structure */
3015 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3019 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3022 /* Now apply the translation and rotation to the ED structure */
3023 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3025 /* Find out how well we fit to the reference (just for output steps) */
3026 if (do_per_step(step, edi->outfrq) && MASTER(cr))
3030 /* Indices of reference and average structures are identical,
3031 * thus we can calculate the rmsd to SREF using xcoll */
3032 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3036 /* We have to translate & rotate the reference atoms first */
3037 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3038 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3042 /* update radsam references, when required */
3043 if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3045 project(buf->xcoll, edi);
3046 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3047 rad_project(edi, buf->xcoll, &edi->vecs.radfix);
3048 buf->oldrad = -1.e5;
3051 /* update radacc references, when required */
3052 if (do_per_step(step, iupdate) && step >= edi->presteps)
3054 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
3055 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3057 project(buf->xcoll, edi);
3058 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3063 buf->oldrad = edi->vecs.radacc.radius;
3067 /* apply the constraints */
3068 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3070 /* ED constraints should be applied already in the first MD step
3071 * (which is step 0), therefore we pass step+1 to the routine */
3072 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3075 /* write to edo, when required */
3076 if (do_per_step(step, edi->outfrq))
3078 project(buf->xcoll, edi);
3079 if (MASTER(cr) && !bSuppress)
3081 write_edo(edi, ed->edo, rmsdev);
3085 /* Copy back the positions unless monitoring only */
3086 if (ed_constraints(ed->eEDtype, edi))
3088 /* remove fitting */
3089 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3091 /* Copy the ED corrected positions into the coordinate array */
3092 /* Each node copies its local part. In the serial case, nat_loc is the
3093 * total number of ED atoms */
3094 for (i = 0; i < edi->sav.nr_loc; i++)
3096 /* Unshift local ED coordinate and store in x_unsh */
3097 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3098 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3100 /* dx is the ED correction to the positions: */
3101 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3105 /* dv is the ED correction to the velocity: */
3106 svmul(dt_1, dx, dv);
3107 /* apply the velocity correction: */
3108 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3110 /* Finally apply the position correction due to ED: */
3111 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3114 } /* END of if ( bNeedDoEdsam(edi) ) */
3116 /* Prepare for the next ED group */
3117 edi = edi->next_edi;
3119 } /* END of loop over ED groups */
3124 void done_ed(gmx_edsam_t *ed)
3128 /* ed->edo is opened sometimes with xvgropen, sometimes with
3129 * gmx_fio_fopen, so we use the least common denominator for
3131 gmx_fio_fclose((*ed)->edo);
3134 /* TODO deallocate ed and set pointer to NULL */