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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
59 #include "mtop_util.h"
61 #include "mpelogging.h"
64 #include "groupcoord.h"
67 /* We use the same defines as in mvdata.c here */
68 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
69 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
70 #define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
72 /* These macros determine the column width in the output file */
73 #define EDcol_sfmt "%17s"
74 #define EDcol_efmt "%17.5e"
75 #define EDcol_ffmt "%17f"
77 /* enum to identify the type of ED: none, normal ED, flooding */
79 eEDnone, eEDedsam, eEDflood, eEDnr
82 /* enum to identify operations on reference, average, origin, target structures */
84 eedREF, eedAV, eedORI, eedTAR, eedNR
90 int neig; /* nr of eigenvectors */
91 int *ieig; /* index nrs of eigenvectors */
92 real *stpsz; /* stepsizes (per eigenvector) */
93 rvec **vec; /* eigenvector components */
94 real *xproj; /* instantaneous x projections */
95 real *fproj; /* instantaneous f projections */
96 real radius; /* instantaneous radius */
97 real *refproj; /* starting or target projecions */
98 /* When using flooding as harmonic restraint: The current reference projection
99 * is at each step calculated from the initial refproj0 and the slope. */
100 real *refproj0, *refprojslope;
106 t_eigvec mon; /* only monitored, no constraints */
107 t_eigvec linfix; /* fixed linear constraints */
108 t_eigvec linacc; /* acceptance linear constraints */
109 t_eigvec radfix; /* fixed radial constraints (exp) */
110 t_eigvec radacc; /* acceptance radial constraints (exp) */
111 t_eigvec radcon; /* acceptance rad. contraction constr. */
118 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
120 gmx_bool bConstForce; /* Do not calculate a flooding potential,
121 instead flood with a constant force */
130 rvec *forces_cartesian;
131 t_eigvec vecs; /* use flooding for these */
135 /* This type is for the average, reference, target, and origin structure */
136 typedef struct gmx_edx
138 int nr; /* number of atoms this structure contains */
139 int nr_loc; /* number of atoms on local node */
140 int *anrs; /* atom index numbers */
141 int *anrs_loc; /* local atom index numbers */
142 int nalloc_loc; /* allocation size of anrs_loc */
143 int *c_ind; /* at which position of the whole anrs
144 * array is a local atom?, i.e.
145 * c_ind[0...nr_loc-1] gives the atom index
146 * with respect to the collective
147 * anrs[0...nr-1] array */
148 rvec *x; /* positions for this structure */
149 rvec *x_old; /* Last positions which have the correct PBC
150 representation of the ED group. In
151 combination with keeping track of the
152 shift vectors, the ED group can always
154 real *m; /* masses */
155 real mtot; /* total mass (only used in sref) */
156 real *sqrtm; /* sqrt of the masses used for mass-
157 * weighting of analysis (only used in sav) */
163 int nini; /* total Nr of atoms */
164 gmx_bool fitmas; /* true if trans fit with cm */
165 gmx_bool pcamas; /* true if mass-weighted PCA */
166 int presteps; /* number of steps to run without any
167 * perturbations ... just monitoring */
168 int outfrq; /* freq (in steps) of writing to edo */
169 int maxedsteps; /* max nr of steps per cycle */
171 /* all gmx_edx datasets are copied to all nodes in the parallel case */
172 struct gmx_edx sref; /* reference positions, to these fitting
174 gmx_bool bRefEqAv; /* If true, reference & average indices
175 * are the same. Used for optimization */
176 struct gmx_edx sav; /* average positions */
177 struct gmx_edx star; /* target positions */
178 struct gmx_edx sori; /* origin positions */
180 t_edvecs vecs; /* eigenvectors */
181 real slope; /* minimal slope in acceptance radexp */
183 gmx_bool bNeedDoEdsam; /* if any of the options mon, linfix, ...
184 * is used (i.e. apart from flooding) */
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, edsamstate_t *EDstate);
236 static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate);
237 static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv);
238 /* End function declarations */
241 /* Multiple ED groups will be labeled with letters instead of numbers
242 * to avoid confusion with eigenvector indices */
243 static char get_EDgroupChar(int nr_edi, int nED)
251 * nr_edi = 2 -> B ...
253 return 'A' + nr_edi - 1;
257 /* Does not subtract average positions, projection on single eigenvector is returned
258 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
259 * Average position is subtracted in ed_apply_constraints prior to calling projectx
261 static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec)
267 for (i = 0; i < edi->sav.nr; i++)
269 proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
276 /* Specialized: projection is stored in vec->refproj
277 * -> used for radacc, radfix, radcon and center of flooding potential
278 * subtracts average positions, projects vector x */
279 static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec)
284 /* Subtract average positions */
285 for (i = 0; i < edi->sav.nr; i++)
287 rvec_dec(x[i], edi->sav.x[i]);
290 for (i = 0; i < vec->neig; i++)
292 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
293 rad += pow((vec->refproj[i]-vec->xproj[i]), 2);
295 vec->radius = sqrt(rad);
297 /* Add average positions */
298 for (i = 0; i < edi->sav.nr; i++)
300 rvec_inc(x[i], edi->sav.x[i]);
305 /* Project vector x, subtract average positions prior to projection and add
306 * them afterwards to retain the unchanged vector. Store in xproj. Mass-weighting
308 static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
309 t_eigvec *vec, /* The eigenvectors */
320 /* Subtract average positions */
321 for (i = 0; i < edi->sav.nr; i++)
323 rvec_dec(x[i], edi->sav.x[i]);
326 for (i = 0; i < vec->neig; i++)
328 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
331 /* Add average positions */
332 for (i = 0; i < edi->sav.nr; i++)
334 rvec_inc(x[i], edi->sav.x[i]);
339 /* Project vector x onto all edi->vecs (mon, linfix,...) */
340 static void project(rvec *x, /* positions to project */
341 t_edpar *edi) /* edi data set */
343 /* It is not more work to subtract the average position in every
344 * subroutine again, because these routines are rarely used simultanely */
345 project_to_eigvectors(x, &edi->vecs.mon, edi);
346 project_to_eigvectors(x, &edi->vecs.linfix, edi);
347 project_to_eigvectors(x, &edi->vecs.linacc, edi);
348 project_to_eigvectors(x, &edi->vecs.radfix, edi);
349 project_to_eigvectors(x, &edi->vecs.radacc, edi);
350 project_to_eigvectors(x, &edi->vecs.radcon, edi);
354 static real calc_radius(t_eigvec *vec)
360 for (i = 0; i < vec->neig; i++)
362 rad += pow((vec->refproj[i]-vec->xproj[i]), 2);
365 return rad = sqrt(rad);
371 static void dump_xcoll(t_edpar *edi, struct t_do_edsam *buf, t_commrec *cr,
378 ivec *shifts, *eshifts;
387 shifts = buf->shifts_xcoll;
388 eshifts = buf->extra_shifts_xcoll;
390 sprintf(fn, "xcolldump_step%d.txt", step);
393 for (i = 0; i < edi->sav.nr; i++)
395 fprintf(fp, "%d %9.5f %9.5f %9.5f %d %d %d %d %d %d\n",
397 xcoll[i][XX], xcoll[i][YY], xcoll[i][ZZ],
398 shifts[i][XX], shifts[i][YY], shifts[i][ZZ],
399 eshifts[i][XX], eshifts[i][YY], eshifts[i][ZZ]);
407 static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
412 fprintf(out, "#%s positions:\n%d\n", name, s->nr);
418 fprintf(out, "#index, x, y, z");
421 fprintf(out, ", sqrt(m)");
423 for (i = 0; i < s->nr; i++)
425 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]);
428 fprintf(out, "%9.3f", s->sqrtm[i]);
436 static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
437 const char name[], int length)
442 fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig);
443 /* Dump the data for every eigenvector: */
444 for (i = 0; i < ev->neig; i++)
446 fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
447 ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
448 for (j = 0; j < length; j++)
450 fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
457 static void dump_edi(t_edpar *edpars, t_commrec *cr, int nr_edi)
463 sprintf(fn, "EDdump_node%d_edi%d", cr->nodeid, nr_edi);
464 out = ffopen(fn, "w");
466 fprintf(out, "#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
467 edpars->nini, edpars->fitmas, edpars->pcamas);
468 fprintf(out, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
469 edpars->outfrq, edpars->maxedsteps, edpars->slope);
470 fprintf(out, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n",
471 edpars->presteps, edpars->flood.deltaF0, edpars->flood.tau,
472 edpars->flood.constEfl, edpars->flood.alpha2);
474 /* Dump reference, average, target, origin positions */
475 dump_edi_positions(out, &edpars->sref, "REFERENCE");
476 dump_edi_positions(out, &edpars->sav, "AVERAGE" );
477 dump_edi_positions(out, &edpars->star, "TARGET" );
478 dump_edi_positions(out, &edpars->sori, "ORIGIN" );
480 /* Dump eigenvectors */
481 dump_edi_eigenvecs(out, &edpars->vecs.mon, "MONITORED", edpars->sav.nr);
482 dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX", edpars->sav.nr);
483 dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC", edpars->sav.nr);
484 dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX", edpars->sav.nr);
485 dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC", edpars->sav.nr);
486 dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON", edpars->sav.nr);
488 /* Dump flooding eigenvectors */
489 dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING", edpars->sav.nr);
491 /* Dump ed local buffer */
492 fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit );
493 fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam );
494 fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon );
501 static void dump_rotmat(FILE* out, matrix rotmat)
503 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[XX][XX], rotmat[XX][YY], rotmat[XX][ZZ]);
504 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[YY][XX], rotmat[YY][YY], rotmat[YY][ZZ]);
505 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[ZZ][XX], rotmat[ZZ][YY], rotmat[ZZ][ZZ]);
510 static void dump_rvec(FILE *out, int dim, rvec *x)
515 for (i = 0; i < dim; i++)
517 fprintf(out, "%4d %f %f %f\n", i, x[i][XX], x[i][YY], x[i][ZZ]);
523 static void dump_mat(FILE* out, int dim, double** mat)
528 fprintf(out, "MATRIX:\n");
529 for (i = 0; i < dim; i++)
531 for (j = 0; j < dim; j++)
533 fprintf(out, "%f ", mat[i][j]);
546 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
548 /* this is a copy of do_fit with some modifications */
549 int c, r, n, j, i, irot;
550 double d[6], xnr, xpc;
555 struct t_do_edfit *loc;
558 if (edi->buf->do_edfit != NULL)
565 snew(edi->buf->do_edfit, 1);
567 loc = edi->buf->do_edfit;
571 snew(loc->omega, 2*DIM);
572 snew(loc->om, 2*DIM);
573 for (i = 0; i < 2*DIM; i++)
575 snew(loc->omega[i], 2*DIM);
576 snew(loc->om[i], 2*DIM);
580 for (i = 0; (i < 6); i++)
583 for (j = 0; (j < 6); j++)
585 loc->omega[i][j] = 0;
590 /* calculate the matrix U */
592 for (n = 0; (n < natoms); n++)
594 for (c = 0; (c < DIM); c++)
597 for (r = 0; (r < DIM); r++)
605 /* construct loc->omega */
606 /* loc->omega is symmetric -> loc->omega==loc->omega' */
607 for (r = 0; (r < 6); r++)
609 for (c = 0; (c <= r); c++)
611 if ((r >= 3) && (c < 3))
613 loc->omega[r][c] = u[r-3][c];
614 loc->omega[c][r] = u[r-3][c];
618 loc->omega[r][c] = 0;
619 loc->omega[c][r] = 0;
624 /* determine h and k */
628 dump_mat(stderr, 2*DIM, loc->omega);
629 for (i = 0; i < 6; i++)
631 fprintf(stderr, "d[%d] = %f\n", i, d[i]);
635 jacobi(loc->omega, 6, d, loc->om, &irot);
639 fprintf(stderr, "IROT=0\n");
642 index = 0; /* For the compiler only */
644 for (j = 0; (j < 3); j++)
647 for (i = 0; (i < 6); i++)
656 for (i = 0; (i < 3); i++)
658 vh[j][i] = M_SQRT2*loc->om[i][index];
659 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
664 for (c = 0; (c < 3); c++)
666 for (r = 0; (r < 3); r++)
668 R[c][r] = vk[0][r]*vh[0][c]+
675 for (c = 0; (c < 3); c++)
677 for (r = 0; (r < 3); r++)
679 R[c][r] = vk[0][r]*vh[0][c]+
688 static void rmfit(int nat, rvec *xcoll, rvec transvec, matrix rotmat)
695 * The inverse rotation is described by the transposed rotation matrix */
696 transpose(rotmat, tmat);
697 rotate_x(xcoll, nat, tmat);
699 /* Remove translation */
700 vec[XX] = -transvec[XX];
701 vec[YY] = -transvec[YY];
702 vec[ZZ] = -transvec[ZZ];
703 translate_x(xcoll, nat, vec);
707 /**********************************************************************************
708 ******************** FLOODING ****************************************************
709 **********************************************************************************
711 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
712 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
713 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
715 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
716 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
718 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
719 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
721 do_flood makes a copy of the positions,
722 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
723 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
724 fit. Then do_flood adds these forces to the forcefield-forces
725 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
727 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
728 structure is projected to the system of eigenvectors and then this position in the subspace is used as
729 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
730 i.e. the average structure as given in the make_edi file.
732 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
733 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
734 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
735 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
736 further adaption is applied, Efl will stay constant at zero.
738 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
739 used as spring constants for the harmonic potential.
740 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
741 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
743 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
744 the routine read_edi_file reads all of theses flooding files.
745 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
746 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
747 edi there is no interdependence whatsoever. The forces are added together.
749 To write energies into the .edr file, call the function
750 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
752 get_flood_energies(real Vfl[],int nnames);
755 - 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.
757 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
758 two edsam files from two peptide chains
761 static void write_edo_flood(t_edpar *edi, FILE *fp, real rmsd)
766 /* Output how well we fit to the reference structure */
767 fprintf(fp, EDcol_ffmt, rmsd);
769 for (i = 0; i < edi->flood.vecs.neig; i++)
771 fprintf(fp, EDcol_efmt, edi->flood.vecs.xproj[i]);
773 /* Check whether the reference projection changes with time (this can happen
774 * in case flooding is used as harmonic restraint). If so, output the
775 * current reference projection */
776 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
778 fprintf(fp, EDcol_efmt, edi->flood.vecs.refproj[i]);
781 /* Output Efl if we are doing adaptive flooding */
782 if (0 != edi->flood.tau)
784 fprintf(fp, EDcol_efmt, edi->flood.Efl);
786 fprintf(fp, EDcol_efmt, edi->flood.Vfl);
788 /* Output deltaF if we are doing adaptive flooding */
789 if (0 != edi->flood.tau)
791 fprintf(fp, EDcol_efmt, edi->flood.deltaF);
793 fprintf(fp, EDcol_efmt, edi->flood.vecs.fproj[i]);
798 /* From flood.xproj compute the Vfl(x) at this point */
799 static real flood_energy(t_edpar *edi, gmx_large_int_t step)
801 /* compute flooding energy Vfl
802 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
803 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
804 it is already computed by make_edi and stored in stpsz[i]
806 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
813 /* Each time this routine is called (i.e. each time step), we add a small
814 * value to the reference projection. This way a harmonic restraint towards
815 * a moving reference is realized. If no value for the additive constant
816 * is provided in the edi file, the reference will not change. */
817 if (edi->flood.bHarmonic)
819 for (i = 0; i < edi->flood.vecs.neig; i++)
821 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
826 /* Compute sum which will be the exponent of the exponential */
827 for (i = 0; i < edi->flood.vecs.neig; i++)
829 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
830 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]);
833 /* Compute the Gauss function*/
834 if (edi->flood.bHarmonic)
836 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
840 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
847 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
848 static void flood_forces(t_edpar *edi)
850 /* compute the forces in the subspace of the flooding eigenvectors
851 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
854 real energy = edi->flood.Vfl;
857 if (edi->flood.bHarmonic)
859 for (i = 0; i < edi->flood.vecs.neig; i++)
861 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
866 for (i = 0; i < edi->flood.vecs.neig; i++)
868 /* if Efl is zero the forces are zero if not use the formula */
869 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;
875 /* Raise forces from subspace into cartesian space */
876 static void flood_blowup(t_edpar *edi, rvec *forces_cart)
878 /* this function lifts the forces from the subspace to the cartesian space
879 all the values not contained in the subspace are assumed to be zero and then
880 a coordinate transformation from eigenvector to cartesian vectors is performed
881 The nonexistent values don't have to be set to zero explicitly, they would occur
882 as zero valued summands, hence we just stop to compute this part of the sum.
884 for every atom we add all the contributions to this atom from all the different eigenvectors.
886 NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
887 field forces_cart prior the computation, but we compute the forces separately
888 to have them accessible for diagnostics
895 forces_sub = edi->flood.vecs.fproj;
898 /* Calculate the cartesian forces for the local atoms */
900 /* Clear forces first */
901 for (j = 0; j < edi->sav.nr_loc; j++)
903 clear_rvec(forces_cart[j]);
906 /* Now compute atomwise */
907 for (j = 0; j < edi->sav.nr_loc; j++)
909 /* Compute forces_cart[edi->sav.anrs[j]] */
910 for (eig = 0; eig < edi->flood.vecs.neig; eig++)
912 /* Force vector is force * eigenvector (compute only atom j) */
913 svmul(forces_sub[eig], edi->flood.vecs.vec[eig][edi->sav.c_ind[j]], dum);
914 /* Add this vector to the cartesian forces */
915 rvec_inc(forces_cart[j], dum);
921 /* Update the values of Efl, deltaF depending on tau and Vfl */
922 static void update_adaption(t_edpar *edi)
924 /* this function updates the parameter Efl and deltaF according to the rules given in
925 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
928 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
930 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
931 /* check if restrain (inverted flooding) -> don't let EFL become positive */
932 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
937 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
942 static void do_single_flood(
947 gmx_large_int_t step,
950 gmx_bool bNS) /* Are we in a neighbor searching step? */
953 matrix rotmat; /* rotation matrix */
954 matrix tmat; /* inverse rotation */
955 rvec transvec; /* translation vector */
957 struct t_do_edsam *buf;
960 buf = edi->buf->do_edsam;
963 /* Broadcast the positions of the AVERAGE structure such that they are known on
964 * every processor. Each node contributes its local positions x and stores them in
965 * the collective ED array buf->xcoll */
966 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
967 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
969 /* Only assembly REFERENCE positions if their indices differ from the average ones */
972 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
973 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
976 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
977 * We do not need to update the shifts until the next NS step */
978 buf->bUpdateShifts = FALSE;
980 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
981 * as well as the indices in edi->sav.anrs */
983 /* Fit the reference indices to the reference structure */
986 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
990 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
993 /* Now apply the translation and rotation to the ED structure */
994 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
996 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
997 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, edi);
999 if (FALSE == edi->flood.bConstForce)
1001 /* Compute Vfl(x) from flood.xproj */
1002 edi->flood.Vfl = flood_energy(edi, step);
1004 update_adaption(edi);
1006 /* Compute the flooding forces */
1010 /* Translate them into cartesian positions */
1011 flood_blowup(edi, edi->flood.forces_cartesian);
1013 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
1014 /* Each node rotates back its local forces */
1015 transpose(rotmat, tmat);
1016 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
1018 /* Finally add forces to the main force variable */
1019 for (i = 0; i < edi->sav.nr_loc; i++)
1021 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
1024 /* Output is written by the master process */
1025 if (do_per_step(step, edi->outfrq) && MASTER(cr))
1027 /* Output how well we fit to the reference */
1030 /* Indices of reference and average structures are identical,
1031 * thus we can calculate the rmsd to SREF using xcoll */
1032 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
1036 /* We have to translate & rotate the reference atoms first */
1037 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1038 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1041 write_edo_flood(edi, edo, rmsdev);
1046 /* Main flooding routine, called from do_force */
1047 extern void do_flood(
1048 t_commrec *cr, /* Communication record */
1049 t_inputrec *ir, /* Input record */
1050 rvec x[], /* Positions on the local processor */
1051 rvec force[], /* forcefield forces, to these the flooding forces are added */
1052 gmx_edsam_t ed, /* ed data structure contains all ED and flooding groups */
1053 matrix box, /* the box */
1054 gmx_large_int_t step, /* The relative time step since ir->init_step is already subtracted */
1055 gmx_bool bNS) /* Are we in a neighbor searching step? */
1062 /* Write time to edo, when required. Output the time anyhow since we need
1063 * it in the output file for ED constraints. */
1064 if (MASTER(cr) && do_per_step(step, edi->outfrq))
1066 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
1069 if (ed->eEDtype != eEDflood)
1076 /* Call flooding for one matrix */
1077 if (edi->flood.vecs.neig)
1079 do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS);
1081 edi = edi->next_edi;
1086 /* Called by init_edi, configure some flooding related variables and structures,
1087 * print headers to output files */
1088 static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt)
1093 edi->flood.Efl = edi->flood.constEfl;
1097 if (edi->flood.vecs.neig)
1099 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1100 ed->eEDtype = eEDflood;
1102 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
1104 if (edi->flood.bConstForce)
1106 /* We have used stpsz as a vehicle to carry the fproj values for constant
1107 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1108 * in const force flooding, fproj is never changed. */
1109 for (i = 0; i < edi->flood.vecs.neig; i++)
1111 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1113 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1114 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1122 /*********** Energy book keeping ******/
1123 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1132 srenew(names, count);
1133 sprintf(buf, "Vfl_%d", count);
1134 names[count-1] = strdup(buf);
1135 actual = actual->next_edi;
1142 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1144 /*fl has to be big enough to capture nnames-many entries*/
1153 Vfl[count-1] = actual->flood.Vfl;
1154 actual = actual->next_edi;
1157 if (nnames != count-1)
1159 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1162 /************* END of FLOODING IMPLEMENTATION ****************************/
1166 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)
1172 /* Allocate space for the ED data structure */
1175 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1176 ed->eEDtype = eEDedsam;
1180 fprintf(stderr, "ED sampling will be performed!\n");
1183 /* Read the edi input file: */
1184 nED = read_edi_file(ftp2fn(efEDI, nfile, fnm), ed->edpar, natoms);
1186 /* Make sure the checkpoint was produced in a run using this .edi file */
1187 if (EDstate->bFromCpt)
1189 crosscheck_edi_file_vs_checkpoint(ed, EDstate);
1195 init_edsamstate(ed, EDstate);
1197 /* The master opens the ED output file */
1198 if (Flags & MD_APPENDFILES)
1200 ed->edo = gmx_fio_fopen(opt2fn("-eo", nfile, fnm), "a+");
1204 ed->edo = xvgropen(opt2fn("-eo", nfile, fnm),
1205 "Essential dynamics / flooding output",
1207 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1209 /* Make a descriptive legend */
1210 write_edo_legend(ed, EDstate->nED, oenv);
1217 /* Broadcasts the structure data */
1218 static void bc_ed_positions(t_commrec *cr, struct gmx_edx *s, int stype)
1220 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1221 snew_bc(cr, s->x, s->nr ); /* Positions */
1222 nblock_bc(cr, s->nr, s->anrs );
1223 nblock_bc(cr, s->nr, s->x );
1225 /* For the average & reference structures we need an array for the collective indices,
1226 * and we need to broadcast the masses as well */
1227 if (stype == eedAV || stype == eedREF)
1229 /* We need these additional variables in the parallel case: */
1230 snew(s->c_ind, s->nr ); /* Collective indices */
1231 /* Local atom indices get assigned in dd_make_local_group_indices.
1232 * There, also memory is allocated */
1233 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1234 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1235 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1238 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1239 if (stype == eedREF)
1241 snew_bc(cr, s->m, s->nr);
1242 nblock_bc(cr, s->nr, s->m);
1245 /* For the average structure we might need the masses for mass-weighting */
1248 snew_bc(cr, s->sqrtm, s->nr);
1249 nblock_bc(cr, s->nr, s->sqrtm);
1250 snew_bc(cr, s->m, s->nr);
1251 nblock_bc(cr, s->nr, s->m);
1256 /* Broadcasts the eigenvector data */
1257 static void bc_ed_vecs(t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
1261 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1262 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1263 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1264 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1265 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1267 nblock_bc(cr, ev->neig, ev->ieig );
1268 nblock_bc(cr, ev->neig, ev->stpsz );
1269 nblock_bc(cr, ev->neig, ev->xproj );
1270 nblock_bc(cr, ev->neig, ev->fproj );
1271 nblock_bc(cr, ev->neig, ev->refproj);
1273 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1274 for (i = 0; i < ev->neig; i++)
1276 snew_bc(cr, ev->vec[i], length);
1277 nblock_bc(cr, length, ev->vec[i]);
1280 /* For harmonic restraints the reference projections can change with time */
1283 snew_bc(cr, ev->refproj0, ev->neig);
1284 snew_bc(cr, ev->refprojslope, ev->neig);
1285 nblock_bc(cr, ev->neig, ev->refproj0 );
1286 nblock_bc(cr, ev->neig, ev->refprojslope);
1291 /* Broadcasts the ED / flooding data to other nodes
1292 * and allocates memory where needed */
1293 static void broadcast_ed_data(t_commrec *cr, gmx_edsam_t ed, int numedis)
1299 /* Master lets the other nodes know if its ED only or also flooding */
1300 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1302 snew_bc(cr, ed->edpar, 1);
1303 /* Now transfer the ED data set(s) */
1305 for (nr = 0; nr < numedis; nr++)
1307 /* Broadcast a single ED data set */
1310 /* Broadcast positions */
1311 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1312 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1313 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1314 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1316 /* Broadcast eigenvectors */
1317 bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr, FALSE);
1318 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
1319 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
1320 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
1321 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
1322 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
1323 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1324 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
1326 /* Set the pointer to the next ED group */
1329 snew_bc(cr, edi->next_edi, 1);
1330 edi = edi->next_edi;
1336 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1337 static void init_edi(gmx_mtop_t *mtop, t_edpar *edi)
1340 real totalmass = 0.0;
1342 gmx_mtop_atomlookup_t alook = NULL;
1345 /* NOTE Init_edi is executed on the master process only
1346 * The initialized data sets are then transmitted to the
1347 * other nodes in broadcast_ed_data */
1349 alook = gmx_mtop_atomlookup_init(mtop);
1351 /* evaluate masses (reference structure) */
1352 snew(edi->sref.m, edi->sref.nr);
1353 for (i = 0; i < edi->sref.nr; i++)
1357 gmx_mtop_atomnr_to_atom(alook, edi->sref.anrs[i], &atom);
1358 edi->sref.m[i] = atom->m;
1362 edi->sref.m[i] = 1.0;
1365 /* Check that every m > 0. Bad things will happen otherwise. */
1366 if (edi->sref.m[i] <= 0.0)
1368 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1369 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1370 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1371 "atoms from the reference structure by creating a proper index group.\n",
1372 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1375 totalmass += edi->sref.m[i];
1377 edi->sref.mtot = totalmass;
1379 /* Masses m and sqrt(m) for the average structure. Note that m
1380 * is needed if forces have to be evaluated in do_edsam */
1381 snew(edi->sav.sqrtm, edi->sav.nr );
1382 snew(edi->sav.m, edi->sav.nr );
1383 for (i = 0; i < edi->sav.nr; i++)
1385 gmx_mtop_atomnr_to_atom(alook, edi->sav.anrs[i], &atom);
1386 edi->sav.m[i] = atom->m;
1389 edi->sav.sqrtm[i] = sqrt(atom->m);
1393 edi->sav.sqrtm[i] = 1.0;
1396 /* Check that every m > 0. Bad things will happen otherwise. */
1397 if (edi->sav.sqrtm[i] <= 0.0)
1399 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1400 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1401 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1402 "atoms from the average structure by creating a proper index group.\n",
1403 i, edi->sav.anrs[i]+1, atom->m);
1407 gmx_mtop_atomlookup_destroy(alook);
1409 /* put reference structure in origin */
1410 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1414 translate_x(edi->sref.x, edi->sref.nr, com);
1416 /* Init ED buffer */
1421 static void check(const char *line, const char *label)
1423 if (!strstr(line, label))
1425 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1430 static int read_checked_edint(FILE *file, const char *label)
1432 char line[STRLEN+1];
1436 fgets2 (line, STRLEN, file);
1438 fgets2 (line, STRLEN, file);
1439 sscanf (line, "%d", &idum);
1444 static int read_edint(FILE *file, gmx_bool *bEOF)
1446 char line[STRLEN+1];
1451 eof = fgets2 (line, STRLEN, file);
1457 eof = fgets2 (line, STRLEN, file);
1463 sscanf (line, "%d", &idum);
1469 static real read_checked_edreal(FILE *file, const char *label)
1471 char line[STRLEN+1];
1475 fgets2 (line, STRLEN, file);
1477 fgets2 (line, STRLEN, file);
1478 sscanf (line, "%lf", &rdum);
1479 return (real) rdum; /* always read as double and convert to single */
1483 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1486 char line[STRLEN+1];
1490 for (i = 0; i < number; i++)
1492 fgets2 (line, STRLEN, file);
1493 sscanf (line, "%d%lf%lf%lf", &anrs[i], &d[0], &d[1], &d[2]);
1494 anrs[i]--; /* we are reading FORTRAN indices */
1495 for (j = 0; j < 3; j++)
1497 x[i][j] = d[j]; /* always read as double and convert to single */
1503 static void scan_edvec(FILE *in, int nr, rvec *vec)
1505 char line[STRLEN+1];
1510 for (i = 0; (i < nr); i++)
1512 fgets2 (line, STRLEN, in);
1513 sscanf (line, "%le%le%le", &x, &y, &z);
1521 static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1524 double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0;
1525 char line[STRLEN+1];
1528 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1531 snew(tvec->ieig, tvec->neig);
1532 snew(tvec->stpsz, tvec->neig);
1533 snew(tvec->vec, tvec->neig);
1534 snew(tvec->xproj, tvec->neig);
1535 snew(tvec->fproj, tvec->neig);
1536 snew(tvec->refproj, tvec->neig);
1539 snew(tvec->refproj0, tvec->neig);
1540 snew(tvec->refprojslope, tvec->neig);
1543 for (i = 0; (i < tvec->neig); i++)
1545 fgets2 (line, STRLEN, in);
1546 if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1548 nscan = sscanf(line, "%d%lf%lf%lf", &idum, &rdum, &refproj_dum, &refprojslope_dum);
1549 /* Zero out values which were not scanned */
1553 /* Every 4 values read, including reference position */
1554 *bHaveReference = TRUE;
1557 /* A reference position is provided */
1558 *bHaveReference = TRUE;
1559 /* No value for slope, set to 0 */
1560 refprojslope_dum = 0.0;
1563 /* No values for reference projection and slope, set to 0 */
1565 refprojslope_dum = 0.0;
1568 gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1571 tvec->refproj[i] = refproj_dum;
1572 tvec->refproj0[i] = refproj_dum;
1573 tvec->refprojslope[i] = refprojslope_dum;
1575 else /* Normal flooding */
1577 nscan = sscanf(line, "%d%lf", &idum, &rdum);
1580 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1583 tvec->ieig[i] = idum;
1584 tvec->stpsz[i] = rdum;
1585 } /* end of loop over eigenvectors */
1587 for (i = 0; (i < tvec->neig); i++)
1589 snew(tvec->vec[i], nr);
1590 scan_edvec(in, nr, tvec->vec[i]);
1596 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1597 static void read_edvecs(FILE *in, int nr, t_edvecs *vecs)
1599 gmx_bool bHaveReference = FALSE;
1602 read_edvec(in, nr, &vecs->mon, FALSE, &bHaveReference);
1603 read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1604 read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1605 read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1606 read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1607 read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1611 /* Check if the same atom indices are used for reference and average positions */
1612 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1617 /* If the number of atoms differs between the two structures,
1618 * they cannot be identical */
1619 if (sref.nr != sav.nr)
1624 /* Now that we know that both stuctures have the same number of atoms,
1625 * check if also the indices are identical */
1626 for (i = 0; i < sav.nr; i++)
1628 if (sref.anrs[i] != sav.anrs[i])
1633 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1639 static int read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, const char *fn)
1642 const int magic = 670;
1645 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1646 gmx_bool bHaveReference = FALSE;
1649 /* the edi file is not free format, so expect problems if the input is corrupt. */
1651 /* check the magic number */
1652 readmagic = read_edint(in, &bEOF);
1653 /* Check whether we have reached the end of the input file */
1659 if (readmagic != magic)
1661 if (readmagic == 666 || readmagic == 667 || readmagic == 668)
1663 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1665 else if (readmagic != 669)
1667 gmx_fatal(FARGS, "Wrong magic number %d in %s", readmagic, fn);
1671 /* check the number of atoms */
1672 edi->nini = read_edint(in, &bEOF);
1673 if (edi->nini != nr_mdatoms)
1675 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1678 /* Done checking. For the rest we blindly trust the input */
1679 edi->fitmas = read_checked_edint(in, "FITMAS");
1680 edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS");
1681 edi->outfrq = read_checked_edint(in, "OUTFRQ");
1682 edi->maxedsteps = read_checked_edint(in, "MAXLEN");
1683 edi->slope = read_checked_edreal(in, "SLOPECRIT");
1685 edi->presteps = read_checked_edint(in, "PRESTEPS");
1686 edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1687 edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1688 edi->flood.tau = read_checked_edreal(in, "TAU");
1689 edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1690 edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1691 edi->flood.kT = read_checked_edreal(in, "KT");
1692 edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1693 if (readmagic > 669)
1695 edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1699 edi->flood.bConstForce = FALSE;
1701 edi->sref.nr = read_checked_edint(in, "NREF");
1703 /* allocate space for reference positions and read them */
1704 snew(edi->sref.anrs, edi->sref.nr);
1705 snew(edi->sref.x, edi->sref.nr);
1706 snew(edi->sref.x_old, edi->sref.nr);
1707 edi->sref.sqrtm = NULL;
1708 read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1710 /* average positions. they define which atoms will be used for ED sampling */
1711 edi->sav.nr = read_checked_edint(in, "NAV");
1712 snew(edi->sav.anrs, edi->sav.nr);
1713 snew(edi->sav.x, edi->sav.nr);
1714 snew(edi->sav.x_old, edi->sav.nr);
1715 read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1717 /* Check if the same atom indices are used for reference and average positions */
1718 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1721 read_edvecs(in, edi->sav.nr, &edi->vecs);
1722 read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
1724 /* target positions */
1725 edi->star.nr = read_edint(in, &bEOF);
1726 if (edi->star.nr > 0)
1728 snew(edi->star.anrs, edi->star.nr);
1729 snew(edi->star.x, edi->star.nr);
1730 edi->star.sqrtm = NULL;
1731 read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1734 /* positions defining origin of expansion circle */
1735 edi->sori.nr = read_edint(in, &bEOF);
1736 if (edi->sori.nr > 0)
1740 /* Both an -ori structure and a at least one manual reference point have been
1741 * specified. That's ambiguous and probably not intentional. */
1742 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1743 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1745 snew(edi->sori.anrs, edi->sori.nr);
1746 snew(edi->sori.x, edi->sori.nr);
1747 edi->sori.sqrtm = NULL;
1748 read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1757 /* Read in the edi input file. Note that it may contain several ED data sets which were
1758 * achieved by concatenating multiple edi files. The standard case would be a single ED
1759 * data set, though. */
1760 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1763 t_edpar *curr_edi, *last_edi;
1768 /* This routine is executed on the master only */
1770 /* Open the .edi parameter input file */
1771 in = gmx_fio_fopen(fn, "r");
1772 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1774 /* Now read a sequence of ED input parameter sets from the edi file */
1777 while (read_edi(in, curr_edi, nr_mdatoms, fn) )
1781 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1782 /* We need to allocate space for the data: */
1784 /* Point the 'next_edi' entry to the next edi: */
1785 curr_edi->next_edi = edi_read;
1786 /* Keep the curr_edi pointer for the case that the next group is empty: */
1787 last_edi = curr_edi;
1788 /* Let's prepare to read in the next edi data set: */
1789 curr_edi = edi_read;
1793 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1796 /* Terminate the edi group list with a NULL pointer: */
1797 last_edi->next_edi = NULL;
1799 fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1801 /* Close the .edi file again */
1808 struct t_fit_to_ref {
1809 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1812 /* Fit the current positions to the reference positions
1813 * Do not actually do the fit, just return rotation and translation.
1814 * Note that the COM of the reference structure was already put into
1815 * the origin by init_edi. */
1816 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1817 rvec transvec, /* The translation vector */
1818 matrix rotmat, /* The rotation matrix */
1819 t_edpar *edi) /* Just needed for do_edfit */
1821 rvec com; /* center of mass */
1823 struct t_fit_to_ref *loc;
1826 GMX_MPE_LOG(ev_fit_to_reference_start);
1828 /* Allocate memory the first time this routine is called for each edi group */
1829 if (NULL == edi->buf->fit_to_ref)
1831 snew(edi->buf->fit_to_ref, 1);
1832 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1834 loc = edi->buf->fit_to_ref;
1836 /* We do not touch the original positions but work on a copy. */
1837 for (i = 0; i < edi->sref.nr; i++)
1839 copy_rvec(xcoll[i], loc->xcopy[i]);
1842 /* Calculate the center of mass */
1843 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1845 transvec[XX] = -com[XX];
1846 transvec[YY] = -com[YY];
1847 transvec[ZZ] = -com[ZZ];
1849 /* Subtract the center of mass from the copy */
1850 translate_x(loc->xcopy, edi->sref.nr, transvec);
1852 /* Determine the rotation matrix */
1853 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1855 GMX_MPE_LOG(ev_fit_to_reference_finish);
1859 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1860 int nat, /* How many positions are there? */
1861 rvec transvec, /* The translation vector */
1862 matrix rotmat) /* The rotation matrix */
1865 translate_x(x, nat, transvec);
1868 rotate_x(x, nat, rotmat);
1872 /* Gets the rms deviation of the positions to the structure s */
1873 /* fit_to_structure has to be called before calling this routine! */
1874 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1875 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1881 for (i = 0; i < s->nr; i++)
1883 rmsd += distance2(s->x[i], x[i]);
1886 rmsd /= (real) s->nr;
1893 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1898 if (ed->eEDtype != eEDnone)
1900 /* Loop over ED groups */
1904 /* Local atoms of the reference structure (for fitting), need only be assembled
1905 * if their indices differ from the average ones */
1908 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1909 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1912 /* Local atoms of the average structure (on these ED will be performed) */
1913 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1914 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1916 /* Indicate that the ED shift vectors for this structure need to be updated
1917 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1918 edi->buf->do_edsam->bUpdateShifts = TRUE;
1920 /* Set the pointer to the next ED group (if any) */
1921 edi = edi->next_edi;
1927 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1932 GMX_MPE_LOG(ev_unshift_start);
1940 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1941 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1942 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1946 xu[XX] = x[XX]-tx*box[XX][XX];
1947 xu[YY] = x[YY]-ty*box[YY][YY];
1948 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1951 GMX_MPE_LOG(ev_unshift_finish);
1955 static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_large_int_t step)
1962 /* loop over linfix vectors */
1963 for (i = 0; i < edi->vecs.linfix.neig; i++)
1965 /* calculate the projection */
1966 proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1968 /* calculate the correction */
1969 add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1971 /* apply the correction */
1972 add /= edi->sav.sqrtm[i];
1973 for (j = 0; j < edi->sav.nr; j++)
1975 svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1976 rvec_inc(xcoll[j], vec_dum);
1982 static void do_linacc(rvec *xcoll, t_edpar *edi)
1989 /* loop over linacc vectors */
1990 for (i = 0; i < edi->vecs.linacc.neig; i++)
1992 /* calculate the projection */
1993 proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1995 /* calculate the correction */
1997 if (edi->vecs.linacc.stpsz[i] > 0.0)
1999 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
2001 add = edi->vecs.linacc.refproj[i] - proj;
2004 if (edi->vecs.linacc.stpsz[i] < 0.0)
2006 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
2008 add = edi->vecs.linacc.refproj[i] - proj;
2012 /* apply the correction */
2013 add /= edi->sav.sqrtm[i];
2014 for (j = 0; j < edi->sav.nr; j++)
2016 svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
2017 rvec_inc(xcoll[j], vec_dum);
2020 /* new positions will act as reference */
2021 edi->vecs.linacc.refproj[i] = proj + add;
2026 static void do_radfix(rvec *xcoll, t_edpar *edi)
2029 real *proj, rad = 0.0, ratio;
2033 if (edi->vecs.radfix.neig == 0)
2038 snew(proj, edi->vecs.radfix.neig);
2040 /* loop over radfix vectors */
2041 for (i = 0; i < edi->vecs.radfix.neig; i++)
2043 /* calculate the projections, radius */
2044 proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
2045 rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
2049 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2050 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2052 /* loop over radfix vectors */
2053 for (i = 0; i < edi->vecs.radfix.neig; i++)
2055 proj[i] -= edi->vecs.radfix.refproj[i];
2057 /* apply the correction */
2058 proj[i] /= edi->sav.sqrtm[i];
2060 for (j = 0; j < edi->sav.nr; j++)
2062 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2063 rvec_inc(xcoll[j], vec_dum);
2071 static void do_radacc(rvec *xcoll, t_edpar *edi)
2074 real *proj, rad = 0.0, ratio = 0.0;
2078 if (edi->vecs.radacc.neig == 0)
2083 snew(proj, edi->vecs.radacc.neig);
2085 /* loop over radacc vectors */
2086 for (i = 0; i < edi->vecs.radacc.neig; i++)
2088 /* calculate the projections, radius */
2089 proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2090 rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
2094 /* only correct when radius decreased */
2095 if (rad < edi->vecs.radacc.radius)
2097 ratio = edi->vecs.radacc.radius/rad - 1.0;
2098 rad = edi->vecs.radacc.radius;
2102 edi->vecs.radacc.radius = rad;
2105 /* loop over radacc vectors */
2106 for (i = 0; i < edi->vecs.radacc.neig; i++)
2108 proj[i] -= edi->vecs.radacc.refproj[i];
2110 /* apply the correction */
2111 proj[i] /= edi->sav.sqrtm[i];
2113 for (j = 0; j < edi->sav.nr; j++)
2115 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2116 rvec_inc(xcoll[j], vec_dum);
2123 struct t_do_radcon {
2127 static void do_radcon(rvec *xcoll, t_edpar *edi)
2130 real rad = 0.0, ratio = 0.0;
2131 struct t_do_radcon *loc;
2136 if (edi->buf->do_radcon != NULL)
2139 loc = edi->buf->do_radcon;
2144 snew(edi->buf->do_radcon, 1);
2146 loc = edi->buf->do_radcon;
2148 if (edi->vecs.radcon.neig == 0)
2155 snew(loc->proj, edi->vecs.radcon.neig);
2158 /* loop over radcon vectors */
2159 for (i = 0; i < edi->vecs.radcon.neig; i++)
2161 /* calculate the projections, radius */
2162 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2163 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2166 /* only correct when radius increased */
2167 if (rad > edi->vecs.radcon.radius)
2169 ratio = edi->vecs.radcon.radius/rad - 1.0;
2171 /* loop over radcon vectors */
2172 for (i = 0; i < edi->vecs.radcon.neig; i++)
2174 /* apply the correction */
2175 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2176 loc->proj[i] /= edi->sav.sqrtm[i];
2177 loc->proj[i] *= ratio;
2179 for (j = 0; j < edi->sav.nr; j++)
2181 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2182 rvec_inc(xcoll[j], vec_dum);
2188 edi->vecs.radcon.radius = rad;
2191 if (rad != edi->vecs.radcon.radius)
2194 for (i = 0; i < edi->vecs.radcon.neig; i++)
2196 /* calculate the projections, radius */
2197 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2198 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2205 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step)
2210 GMX_MPE_LOG(ev_ed_apply_cons_start);
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]);
2237 GMX_MPE_LOG(ev_ed_apply_cons_finish);
2241 /* Write out the projections onto the eigenvectors. The order of output
2242 * corresponds to ed_output_legend() */
2243 static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
2248 /* Output how well we fit to the reference structure */
2249 fprintf(fp, EDcol_ffmt, rmsd);
2251 for (i = 0; i < edi->vecs.mon.neig; i++)
2253 fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
2256 for (i = 0; i < edi->vecs.linfix.neig; i++)
2258 fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
2261 for (i = 0; i < edi->vecs.linacc.neig; i++)
2263 fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
2266 for (i = 0; i < edi->vecs.radfix.neig; i++)
2268 fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
2270 if (edi->vecs.radfix.neig)
2272 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
2275 for (i = 0; i < edi->vecs.radacc.neig; i++)
2277 fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
2279 if (edi->vecs.radacc.neig)
2281 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
2284 for (i = 0; i < edi->vecs.radcon.neig; i++)
2286 fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
2288 if (edi->vecs.radcon.neig)
2290 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
2294 /* Returns if any constraints are switched on */
2295 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2297 if (edtype == eEDedsam || edtype == eEDflood)
2299 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2300 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2301 edi->vecs.radcon.neig);
2307 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2308 * umbrella sampling simulations. */
2309 static void copyEvecReference(t_eigvec* floodvecs)
2314 if (NULL == floodvecs->refproj0)
2316 snew(floodvecs->refproj0, floodvecs->neig);
2319 for (i = 0; i < floodvecs->neig; i++)
2321 floodvecs->refproj0[i] = floodvecs->refproj[i];
2326 /* Call on MASTER only. Check whether the essential dynamics / flooding
2327 * groups of the checkpoint file are consistent with the provided .edi file. */
2328 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate)
2330 t_edpar *edi = NULL; /* points to a single edi data set */
2334 if (NULL == EDstate->nref || NULL == EDstate->nav)
2336 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2337 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2338 "it must also continue with/without ED constraints when checkpointing.\n"
2339 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2340 "from without a checkpoint.\n");
2347 /* Check number of atoms in the reference and average structures */
2348 if (EDstate->nref[edinum] != edi->sref.nr)
2350 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2351 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2352 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2354 if (EDstate->nav[edinum] != edi->sav.nr)
2356 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2357 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2358 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2360 edi = edi->next_edi;
2364 if (edinum != EDstate->nED)
2366 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2367 "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2368 "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2373 /* The edsamstate struct stores the information we need to make the ED group
2374 * whole again after restarts from a checkpoint file. Here we do the following:
2375 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2376 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2377 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2378 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2379 * all ED structures at the last time step. */
2380 static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate)
2386 snew(EDstate->old_sref_p, EDstate->nED);
2387 snew(EDstate->old_sav_p, EDstate->nED);
2389 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2390 if (!EDstate->bFromCpt)
2392 snew(EDstate->nref, EDstate->nED);
2393 snew(EDstate->nav, EDstate->nED);
2396 /* Loop over all ED/flooding data sets (usually only one, though) */
2398 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2400 /* We always need the last reference and average positions such that
2401 * in the next time step we can make the ED group whole again
2402 * if the atoms do not have the correct PBC representation */
2403 if (EDstate->bFromCpt)
2405 /* Copy the last whole positions of reference and average group from .cpt */
2406 for (i = 0; i < edi->sref.nr; i++)
2408 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2410 for (i = 0; i < edi->sav.nr; i++)
2412 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2417 EDstate->nref[nr_edi-1] = edi->sref.nr;
2418 EDstate->nav [nr_edi-1] = edi->sav.nr;
2421 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2422 EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2423 EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2425 edi = edi->next_edi;
2430 /* Adds 'buf' to 'str' */
2431 static void add_to_string(char **str, char *buf)
2436 len = strlen(*str) + strlen(buf) + 1;
2442 static void add_to_string_aligned(char **str, char *buf)
2444 char buf_aligned[STRLEN];
2446 sprintf(buf_aligned, EDcol_sfmt, buf);
2447 add_to_string(str, buf_aligned);
2451 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, char *value, char *unit, char EDgroupchar)
2453 char tmp[STRLEN], tmp2[STRLEN];
2456 sprintf(tmp, "%c %s", EDgroupchar, value);
2457 add_to_string_aligned(LegendStr, tmp);
2458 sprintf(tmp2, "%s (%s)", tmp, unit);
2459 (*setname)[*nsets] = strdup(tmp2);
2464 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2470 for (i = 0; i < evec->neig; i++)
2472 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2473 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2478 /* Makes a legend for the xvg output file. Call on MASTER only! */
2479 static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv)
2481 t_edpar *edi = NULL;
2483 int nr_edi, nsets, n_flood, n_edsam;
2484 const char **setname;
2486 char *LegendStr = NULL;
2491 fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2493 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2495 /* Remember for each ED group whether we have to do essential dynamics
2496 * constraints or possibly only flooding */
2497 edi->bNeedDoEdsam = edi->vecs.mon.neig
2498 || edi->vecs.linfix.neig
2499 || edi->vecs.linacc.neig
2500 || edi->vecs.radfix.neig
2501 || edi->vecs.radacc.neig
2502 || edi->vecs.radcon.neig;
2504 fprintf(ed->edo, "#\n");
2505 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2506 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2507 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2508 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2509 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2510 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2511 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2512 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2513 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2515 if (edi->flood.vecs.neig)
2517 /* If in any of the groups we find a flooding vector, flooding is turned on */
2518 ed->eEDtype = eEDflood;
2520 /* Print what flavor of flooding we will do */
2521 if (0 == edi->flood.tau) /* constant flooding strength */
2523 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2524 if (edi->flood.bHarmonic)
2526 fprintf(ed->edo, ", harmonic");
2529 else /* adaptive flooding */
2531 fprintf(ed->edo, ", adaptive");
2534 fprintf(ed->edo, "\n");
2536 edi = edi->next_edi;
2539 /* Print a nice legend */
2541 LegendStr[0] = '\0';
2542 sprintf(buf, "# %6s", "time");
2543 add_to_string(&LegendStr, buf);
2545 /* Calculate the maximum number of columns we could end up with */
2548 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2550 nsets += 5 +edi->vecs.mon.neig
2551 +edi->vecs.linfix.neig
2552 +edi->vecs.linacc.neig
2553 +edi->vecs.radfix.neig
2554 +edi->vecs.radacc.neig
2555 +edi->vecs.radcon.neig
2556 + 6*edi->flood.vecs.neig;
2557 edi = edi->next_edi;
2559 snew(setname, nsets);
2561 /* In the mdrun time step in a first function call (do_flood()) the flooding
2562 * forces are calculated and in a second function call (do_edsam()) the
2563 * ED constraints. To get a corresponding legend, we need to loop twice
2564 * over the edi groups and output first the flooding, then the ED part */
2566 /* The flooding-related legend entries, if flooding is done */
2568 if (eEDflood == ed->eEDtype)
2571 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2573 /* Always write out the projection on the flooding EVs. Of course, this can also
2574 * be achieved with the monitoring option in do_edsam() (if switched on by the
2575 * user), but in that case the positions need to be communicated in do_edsam(),
2576 * which is not necessary when doing flooding only. */
2577 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2579 for (i = 0; i < edi->flood.vecs.neig; i++)
2581 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2582 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2584 /* Output the current reference projection if it changes with time;
2585 * this can happen when flooding is used as harmonic restraint */
2586 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
2588 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2589 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2592 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2593 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2595 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2596 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2599 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2600 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2602 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2604 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2605 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2608 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2609 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2612 edi = edi->next_edi;
2613 } /* End of flooding-related legend entries */
2617 /* Now the ED-related entries, if essential dynamics is done */
2619 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2621 if (edi->bNeedDoEdsam) /* Only print ED legend if at least one ED option is on */
2623 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2625 /* Essential dynamics, projections on eigenvectors */
2626 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2627 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2628 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2629 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2630 if (edi->vecs.radfix.neig)
2632 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2634 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2635 if (edi->vecs.radacc.neig)
2637 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2639 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2640 if (edi->vecs.radcon.neig)
2642 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2645 edi = edi->next_edi;
2646 } /* end of 'pure' essential dynamics legend entries */
2647 n_edsam = nsets - n_flood;
2649 xvgr_legend(ed->edo, nsets, setname, oenv);
2652 fprintf(ed->edo, "#\n"
2653 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2654 n_flood, 1 == n_flood ? "" : "s",
2655 n_edsam, 1 == n_edsam ? "" : "s");
2656 fprintf(ed->edo, "%s", LegendStr);
2663 void init_edsam(gmx_mtop_t *mtop, /* global topology */
2664 t_inputrec *ir, /* input record */
2665 t_commrec *cr, /* communication record */
2666 gmx_edsam_t ed, /* contains all ED data */
2667 rvec x[], /* positions of the whole MD system */
2668 matrix box, /* the box */
2669 edsamstate_t *EDstate)
2671 t_edpar *edi = NULL; /* points to a single edi data set */
2672 int i, nr_edi, avindex;
2673 rvec *x_pbc = NULL; /* positions of the whole MD system with pbc removed */
2674 rvec *xfit = NULL, *xstart = NULL; /* dummy arrays to determine initial RMSDs */
2675 rvec fit_transvec; /* translation ... */
2676 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2679 if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
2681 gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
2684 GMX_MPE_LOG(ev_edsam_start);
2688 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2692 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n"
2693 "flooding simulation. Please also provide the correct .edi file with -ei.\n");
2697 /* Needed for initializing radacc radius in do_edsam */
2700 /* The input file is read by the master and the edi structures are
2701 * initialized here. Input is stored in ed->edpar. Then the edi
2702 * structures are transferred to the other nodes */
2705 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2706 * flooding vector, Essential dynamics can be applied to more than one structure
2707 * as well, but will be done in the order given in the edi file, so
2708 * expect different results for different order of edi file concatenation! */
2712 init_edi(mtop, edi);
2713 init_flood(edi, ed, ir->delta_t);
2714 edi = edi->next_edi;
2718 /* The master does the work here. The other nodes get the positions
2719 * not before dd_partition_system which is called after init_edsam */
2722 /* Remove pbc, make molecule whole.
2723 * When ir->bContinuation=TRUE this has already been done, but ok.
2725 snew(x_pbc, mtop->natoms);
2726 m_rveccopy(mtop->natoms, x, x_pbc);
2727 do_pbc_first_mtop(NULL, ir->ePBC, box, mtop, x_pbc);
2729 /* Reset pointer to first ED data set which contains the actual ED data */
2731 /* Loop over all ED/flooding data sets (usually only one, though) */
2732 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2734 /* For multiple ED groups we use the output frequency that was specified
2735 * in the first set */
2738 edi->outfrq = ed->edpar->outfrq;
2741 /* Extract the initial reference and average positions. When starting
2742 * from .cpt, these have already been read into sref.x_old
2743 * in init_edsamstate() */
2744 if (!EDstate->bFromCpt)
2746 /* If this is the first run (i.e. no checkpoint present) we assume
2747 * that the starting positions give us the correct PBC representation */
2748 for (i = 0; i < edi->sref.nr; i++)
2750 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2753 for (i = 0; i < edi->sav.nr; i++)
2755 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2759 /* Now we have the PBC-correct start positions of the reference and
2760 average structure. We copy that over to dummy arrays on which we
2761 can apply fitting to print out the RMSD. We srenew the memory since
2762 the size of the buffers is likely different for every ED group */
2763 srenew(xfit, edi->sref.nr );
2764 srenew(xstart, edi->sav.nr );
2765 copy_rvecn(edi->sref.x_old, xfit, 0, edi->sref.nr);
2766 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2768 /* Make the fit to the REFERENCE structure, get translation and rotation */
2769 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2771 /* Output how well we fit to the reference at the start */
2772 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2773 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2774 rmsd_from_structure(xfit, &edi->sref));
2775 if (EDstate->nED > 1)
2777 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2779 fprintf(stderr, "\n");
2781 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2782 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2784 /* calculate initial projections */
2785 project(xstart, edi);
2787 /* For the target and origin structure both a reference (fit) and an
2788 * average structure can be provided in make_edi. If both structures
2789 * are the same, make_edi only stores one of them in the .edi file.
2790 * If they differ, first the fit and then the average structure is stored
2791 * in star (or sor), thus the number of entries in star/sor is
2792 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2793 * the size of the average group. */
2795 /* process target structure, if required */
2796 if (edi->star.nr > 0)
2798 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2800 /* get translation & rotation for fit of target structure to reference structure */
2801 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2803 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2804 if (edi->star.nr == edi->sav.nr)
2808 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2810 /* The last sav.nr indices of the target structure correspond to
2811 * the average structure, which must be projected */
2812 avindex = edi->star.nr - edi->sav.nr;
2814 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
2818 rad_project(edi, xstart, &edi->vecs.radcon);
2821 /* process structure that will serve as origin of expansion circle */
2822 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2824 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2827 if (edi->sori.nr > 0)
2829 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2831 /* fit this structure to reference structure */
2832 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2834 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2835 if (edi->sori.nr == edi->sav.nr)
2839 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2841 /* For the projection, we need the last sav.nr indices of sori */
2842 avindex = edi->sori.nr - edi->sav.nr;
2845 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2846 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2847 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2849 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2850 /* Set center of flooding potential to the ORIGIN structure */
2851 rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
2852 /* We already know that no (moving) reference position was provided,
2853 * therefore we can overwrite refproj[0]*/
2854 copyEvecReference(&edi->flood.vecs);
2857 else /* No origin structure given */
2859 rad_project(edi, xstart, &edi->vecs.radacc);
2860 rad_project(edi, xstart, &edi->vecs.radfix);
2861 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2863 if (edi->flood.bHarmonic)
2865 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2866 for (i = 0; i < edi->flood.vecs.neig; i++)
2868 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2873 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2874 /* Set center of flooding potential to the center of the covariance matrix,
2875 * i.e. the average structure, i.e. zero in the projected system */
2876 for (i = 0; i < edi->flood.vecs.neig; i++)
2878 edi->flood.vecs.refproj[i] = 0.0;
2883 /* For convenience, output the center of the flooding potential for the eigenvectors */
2884 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2886 for (i = 0; i < edi->flood.vecs.neig; i++)
2888 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2889 if (edi->flood.bHarmonic)
2891 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2893 fprintf(stdout, "\n");
2897 /* set starting projections for linsam */
2898 rad_project(edi, xstart, &edi->vecs.linacc);
2899 rad_project(edi, xstart, &edi->vecs.linfix);
2901 /* Prepare for the next edi data set: */
2902 edi = edi->next_edi;
2904 /* Cleaning up on the master node: */
2909 } /* end of MASTER only section */
2913 /* First let everybody know how many ED data sets to expect */
2914 gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr);
2915 /* Broadcast the essential dynamics / flooding data to all nodes */
2916 broadcast_ed_data(cr, ed, EDstate->nED);
2920 /* In the single-CPU case, point the local atom numbers pointers to the global
2921 * one, so that we can use the same notation in serial and parallel case: */
2923 /* Loop over all ED data sets (usually only one, though) */
2925 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2927 edi->sref.anrs_loc = edi->sref.anrs;
2928 edi->sav.anrs_loc = edi->sav.anrs;
2929 edi->star.anrs_loc = edi->star.anrs;
2930 edi->sori.anrs_loc = edi->sori.anrs;
2931 /* For the same reason as above, make a dummy c_ind array: */
2932 snew(edi->sav.c_ind, edi->sav.nr);
2933 /* Initialize the array */
2934 for (i = 0; i < edi->sav.nr; i++)
2936 edi->sav.c_ind[i] = i;
2938 /* In the general case we will need a different-sized array for the reference indices: */
2941 snew(edi->sref.c_ind, edi->sref.nr);
2942 for (i = 0; i < edi->sref.nr; i++)
2944 edi->sref.c_ind[i] = i;
2947 /* Point to the very same array in case of other structures: */
2948 edi->star.c_ind = edi->sav.c_ind;
2949 edi->sori.c_ind = edi->sav.c_ind;
2950 /* In the serial case, the local number of atoms is the global one: */
2951 edi->sref.nr_loc = edi->sref.nr;
2952 edi->sav.nr_loc = edi->sav.nr;
2953 edi->star.nr_loc = edi->star.nr;
2954 edi->sori.nr_loc = edi->sori.nr;
2956 /* An on we go to the next ED group */
2957 edi = edi->next_edi;
2961 /* Allocate space for ED buffer variables */
2962 /* Again, loop over ED data sets */
2964 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2966 /* Allocate space for ED buffer */
2968 snew(edi->buf->do_edsam, 1);
2970 /* Space for collective ED buffer variables */
2972 /* Collective positions of atoms with the average indices */
2973 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2974 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
2975 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2976 /* Collective positions of atoms with the reference indices */
2979 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2980 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
2981 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2984 /* Get memory for flooding forces */
2985 snew(edi->flood.forces_cartesian, edi->sav.nr);
2988 /* Dump it all into one file per process */
2989 dump_edi(edi, cr, nr_edi);
2993 edi = edi->next_edi;
2996 /* Flush the edo file so that the user can check some things
2997 * when the simulation has started */
3003 GMX_MPE_LOG(ev_edsam_finish);
3007 void do_edsam(t_inputrec *ir,
3008 gmx_large_int_t step,
3010 rvec xs[], /* The local current positions on this processor */
3011 rvec v[], /* The velocities */
3015 int i, edinr, iupdate = 500;
3016 matrix rotmat; /* rotation matrix */
3017 rvec transvec; /* translation vector */
3018 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3019 real dt_1; /* 1/dt */
3020 struct t_do_edsam *buf;
3022 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3023 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3026 /* Check if ED sampling has to be performed */
3027 if (ed->eEDtype == eEDnone)
3032 /* Suppress output on first call of do_edsam if
3033 * two-step sd2 integrator is used */
3034 if ( (ir->eI == eiSD2) && (v != NULL) )
3039 dt_1 = 1.0/ir->delta_t;
3041 /* Loop over all ED groups (usually one) */
3047 if (edi->bNeedDoEdsam)
3050 buf = edi->buf->do_edsam;
3054 /* initialise radacc radius for slope criterion */
3055 buf->oldrad = calc_radius(&edi->vecs.radacc);
3058 /* Copy the positions into buf->xc* arrays and after ED
3059 * feed back corrections to the official positions */
3061 /* Broadcast the ED positions such that every node has all of them
3062 * Every node contributes its local positions xs and stores it in
3063 * the collective buf->xcoll array. Note that for edinr > 1
3064 * xs could already have been modified by an earlier ED */
3066 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3067 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
3069 /* Only assembly reference positions if their indices differ from the average ones */
3072 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3073 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3076 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3077 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3078 * set bUpdateShifts=TRUE in the parallel case. */
3079 buf->bUpdateShifts = FALSE;
3081 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3082 * as well as the indices in edi->sav.anrs */
3084 /* Fit the reference indices to the reference structure */
3087 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3091 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3094 /* Now apply the translation and rotation to the ED structure */
3095 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3097 /* Find out how well we fit to the reference (just for output steps) */
3098 if (do_per_step(step, edi->outfrq) && MASTER(cr))
3102 /* Indices of reference and average structures are identical,
3103 * thus we can calculate the rmsd to SREF using xcoll */
3104 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3108 /* We have to translate & rotate the reference atoms first */
3109 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3110 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3114 /* update radsam references, when required */
3115 if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3117 project(buf->xcoll, edi);
3118 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3119 rad_project(edi, buf->xcoll, &edi->vecs.radfix);
3120 buf->oldrad = -1.e5;
3123 /* update radacc references, when required */
3124 if (do_per_step(step, iupdate) && step >= edi->presteps)
3126 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
3127 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3129 project(buf->xcoll, edi);
3130 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3135 buf->oldrad = edi->vecs.radacc.radius;
3139 /* apply the constraints */
3140 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3142 /* ED constraints should be applied already in the first MD step
3143 * (which is step 0), therefore we pass step+1 to the routine */
3144 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3147 /* write to edo, when required */
3148 if (do_per_step(step, edi->outfrq))
3150 project(buf->xcoll, edi);
3151 if (MASTER(cr) && !bSuppress)
3153 write_edo(edi, ed->edo, rmsdev);
3157 /* Copy back the positions unless monitoring only */
3158 if (ed_constraints(ed->eEDtype, edi))
3160 /* remove fitting */
3161 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3163 /* Copy the ED corrected positions into the coordinate array */
3164 /* Each node copies its local part. In the serial case, nat_loc is the
3165 * total number of ED atoms */
3166 for (i = 0; i < edi->sav.nr_loc; i++)
3168 /* Unshift local ED coordinate and store in x_unsh */
3169 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3170 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3172 /* dx is the ED correction to the positions: */
3173 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3177 /* dv is the ED correction to the velocity: */
3178 svmul(dt_1, dx, dv);
3179 /* apply the velocity correction: */
3180 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3182 /* Finally apply the position correction due to ED: */
3183 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3186 } /* END of if (edi->bNeedDoEdsam) */
3188 /* Prepare for the next ED group */
3189 edi = edi->next_edi;
3191 } /* END of loop over ED groups */
3195 GMX_MPE_LOG(ev_edsam_finish);