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, 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.
45 #include "gromacs/legacyheaders/typedefs.h"
46 #include "gromacs/utility/cstringutil.h"
47 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/legacyheaders/names.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/legacyheaders/txtdump.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/legacyheaders/nrnb.h"
53 #include "gromacs/legacyheaders/mdrun.h"
54 #include "gromacs/legacyheaders/update.h"
55 #include "gromacs/topology/mtop_util.h"
56 #include "gromacs/essentialdynamics/edsam.h"
57 #include "gromacs/fileio/gmxfio.h"
58 #include "gromacs/fileio/xvgr.h"
59 #include "gromacs/mdlib/groupcoord.h"
61 #include "gromacs/linearalgebra/nrjac.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/utility/fatalerror.h"
65 /* We use the same defines as in mvdata.c here */
66 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
67 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
68 #define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
70 /* These macros determine the column width in the output file */
71 #define EDcol_sfmt "%17s"
72 #define EDcol_efmt "%17.5e"
73 #define EDcol_ffmt "%17f"
75 /* enum to identify the type of ED: none, normal ED, flooding */
77 eEDnone, eEDedsam, eEDflood, eEDnr
80 /* enum to identify operations on reference, average, origin, target structures */
82 eedREF, eedAV, eedORI, eedTAR, eedNR
88 int neig; /* nr of eigenvectors */
89 int *ieig; /* index nrs of eigenvectors */
90 real *stpsz; /* stepsizes (per eigenvector) */
91 rvec **vec; /* eigenvector components */
92 real *xproj; /* instantaneous x projections */
93 real *fproj; /* instantaneous f projections */
94 real radius; /* instantaneous radius */
95 real *refproj; /* starting or target projections */
96 /* When using flooding as harmonic restraint: The current reference projection
97 * is at each step calculated from the initial refproj0 and the slope. */
98 real *refproj0, *refprojslope;
104 t_eigvec mon; /* only monitored, no constraints */
105 t_eigvec linfix; /* fixed linear constraints */
106 t_eigvec linacc; /* acceptance linear constraints */
107 t_eigvec radfix; /* fixed radial constraints (exp) */
108 t_eigvec radacc; /* acceptance radial constraints (exp) */
109 t_eigvec radcon; /* acceptance rad. contraction constr. */
116 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
118 gmx_bool bConstForce; /* Do not calculate a flooding potential,
119 instead flood with a constant force */
128 rvec *forces_cartesian;
129 t_eigvec vecs; /* use flooding for these */
133 /* This type is for the average, reference, target, and origin structure */
134 typedef struct gmx_edx
136 int nr; /* number of atoms this structure contains */
137 int nr_loc; /* number of atoms on local node */
138 int *anrs; /* atom index numbers */
139 int *anrs_loc; /* local atom index numbers */
140 int nalloc_loc; /* allocation size of anrs_loc */
141 int *c_ind; /* at which position of the whole anrs
142 * array is a local atom?, i.e.
143 * c_ind[0...nr_loc-1] gives the atom index
144 * with respect to the collective
145 * anrs[0...nr-1] array */
146 rvec *x; /* positions for this structure */
147 rvec *x_old; /* Last positions which have the correct PBC
148 representation of the ED group. In
149 combination with keeping track of the
150 shift vectors, the ED group can always
152 real *m; /* masses */
153 real mtot; /* total mass (only used in sref) */
154 real *sqrtm; /* sqrt of the masses used for mass-
155 * weighting of analysis (only used in sav) */
161 int nini; /* total Nr of atoms */
162 gmx_bool fitmas; /* true if trans fit with cm */
163 gmx_bool pcamas; /* true if mass-weighted PCA */
164 int presteps; /* number of steps to run without any
165 * perturbations ... just monitoring */
166 int outfrq; /* freq (in steps) of writing to edo */
167 int maxedsteps; /* max nr of steps per cycle */
169 /* all gmx_edx datasets are copied to all nodes in the parallel case */
170 struct gmx_edx sref; /* reference positions, to these fitting
172 gmx_bool bRefEqAv; /* If true, reference & average indices
173 * are the same. Used for optimization */
174 struct gmx_edx sav; /* average positions */
175 struct gmx_edx star; /* target positions */
176 struct gmx_edx sori; /* origin positions */
178 t_edvecs vecs; /* eigenvectors */
179 real slope; /* minimal slope in acceptance radexp */
181 t_edflood flood; /* parameters especially for flooding */
182 struct t_ed_buffer *buf; /* handle to local buffers */
183 struct edpar *next_edi; /* Pointer to another ED group */
187 typedef struct gmx_edsam
189 int eEDtype; /* Type of ED: see enums above */
190 FILE *edo; /* output file pointer */
200 rvec old_transvec, older_transvec, transvec_compact;
201 rvec *xcoll; /* Positions from all nodes, this is the
202 collective set we work on.
203 These are the positions of atoms with
204 average structure indices */
205 rvec *xc_ref; /* same but with reference structure indices */
206 ivec *shifts_xcoll; /* Shifts for xcoll */
207 ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
208 ivec *shifts_xc_ref; /* Shifts for xc_ref */
209 ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
210 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
211 ED shifts for this ED group need to
216 /* definition of ED buffer structure */
219 struct t_fit_to_ref * fit_to_ref;
220 struct t_do_edfit * do_edfit;
221 struct t_do_edsam * do_edsam;
222 struct t_do_radcon * do_radcon;
226 /* Function declarations */
227 static void fit_to_reference(rvec *xcoll, rvec transvec, matrix rotmat, t_edpar *edi);
228 static void translate_and_rotate(rvec *x, int nat, rvec transvec, matrix rotmat);
229 static real rmsd_from_structure(rvec *x, struct gmx_edx *s);
230 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms);
231 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate);
232 static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate);
233 static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv);
234 /* End function declarations */
237 /* Do we have to perform essential dynamics constraints or possibly only flooding
238 * for any of the ED groups? */
239 static gmx_bool bNeedDoEdsam(t_edpar *edi)
241 return edi->vecs.mon.neig
242 || edi->vecs.linfix.neig
243 || edi->vecs.linacc.neig
244 || edi->vecs.radfix.neig
245 || edi->vecs.radacc.neig
246 || edi->vecs.radcon.neig;
250 /* Multiple ED groups will be labeled with letters instead of numbers
251 * to avoid confusion with eigenvector indices */
252 static char get_EDgroupChar(int nr_edi, int nED)
260 * nr_edi = 2 -> B ...
262 return 'A' + nr_edi - 1;
266 /* Does not subtract average positions, projection on single eigenvector is returned
267 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
268 * Average position is subtracted in ed_apply_constraints prior to calling projectx
270 static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec)
276 for (i = 0; i < edi->sav.nr; i++)
278 proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
285 /* Specialized: projection is stored in vec->refproj
286 * -> used for radacc, radfix, radcon and center of flooding potential
287 * subtracts average positions, projects vector x */
288 static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec)
293 /* Subtract average positions */
294 for (i = 0; i < edi->sav.nr; i++)
296 rvec_dec(x[i], edi->sav.x[i]);
299 for (i = 0; i < vec->neig; i++)
301 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
302 rad += pow((vec->refproj[i]-vec->xproj[i]), 2);
304 vec->radius = sqrt(rad);
306 /* Add average positions */
307 for (i = 0; i < edi->sav.nr; i++)
309 rvec_inc(x[i], edi->sav.x[i]);
314 /* Project vector x, subtract average positions prior to projection and add
315 * them afterwards to retain the unchanged vector. Store in xproj. Mass-weighting
317 static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
318 t_eigvec *vec, /* The eigenvectors */
329 /* Subtract average positions */
330 for (i = 0; i < edi->sav.nr; i++)
332 rvec_dec(x[i], edi->sav.x[i]);
335 for (i = 0; i < vec->neig; i++)
337 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
340 /* Add average positions */
341 for (i = 0; i < edi->sav.nr; i++)
343 rvec_inc(x[i], edi->sav.x[i]);
348 /* Project vector x onto all edi->vecs (mon, linfix,...) */
349 static void project(rvec *x, /* positions to project */
350 t_edpar *edi) /* edi data set */
352 /* It is not more work to subtract the average position in every
353 * subroutine again, because these routines are rarely used simultaneously */
354 project_to_eigvectors(x, &edi->vecs.mon, edi);
355 project_to_eigvectors(x, &edi->vecs.linfix, edi);
356 project_to_eigvectors(x, &edi->vecs.linacc, edi);
357 project_to_eigvectors(x, &edi->vecs.radfix, edi);
358 project_to_eigvectors(x, &edi->vecs.radacc, edi);
359 project_to_eigvectors(x, &edi->vecs.radcon, edi);
363 static real calc_radius(t_eigvec *vec)
369 for (i = 0; i < vec->neig; i++)
371 rad += pow((vec->refproj[i]-vec->xproj[i]), 2);
374 return rad = sqrt(rad);
380 static void dump_xcoll(t_edpar *edi, struct t_do_edsam *buf, t_commrec *cr,
387 ivec *shifts, *eshifts;
396 shifts = buf->shifts_xcoll;
397 eshifts = buf->extra_shifts_xcoll;
399 sprintf(fn, "xcolldump_step%d.txt", step);
402 for (i = 0; i < edi->sav.nr; i++)
404 fprintf(fp, "%d %9.5f %9.5f %9.5f %d %d %d %d %d %d\n",
406 xcoll[i][XX], xcoll[i][YY], xcoll[i][ZZ],
407 shifts[i][XX], shifts[i][YY], shifts[i][ZZ],
408 eshifts[i][XX], eshifts[i][YY], eshifts[i][ZZ]);
416 static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
421 fprintf(out, "#%s positions:\n%d\n", name, s->nr);
427 fprintf(out, "#index, x, y, z");
430 fprintf(out, ", sqrt(m)");
432 for (i = 0; i < s->nr; i++)
434 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]);
437 fprintf(out, "%9.3f", s->sqrtm[i]);
445 static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
446 const char name[], int length)
451 fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig);
452 /* Dump the data for every eigenvector: */
453 for (i = 0; i < ev->neig; i++)
455 fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
456 ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
457 for (j = 0; j < length; j++)
459 fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
466 static void dump_edi(t_edpar *edpars, t_commrec *cr, int nr_edi)
472 sprintf(fn, "EDdump_rank%d_edi%d", cr->nodeid, nr_edi);
473 out = gmx_ffopen(fn, "w");
475 fprintf(out, "#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
476 edpars->nini, edpars->fitmas, edpars->pcamas);
477 fprintf(out, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
478 edpars->outfrq, edpars->maxedsteps, edpars->slope);
479 fprintf(out, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n",
480 edpars->presteps, edpars->flood.deltaF0, edpars->flood.tau,
481 edpars->flood.constEfl, edpars->flood.alpha2);
483 /* Dump reference, average, target, origin positions */
484 dump_edi_positions(out, &edpars->sref, "REFERENCE");
485 dump_edi_positions(out, &edpars->sav, "AVERAGE" );
486 dump_edi_positions(out, &edpars->star, "TARGET" );
487 dump_edi_positions(out, &edpars->sori, "ORIGIN" );
489 /* Dump eigenvectors */
490 dump_edi_eigenvecs(out, &edpars->vecs.mon, "MONITORED", edpars->sav.nr);
491 dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX", edpars->sav.nr);
492 dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC", edpars->sav.nr);
493 dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX", edpars->sav.nr);
494 dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC", edpars->sav.nr);
495 dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON", edpars->sav.nr);
497 /* Dump flooding eigenvectors */
498 dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING", edpars->sav.nr);
500 /* Dump ed local buffer */
501 fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit );
502 fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam );
503 fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon );
510 static void dump_rotmat(FILE* out, matrix rotmat)
512 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[XX][XX], rotmat[XX][YY], rotmat[XX][ZZ]);
513 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[YY][XX], rotmat[YY][YY], rotmat[YY][ZZ]);
514 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[ZZ][XX], rotmat[ZZ][YY], rotmat[ZZ][ZZ]);
519 static void dump_rvec(FILE *out, int dim, rvec *x)
524 for (i = 0; i < dim; i++)
526 fprintf(out, "%4d %f %f %f\n", i, x[i][XX], x[i][YY], x[i][ZZ]);
532 static void dump_mat(FILE* out, int dim, double** mat)
537 fprintf(out, "MATRIX:\n");
538 for (i = 0; i < dim; i++)
540 for (j = 0; j < dim; j++)
542 fprintf(out, "%f ", mat[i][j]);
555 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
557 /* this is a copy of do_fit with some modifications */
558 int c, r, n, j, i, irot;
559 double d[6], xnr, xpc;
564 struct t_do_edfit *loc;
567 if (edi->buf->do_edfit != NULL)
574 snew(edi->buf->do_edfit, 1);
576 loc = edi->buf->do_edfit;
580 snew(loc->omega, 2*DIM);
581 snew(loc->om, 2*DIM);
582 for (i = 0; i < 2*DIM; i++)
584 snew(loc->omega[i], 2*DIM);
585 snew(loc->om[i], 2*DIM);
589 for (i = 0; (i < 6); i++)
592 for (j = 0; (j < 6); j++)
594 loc->omega[i][j] = 0;
599 /* calculate the matrix U */
601 for (n = 0; (n < natoms); n++)
603 for (c = 0; (c < DIM); c++)
606 for (r = 0; (r < DIM); r++)
614 /* construct loc->omega */
615 /* loc->omega is symmetric -> loc->omega==loc->omega' */
616 for (r = 0; (r < 6); r++)
618 for (c = 0; (c <= r); c++)
620 if ((r >= 3) && (c < 3))
622 loc->omega[r][c] = u[r-3][c];
623 loc->omega[c][r] = u[r-3][c];
627 loc->omega[r][c] = 0;
628 loc->omega[c][r] = 0;
633 /* determine h and k */
637 dump_mat(stderr, 2*DIM, loc->omega);
638 for (i = 0; i < 6; i++)
640 fprintf(stderr, "d[%d] = %f\n", i, d[i]);
644 jacobi(loc->omega, 6, d, loc->om, &irot);
648 fprintf(stderr, "IROT=0\n");
651 index = 0; /* For the compiler only */
653 for (j = 0; (j < 3); j++)
656 for (i = 0; (i < 6); i++)
665 for (i = 0; (i < 3); i++)
667 vh[j][i] = M_SQRT2*loc->om[i][index];
668 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
673 for (c = 0; (c < 3); c++)
675 for (r = 0; (r < 3); r++)
677 R[c][r] = vk[0][r]*vh[0][c]+
684 for (c = 0; (c < 3); c++)
686 for (r = 0; (r < 3); r++)
688 R[c][r] = vk[0][r]*vh[0][c]+
697 static void rmfit(int nat, rvec *xcoll, rvec transvec, matrix rotmat)
704 * The inverse rotation is described by the transposed rotation matrix */
705 transpose(rotmat, tmat);
706 rotate_x(xcoll, nat, tmat);
708 /* Remove translation */
709 vec[XX] = -transvec[XX];
710 vec[YY] = -transvec[YY];
711 vec[ZZ] = -transvec[ZZ];
712 translate_x(xcoll, nat, vec);
716 /**********************************************************************************
717 ******************** FLOODING ****************************************************
718 **********************************************************************************
720 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
721 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
722 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
724 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
725 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
727 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
728 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
730 do_flood makes a copy of the positions,
731 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
732 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
733 fit. Then do_flood adds these forces to the forcefield-forces
734 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
736 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
737 structure is projected to the system of eigenvectors and then this position in the subspace is used as
738 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
739 i.e. the average structure as given in the make_edi file.
741 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
742 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
743 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
744 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
745 further adaption is applied, Efl will stay constant at zero.
747 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
748 used as spring constants for the harmonic potential.
749 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
750 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
752 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
753 the routine read_edi_file reads all of theses flooding files.
754 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
755 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
756 edi there is no interdependence whatsoever. The forces are added together.
758 To write energies into the .edr file, call the function
759 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
761 get_flood_energies(real Vfl[],int nnames);
764 - 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.
766 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
767 two edsam files from two peptide chains
770 static void write_edo_flood(t_edpar *edi, FILE *fp, real rmsd)
775 /* Output how well we fit to the reference structure */
776 fprintf(fp, EDcol_ffmt, rmsd);
778 for (i = 0; i < edi->flood.vecs.neig; i++)
780 fprintf(fp, EDcol_efmt, edi->flood.vecs.xproj[i]);
782 /* Check whether the reference projection changes with time (this can happen
783 * in case flooding is used as harmonic restraint). If so, output the
784 * current reference projection */
785 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
787 fprintf(fp, EDcol_efmt, edi->flood.vecs.refproj[i]);
790 /* Output Efl if we are doing adaptive flooding */
791 if (0 != edi->flood.tau)
793 fprintf(fp, EDcol_efmt, edi->flood.Efl);
795 fprintf(fp, EDcol_efmt, edi->flood.Vfl);
797 /* Output deltaF if we are doing adaptive flooding */
798 if (0 != edi->flood.tau)
800 fprintf(fp, EDcol_efmt, edi->flood.deltaF);
802 fprintf(fp, EDcol_efmt, edi->flood.vecs.fproj[i]);
807 /* From flood.xproj compute the Vfl(x) at this point */
808 static real flood_energy(t_edpar *edi, gmx_int64_t step)
810 /* compute flooding energy Vfl
811 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
812 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
813 it is already computed by make_edi and stored in stpsz[i]
815 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
822 /* Each time this routine is called (i.e. each time step), we add a small
823 * value to the reference projection. This way a harmonic restraint towards
824 * a moving reference is realized. If no value for the additive constant
825 * is provided in the edi file, the reference will not change. */
826 if (edi->flood.bHarmonic)
828 for (i = 0; i < edi->flood.vecs.neig; i++)
830 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
835 /* Compute sum which will be the exponent of the exponential */
836 for (i = 0; i < edi->flood.vecs.neig; i++)
838 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
839 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]);
842 /* Compute the Gauss function*/
843 if (edi->flood.bHarmonic)
845 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
849 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
856 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
857 static void flood_forces(t_edpar *edi)
859 /* compute the forces in the subspace of the flooding eigenvectors
860 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
863 real energy = edi->flood.Vfl;
866 if (edi->flood.bHarmonic)
868 for (i = 0; i < edi->flood.vecs.neig; i++)
870 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
875 for (i = 0; i < edi->flood.vecs.neig; i++)
877 /* if Efl is zero the forces are zero if not use the formula */
878 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;
884 /* Raise forces from subspace into cartesian space */
885 static void flood_blowup(t_edpar *edi, rvec *forces_cart)
887 /* this function lifts the forces from the subspace to the cartesian space
888 all the values not contained in the subspace are assumed to be zero and then
889 a coordinate transformation from eigenvector to cartesian vectors is performed
890 The nonexistent values don't have to be set to zero explicitly, they would occur
891 as zero valued summands, hence we just stop to compute this part of the sum.
893 for every atom we add all the contributions to this atom from all the different eigenvectors.
895 NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
896 field forces_cart prior the computation, but we compute the forces separately
897 to have them accessible for diagnostics
904 forces_sub = edi->flood.vecs.fproj;
907 /* Calculate the cartesian forces for the local atoms */
909 /* Clear forces first */
910 for (j = 0; j < edi->sav.nr_loc; j++)
912 clear_rvec(forces_cart[j]);
915 /* Now compute atomwise */
916 for (j = 0; j < edi->sav.nr_loc; j++)
918 /* Compute forces_cart[edi->sav.anrs[j]] */
919 for (eig = 0; eig < edi->flood.vecs.neig; eig++)
921 /* Force vector is force * eigenvector (compute only atom j) */
922 svmul(forces_sub[eig], edi->flood.vecs.vec[eig][edi->sav.c_ind[j]], dum);
923 /* Add this vector to the cartesian forces */
924 rvec_inc(forces_cart[j], dum);
930 /* Update the values of Efl, deltaF depending on tau and Vfl */
931 static void update_adaption(t_edpar *edi)
933 /* this function updates the parameter Efl and deltaF according to the rules given in
934 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
937 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
939 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
940 /* check if restrain (inverted flooding) -> don't let EFL become positive */
941 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
946 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
951 static void do_single_flood(
959 gmx_bool bNS) /* Are we in a neighbor searching step? */
962 matrix rotmat; /* rotation matrix */
963 matrix tmat; /* inverse rotation */
964 rvec transvec; /* translation vector */
966 struct t_do_edsam *buf;
969 buf = edi->buf->do_edsam;
972 /* Broadcast the positions of the AVERAGE structure such that they are known on
973 * every processor. Each node contributes its local positions x and stores them in
974 * the collective ED array buf->xcoll */
975 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
976 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
978 /* Only assembly REFERENCE positions if their indices differ from the average ones */
981 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
982 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
985 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
986 * We do not need to update the shifts until the next NS step */
987 buf->bUpdateShifts = FALSE;
989 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
990 * as well as the indices in edi->sav.anrs */
992 /* Fit the reference indices to the reference structure */
995 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
999 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
1002 /* Now apply the translation and rotation to the ED structure */
1003 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
1005 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
1006 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, edi);
1008 if (FALSE == edi->flood.bConstForce)
1010 /* Compute Vfl(x) from flood.xproj */
1011 edi->flood.Vfl = flood_energy(edi, step);
1013 update_adaption(edi);
1015 /* Compute the flooding forces */
1019 /* Translate them into cartesian positions */
1020 flood_blowup(edi, edi->flood.forces_cartesian);
1022 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
1023 /* Each node rotates back its local forces */
1024 transpose(rotmat, tmat);
1025 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
1027 /* Finally add forces to the main force variable */
1028 for (i = 0; i < edi->sav.nr_loc; i++)
1030 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
1033 /* Output is written by the master process */
1034 if (do_per_step(step, edi->outfrq) && MASTER(cr))
1036 /* Output how well we fit to the reference */
1039 /* Indices of reference and average structures are identical,
1040 * thus we can calculate the rmsd to SREF using xcoll */
1041 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
1045 /* We have to translate & rotate the reference atoms first */
1046 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1047 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1050 write_edo_flood(edi, edo, rmsdev);
1055 /* Main flooding routine, called from do_force */
1056 extern void do_flood(
1071 /* Write time to edo, when required. Output the time anyhow since we need
1072 * it in the output file for ED constraints. */
1073 if (MASTER(cr) && do_per_step(step, edi->outfrq))
1075 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
1078 if (ed->eEDtype != eEDflood)
1085 /* Call flooding for one matrix */
1086 if (edi->flood.vecs.neig)
1088 do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS);
1090 edi = edi->next_edi;
1095 /* Called by init_edi, configure some flooding related variables and structures,
1096 * print headers to output files */
1097 static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt)
1102 edi->flood.Efl = edi->flood.constEfl;
1106 if (edi->flood.vecs.neig)
1108 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1109 ed->eEDtype = eEDflood;
1111 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
1113 if (edi->flood.bConstForce)
1115 /* We have used stpsz as a vehicle to carry the fproj values for constant
1116 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1117 * in const force flooding, fproj is never changed. */
1118 for (i = 0; i < edi->flood.vecs.neig; i++)
1120 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1122 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1123 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1131 /*********** Energy book keeping ******/
1132 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1141 srenew(names, count);
1142 sprintf(buf, "Vfl_%d", count);
1143 names[count-1] = gmx_strdup(buf);
1144 actual = actual->next_edi;
1151 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1153 /*fl has to be big enough to capture nnames-many entries*/
1162 Vfl[count-1] = actual->flood.Vfl;
1163 actual = actual->next_edi;
1166 if (nnames != count-1)
1168 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1171 /************* END of FLOODING IMPLEMENTATION ****************************/
1175 gmx_edsam_t ed_open(int natoms, edsamstate_t *EDstate, int nfile, const t_filenm fnm[], unsigned long Flags, const output_env_t oenv, t_commrec *cr)
1181 /* Allocate space for the ED data structure */
1184 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1185 ed->eEDtype = eEDedsam;
1189 fprintf(stderr, "ED sampling will be performed!\n");
1192 /* Read the edi input file: */
1193 nED = read_edi_file(ftp2fn(efEDI, nfile, fnm), ed->edpar, natoms);
1195 /* Make sure the checkpoint was produced in a run using this .edi file */
1196 if (EDstate->bFromCpt)
1198 crosscheck_edi_file_vs_checkpoint(ed, EDstate);
1204 init_edsamstate(ed, EDstate);
1206 /* The master opens the ED output file */
1207 if (Flags & MD_APPENDFILES)
1209 ed->edo = gmx_fio_fopen(opt2fn("-eo", nfile, fnm), "a+");
1213 ed->edo = xvgropen(opt2fn("-eo", nfile, fnm),
1214 "Essential dynamics / flooding output",
1216 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1218 /* Make a descriptive legend */
1219 write_edo_legend(ed, EDstate->nED, oenv);
1226 /* Broadcasts the structure data */
1227 static void bc_ed_positions(t_commrec *cr, struct gmx_edx *s, int stype)
1229 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1230 snew_bc(cr, s->x, s->nr ); /* Positions */
1231 nblock_bc(cr, s->nr, s->anrs );
1232 nblock_bc(cr, s->nr, s->x );
1234 /* For the average & reference structures we need an array for the collective indices,
1235 * and we need to broadcast the masses as well */
1236 if (stype == eedAV || stype == eedREF)
1238 /* We need these additional variables in the parallel case: */
1239 snew(s->c_ind, s->nr ); /* Collective indices */
1240 /* Local atom indices get assigned in dd_make_local_group_indices.
1241 * There, also memory is allocated */
1242 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1243 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1244 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1247 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1248 if (stype == eedREF)
1250 snew_bc(cr, s->m, s->nr);
1251 nblock_bc(cr, s->nr, s->m);
1254 /* For the average structure we might need the masses for mass-weighting */
1257 snew_bc(cr, s->sqrtm, s->nr);
1258 nblock_bc(cr, s->nr, s->sqrtm);
1259 snew_bc(cr, s->m, s->nr);
1260 nblock_bc(cr, s->nr, s->m);
1265 /* Broadcasts the eigenvector data */
1266 static void bc_ed_vecs(t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
1270 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1271 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1272 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1273 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1274 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1276 nblock_bc(cr, ev->neig, ev->ieig );
1277 nblock_bc(cr, ev->neig, ev->stpsz );
1278 nblock_bc(cr, ev->neig, ev->xproj );
1279 nblock_bc(cr, ev->neig, ev->fproj );
1280 nblock_bc(cr, ev->neig, ev->refproj);
1282 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1283 for (i = 0; i < ev->neig; i++)
1285 snew_bc(cr, ev->vec[i], length);
1286 nblock_bc(cr, length, ev->vec[i]);
1289 /* For harmonic restraints the reference projections can change with time */
1292 snew_bc(cr, ev->refproj0, ev->neig);
1293 snew_bc(cr, ev->refprojslope, ev->neig);
1294 nblock_bc(cr, ev->neig, ev->refproj0 );
1295 nblock_bc(cr, ev->neig, ev->refprojslope);
1300 /* Broadcasts the ED / flooding data to other nodes
1301 * and allocates memory where needed */
1302 static void broadcast_ed_data(t_commrec *cr, gmx_edsam_t ed, int numedis)
1308 /* Master lets the other nodes know if its ED only or also flooding */
1309 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1311 snew_bc(cr, ed->edpar, 1);
1312 /* Now transfer the ED data set(s) */
1314 for (nr = 0; nr < numedis; nr++)
1316 /* Broadcast a single ED data set */
1319 /* Broadcast positions */
1320 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1321 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1322 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1323 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1325 /* Broadcast eigenvectors */
1326 bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr, FALSE);
1327 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
1328 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
1329 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
1330 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
1331 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
1332 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1333 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
1335 /* Set the pointer to the next ED group */
1338 snew_bc(cr, edi->next_edi, 1);
1339 edi = edi->next_edi;
1345 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1346 static void init_edi(gmx_mtop_t *mtop, t_edpar *edi)
1349 real totalmass = 0.0;
1351 gmx_mtop_atomlookup_t alook = NULL;
1354 /* NOTE Init_edi is executed on the master process only
1355 * The initialized data sets are then transmitted to the
1356 * other nodes in broadcast_ed_data */
1358 alook = gmx_mtop_atomlookup_init(mtop);
1360 /* evaluate masses (reference structure) */
1361 snew(edi->sref.m, edi->sref.nr);
1362 for (i = 0; i < edi->sref.nr; i++)
1366 gmx_mtop_atomnr_to_atom(alook, edi->sref.anrs[i], &atom);
1367 edi->sref.m[i] = atom->m;
1371 edi->sref.m[i] = 1.0;
1374 /* Check that every m > 0. Bad things will happen otherwise. */
1375 if (edi->sref.m[i] <= 0.0)
1377 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1378 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1379 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1380 "atoms from the reference structure by creating a proper index group.\n",
1381 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1384 totalmass += edi->sref.m[i];
1386 edi->sref.mtot = totalmass;
1388 /* Masses m and sqrt(m) for the average structure. Note that m
1389 * is needed if forces have to be evaluated in do_edsam */
1390 snew(edi->sav.sqrtm, edi->sav.nr );
1391 snew(edi->sav.m, edi->sav.nr );
1392 for (i = 0; i < edi->sav.nr; i++)
1394 gmx_mtop_atomnr_to_atom(alook, edi->sav.anrs[i], &atom);
1395 edi->sav.m[i] = atom->m;
1398 edi->sav.sqrtm[i] = sqrt(atom->m);
1402 edi->sav.sqrtm[i] = 1.0;
1405 /* Check that every m > 0. Bad things will happen otherwise. */
1406 if (edi->sav.sqrtm[i] <= 0.0)
1408 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1409 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1410 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1411 "atoms from the average structure by creating a proper index group.\n",
1412 i, edi->sav.anrs[i]+1, atom->m);
1416 gmx_mtop_atomlookup_destroy(alook);
1418 /* put reference structure in origin */
1419 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1423 translate_x(edi->sref.x, edi->sref.nr, com);
1425 /* Init ED buffer */
1430 static void check(const char *line, const char *label)
1432 if (!strstr(line, label))
1434 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1439 static int read_checked_edint(FILE *file, const char *label)
1441 char line[STRLEN+1];
1445 fgets2 (line, STRLEN, file);
1447 fgets2 (line, STRLEN, file);
1448 sscanf (line, "%d", &idum);
1453 static int read_edint(FILE *file, gmx_bool *bEOF)
1455 char line[STRLEN+1];
1460 eof = fgets2 (line, STRLEN, file);
1466 eof = fgets2 (line, STRLEN, file);
1472 sscanf (line, "%d", &idum);
1478 static real read_checked_edreal(FILE *file, const char *label)
1480 char line[STRLEN+1];
1484 fgets2 (line, STRLEN, file);
1486 fgets2 (line, STRLEN, file);
1487 sscanf (line, "%lf", &rdum);
1488 return (real) rdum; /* always read as double and convert to single */
1492 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1495 char line[STRLEN+1];
1499 for (i = 0; i < number; i++)
1501 fgets2 (line, STRLEN, file);
1502 sscanf (line, "%d%lf%lf%lf", &anrs[i], &d[0], &d[1], &d[2]);
1503 anrs[i]--; /* we are reading FORTRAN indices */
1504 for (j = 0; j < 3; j++)
1506 x[i][j] = d[j]; /* always read as double and convert to single */
1512 static void scan_edvec(FILE *in, int nr, rvec *vec)
1514 char line[STRLEN+1];
1519 for (i = 0; (i < nr); i++)
1521 fgets2 (line, STRLEN, in);
1522 sscanf (line, "%le%le%le", &x, &y, &z);
1530 static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1533 double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0;
1534 char line[STRLEN+1];
1537 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1540 snew(tvec->ieig, tvec->neig);
1541 snew(tvec->stpsz, tvec->neig);
1542 snew(tvec->vec, tvec->neig);
1543 snew(tvec->xproj, tvec->neig);
1544 snew(tvec->fproj, tvec->neig);
1545 snew(tvec->refproj, tvec->neig);
1548 snew(tvec->refproj0, tvec->neig);
1549 snew(tvec->refprojslope, tvec->neig);
1552 for (i = 0; (i < tvec->neig); i++)
1554 fgets2 (line, STRLEN, in);
1555 if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1557 nscan = sscanf(line, "%d%lf%lf%lf", &idum, &rdum, &refproj_dum, &refprojslope_dum);
1558 /* Zero out values which were not scanned */
1562 /* Every 4 values read, including reference position */
1563 *bHaveReference = TRUE;
1566 /* A reference position is provided */
1567 *bHaveReference = TRUE;
1568 /* No value for slope, set to 0 */
1569 refprojslope_dum = 0.0;
1572 /* No values for reference projection and slope, set to 0 */
1574 refprojslope_dum = 0.0;
1577 gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1580 tvec->refproj[i] = refproj_dum;
1581 tvec->refproj0[i] = refproj_dum;
1582 tvec->refprojslope[i] = refprojslope_dum;
1584 else /* Normal flooding */
1586 nscan = sscanf(line, "%d%lf", &idum, &rdum);
1589 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1592 tvec->ieig[i] = idum;
1593 tvec->stpsz[i] = rdum;
1594 } /* end of loop over eigenvectors */
1596 for (i = 0; (i < tvec->neig); i++)
1598 snew(tvec->vec[i], nr);
1599 scan_edvec(in, nr, tvec->vec[i]);
1605 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1606 static void read_edvecs(FILE *in, int nr, t_edvecs *vecs)
1608 gmx_bool bHaveReference = FALSE;
1611 read_edvec(in, nr, &vecs->mon, FALSE, &bHaveReference);
1612 read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1613 read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1614 read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1615 read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1616 read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1620 /* Check if the same atom indices are used for reference and average positions */
1621 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1626 /* If the number of atoms differs between the two structures,
1627 * they cannot be identical */
1628 if (sref.nr != sav.nr)
1633 /* Now that we know that both stuctures have the same number of atoms,
1634 * check if also the indices are identical */
1635 for (i = 0; i < sav.nr; i++)
1637 if (sref.anrs[i] != sav.anrs[i])
1642 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1648 static int read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, const char *fn)
1651 const int magic = 670;
1654 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1655 gmx_bool bHaveReference = FALSE;
1658 /* the edi file is not free format, so expect problems if the input is corrupt. */
1660 /* check the magic number */
1661 readmagic = read_edint(in, &bEOF);
1662 /* Check whether we have reached the end of the input file */
1668 if (readmagic != magic)
1670 if (readmagic == 666 || readmagic == 667 || readmagic == 668)
1672 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1674 else if (readmagic != 669)
1676 gmx_fatal(FARGS, "Wrong magic number %d in %s", readmagic, fn);
1680 /* check the number of atoms */
1681 edi->nini = read_edint(in, &bEOF);
1682 if (edi->nini != nr_mdatoms)
1684 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1687 /* Done checking. For the rest we blindly trust the input */
1688 edi->fitmas = read_checked_edint(in, "FITMAS");
1689 edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS");
1690 edi->outfrq = read_checked_edint(in, "OUTFRQ");
1691 edi->maxedsteps = read_checked_edint(in, "MAXLEN");
1692 edi->slope = read_checked_edreal(in, "SLOPECRIT");
1694 edi->presteps = read_checked_edint(in, "PRESTEPS");
1695 edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1696 edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1697 edi->flood.tau = read_checked_edreal(in, "TAU");
1698 edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1699 edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1700 edi->flood.kT = read_checked_edreal(in, "KT");
1701 edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1702 if (readmagic > 669)
1704 edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1708 edi->flood.bConstForce = FALSE;
1710 edi->sref.nr = read_checked_edint(in, "NREF");
1712 /* allocate space for reference positions and read them */
1713 snew(edi->sref.anrs, edi->sref.nr);
1714 snew(edi->sref.x, edi->sref.nr);
1715 snew(edi->sref.x_old, edi->sref.nr);
1716 edi->sref.sqrtm = NULL;
1717 read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1719 /* average positions. they define which atoms will be used for ED sampling */
1720 edi->sav.nr = read_checked_edint(in, "NAV");
1721 snew(edi->sav.anrs, edi->sav.nr);
1722 snew(edi->sav.x, edi->sav.nr);
1723 snew(edi->sav.x_old, edi->sav.nr);
1724 read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1726 /* Check if the same atom indices are used for reference and average positions */
1727 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1730 read_edvecs(in, edi->sav.nr, &edi->vecs);
1731 read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
1733 /* target positions */
1734 edi->star.nr = read_edint(in, &bEOF);
1735 if (edi->star.nr > 0)
1737 snew(edi->star.anrs, edi->star.nr);
1738 snew(edi->star.x, edi->star.nr);
1739 edi->star.sqrtm = NULL;
1740 read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1743 /* positions defining origin of expansion circle */
1744 edi->sori.nr = read_edint(in, &bEOF);
1745 if (edi->sori.nr > 0)
1749 /* Both an -ori structure and a at least one manual reference point have been
1750 * specified. That's ambiguous and probably not intentional. */
1751 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1752 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1754 snew(edi->sori.anrs, edi->sori.nr);
1755 snew(edi->sori.x, edi->sori.nr);
1756 edi->sori.sqrtm = NULL;
1757 read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1766 /* Read in the edi input file. Note that it may contain several ED data sets which were
1767 * achieved by concatenating multiple edi files. The standard case would be a single ED
1768 * data set, though. */
1769 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1772 t_edpar *curr_edi, *last_edi;
1777 /* This routine is executed on the master only */
1779 /* Open the .edi parameter input file */
1780 in = gmx_fio_fopen(fn, "r");
1781 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1783 /* Now read a sequence of ED input parameter sets from the edi file */
1786 while (read_edi(in, curr_edi, nr_mdatoms, fn) )
1790 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1791 /* We need to allocate space for the data: */
1793 /* Point the 'next_edi' entry to the next edi: */
1794 curr_edi->next_edi = edi_read;
1795 /* Keep the curr_edi pointer for the case that the next group is empty: */
1796 last_edi = curr_edi;
1797 /* Let's prepare to read in the next edi data set: */
1798 /* cppcheck-suppress uninitvar Fixed in cppcheck 1.65 */
1799 curr_edi = edi_read;
1803 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1806 /* Terminate the edi group list with a NULL pointer: */
1807 last_edi->next_edi = NULL;
1809 fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1811 /* Close the .edi file again */
1818 struct t_fit_to_ref {
1819 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1822 /* Fit the current positions to the reference positions
1823 * Do not actually do the fit, just return rotation and translation.
1824 * Note that the COM of the reference structure was already put into
1825 * the origin by init_edi. */
1826 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1827 rvec transvec, /* The translation vector */
1828 matrix rotmat, /* The rotation matrix */
1829 t_edpar *edi) /* Just needed for do_edfit */
1831 rvec com; /* center of mass */
1833 struct t_fit_to_ref *loc;
1836 /* Allocate memory the first time this routine is called for each edi group */
1837 if (NULL == edi->buf->fit_to_ref)
1839 snew(edi->buf->fit_to_ref, 1);
1840 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1842 loc = edi->buf->fit_to_ref;
1844 /* We do not touch the original positions but work on a copy. */
1845 for (i = 0; i < edi->sref.nr; i++)
1847 copy_rvec(xcoll[i], loc->xcopy[i]);
1850 /* Calculate the center of mass */
1851 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1853 transvec[XX] = -com[XX];
1854 transvec[YY] = -com[YY];
1855 transvec[ZZ] = -com[ZZ];
1857 /* Subtract the center of mass from the copy */
1858 translate_x(loc->xcopy, edi->sref.nr, transvec);
1860 /* Determine the rotation matrix */
1861 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1865 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1866 int nat, /* How many positions are there? */
1867 rvec transvec, /* The translation vector */
1868 matrix rotmat) /* The rotation matrix */
1871 translate_x(x, nat, transvec);
1874 rotate_x(x, nat, rotmat);
1878 /* Gets the rms deviation of the positions to the structure s */
1879 /* fit_to_structure has to be called before calling this routine! */
1880 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1881 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1887 for (i = 0; i < s->nr; i++)
1889 rmsd += distance2(s->x[i], x[i]);
1892 rmsd /= (real) s->nr;
1899 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1904 if (ed->eEDtype != eEDnone)
1906 /* Loop over ED groups */
1910 /* Local atoms of the reference structure (for fitting), need only be assembled
1911 * if their indices differ from the average ones */
1914 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1915 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1918 /* Local atoms of the average structure (on these ED will be performed) */
1919 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1920 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1922 /* Indicate that the ED shift vectors for this structure need to be updated
1923 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1924 edi->buf->do_edsam->bUpdateShifts = TRUE;
1926 /* Set the pointer to the next ED group (if any) */
1927 edi = edi->next_edi;
1933 static gmx_inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1944 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1945 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1946 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1950 xu[XX] = x[XX]-tx*box[XX][XX];
1951 xu[YY] = x[YY]-ty*box[YY][YY];
1952 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1957 static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
1964 /* loop over linfix vectors */
1965 for (i = 0; i < edi->vecs.linfix.neig; i++)
1967 /* calculate the projection */
1968 proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1970 /* calculate the correction */
1971 add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1973 /* apply the correction */
1974 add /= edi->sav.sqrtm[i];
1975 for (j = 0; j < edi->sav.nr; j++)
1977 svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1978 rvec_inc(xcoll[j], vec_dum);
1984 static void do_linacc(rvec *xcoll, t_edpar *edi)
1991 /* loop over linacc vectors */
1992 for (i = 0; i < edi->vecs.linacc.neig; i++)
1994 /* calculate the projection */
1995 proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1997 /* calculate the correction */
1999 if (edi->vecs.linacc.stpsz[i] > 0.0)
2001 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
2003 add = edi->vecs.linacc.refproj[i] - proj;
2006 if (edi->vecs.linacc.stpsz[i] < 0.0)
2008 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
2010 add = edi->vecs.linacc.refproj[i] - proj;
2014 /* apply the correction */
2015 add /= edi->sav.sqrtm[i];
2016 for (j = 0; j < edi->sav.nr; j++)
2018 svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
2019 rvec_inc(xcoll[j], vec_dum);
2022 /* new positions will act as reference */
2023 edi->vecs.linacc.refproj[i] = proj + add;
2028 static void do_radfix(rvec *xcoll, t_edpar *edi)
2031 real *proj, rad = 0.0, ratio;
2035 if (edi->vecs.radfix.neig == 0)
2040 snew(proj, edi->vecs.radfix.neig);
2042 /* loop over radfix vectors */
2043 for (i = 0; i < edi->vecs.radfix.neig; i++)
2045 /* calculate the projections, radius */
2046 proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
2047 rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
2051 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2052 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2054 /* loop over radfix vectors */
2055 for (i = 0; i < edi->vecs.radfix.neig; i++)
2057 proj[i] -= edi->vecs.radfix.refproj[i];
2059 /* apply the correction */
2060 proj[i] /= edi->sav.sqrtm[i];
2062 for (j = 0; j < edi->sav.nr; j++)
2064 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2065 rvec_inc(xcoll[j], vec_dum);
2073 static void do_radacc(rvec *xcoll, t_edpar *edi)
2076 real *proj, rad = 0.0, ratio = 0.0;
2080 if (edi->vecs.radacc.neig == 0)
2085 snew(proj, edi->vecs.radacc.neig);
2087 /* loop over radacc vectors */
2088 for (i = 0; i < edi->vecs.radacc.neig; i++)
2090 /* calculate the projections, radius */
2091 proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2092 rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
2096 /* only correct when radius decreased */
2097 if (rad < edi->vecs.radacc.radius)
2099 ratio = edi->vecs.radacc.radius/rad - 1.0;
2100 rad = edi->vecs.radacc.radius;
2104 edi->vecs.radacc.radius = rad;
2107 /* loop over radacc vectors */
2108 for (i = 0; i < edi->vecs.radacc.neig; i++)
2110 proj[i] -= edi->vecs.radacc.refproj[i];
2112 /* apply the correction */
2113 proj[i] /= edi->sav.sqrtm[i];
2115 for (j = 0; j < edi->sav.nr; j++)
2117 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2118 rvec_inc(xcoll[j], vec_dum);
2125 struct t_do_radcon {
2129 static void do_radcon(rvec *xcoll, t_edpar *edi)
2132 real rad = 0.0, ratio = 0.0;
2133 struct t_do_radcon *loc;
2138 if (edi->buf->do_radcon != NULL)
2141 loc = edi->buf->do_radcon;
2146 snew(edi->buf->do_radcon, 1);
2148 loc = edi->buf->do_radcon;
2150 if (edi->vecs.radcon.neig == 0)
2157 snew(loc->proj, edi->vecs.radcon.neig);
2160 /* loop over radcon vectors */
2161 for (i = 0; i < edi->vecs.radcon.neig; i++)
2163 /* calculate the projections, radius */
2164 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2165 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2168 /* only correct when radius increased */
2169 if (rad > edi->vecs.radcon.radius)
2171 ratio = edi->vecs.radcon.radius/rad - 1.0;
2173 /* loop over radcon vectors */
2174 for (i = 0; i < edi->vecs.radcon.neig; i++)
2176 /* apply the correction */
2177 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2178 loc->proj[i] /= edi->sav.sqrtm[i];
2179 loc->proj[i] *= ratio;
2181 for (j = 0; j < edi->sav.nr; j++)
2183 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2184 rvec_inc(xcoll[j], vec_dum);
2190 edi->vecs.radcon.radius = rad;
2193 if (rad != edi->vecs.radcon.radius)
2196 for (i = 0; i < edi->vecs.radcon.neig; i++)
2198 /* calculate the projections, radius */
2199 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2200 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2207 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
2212 /* subtract the average positions */
2213 for (i = 0; i < edi->sav.nr; i++)
2215 rvec_dec(xcoll[i], edi->sav.x[i]);
2218 /* apply the constraints */
2221 do_linfix(xcoll, edi, step);
2223 do_linacc(xcoll, edi);
2226 do_radfix(xcoll, edi);
2228 do_radacc(xcoll, edi);
2229 do_radcon(xcoll, edi);
2231 /* add back the average positions */
2232 for (i = 0; i < edi->sav.nr; i++)
2234 rvec_inc(xcoll[i], edi->sav.x[i]);
2239 /* Write out the projections onto the eigenvectors. The order of output
2240 * corresponds to ed_output_legend() */
2241 static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
2246 /* Output how well we fit to the reference structure */
2247 fprintf(fp, EDcol_ffmt, rmsd);
2249 for (i = 0; i < edi->vecs.mon.neig; i++)
2251 fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
2254 for (i = 0; i < edi->vecs.linfix.neig; i++)
2256 fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
2259 for (i = 0; i < edi->vecs.linacc.neig; i++)
2261 fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
2264 for (i = 0; i < edi->vecs.radfix.neig; i++)
2266 fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
2268 if (edi->vecs.radfix.neig)
2270 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
2273 for (i = 0; i < edi->vecs.radacc.neig; i++)
2275 fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
2277 if (edi->vecs.radacc.neig)
2279 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
2282 for (i = 0; i < edi->vecs.radcon.neig; i++)
2284 fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
2286 if (edi->vecs.radcon.neig)
2288 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
2292 /* Returns if any constraints are switched on */
2293 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2295 if (edtype == eEDedsam || edtype == eEDflood)
2297 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2298 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2299 edi->vecs.radcon.neig);
2305 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2306 * umbrella sampling simulations. */
2307 static void copyEvecReference(t_eigvec* floodvecs)
2312 if (NULL == floodvecs->refproj0)
2314 snew(floodvecs->refproj0, floodvecs->neig);
2317 for (i = 0; i < floodvecs->neig; i++)
2319 floodvecs->refproj0[i] = floodvecs->refproj[i];
2324 /* Call on MASTER only. Check whether the essential dynamics / flooding
2325 * groups of the checkpoint file are consistent with the provided .edi file. */
2326 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate)
2328 t_edpar *edi = NULL; /* points to a single edi data set */
2332 if (NULL == EDstate->nref || NULL == EDstate->nav)
2334 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2335 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2336 "it must also continue with/without ED constraints when checkpointing.\n"
2337 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2338 "from without a checkpoint.\n");
2345 /* Check number of atoms in the reference and average structures */
2346 if (EDstate->nref[edinum] != edi->sref.nr)
2348 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2349 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2350 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2352 if (EDstate->nav[edinum] != edi->sav.nr)
2354 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2355 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2356 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2358 edi = edi->next_edi;
2362 if (edinum != EDstate->nED)
2364 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2365 "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2366 "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2371 /* The edsamstate struct stores the information we need to make the ED group
2372 * whole again after restarts from a checkpoint file. Here we do the following:
2373 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2374 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2375 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2376 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2377 * all ED structures at the last time step. */
2378 static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate)
2384 snew(EDstate->old_sref_p, EDstate->nED);
2385 snew(EDstate->old_sav_p, EDstate->nED);
2387 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2388 if (!EDstate->bFromCpt)
2390 snew(EDstate->nref, EDstate->nED);
2391 snew(EDstate->nav, EDstate->nED);
2394 /* Loop over all ED/flooding data sets (usually only one, though) */
2396 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2398 /* We always need the last reference and average positions such that
2399 * in the next time step we can make the ED group whole again
2400 * if the atoms do not have the correct PBC representation */
2401 if (EDstate->bFromCpt)
2403 /* Copy the last whole positions of reference and average group from .cpt */
2404 for (i = 0; i < edi->sref.nr; i++)
2406 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2408 for (i = 0; i < edi->sav.nr; i++)
2410 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2415 EDstate->nref[nr_edi-1] = edi->sref.nr;
2416 EDstate->nav [nr_edi-1] = edi->sav.nr;
2419 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2420 EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2421 EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2423 edi = edi->next_edi;
2428 /* Adds 'buf' to 'str' */
2429 static void add_to_string(char **str, char *buf)
2434 len = strlen(*str) + strlen(buf) + 1;
2440 static void add_to_string_aligned(char **str, char *buf)
2442 char buf_aligned[STRLEN];
2444 sprintf(buf_aligned, EDcol_sfmt, buf);
2445 add_to_string(str, buf_aligned);
2449 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, char *value, char *unit, char EDgroupchar)
2451 char tmp[STRLEN], tmp2[STRLEN];
2454 sprintf(tmp, "%c %s", EDgroupchar, value);
2455 add_to_string_aligned(LegendStr, tmp);
2456 sprintf(tmp2, "%s (%s)", tmp, unit);
2457 (*setname)[*nsets] = gmx_strdup(tmp2);
2462 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2468 for (i = 0; i < evec->neig; i++)
2470 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2471 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2476 /* Makes a legend for the xvg output file. Call on MASTER only! */
2477 static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv)
2479 t_edpar *edi = NULL;
2481 int nr_edi, nsets, n_flood, n_edsam;
2482 const char **setname;
2484 char *LegendStr = NULL;
2489 fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2491 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2493 fprintf(ed->edo, "#\n");
2494 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2495 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2496 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2497 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2498 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2499 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2500 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2501 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2502 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2504 if (edi->flood.vecs.neig)
2506 /* If in any of the groups we find a flooding vector, flooding is turned on */
2507 ed->eEDtype = eEDflood;
2509 /* Print what flavor of flooding we will do */
2510 if (0 == edi->flood.tau) /* constant flooding strength */
2512 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2513 if (edi->flood.bHarmonic)
2515 fprintf(ed->edo, ", harmonic");
2518 else /* adaptive flooding */
2520 fprintf(ed->edo, ", adaptive");
2523 fprintf(ed->edo, "\n");
2525 edi = edi->next_edi;
2528 /* Print a nice legend */
2530 LegendStr[0] = '\0';
2531 sprintf(buf, "# %6s", "time");
2532 add_to_string(&LegendStr, buf);
2534 /* Calculate the maximum number of columns we could end up with */
2537 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2539 nsets += 5 +edi->vecs.mon.neig
2540 +edi->vecs.linfix.neig
2541 +edi->vecs.linacc.neig
2542 +edi->vecs.radfix.neig
2543 +edi->vecs.radacc.neig
2544 +edi->vecs.radcon.neig
2545 + 6*edi->flood.vecs.neig;
2546 edi = edi->next_edi;
2548 snew(setname, nsets);
2550 /* In the mdrun time step in a first function call (do_flood()) the flooding
2551 * forces are calculated and in a second function call (do_edsam()) the
2552 * ED constraints. To get a corresponding legend, we need to loop twice
2553 * over the edi groups and output first the flooding, then the ED part */
2555 /* The flooding-related legend entries, if flooding is done */
2557 if (eEDflood == ed->eEDtype)
2560 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2562 /* Always write out the projection on the flooding EVs. Of course, this can also
2563 * be achieved with the monitoring option in do_edsam() (if switched on by the
2564 * user), but in that case the positions need to be communicated in do_edsam(),
2565 * which is not necessary when doing flooding only. */
2566 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2568 for (i = 0; i < edi->flood.vecs.neig; i++)
2570 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2571 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2573 /* Output the current reference projection if it changes with time;
2574 * this can happen when flooding is used as harmonic restraint */
2575 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
2577 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2578 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2581 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2582 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2584 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2585 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2588 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2589 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2591 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2593 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2594 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2597 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2598 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2601 edi = edi->next_edi;
2602 } /* End of flooding-related legend entries */
2606 /* Now the ED-related entries, if essential dynamics is done */
2608 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2610 if (bNeedDoEdsam(edi)) /* Only print ED legend if at least one ED option is on */
2612 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2614 /* Essential dynamics, projections on eigenvectors */
2615 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2616 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2617 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2618 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2619 if (edi->vecs.radfix.neig)
2621 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2623 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2624 if (edi->vecs.radacc.neig)
2626 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2628 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2629 if (edi->vecs.radcon.neig)
2631 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2634 edi = edi->next_edi;
2635 } /* end of 'pure' essential dynamics legend entries */
2636 n_edsam = nsets - n_flood;
2638 xvgr_legend(ed->edo, nsets, setname, oenv);
2641 fprintf(ed->edo, "#\n"
2642 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2643 n_flood, 1 == n_flood ? "" : "s",
2644 n_edsam, 1 == n_edsam ? "" : "s");
2645 fprintf(ed->edo, "%s", LegendStr);
2652 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2653 * contained in the input file, creates a NULL terminated list of t_edpar structures */
2654 void init_edsam(gmx_mtop_t *mtop,
2660 edsamstate_t *EDstate)
2662 t_edpar *edi = NULL; /* points to a single edi data set */
2663 int i, nr_edi, avindex;
2664 rvec *x_pbc = NULL; /* positions of the whole MD system with pbc removed */
2665 rvec *xfit = NULL, *xstart = NULL; /* dummy arrays to determine initial RMSDs */
2666 rvec fit_transvec; /* translation ... */
2667 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2668 rvec *ref_x_old = NULL; /* helper pointer */
2672 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2676 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n"
2677 "flooding simulation. Please also provide the correct .edi file with -ei.\n");
2681 /* Needed for initializing radacc radius in do_edsam */
2684 /* The input file is read by the master and the edi structures are
2685 * initialized here. Input is stored in ed->edpar. Then the edi
2686 * structures are transferred to the other nodes */
2689 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2690 * flooding vector, Essential dynamics can be applied to more than one structure
2691 * as well, but will be done in the order given in the edi file, so
2692 * expect different results for different order of edi file concatenation! */
2696 init_edi(mtop, edi);
2697 init_flood(edi, ed, ir->delta_t);
2698 edi = edi->next_edi;
2702 /* The master does the work here. The other nodes get the positions
2703 * not before dd_partition_system which is called after init_edsam */
2706 if (!EDstate->bFromCpt)
2708 /* Remove PBC, make molecule(s) subject to ED whole. */
2709 snew(x_pbc, mtop->natoms);
2710 m_rveccopy(mtop->natoms, x, x_pbc);
2711 do_pbc_first_mtop(NULL, ir->ePBC, box, mtop, x_pbc);
2713 /* Reset pointer to first ED data set which contains the actual ED data */
2715 /* Loop over all ED/flooding data sets (usually only one, though) */
2716 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2718 /* For multiple ED groups we use the output frequency that was specified
2719 * in the first set */
2722 edi->outfrq = ed->edpar->outfrq;
2725 /* Extract the initial reference and average positions. When starting
2726 * from .cpt, these have already been read into sref.x_old
2727 * in init_edsamstate() */
2728 if (!EDstate->bFromCpt)
2730 /* If this is the first run (i.e. no checkpoint present) we assume
2731 * that the starting positions give us the correct PBC representation */
2732 for (i = 0; i < edi->sref.nr; i++)
2734 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2737 for (i = 0; i < edi->sav.nr; i++)
2739 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2743 /* Now we have the PBC-correct start positions of the reference and
2744 average structure. We copy that over to dummy arrays on which we
2745 can apply fitting to print out the RMSD. We srenew the memory since
2746 the size of the buffers is likely different for every ED group */
2747 srenew(xfit, edi->sref.nr );
2748 srenew(xstart, edi->sav.nr );
2751 /* Reference indices are the same as average indices */
2752 ref_x_old = edi->sav.x_old;
2756 ref_x_old = edi->sref.x_old;
2758 copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2759 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2761 /* Make the fit to the REFERENCE structure, get translation and rotation */
2762 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2764 /* Output how well we fit to the reference at the start */
2765 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2766 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2767 rmsd_from_structure(xfit, &edi->sref));
2768 if (EDstate->nED > 1)
2770 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2772 fprintf(stderr, "\n");
2774 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2775 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2777 /* calculate initial projections */
2778 project(xstart, edi);
2780 /* For the target and origin structure both a reference (fit) and an
2781 * average structure can be provided in make_edi. If both structures
2782 * are the same, make_edi only stores one of them in the .edi file.
2783 * If they differ, first the fit and then the average structure is stored
2784 * in star (or sor), thus the number of entries in star/sor is
2785 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2786 * the size of the average group. */
2788 /* process target structure, if required */
2789 if (edi->star.nr > 0)
2791 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2793 /* get translation & rotation for fit of target structure to reference structure */
2794 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2796 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2797 if (edi->star.nr == edi->sav.nr)
2801 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2803 /* The last sav.nr indices of the target structure correspond to
2804 * the average structure, which must be projected */
2805 avindex = edi->star.nr - edi->sav.nr;
2807 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
2811 rad_project(edi, xstart, &edi->vecs.radcon);
2814 /* process structure that will serve as origin of expansion circle */
2815 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2817 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2820 if (edi->sori.nr > 0)
2822 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2824 /* fit this structure to reference structure */
2825 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2827 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2828 if (edi->sori.nr == edi->sav.nr)
2832 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2834 /* For the projection, we need the last sav.nr indices of sori */
2835 avindex = edi->sori.nr - edi->sav.nr;
2838 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2839 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2840 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2842 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2843 /* Set center of flooding potential to the ORIGIN structure */
2844 rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
2845 /* We already know that no (moving) reference position was provided,
2846 * therefore we can overwrite refproj[0]*/
2847 copyEvecReference(&edi->flood.vecs);
2850 else /* No origin structure given */
2852 rad_project(edi, xstart, &edi->vecs.radacc);
2853 rad_project(edi, xstart, &edi->vecs.radfix);
2854 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2856 if (edi->flood.bHarmonic)
2858 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2859 for (i = 0; i < edi->flood.vecs.neig; i++)
2861 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2866 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2867 /* Set center of flooding potential to the center of the covariance matrix,
2868 * i.e. the average structure, i.e. zero in the projected system */
2869 for (i = 0; i < edi->flood.vecs.neig; i++)
2871 edi->flood.vecs.refproj[i] = 0.0;
2876 /* For convenience, output the center of the flooding potential for the eigenvectors */
2877 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2879 for (i = 0; i < edi->flood.vecs.neig; i++)
2881 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2882 if (edi->flood.bHarmonic)
2884 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2886 fprintf(stdout, "\n");
2890 /* set starting projections for linsam */
2891 rad_project(edi, xstart, &edi->vecs.linacc);
2892 rad_project(edi, xstart, &edi->vecs.linfix);
2894 /* Prepare for the next edi data set: */
2895 edi = edi->next_edi;
2897 /* Cleaning up on the master node: */
2898 if (!EDstate->bFromCpt)
2905 } /* end of MASTER only section */
2909 /* First let everybody know how many ED data sets to expect */
2910 gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr);
2911 /* Broadcast the essential dynamics / flooding data to all nodes */
2912 broadcast_ed_data(cr, ed, EDstate->nED);
2916 /* In the single-CPU case, point the local atom numbers pointers to the global
2917 * one, so that we can use the same notation in serial and parallel case: */
2919 /* Loop over all ED data sets (usually only one, though) */
2921 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2923 edi->sref.anrs_loc = edi->sref.anrs;
2924 edi->sav.anrs_loc = edi->sav.anrs;
2925 edi->star.anrs_loc = edi->star.anrs;
2926 edi->sori.anrs_loc = edi->sori.anrs;
2927 /* For the same reason as above, make a dummy c_ind array: */
2928 snew(edi->sav.c_ind, edi->sav.nr);
2929 /* Initialize the array */
2930 for (i = 0; i < edi->sav.nr; i++)
2932 edi->sav.c_ind[i] = i;
2934 /* In the general case we will need a different-sized array for the reference indices: */
2937 snew(edi->sref.c_ind, edi->sref.nr);
2938 for (i = 0; i < edi->sref.nr; i++)
2940 edi->sref.c_ind[i] = i;
2943 /* Point to the very same array in case of other structures: */
2944 edi->star.c_ind = edi->sav.c_ind;
2945 edi->sori.c_ind = edi->sav.c_ind;
2946 /* In the serial case, the local number of atoms is the global one: */
2947 edi->sref.nr_loc = edi->sref.nr;
2948 edi->sav.nr_loc = edi->sav.nr;
2949 edi->star.nr_loc = edi->star.nr;
2950 edi->sori.nr_loc = edi->sori.nr;
2952 /* An on we go to the next ED group */
2953 edi = edi->next_edi;
2957 /* Allocate space for ED buffer variables */
2958 /* Again, loop over ED data sets */
2960 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2962 /* Allocate space for ED buffer variables */
2963 snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
2964 snew(edi->buf->do_edsam, 1);
2966 /* Space for collective ED buffer variables */
2968 /* Collective positions of atoms with the average indices */
2969 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2970 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
2971 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2972 /* Collective positions of atoms with the reference indices */
2975 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2976 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
2977 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2980 /* Get memory for flooding forces */
2981 snew(edi->flood.forces_cartesian, edi->sav.nr);
2984 /* Dump it all into one file per process */
2985 dump_edi(edi, cr, nr_edi);
2989 edi = edi->next_edi;
2992 /* Flush the edo file so that the user can check some things
2993 * when the simulation has started */
3001 void do_edsam(t_inputrec *ir,
3009 int i, edinr, iupdate = 500;
3010 matrix rotmat; /* rotation matrix */
3011 rvec transvec; /* translation vector */
3012 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3013 real dt_1; /* 1/dt */
3014 struct t_do_edsam *buf;
3016 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3017 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3020 /* Check if ED sampling has to be performed */
3021 if (ed->eEDtype == eEDnone)
3026 /* Suppress output on first call of do_edsam if
3027 * two-step sd2 integrator is used */
3028 if ( (ir->eI == eiSD2) && (v != NULL) )
3033 dt_1 = 1.0/ir->delta_t;
3035 /* Loop over all ED groups (usually one) */
3041 if (bNeedDoEdsam(edi))
3044 buf = edi->buf->do_edsam;
3048 /* initialize radacc radius for slope criterion */
3049 buf->oldrad = calc_radius(&edi->vecs.radacc);
3052 /* Copy the positions into buf->xc* arrays and after ED
3053 * feed back corrections to the official positions */
3055 /* Broadcast the ED positions such that every node has all of them
3056 * Every node contributes its local positions xs and stores it in
3057 * the collective buf->xcoll array. Note that for edinr > 1
3058 * xs could already have been modified by an earlier ED */
3060 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3061 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
3063 /* Only assembly reference positions if their indices differ from the average ones */
3066 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3067 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3070 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3071 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3072 * set bUpdateShifts=TRUE in the parallel case. */
3073 buf->bUpdateShifts = FALSE;
3075 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3076 * as well as the indices in edi->sav.anrs */
3078 /* Fit the reference indices to the reference structure */
3081 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3085 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3088 /* Now apply the translation and rotation to the ED structure */
3089 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3091 /* Find out how well we fit to the reference (just for output steps) */
3092 if (do_per_step(step, edi->outfrq) && MASTER(cr))
3096 /* Indices of reference and average structures are identical,
3097 * thus we can calculate the rmsd to SREF using xcoll */
3098 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3102 /* We have to translate & rotate the reference atoms first */
3103 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3104 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3108 /* update radsam references, when required */
3109 if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3111 project(buf->xcoll, edi);
3112 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3113 rad_project(edi, buf->xcoll, &edi->vecs.radfix);
3114 buf->oldrad = -1.e5;
3117 /* update radacc references, when required */
3118 if (do_per_step(step, iupdate) && step >= edi->presteps)
3120 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
3121 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3123 project(buf->xcoll, edi);
3124 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3129 buf->oldrad = edi->vecs.radacc.radius;
3133 /* apply the constraints */
3134 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3136 /* ED constraints should be applied already in the first MD step
3137 * (which is step 0), therefore we pass step+1 to the routine */
3138 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3141 /* write to edo, when required */
3142 if (do_per_step(step, edi->outfrq))
3144 project(buf->xcoll, edi);
3145 if (MASTER(cr) && !bSuppress)
3147 write_edo(edi, ed->edo, rmsdev);
3151 /* Copy back the positions unless monitoring only */
3152 if (ed_constraints(ed->eEDtype, edi))
3154 /* remove fitting */
3155 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3157 /* Copy the ED corrected positions into the coordinate array */
3158 /* Each node copies its local part. In the serial case, nat_loc is the
3159 * total number of ED atoms */
3160 for (i = 0; i < edi->sav.nr_loc; i++)
3162 /* Unshift local ED coordinate and store in x_unsh */
3163 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3164 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3166 /* dx is the ED correction to the positions: */
3167 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3171 /* dv is the ED correction to the velocity: */
3172 svmul(dt_1, dx, dv);
3173 /* apply the velocity correction: */
3174 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3176 /* Finally apply the position correction due to ED: */
3177 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3180 } /* END of if ( bNeedDoEdsam(edi) ) */
3182 /* Prepare for the next ED group */
3183 edi = edi->next_edi;
3185 } /* END of loop over ED groups */