1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
56 #include "mtop_util.h"
60 #include "groupcoord.h"
63 /* We use the same defines as in mvdata.c here */
64 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
65 #define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
66 #define snew_bc(cr, d, nr) { if (!MASTER(cr)) {snew((d), (nr)); }}
68 /* These macros determine the column width in the output file */
69 #define EDcol_sfmt "%17s"
70 #define EDcol_efmt "%17.5e"
71 #define EDcol_ffmt "%17f"
73 /* enum to identify the type of ED: none, normal ED, flooding */
75 eEDnone, eEDedsam, eEDflood, eEDnr
78 /* enum to identify operations on reference, average, origin, target structures */
80 eedREF, eedAV, eedORI, eedTAR, eedNR
86 int neig; /* nr of eigenvectors */
87 int *ieig; /* index nrs of eigenvectors */
88 real *stpsz; /* stepsizes (per eigenvector) */
89 rvec **vec; /* eigenvector components */
90 real *xproj; /* instantaneous x projections */
91 real *fproj; /* instantaneous f projections */
92 real radius; /* instantaneous radius */
93 real *refproj; /* starting or target projecions */
94 /* When using flooding as harmonic restraint: The current reference projection
95 * is at each step calculated from the initial refproj0 and the slope. */
96 real *refproj0, *refprojslope;
102 t_eigvec mon; /* only monitored, no constraints */
103 t_eigvec linfix; /* fixed linear constraints */
104 t_eigvec linacc; /* acceptance linear constraints */
105 t_eigvec radfix; /* fixed radial constraints (exp) */
106 t_eigvec radacc; /* acceptance radial constraints (exp) */
107 t_eigvec radcon; /* acceptance rad. contraction constr. */
114 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
116 gmx_bool bConstForce; /* Do not calculate a flooding potential,
117 instead flood with a constant force */
126 rvec *forces_cartesian;
127 t_eigvec vecs; /* use flooding for these */
131 /* This type is for the average, reference, target, and origin structure */
132 typedef struct gmx_edx
134 int nr; /* number of atoms this structure contains */
135 int nr_loc; /* number of atoms on local node */
136 int *anrs; /* atom index numbers */
137 int *anrs_loc; /* local atom index numbers */
138 int nalloc_loc; /* allocation size of anrs_loc */
139 int *c_ind; /* at which position of the whole anrs
140 * array is a local atom?, i.e.
141 * c_ind[0...nr_loc-1] gives the atom index
142 * with respect to the collective
143 * anrs[0...nr-1] array */
144 rvec *x; /* positions for this structure */
145 rvec *x_old; /* Last positions which have the correct PBC
146 representation of the ED group. In
147 combination with keeping track of the
148 shift vectors, the ED group can always
150 real *m; /* masses */
151 real mtot; /* total mass (only used in sref) */
152 real *sqrtm; /* sqrt of the masses used for mass-
153 * weighting of analysis (only used in sav) */
159 int nini; /* total Nr of atoms */
160 gmx_bool fitmas; /* true if trans fit with cm */
161 gmx_bool pcamas; /* true if mass-weighted PCA */
162 int presteps; /* number of steps to run without any
163 * perturbations ... just monitoring */
164 int outfrq; /* freq (in steps) of writing to edo */
165 int maxedsteps; /* max nr of steps per cycle */
167 /* all gmx_edx datasets are copied to all nodes in the parallel case */
168 struct gmx_edx sref; /* reference positions, to these fitting
170 gmx_bool bRefEqAv; /* If true, reference & average indices
171 * are the same. Used for optimization */
172 struct gmx_edx sav; /* average positions */
173 struct gmx_edx star; /* target positions */
174 struct gmx_edx sori; /* origin positions */
176 t_edvecs vecs; /* eigenvectors */
177 real slope; /* minimal slope in acceptance radexp */
179 gmx_bool bNeedDoEdsam; /* if any of the options mon, linfix, ...
180 * is used (i.e. apart from flooding) */
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 /* Multiple ED groups will be labeled with letters instead of numbers
238 * to avoid confusion with eigenvector indices */
239 static char get_EDgroupChar(int nr_edi, int nED)
247 * nr_edi = 2 -> B ...
249 return 'A' + nr_edi - 1;
253 /* Does not subtract average positions, projection on single eigenvector is returned
254 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
255 * Average position is subtracted in ed_apply_constraints prior to calling projectx
257 static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec)
263 for (i = 0; i < edi->sav.nr; i++)
265 proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
272 /* Specialized: projection is stored in vec->refproj
273 * -> used for radacc, radfix, radcon and center of flooding potential
274 * subtracts average positions, projects vector x */
275 static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec)
280 /* Subtract average positions */
281 for (i = 0; i < edi->sav.nr; i++)
283 rvec_dec(x[i], edi->sav.x[i]);
286 for (i = 0; i < vec->neig; i++)
288 vec->refproj[i] = projectx(edi, x, vec->vec[i]);
289 rad += pow((vec->refproj[i]-vec->xproj[i]), 2);
291 vec->radius = sqrt(rad);
293 /* Add average positions */
294 for (i = 0; i < edi->sav.nr; i++)
296 rvec_inc(x[i], edi->sav.x[i]);
301 /* Project vector x, subtract average positions prior to projection and add
302 * them afterwards to retain the unchanged vector. Store in xproj. Mass-weighting
304 static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
305 t_eigvec *vec, /* The eigenvectors */
316 /* Subtract average positions */
317 for (i = 0; i < edi->sav.nr; i++)
319 rvec_dec(x[i], edi->sav.x[i]);
322 for (i = 0; i < vec->neig; i++)
324 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
327 /* Add average positions */
328 for (i = 0; i < edi->sav.nr; i++)
330 rvec_inc(x[i], edi->sav.x[i]);
335 /* Project vector x onto all edi->vecs (mon, linfix,...) */
336 static void project(rvec *x, /* positions to project */
337 t_edpar *edi) /* edi data set */
339 /* It is not more work to subtract the average position in every
340 * subroutine again, because these routines are rarely used simultanely */
341 project_to_eigvectors(x, &edi->vecs.mon, edi);
342 project_to_eigvectors(x, &edi->vecs.linfix, edi);
343 project_to_eigvectors(x, &edi->vecs.linacc, edi);
344 project_to_eigvectors(x, &edi->vecs.radfix, edi);
345 project_to_eigvectors(x, &edi->vecs.radacc, edi);
346 project_to_eigvectors(x, &edi->vecs.radcon, edi);
350 static real calc_radius(t_eigvec *vec)
356 for (i = 0; i < vec->neig; i++)
358 rad += pow((vec->refproj[i]-vec->xproj[i]), 2);
361 return rad = sqrt(rad);
367 static void dump_xcoll(t_edpar *edi, struct t_do_edsam *buf, t_commrec *cr,
374 ivec *shifts, *eshifts;
383 shifts = buf->shifts_xcoll;
384 eshifts = buf->extra_shifts_xcoll;
386 sprintf(fn, "xcolldump_step%d.txt", step);
389 for (i = 0; i < edi->sav.nr; i++)
391 fprintf(fp, "%d %9.5f %9.5f %9.5f %d %d %d %d %d %d\n",
393 xcoll[i][XX], xcoll[i][YY], xcoll[i][ZZ],
394 shifts[i][XX], shifts[i][YY], shifts[i][ZZ],
395 eshifts[i][XX], eshifts[i][YY], eshifts[i][ZZ]);
403 static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
408 fprintf(out, "#%s positions:\n%d\n", name, s->nr);
414 fprintf(out, "#index, x, y, z");
417 fprintf(out, ", sqrt(m)");
419 for (i = 0; i < s->nr; i++)
421 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]);
424 fprintf(out, "%9.3f", s->sqrtm[i]);
432 static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
433 const char name[], int length)
438 fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig);
439 /* Dump the data for every eigenvector: */
440 for (i = 0; i < ev->neig; i++)
442 fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
443 ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
444 for (j = 0; j < length; j++)
446 fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
453 static void dump_edi(t_edpar *edpars, t_commrec *cr, int nr_edi)
459 sprintf(fn, "EDdump_node%d_edi%d", cr->nodeid, nr_edi);
460 out = ffopen(fn, "w");
462 fprintf(out, "#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
463 edpars->nini, edpars->fitmas, edpars->pcamas);
464 fprintf(out, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
465 edpars->outfrq, edpars->maxedsteps, edpars->slope);
466 fprintf(out, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n",
467 edpars->presteps, edpars->flood.deltaF0, edpars->flood.tau,
468 edpars->flood.constEfl, edpars->flood.alpha2);
470 /* Dump reference, average, target, origin positions */
471 dump_edi_positions(out, &edpars->sref, "REFERENCE");
472 dump_edi_positions(out, &edpars->sav, "AVERAGE" );
473 dump_edi_positions(out, &edpars->star, "TARGET" );
474 dump_edi_positions(out, &edpars->sori, "ORIGIN" );
476 /* Dump eigenvectors */
477 dump_edi_eigenvecs(out, &edpars->vecs.mon, "MONITORED", edpars->sav.nr);
478 dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX", edpars->sav.nr);
479 dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC", edpars->sav.nr);
480 dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX", edpars->sav.nr);
481 dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC", edpars->sav.nr);
482 dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON", edpars->sav.nr);
484 /* Dump flooding eigenvectors */
485 dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING", edpars->sav.nr);
487 /* Dump ed local buffer */
488 fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit );
489 fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam );
490 fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon );
497 static void dump_rotmat(FILE* out, matrix rotmat)
499 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[XX][XX], rotmat[XX][YY], rotmat[XX][ZZ]);
500 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[YY][XX], rotmat[YY][YY], rotmat[YY][ZZ]);
501 fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[ZZ][XX], rotmat[ZZ][YY], rotmat[ZZ][ZZ]);
506 static void dump_rvec(FILE *out, int dim, rvec *x)
511 for (i = 0; i < dim; i++)
513 fprintf(out, "%4d %f %f %f\n", i, x[i][XX], x[i][YY], x[i][ZZ]);
519 static void dump_mat(FILE* out, int dim, double** mat)
524 fprintf(out, "MATRIX:\n");
525 for (i = 0; i < dim; i++)
527 for (j = 0; j < dim; j++)
529 fprintf(out, "%f ", mat[i][j]);
542 static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi)
544 /* this is a copy of do_fit with some modifications */
545 int c, r, n, j, i, irot;
546 double d[6], xnr, xpc;
551 struct t_do_edfit *loc;
554 if (edi->buf->do_edfit != NULL)
561 snew(edi->buf->do_edfit, 1);
563 loc = edi->buf->do_edfit;
567 snew(loc->omega, 2*DIM);
568 snew(loc->om, 2*DIM);
569 for (i = 0; i < 2*DIM; i++)
571 snew(loc->omega[i], 2*DIM);
572 snew(loc->om[i], 2*DIM);
576 for (i = 0; (i < 6); i++)
579 for (j = 0; (j < 6); j++)
581 loc->omega[i][j] = 0;
586 /* calculate the matrix U */
588 for (n = 0; (n < natoms); n++)
590 for (c = 0; (c < DIM); c++)
593 for (r = 0; (r < DIM); r++)
601 /* construct loc->omega */
602 /* loc->omega is symmetric -> loc->omega==loc->omega' */
603 for (r = 0; (r < 6); r++)
605 for (c = 0; (c <= r); c++)
607 if ((r >= 3) && (c < 3))
609 loc->omega[r][c] = u[r-3][c];
610 loc->omega[c][r] = u[r-3][c];
614 loc->omega[r][c] = 0;
615 loc->omega[c][r] = 0;
620 /* determine h and k */
624 dump_mat(stderr, 2*DIM, loc->omega);
625 for (i = 0; i < 6; i++)
627 fprintf(stderr, "d[%d] = %f\n", i, d[i]);
631 jacobi(loc->omega, 6, d, loc->om, &irot);
635 fprintf(stderr, "IROT=0\n");
638 index = 0; /* For the compiler only */
640 for (j = 0; (j < 3); j++)
643 for (i = 0; (i < 6); i++)
652 for (i = 0; (i < 3); i++)
654 vh[j][i] = M_SQRT2*loc->om[i][index];
655 vk[j][i] = M_SQRT2*loc->om[i+DIM][index];
660 for (c = 0; (c < 3); c++)
662 for (r = 0; (r < 3); r++)
664 R[c][r] = vk[0][r]*vh[0][c]+
671 for (c = 0; (c < 3); c++)
673 for (r = 0; (r < 3); r++)
675 R[c][r] = vk[0][r]*vh[0][c]+
684 static void rmfit(int nat, rvec *xcoll, rvec transvec, matrix rotmat)
691 * The inverse rotation is described by the transposed rotation matrix */
692 transpose(rotmat, tmat);
693 rotate_x(xcoll, nat, tmat);
695 /* Remove translation */
696 vec[XX] = -transvec[XX];
697 vec[YY] = -transvec[YY];
698 vec[ZZ] = -transvec[ZZ];
699 translate_x(xcoll, nat, vec);
703 /**********************************************************************************
704 ******************** FLOODING ****************************************************
705 **********************************************************************************
707 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
708 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
709 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
711 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
712 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
714 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
715 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
717 do_flood makes a copy of the positions,
718 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
719 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
720 fit. Then do_flood adds these forces to the forcefield-forces
721 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
723 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
724 structure is projected to the system of eigenvectors and then this position in the subspace is used as
725 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
726 i.e. the average structure as given in the make_edi file.
728 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
729 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
730 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
731 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
732 further adaption is applied, Efl will stay constant at zero.
734 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
735 used as spring constants for the harmonic potential.
736 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
737 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
739 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
740 the routine read_edi_file reads all of theses flooding files.
741 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
742 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
743 edi there is no interdependence whatsoever. The forces are added together.
745 To write energies into the .edr file, call the function
746 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
748 get_flood_energies(real Vfl[],int nnames);
751 - 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.
753 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
754 two edsam files from two peptide chains
757 static void write_edo_flood(t_edpar *edi, FILE *fp, real rmsd)
762 /* Output how well we fit to the reference structure */
763 fprintf(fp, EDcol_ffmt, rmsd);
765 for (i = 0; i < edi->flood.vecs.neig; i++)
767 fprintf(fp, EDcol_efmt, edi->flood.vecs.xproj[i]);
769 /* Check whether the reference projection changes with time (this can happen
770 * in case flooding is used as harmonic restraint). If so, output the
771 * current reference projection */
772 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
774 fprintf(fp, EDcol_efmt, edi->flood.vecs.refproj[i]);
777 /* Output Efl if we are doing adaptive flooding */
778 if (0 != edi->flood.tau)
780 fprintf(fp, EDcol_efmt, edi->flood.Efl);
782 fprintf(fp, EDcol_efmt, edi->flood.Vfl);
784 /* Output deltaF if we are doing adaptive flooding */
785 if (0 != edi->flood.tau)
787 fprintf(fp, EDcol_efmt, edi->flood.deltaF);
789 fprintf(fp, EDcol_efmt, edi->flood.vecs.fproj[i]);
794 /* From flood.xproj compute the Vfl(x) at this point */
795 static real flood_energy(t_edpar *edi, gmx_large_int_t step)
797 /* compute flooding energy Vfl
798 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
799 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
800 it is already computed by make_edi and stored in stpsz[i]
802 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
809 /* Each time this routine is called (i.e. each time step), we add a small
810 * value to the reference projection. This way a harmonic restraint towards
811 * a moving reference is realized. If no value for the additive constant
812 * is provided in the edi file, the reference will not change. */
813 if (edi->flood.bHarmonic)
815 for (i = 0; i < edi->flood.vecs.neig; i++)
817 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
822 /* Compute sum which will be the exponent of the exponential */
823 for (i = 0; i < edi->flood.vecs.neig; i++)
825 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
826 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]);
829 /* Compute the Gauss function*/
830 if (edi->flood.bHarmonic)
832 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
836 Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0;
843 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
844 static void flood_forces(t_edpar *edi)
846 /* compute the forces in the subspace of the flooding eigenvectors
847 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
850 real energy = edi->flood.Vfl;
853 if (edi->flood.bHarmonic)
855 for (i = 0; i < edi->flood.vecs.neig; i++)
857 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
862 for (i = 0; i < edi->flood.vecs.neig; i++)
864 /* if Efl is zero the forces are zero if not use the formula */
865 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;
871 /* Raise forces from subspace into cartesian space */
872 static void flood_blowup(t_edpar *edi, rvec *forces_cart)
874 /* this function lifts the forces from the subspace to the cartesian space
875 all the values not contained in the subspace are assumed to be zero and then
876 a coordinate transformation from eigenvector to cartesian vectors is performed
877 The nonexistent values don't have to be set to zero explicitly, they would occur
878 as zero valued summands, hence we just stop to compute this part of the sum.
880 for every atom we add all the contributions to this atom from all the different eigenvectors.
882 NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
883 field forces_cart prior the computation, but we compute the forces separately
884 to have them accessible for diagnostics
891 forces_sub = edi->flood.vecs.fproj;
894 /* Calculate the cartesian forces for the local atoms */
896 /* Clear forces first */
897 for (j = 0; j < edi->sav.nr_loc; j++)
899 clear_rvec(forces_cart[j]);
902 /* Now compute atomwise */
903 for (j = 0; j < edi->sav.nr_loc; j++)
905 /* Compute forces_cart[edi->sav.anrs[j]] */
906 for (eig = 0; eig < edi->flood.vecs.neig; eig++)
908 /* Force vector is force * eigenvector (compute only atom j) */
909 svmul(forces_sub[eig], edi->flood.vecs.vec[eig][edi->sav.c_ind[j]], dum);
910 /* Add this vector to the cartesian forces */
911 rvec_inc(forces_cart[j], dum);
917 /* Update the values of Efl, deltaF depending on tau and Vfl */
918 static void update_adaption(t_edpar *edi)
920 /* this function updates the parameter Efl and deltaF according to the rules given in
921 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
924 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
926 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
927 /* check if restrain (inverted flooding) -> don't let EFL become positive */
928 if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
933 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
938 static void do_single_flood(
943 gmx_large_int_t step,
946 gmx_bool bNS) /* Are we in a neighbor searching step? */
949 matrix rotmat; /* rotation matrix */
950 matrix tmat; /* inverse rotation */
951 rvec transvec; /* translation vector */
953 struct t_do_edsam *buf;
956 buf = edi->buf->do_edsam;
959 /* Broadcast the positions of the AVERAGE structure such that they are known on
960 * every processor. Each node contributes its local positions x and stores them in
961 * the collective ED array buf->xcoll */
962 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
963 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
965 /* Only assembly REFERENCE positions if their indices differ from the average ones */
968 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
969 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
972 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
973 * We do not need to update the shifts until the next NS step */
974 buf->bUpdateShifts = FALSE;
976 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
977 * as well as the indices in edi->sav.anrs */
979 /* Fit the reference indices to the reference structure */
982 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
986 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
989 /* Now apply the translation and rotation to the ED structure */
990 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
992 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
993 project_to_eigvectors(buf->xcoll, &edi->flood.vecs, edi);
995 if (FALSE == edi->flood.bConstForce)
997 /* Compute Vfl(x) from flood.xproj */
998 edi->flood.Vfl = flood_energy(edi, step);
1000 update_adaption(edi);
1002 /* Compute the flooding forces */
1006 /* Translate them into cartesian positions */
1007 flood_blowup(edi, edi->flood.forces_cartesian);
1009 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
1010 /* Each node rotates back its local forces */
1011 transpose(rotmat, tmat);
1012 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
1014 /* Finally add forces to the main force variable */
1015 for (i = 0; i < edi->sav.nr_loc; i++)
1017 rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
1020 /* Output is written by the master process */
1021 if (do_per_step(step, edi->outfrq) && MASTER(cr))
1023 /* Output how well we fit to the reference */
1026 /* Indices of reference and average structures are identical,
1027 * thus we can calculate the rmsd to SREF using xcoll */
1028 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
1032 /* We have to translate & rotate the reference atoms first */
1033 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1034 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1037 write_edo_flood(edi, edo, rmsdev);
1042 /* Main flooding routine, called from do_force */
1043 extern void do_flood(
1044 t_commrec *cr, /* Communication record */
1045 t_inputrec *ir, /* Input record */
1046 rvec x[], /* Positions on the local processor */
1047 rvec force[], /* forcefield forces, to these the flooding forces are added */
1048 gmx_edsam_t ed, /* ed data structure contains all ED and flooding groups */
1049 matrix box, /* the box */
1050 gmx_large_int_t step, /* The relative time step since ir->init_step is already subtracted */
1051 gmx_bool bNS) /* Are we in a neighbor searching step? */
1058 /* Write time to edo, when required. Output the time anyhow since we need
1059 * it in the output file for ED constraints. */
1060 if (MASTER(cr) && do_per_step(step, edi->outfrq))
1062 fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
1065 if (ed->eEDtype != eEDflood)
1072 /* Call flooding for one matrix */
1073 if (edi->flood.vecs.neig)
1075 do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS);
1077 edi = edi->next_edi;
1082 /* Called by init_edi, configure some flooding related variables and structures,
1083 * print headers to output files */
1084 static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt)
1089 edi->flood.Efl = edi->flood.constEfl;
1093 if (edi->flood.vecs.neig)
1095 /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1096 ed->eEDtype = eEDflood;
1098 fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : "");
1100 if (edi->flood.bConstForce)
1102 /* We have used stpsz as a vehicle to carry the fproj values for constant
1103 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1104 * in const force flooding, fproj is never changed. */
1105 for (i = 0; i < edi->flood.vecs.neig; i++)
1107 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1109 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1110 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1118 /*********** Energy book keeping ******/
1119 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1128 srenew(names, count);
1129 sprintf(buf, "Vfl_%d", count);
1130 names[count-1] = strdup(buf);
1131 actual = actual->next_edi;
1138 static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames)
1140 /*fl has to be big enough to capture nnames-many entries*/
1149 Vfl[count-1] = actual->flood.Vfl;
1150 actual = actual->next_edi;
1153 if (nnames != count-1)
1155 gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1158 /************* END of FLOODING IMPLEMENTATION ****************************/
1162 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)
1168 /* Allocate space for the ED data structure */
1171 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1172 ed->eEDtype = eEDedsam;
1176 fprintf(stderr, "ED sampling will be performed!\n");
1179 /* Read the edi input file: */
1180 nED = read_edi_file(ftp2fn(efEDI, nfile, fnm), ed->edpar, natoms);
1182 /* Make sure the checkpoint was produced in a run using this .edi file */
1183 if (EDstate->bFromCpt)
1185 crosscheck_edi_file_vs_checkpoint(ed, EDstate);
1191 init_edsamstate(ed, EDstate);
1193 /* The master opens the ED output file */
1194 if (Flags & MD_APPENDFILES)
1196 ed->edo = gmx_fio_fopen(opt2fn("-eo", nfile, fnm), "a+");
1200 ed->edo = xvgropen(opt2fn("-eo", nfile, fnm),
1201 "Essential dynamics / flooding output",
1203 "RMSDs (nm), projections on EVs (nm), ...", oenv);
1205 /* Make a descriptive legend */
1206 write_edo_legend(ed, EDstate->nED, oenv);
1213 /* Broadcasts the structure data */
1214 static void bc_ed_positions(t_commrec *cr, struct gmx_edx *s, int stype)
1216 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1217 snew_bc(cr, s->x, s->nr ); /* Positions */
1218 nblock_bc(cr, s->nr, s->anrs );
1219 nblock_bc(cr, s->nr, s->x );
1221 /* For the average & reference structures we need an array for the collective indices,
1222 * and we need to broadcast the masses as well */
1223 if (stype == eedAV || stype == eedREF)
1225 /* We need these additional variables in the parallel case: */
1226 snew(s->c_ind, s->nr ); /* Collective indices */
1227 /* Local atom indices get assigned in dd_make_local_group_indices.
1228 * There, also memory is allocated */
1229 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1230 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1231 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1234 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1235 if (stype == eedREF)
1237 snew_bc(cr, s->m, s->nr);
1238 nblock_bc(cr, s->nr, s->m);
1241 /* For the average structure we might need the masses for mass-weighting */
1244 snew_bc(cr, s->sqrtm, s->nr);
1245 nblock_bc(cr, s->nr, s->sqrtm);
1246 snew_bc(cr, s->m, s->nr);
1247 nblock_bc(cr, s->nr, s->m);
1252 /* Broadcasts the eigenvector data */
1253 static void bc_ed_vecs(t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
1257 snew_bc(cr, ev->ieig, ev->neig); /* index numbers of eigenvector */
1258 snew_bc(cr, ev->stpsz, ev->neig); /* stepsizes per eigenvector */
1259 snew_bc(cr, ev->xproj, ev->neig); /* instantaneous x projection */
1260 snew_bc(cr, ev->fproj, ev->neig); /* instantaneous f projection */
1261 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1263 nblock_bc(cr, ev->neig, ev->ieig );
1264 nblock_bc(cr, ev->neig, ev->stpsz );
1265 nblock_bc(cr, ev->neig, ev->xproj );
1266 nblock_bc(cr, ev->neig, ev->fproj );
1267 nblock_bc(cr, ev->neig, ev->refproj);
1269 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1270 for (i = 0; i < ev->neig; i++)
1272 snew_bc(cr, ev->vec[i], length);
1273 nblock_bc(cr, length, ev->vec[i]);
1276 /* For harmonic restraints the reference projections can change with time */
1279 snew_bc(cr, ev->refproj0, ev->neig);
1280 snew_bc(cr, ev->refprojslope, ev->neig);
1281 nblock_bc(cr, ev->neig, ev->refproj0 );
1282 nblock_bc(cr, ev->neig, ev->refprojslope);
1287 /* Broadcasts the ED / flooding data to other nodes
1288 * and allocates memory where needed */
1289 static void broadcast_ed_data(t_commrec *cr, gmx_edsam_t ed, int numedis)
1295 /* Master lets the other nodes know if its ED only or also flooding */
1296 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1298 snew_bc(cr, ed->edpar, 1);
1299 /* Now transfer the ED data set(s) */
1301 for (nr = 0; nr < numedis; nr++)
1303 /* Broadcast a single ED data set */
1306 /* Broadcast positions */
1307 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1308 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1309 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1310 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1312 /* Broadcast eigenvectors */
1313 bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr, FALSE);
1314 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
1315 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
1316 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
1317 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
1318 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
1319 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1320 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
1322 /* Set the pointer to the next ED group */
1325 snew_bc(cr, edi->next_edi, 1);
1326 edi = edi->next_edi;
1332 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1333 static void init_edi(gmx_mtop_t *mtop, t_edpar *edi)
1336 real totalmass = 0.0;
1338 gmx_mtop_atomlookup_t alook = NULL;
1341 /* NOTE Init_edi is executed on the master process only
1342 * The initialized data sets are then transmitted to the
1343 * other nodes in broadcast_ed_data */
1345 alook = gmx_mtop_atomlookup_init(mtop);
1347 /* evaluate masses (reference structure) */
1348 snew(edi->sref.m, edi->sref.nr);
1349 for (i = 0; i < edi->sref.nr; i++)
1353 gmx_mtop_atomnr_to_atom(alook, edi->sref.anrs[i], &atom);
1354 edi->sref.m[i] = atom->m;
1358 edi->sref.m[i] = 1.0;
1361 /* Check that every m > 0. Bad things will happen otherwise. */
1362 if (edi->sref.m[i] <= 0.0)
1364 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1365 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1366 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1367 "atoms from the reference structure by creating a proper index group.\n",
1368 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1371 totalmass += edi->sref.m[i];
1373 edi->sref.mtot = totalmass;
1375 /* Masses m and sqrt(m) for the average structure. Note that m
1376 * is needed if forces have to be evaluated in do_edsam */
1377 snew(edi->sav.sqrtm, edi->sav.nr );
1378 snew(edi->sav.m, edi->sav.nr );
1379 for (i = 0; i < edi->sav.nr; i++)
1381 gmx_mtop_atomnr_to_atom(alook, edi->sav.anrs[i], &atom);
1382 edi->sav.m[i] = atom->m;
1385 edi->sav.sqrtm[i] = sqrt(atom->m);
1389 edi->sav.sqrtm[i] = 1.0;
1392 /* Check that every m > 0. Bad things will happen otherwise. */
1393 if (edi->sav.sqrtm[i] <= 0.0)
1395 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1396 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1397 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1398 "atoms from the average structure by creating a proper index group.\n",
1399 i, edi->sav.anrs[i]+1, atom->m);
1403 gmx_mtop_atomlookup_destroy(alook);
1405 /* put reference structure in origin */
1406 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1410 translate_x(edi->sref.x, edi->sref.nr, com);
1412 /* Init ED buffer */
1417 static void check(const char *line, const char *label)
1419 if (!strstr(line, label))
1421 gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1426 static int read_checked_edint(FILE *file, const char *label)
1428 char line[STRLEN+1];
1432 fgets2 (line, STRLEN, file);
1434 fgets2 (line, STRLEN, file);
1435 sscanf (line, "%d", &idum);
1440 static int read_edint(FILE *file, gmx_bool *bEOF)
1442 char line[STRLEN+1];
1447 eof = fgets2 (line, STRLEN, file);
1453 eof = fgets2 (line, STRLEN, file);
1459 sscanf (line, "%d", &idum);
1465 static real read_checked_edreal(FILE *file, const char *label)
1467 char line[STRLEN+1];
1471 fgets2 (line, STRLEN, file);
1473 fgets2 (line, STRLEN, file);
1474 sscanf (line, "%lf", &rdum);
1475 return (real) rdum; /* always read as double and convert to single */
1479 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1482 char line[STRLEN+1];
1486 for (i = 0; i < number; i++)
1488 fgets2 (line, STRLEN, file);
1489 sscanf (line, "%d%lf%lf%lf", &anrs[i], &d[0], &d[1], &d[2]);
1490 anrs[i]--; /* we are reading FORTRAN indices */
1491 for (j = 0; j < 3; j++)
1493 x[i][j] = d[j]; /* always read as double and convert to single */
1499 static void scan_edvec(FILE *in, int nr, rvec *vec)
1501 char line[STRLEN+1];
1506 for (i = 0; (i < nr); i++)
1508 fgets2 (line, STRLEN, in);
1509 sscanf (line, "%le%le%le", &x, &y, &z);
1517 static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1520 double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0;
1521 char line[STRLEN+1];
1524 tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1527 snew(tvec->ieig, tvec->neig);
1528 snew(tvec->stpsz, tvec->neig);
1529 snew(tvec->vec, tvec->neig);
1530 snew(tvec->xproj, tvec->neig);
1531 snew(tvec->fproj, tvec->neig);
1532 snew(tvec->refproj, tvec->neig);
1535 snew(tvec->refproj0, tvec->neig);
1536 snew(tvec->refprojslope, tvec->neig);
1539 for (i = 0; (i < tvec->neig); i++)
1541 fgets2 (line, STRLEN, in);
1542 if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1544 nscan = sscanf(line, "%d%lf%lf%lf", &idum, &rdum, &refproj_dum, &refprojslope_dum);
1545 /* Zero out values which were not scanned */
1549 /* Every 4 values read, including reference position */
1550 *bHaveReference = TRUE;
1553 /* A reference position is provided */
1554 *bHaveReference = TRUE;
1555 /* No value for slope, set to 0 */
1556 refprojslope_dum = 0.0;
1559 /* No values for reference projection and slope, set to 0 */
1561 refprojslope_dum = 0.0;
1564 gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1567 tvec->refproj[i] = refproj_dum;
1568 tvec->refproj0[i] = refproj_dum;
1569 tvec->refprojslope[i] = refprojslope_dum;
1571 else /* Normal flooding */
1573 nscan = sscanf(line, "%d%lf", &idum, &rdum);
1576 gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1579 tvec->ieig[i] = idum;
1580 tvec->stpsz[i] = rdum;
1581 } /* end of loop over eigenvectors */
1583 for (i = 0; (i < tvec->neig); i++)
1585 snew(tvec->vec[i], nr);
1586 scan_edvec(in, nr, tvec->vec[i]);
1592 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1593 static void read_edvecs(FILE *in, int nr, t_edvecs *vecs)
1595 gmx_bool bHaveReference = FALSE;
1598 read_edvec(in, nr, &vecs->mon, FALSE, &bHaveReference);
1599 read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1600 read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1601 read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1602 read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1603 read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1607 /* Check if the same atom indices are used for reference and average positions */
1608 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1613 /* If the number of atoms differs between the two structures,
1614 * they cannot be identical */
1615 if (sref.nr != sav.nr)
1620 /* Now that we know that both stuctures have the same number of atoms,
1621 * check if also the indices are identical */
1622 for (i = 0; i < sav.nr; i++)
1624 if (sref.anrs[i] != sav.anrs[i])
1629 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1635 static int read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, const char *fn)
1638 const int magic = 670;
1641 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1642 gmx_bool bHaveReference = FALSE;
1645 /* the edi file is not free format, so expect problems if the input is corrupt. */
1647 /* check the magic number */
1648 readmagic = read_edint(in, &bEOF);
1649 /* Check whether we have reached the end of the input file */
1655 if (readmagic != magic)
1657 if (readmagic == 666 || readmagic == 667 || readmagic == 668)
1659 gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1661 else if (readmagic != 669)
1663 gmx_fatal(FARGS, "Wrong magic number %d in %s", readmagic, fn);
1667 /* check the number of atoms */
1668 edi->nini = read_edint(in, &bEOF);
1669 if (edi->nini != nr_mdatoms)
1671 gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1674 /* Done checking. For the rest we blindly trust the input */
1675 edi->fitmas = read_checked_edint(in, "FITMAS");
1676 edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS");
1677 edi->outfrq = read_checked_edint(in, "OUTFRQ");
1678 edi->maxedsteps = read_checked_edint(in, "MAXLEN");
1679 edi->slope = read_checked_edreal(in, "SLOPECRIT");
1681 edi->presteps = read_checked_edint(in, "PRESTEPS");
1682 edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0");
1683 edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F");
1684 edi->flood.tau = read_checked_edreal(in, "TAU");
1685 edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL");
1686 edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2");
1687 edi->flood.kT = read_checked_edreal(in, "KT");
1688 edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1689 if (readmagic > 669)
1691 edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1695 edi->flood.bConstForce = FALSE;
1697 edi->sref.nr = read_checked_edint(in, "NREF");
1699 /* allocate space for reference positions and read them */
1700 snew(edi->sref.anrs, edi->sref.nr);
1701 snew(edi->sref.x, edi->sref.nr);
1702 snew(edi->sref.x_old, edi->sref.nr);
1703 edi->sref.sqrtm = NULL;
1704 read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1706 /* average positions. they define which atoms will be used for ED sampling */
1707 edi->sav.nr = read_checked_edint(in, "NAV");
1708 snew(edi->sav.anrs, edi->sav.nr);
1709 snew(edi->sav.x, edi->sav.nr);
1710 snew(edi->sav.x_old, edi->sav.nr);
1711 read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1713 /* Check if the same atom indices are used for reference and average positions */
1714 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1717 read_edvecs(in, edi->sav.nr, &edi->vecs);
1718 read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
1720 /* target positions */
1721 edi->star.nr = read_edint(in, &bEOF);
1722 if (edi->star.nr > 0)
1724 snew(edi->star.anrs, edi->star.nr);
1725 snew(edi->star.x, edi->star.nr);
1726 edi->star.sqrtm = NULL;
1727 read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1730 /* positions defining origin of expansion circle */
1731 edi->sori.nr = read_edint(in, &bEOF);
1732 if (edi->sori.nr > 0)
1736 /* Both an -ori structure and a at least one manual reference point have been
1737 * specified. That's ambiguous and probably not intentional. */
1738 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1739 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1741 snew(edi->sori.anrs, edi->sori.nr);
1742 snew(edi->sori.x, edi->sori.nr);
1743 edi->sori.sqrtm = NULL;
1744 read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1753 /* Read in the edi input file. Note that it may contain several ED data sets which were
1754 * achieved by concatenating multiple edi files. The standard case would be a single ED
1755 * data set, though. */
1756 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1759 t_edpar *curr_edi, *last_edi;
1764 /* This routine is executed on the master only */
1766 /* Open the .edi parameter input file */
1767 in = gmx_fio_fopen(fn, "r");
1768 fprintf(stderr, "ED: Reading edi file %s\n", fn);
1770 /* Now read a sequence of ED input parameter sets from the edi file */
1773 while (read_edi(in, curr_edi, nr_mdatoms, fn) )
1777 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1778 /* We need to allocate space for the data: */
1780 /* Point the 'next_edi' entry to the next edi: */
1781 curr_edi->next_edi = edi_read;
1782 /* Keep the curr_edi pointer for the case that the next group is empty: */
1783 last_edi = curr_edi;
1784 /* Let's prepare to read in the next edi data set: */
1785 curr_edi = edi_read;
1789 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1792 /* Terminate the edi group list with a NULL pointer: */
1793 last_edi->next_edi = NULL;
1795 fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1797 /* Close the .edi file again */
1804 struct t_fit_to_ref {
1805 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1808 /* Fit the current positions to the reference positions
1809 * Do not actually do the fit, just return rotation and translation.
1810 * Note that the COM of the reference structure was already put into
1811 * the origin by init_edi. */
1812 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1813 rvec transvec, /* The translation vector */
1814 matrix rotmat, /* The rotation matrix */
1815 t_edpar *edi) /* Just needed for do_edfit */
1817 rvec com; /* center of mass */
1819 struct t_fit_to_ref *loc;
1822 /* Allocate memory the first time this routine is called for each edi group */
1823 if (NULL == edi->buf->fit_to_ref)
1825 snew(edi->buf->fit_to_ref, 1);
1826 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1828 loc = edi->buf->fit_to_ref;
1830 /* We do not touch the original positions but work on a copy. */
1831 for (i = 0; i < edi->sref.nr; i++)
1833 copy_rvec(xcoll[i], loc->xcopy[i]);
1836 /* Calculate the center of mass */
1837 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1839 transvec[XX] = -com[XX];
1840 transvec[YY] = -com[YY];
1841 transvec[ZZ] = -com[ZZ];
1843 /* Subtract the center of mass from the copy */
1844 translate_x(loc->xcopy, edi->sref.nr, transvec);
1846 /* Determine the rotation matrix */
1847 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1851 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1852 int nat, /* How many positions are there? */
1853 rvec transvec, /* The translation vector */
1854 matrix rotmat) /* The rotation matrix */
1857 translate_x(x, nat, transvec);
1860 rotate_x(x, nat, rotmat);
1864 /* Gets the rms deviation of the positions to the structure s */
1865 /* fit_to_structure has to be called before calling this routine! */
1866 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1867 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1873 for (i = 0; i < s->nr; i++)
1875 rmsd += distance2(s->x[i], x[i]);
1878 rmsd /= (real) s->nr;
1885 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1890 if (ed->eEDtype != eEDnone)
1892 /* Loop over ED groups */
1896 /* Local atoms of the reference structure (for fitting), need only be assembled
1897 * if their indices differ from the average ones */
1900 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1901 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1904 /* Local atoms of the average structure (on these ED will be performed) */
1905 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1906 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1908 /* Indicate that the ED shift vectors for this structure need to be updated
1909 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1910 edi->buf->do_edsam->bUpdateShifts = TRUE;
1912 /* Set the pointer to the next ED group (if any) */
1913 edi = edi->next_edi;
1919 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1930 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1931 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1932 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1936 xu[XX] = x[XX]-tx*box[XX][XX];
1937 xu[YY] = x[YY]-ty*box[YY][YY];
1938 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1943 static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_large_int_t step)
1950 /* loop over linfix vectors */
1951 for (i = 0; i < edi->vecs.linfix.neig; i++)
1953 /* calculate the projection */
1954 proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1956 /* calculate the correction */
1957 add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1959 /* apply the correction */
1960 add /= edi->sav.sqrtm[i];
1961 for (j = 0; j < edi->sav.nr; j++)
1963 svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1964 rvec_inc(xcoll[j], vec_dum);
1970 static void do_linacc(rvec *xcoll, t_edpar *edi)
1977 /* loop over linacc vectors */
1978 for (i = 0; i < edi->vecs.linacc.neig; i++)
1980 /* calculate the projection */
1981 proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1983 /* calculate the correction */
1985 if (edi->vecs.linacc.stpsz[i] > 0.0)
1987 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1989 add = edi->vecs.linacc.refproj[i] - proj;
1992 if (edi->vecs.linacc.stpsz[i] < 0.0)
1994 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1996 add = edi->vecs.linacc.refproj[i] - proj;
2000 /* apply the correction */
2001 add /= edi->sav.sqrtm[i];
2002 for (j = 0; j < edi->sav.nr; j++)
2004 svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
2005 rvec_inc(xcoll[j], vec_dum);
2008 /* new positions will act as reference */
2009 edi->vecs.linacc.refproj[i] = proj + add;
2014 static void do_radfix(rvec *xcoll, t_edpar *edi)
2017 real *proj, rad = 0.0, ratio;
2021 if (edi->vecs.radfix.neig == 0)
2026 snew(proj, edi->vecs.radfix.neig);
2028 /* loop over radfix vectors */
2029 for (i = 0; i < edi->vecs.radfix.neig; i++)
2031 /* calculate the projections, radius */
2032 proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
2033 rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
2037 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2038 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2040 /* loop over radfix vectors */
2041 for (i = 0; i < edi->vecs.radfix.neig; i++)
2043 proj[i] -= edi->vecs.radfix.refproj[i];
2045 /* apply the correction */
2046 proj[i] /= edi->sav.sqrtm[i];
2048 for (j = 0; j < edi->sav.nr; j++)
2050 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2051 rvec_inc(xcoll[j], vec_dum);
2059 static void do_radacc(rvec *xcoll, t_edpar *edi)
2062 real *proj, rad = 0.0, ratio = 0.0;
2066 if (edi->vecs.radacc.neig == 0)
2071 snew(proj, edi->vecs.radacc.neig);
2073 /* loop over radacc vectors */
2074 for (i = 0; i < edi->vecs.radacc.neig; i++)
2076 /* calculate the projections, radius */
2077 proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2078 rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
2082 /* only correct when radius decreased */
2083 if (rad < edi->vecs.radacc.radius)
2085 ratio = edi->vecs.radacc.radius/rad - 1.0;
2086 rad = edi->vecs.radacc.radius;
2090 edi->vecs.radacc.radius = rad;
2093 /* loop over radacc vectors */
2094 for (i = 0; i < edi->vecs.radacc.neig; i++)
2096 proj[i] -= edi->vecs.radacc.refproj[i];
2098 /* apply the correction */
2099 proj[i] /= edi->sav.sqrtm[i];
2101 for (j = 0; j < edi->sav.nr; j++)
2103 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2104 rvec_inc(xcoll[j], vec_dum);
2111 struct t_do_radcon {
2115 static void do_radcon(rvec *xcoll, t_edpar *edi)
2118 real rad = 0.0, ratio = 0.0;
2119 struct t_do_radcon *loc;
2124 if (edi->buf->do_radcon != NULL)
2127 loc = edi->buf->do_radcon;
2132 snew(edi->buf->do_radcon, 1);
2134 loc = edi->buf->do_radcon;
2136 if (edi->vecs.radcon.neig == 0)
2143 snew(loc->proj, edi->vecs.radcon.neig);
2146 /* loop over radcon vectors */
2147 for (i = 0; i < edi->vecs.radcon.neig; i++)
2149 /* calculate the projections, radius */
2150 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2151 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2154 /* only correct when radius increased */
2155 if (rad > edi->vecs.radcon.radius)
2157 ratio = edi->vecs.radcon.radius/rad - 1.0;
2159 /* loop over radcon vectors */
2160 for (i = 0; i < edi->vecs.radcon.neig; i++)
2162 /* apply the correction */
2163 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2164 loc->proj[i] /= edi->sav.sqrtm[i];
2165 loc->proj[i] *= ratio;
2167 for (j = 0; j < edi->sav.nr; j++)
2169 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2170 rvec_inc(xcoll[j], vec_dum);
2176 edi->vecs.radcon.radius = rad;
2179 if (rad != edi->vecs.radcon.radius)
2182 for (i = 0; i < edi->vecs.radcon.neig; i++)
2184 /* calculate the projections, radius */
2185 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2186 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2193 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step)
2198 /* subtract the average positions */
2199 for (i = 0; i < edi->sav.nr; i++)
2201 rvec_dec(xcoll[i], edi->sav.x[i]);
2204 /* apply the constraints */
2207 do_linfix(xcoll, edi, step);
2209 do_linacc(xcoll, edi);
2212 do_radfix(xcoll, edi);
2214 do_radacc(xcoll, edi);
2215 do_radcon(xcoll, edi);
2217 /* add back the average positions */
2218 for (i = 0; i < edi->sav.nr; i++)
2220 rvec_inc(xcoll[i], edi->sav.x[i]);
2225 /* Write out the projections onto the eigenvectors. The order of output
2226 * corresponds to ed_output_legend() */
2227 static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
2232 /* Output how well we fit to the reference structure */
2233 fprintf(fp, EDcol_ffmt, rmsd);
2235 for (i = 0; i < edi->vecs.mon.neig; i++)
2237 fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
2240 for (i = 0; i < edi->vecs.linfix.neig; i++)
2242 fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
2245 for (i = 0; i < edi->vecs.linacc.neig; i++)
2247 fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
2250 for (i = 0; i < edi->vecs.radfix.neig; i++)
2252 fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
2254 if (edi->vecs.radfix.neig)
2256 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
2259 for (i = 0; i < edi->vecs.radacc.neig; i++)
2261 fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
2263 if (edi->vecs.radacc.neig)
2265 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
2268 for (i = 0; i < edi->vecs.radcon.neig; i++)
2270 fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
2272 if (edi->vecs.radcon.neig)
2274 fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
2278 /* Returns if any constraints are switched on */
2279 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2281 if (edtype == eEDedsam || edtype == eEDflood)
2283 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2284 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2285 edi->vecs.radcon.neig);
2291 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2292 * umbrella sampling simulations. */
2293 static void copyEvecReference(t_eigvec* floodvecs)
2298 if (NULL == floodvecs->refproj0)
2300 snew(floodvecs->refproj0, floodvecs->neig);
2303 for (i = 0; i < floodvecs->neig; i++)
2305 floodvecs->refproj0[i] = floodvecs->refproj[i];
2310 /* Call on MASTER only. Check whether the essential dynamics / flooding
2311 * groups of the checkpoint file are consistent with the provided .edi file. */
2312 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate)
2314 t_edpar *edi = NULL; /* points to a single edi data set */
2318 if (NULL == EDstate->nref || NULL == EDstate->nav)
2320 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2321 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2322 "it must also continue with/without ED constraints when checkpointing.\n"
2323 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2324 "from without a checkpoint.\n");
2331 /* Check number of atoms in the reference and average structures */
2332 if (EDstate->nref[edinum] != edi->sref.nr)
2334 gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2335 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2336 get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2338 if (EDstate->nav[edinum] != edi->sav.nr)
2340 gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2341 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2342 get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2344 edi = edi->next_edi;
2348 if (edinum != EDstate->nED)
2350 gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2351 "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2352 "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2357 /* The edsamstate struct stores the information we need to make the ED group
2358 * whole again after restarts from a checkpoint file. Here we do the following:
2359 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2360 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2361 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2362 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2363 * all ED structures at the last time step. */
2364 static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate)
2370 snew(EDstate->old_sref_p, EDstate->nED);
2371 snew(EDstate->old_sav_p, EDstate->nED);
2373 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2374 if (!EDstate->bFromCpt)
2376 snew(EDstate->nref, EDstate->nED);
2377 snew(EDstate->nav, EDstate->nED);
2380 /* Loop over all ED/flooding data sets (usually only one, though) */
2382 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2384 /* We always need the last reference and average positions such that
2385 * in the next time step we can make the ED group whole again
2386 * if the atoms do not have the correct PBC representation */
2387 if (EDstate->bFromCpt)
2389 /* Copy the last whole positions of reference and average group from .cpt */
2390 for (i = 0; i < edi->sref.nr; i++)
2392 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2394 for (i = 0; i < edi->sav.nr; i++)
2396 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2401 EDstate->nref[nr_edi-1] = edi->sref.nr;
2402 EDstate->nav [nr_edi-1] = edi->sav.nr;
2405 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2406 EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2407 EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2409 edi = edi->next_edi;
2414 /* Adds 'buf' to 'str' */
2415 static void add_to_string(char **str, char *buf)
2420 len = strlen(*str) + strlen(buf) + 1;
2426 static void add_to_string_aligned(char **str, char *buf)
2428 char buf_aligned[STRLEN];
2430 sprintf(buf_aligned, EDcol_sfmt, buf);
2431 add_to_string(str, buf_aligned);
2435 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, char *value, char *unit, char EDgroupchar)
2437 char tmp[STRLEN], tmp2[STRLEN];
2440 sprintf(tmp, "%c %s", EDgroupchar, value);
2441 add_to_string_aligned(LegendStr, tmp);
2442 sprintf(tmp2, "%s (%s)", tmp, unit);
2443 (*setname)[*nsets] = strdup(tmp2);
2448 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2454 for (i = 0; i < evec->neig; i++)
2456 sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2457 nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2462 /* Makes a legend for the xvg output file. Call on MASTER only! */
2463 static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv)
2465 t_edpar *edi = NULL;
2467 int nr_edi, nsets, n_flood, n_edsam;
2468 const char **setname;
2470 char *LegendStr = NULL;
2475 fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2477 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2479 /* Remember for each ED group whether we have to do essential dynamics
2480 * constraints or possibly only flooding */
2481 edi->bNeedDoEdsam = edi->vecs.mon.neig
2482 || edi->vecs.linfix.neig
2483 || edi->vecs.linacc.neig
2484 || edi->vecs.radfix.neig
2485 || edi->vecs.radacc.neig
2486 || edi->vecs.radcon.neig;
2488 fprintf(ed->edo, "#\n");
2489 fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2490 fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2491 fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2492 fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2493 fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2494 fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2495 fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2496 fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2497 fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2499 if (edi->flood.vecs.neig)
2501 /* If in any of the groups we find a flooding vector, flooding is turned on */
2502 ed->eEDtype = eEDflood;
2504 /* Print what flavor of flooding we will do */
2505 if (0 == edi->flood.tau) /* constant flooding strength */
2507 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2508 if (edi->flood.bHarmonic)
2510 fprintf(ed->edo, ", harmonic");
2513 else /* adaptive flooding */
2515 fprintf(ed->edo, ", adaptive");
2518 fprintf(ed->edo, "\n");
2520 edi = edi->next_edi;
2523 /* Print a nice legend */
2525 LegendStr[0] = '\0';
2526 sprintf(buf, "# %6s", "time");
2527 add_to_string(&LegendStr, buf);
2529 /* Calculate the maximum number of columns we could end up with */
2532 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2534 nsets += 5 +edi->vecs.mon.neig
2535 +edi->vecs.linfix.neig
2536 +edi->vecs.linacc.neig
2537 +edi->vecs.radfix.neig
2538 +edi->vecs.radacc.neig
2539 +edi->vecs.radcon.neig
2540 + 6*edi->flood.vecs.neig;
2541 edi = edi->next_edi;
2543 snew(setname, nsets);
2545 /* In the mdrun time step in a first function call (do_flood()) the flooding
2546 * forces are calculated and in a second function call (do_edsam()) the
2547 * ED constraints. To get a corresponding legend, we need to loop twice
2548 * over the edi groups and output first the flooding, then the ED part */
2550 /* The flooding-related legend entries, if flooding is done */
2552 if (eEDflood == ed->eEDtype)
2555 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2557 /* Always write out the projection on the flooding EVs. Of course, this can also
2558 * be achieved with the monitoring option in do_edsam() (if switched on by the
2559 * user), but in that case the positions need to be communicated in do_edsam(),
2560 * which is not necessary when doing flooding only. */
2561 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2563 for (i = 0; i < edi->flood.vecs.neig; i++)
2565 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2566 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2568 /* Output the current reference projection if it changes with time;
2569 * this can happen when flooding is used as harmonic restraint */
2570 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
2572 sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2573 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2576 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2577 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2579 sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2580 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2583 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2584 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2586 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2588 sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2589 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2592 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2593 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2596 edi = edi->next_edi;
2597 } /* End of flooding-related legend entries */
2601 /* Now the ED-related entries, if essential dynamics is done */
2603 for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2605 if (edi->bNeedDoEdsam) /* Only print ED legend if at least one ED option is on */
2607 nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2609 /* Essential dynamics, projections on eigenvectors */
2610 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" );
2611 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2612 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2613 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2614 if (edi->vecs.radfix.neig)
2616 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2618 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2619 if (edi->vecs.radacc.neig)
2621 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2623 nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2624 if (edi->vecs.radcon.neig)
2626 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2629 edi = edi->next_edi;
2630 } /* end of 'pure' essential dynamics legend entries */
2631 n_edsam = nsets - n_flood;
2633 xvgr_legend(ed->edo, nsets, setname, oenv);
2636 fprintf(ed->edo, "#\n"
2637 "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2638 n_flood, 1 == n_flood ? "" : "s",
2639 n_edsam, 1 == n_edsam ? "" : "s");
2640 fprintf(ed->edo, "%s", LegendStr);
2647 void init_edsam(gmx_mtop_t *mtop, /* global topology */
2648 t_inputrec *ir, /* input record */
2649 t_commrec *cr, /* communication record */
2650 gmx_edsam_t ed, /* contains all ED data */
2651 rvec x[], /* positions of the whole MD system */
2652 matrix box, /* the box */
2653 edsamstate_t *EDstate)
2655 t_edpar *edi = NULL; /* points to a single edi data set */
2656 int i, nr_edi, avindex;
2657 rvec *x_pbc = NULL; /* positions of the whole MD system with pbc removed */
2658 rvec *xfit = NULL, *xstart = NULL; /* dummy arrays to determine initial RMSDs */
2659 rvec fit_transvec; /* translation ... */
2660 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2663 if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
2665 gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
2670 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2674 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n"
2675 "flooding simulation. Please also provide the correct .edi file with -ei.\n");
2679 /* Needed for initializing radacc radius in do_edsam */
2682 /* The input file is read by the master and the edi structures are
2683 * initialized here. Input is stored in ed->edpar. Then the edi
2684 * structures are transferred to the other nodes */
2687 /* Initialization for every ED/flooding group. Flooding uses one edi group per
2688 * flooding vector, Essential dynamics can be applied to more than one structure
2689 * as well, but will be done in the order given in the edi file, so
2690 * expect different results for different order of edi file concatenation! */
2694 init_edi(mtop, edi);
2695 init_flood(edi, ed, ir->delta_t);
2696 edi = edi->next_edi;
2700 /* The master does the work here. The other nodes get the positions
2701 * not before dd_partition_system which is called after init_edsam */
2704 /* Remove pbc, make molecule whole.
2705 * When ir->bContinuation=TRUE this has already been done, but ok.
2707 snew(x_pbc, mtop->natoms);
2708 m_rveccopy(mtop->natoms, x, x_pbc);
2709 do_pbc_first_mtop(NULL, ir->ePBC, box, mtop, x_pbc);
2711 /* Reset pointer to first ED data set which contains the actual ED data */
2713 /* Loop over all ED/flooding data sets (usually only one, though) */
2714 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2716 /* For multiple ED groups we use the output frequency that was specified
2717 * in the first set */
2720 edi->outfrq = ed->edpar->outfrq;
2723 /* Extract the initial reference and average positions. When starting
2724 * from .cpt, these have already been read into sref.x_old
2725 * in init_edsamstate() */
2726 if (!EDstate->bFromCpt)
2728 /* If this is the first run (i.e. no checkpoint present) we assume
2729 * that the starting positions give us the correct PBC representation */
2730 for (i = 0; i < edi->sref.nr; i++)
2732 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2735 for (i = 0; i < edi->sav.nr; i++)
2737 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2741 /* Now we have the PBC-correct start positions of the reference and
2742 average structure. We copy that over to dummy arrays on which we
2743 can apply fitting to print out the RMSD. We srenew the memory since
2744 the size of the buffers is likely different for every ED group */
2745 srenew(xfit, edi->sref.nr );
2746 srenew(xstart, edi->sav.nr );
2747 copy_rvecn(edi->sref.x_old, xfit, 0, edi->sref.nr);
2748 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2750 /* Make the fit to the REFERENCE structure, get translation and rotation */
2751 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2753 /* Output how well we fit to the reference at the start */
2754 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2755 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2756 rmsd_from_structure(xfit, &edi->sref));
2757 if (EDstate->nED > 1)
2759 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2761 fprintf(stderr, "\n");
2763 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2764 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2766 /* calculate initial projections */
2767 project(xstart, edi);
2769 /* For the target and origin structure both a reference (fit) and an
2770 * average structure can be provided in make_edi. If both structures
2771 * are the same, make_edi only stores one of them in the .edi file.
2772 * If they differ, first the fit and then the average structure is stored
2773 * in star (or sor), thus the number of entries in star/sor is
2774 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2775 * the size of the average group. */
2777 /* process target structure, if required */
2778 if (edi->star.nr > 0)
2780 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2782 /* get translation & rotation for fit of target structure to reference structure */
2783 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2785 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2786 if (edi->star.nr == edi->sav.nr)
2790 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2792 /* The last sav.nr indices of the target structure correspond to
2793 * the average structure, which must be projected */
2794 avindex = edi->star.nr - edi->sav.nr;
2796 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
2800 rad_project(edi, xstart, &edi->vecs.radcon);
2803 /* process structure that will serve as origin of expansion circle */
2804 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2806 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2809 if (edi->sori.nr > 0)
2811 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2813 /* fit this structure to reference structure */
2814 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2816 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2817 if (edi->sori.nr == edi->sav.nr)
2821 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2823 /* For the projection, we need the last sav.nr indices of sori */
2824 avindex = edi->sori.nr - edi->sav.nr;
2827 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2828 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2829 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2831 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2832 /* Set center of flooding potential to the ORIGIN structure */
2833 rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
2834 /* We already know that no (moving) reference position was provided,
2835 * therefore we can overwrite refproj[0]*/
2836 copyEvecReference(&edi->flood.vecs);
2839 else /* No origin structure given */
2841 rad_project(edi, xstart, &edi->vecs.radacc);
2842 rad_project(edi, xstart, &edi->vecs.radfix);
2843 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2845 if (edi->flood.bHarmonic)
2847 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2848 for (i = 0; i < edi->flood.vecs.neig; i++)
2850 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2855 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2856 /* Set center of flooding potential to the center of the covariance matrix,
2857 * i.e. the average structure, i.e. zero in the projected system */
2858 for (i = 0; i < edi->flood.vecs.neig; i++)
2860 edi->flood.vecs.refproj[i] = 0.0;
2865 /* For convenience, output the center of the flooding potential for the eigenvectors */
2866 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2868 for (i = 0; i < edi->flood.vecs.neig; i++)
2870 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2871 if (edi->flood.bHarmonic)
2873 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2875 fprintf(stdout, "\n");
2879 /* set starting projections for linsam */
2880 rad_project(edi, xstart, &edi->vecs.linacc);
2881 rad_project(edi, xstart, &edi->vecs.linfix);
2883 /* Prepare for the next edi data set: */
2884 edi = edi->next_edi;
2886 /* Cleaning up on the master node: */
2891 } /* end of MASTER only section */
2895 /* First let everybody know how many ED data sets to expect */
2896 gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr);
2897 /* Broadcast the essential dynamics / flooding data to all nodes */
2898 broadcast_ed_data(cr, ed, EDstate->nED);
2902 /* In the single-CPU case, point the local atom numbers pointers to the global
2903 * one, so that we can use the same notation in serial and parallel case: */
2905 /* Loop over all ED data sets (usually only one, though) */
2907 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2909 edi->sref.anrs_loc = edi->sref.anrs;
2910 edi->sav.anrs_loc = edi->sav.anrs;
2911 edi->star.anrs_loc = edi->star.anrs;
2912 edi->sori.anrs_loc = edi->sori.anrs;
2913 /* For the same reason as above, make a dummy c_ind array: */
2914 snew(edi->sav.c_ind, edi->sav.nr);
2915 /* Initialize the array */
2916 for (i = 0; i < edi->sav.nr; i++)
2918 edi->sav.c_ind[i] = i;
2920 /* In the general case we will need a different-sized array for the reference indices: */
2923 snew(edi->sref.c_ind, edi->sref.nr);
2924 for (i = 0; i < edi->sref.nr; i++)
2926 edi->sref.c_ind[i] = i;
2929 /* Point to the very same array in case of other structures: */
2930 edi->star.c_ind = edi->sav.c_ind;
2931 edi->sori.c_ind = edi->sav.c_ind;
2932 /* In the serial case, the local number of atoms is the global one: */
2933 edi->sref.nr_loc = edi->sref.nr;
2934 edi->sav.nr_loc = edi->sav.nr;
2935 edi->star.nr_loc = edi->star.nr;
2936 edi->sori.nr_loc = edi->sori.nr;
2938 /* An on we go to the next ED group */
2939 edi = edi->next_edi;
2943 /* Allocate space for ED buffer variables */
2944 /* Again, loop over ED data sets */
2946 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2948 /* Allocate space for ED buffer */
2950 snew(edi->buf->do_edsam, 1);
2952 /* Space for collective ED buffer variables */
2954 /* Collective positions of atoms with the average indices */
2955 snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2956 snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
2957 snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2958 /* Collective positions of atoms with the reference indices */
2961 snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2962 snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
2963 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2966 /* Get memory for flooding forces */
2967 snew(edi->flood.forces_cartesian, edi->sav.nr);
2970 /* Dump it all into one file per process */
2971 dump_edi(edi, cr, nr_edi);
2975 edi = edi->next_edi;
2978 /* Flush the edo file so that the user can check some things
2979 * when the simulation has started */
2987 void do_edsam(t_inputrec *ir,
2988 gmx_large_int_t step,
2990 rvec xs[], /* The local current positions on this processor */
2991 rvec v[], /* The velocities */
2995 int i, edinr, iupdate = 500;
2996 matrix rotmat; /* rotation matrix */
2997 rvec transvec; /* translation vector */
2998 rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
2999 real dt_1; /* 1/dt */
3000 struct t_do_edsam *buf;
3002 real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */
3003 gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3006 /* Check if ED sampling has to be performed */
3007 if (ed->eEDtype == eEDnone)
3012 /* Suppress output on first call of do_edsam if
3013 * two-step sd2 integrator is used */
3014 if ( (ir->eI == eiSD2) && (v != NULL) )
3019 dt_1 = 1.0/ir->delta_t;
3021 /* Loop over all ED groups (usually one) */
3027 if (edi->bNeedDoEdsam)
3030 buf = edi->buf->do_edsam;
3034 /* initialise radacc radius for slope criterion */
3035 buf->oldrad = calc_radius(&edi->vecs.radacc);
3038 /* Copy the positions into buf->xc* arrays and after ED
3039 * feed back corrections to the official positions */
3041 /* Broadcast the ED positions such that every node has all of them
3042 * Every node contributes its local positions xs and stores it in
3043 * the collective buf->xcoll array. Note that for edinr > 1
3044 * xs could already have been modified by an earlier ED */
3046 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3047 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
3049 /* Only assembly reference positions if their indices differ from the average ones */
3052 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3053 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3056 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3057 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3058 * set bUpdateShifts=TRUE in the parallel case. */
3059 buf->bUpdateShifts = FALSE;
3061 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3062 * as well as the indices in edi->sav.anrs */
3064 /* Fit the reference indices to the reference structure */
3067 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3071 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3074 /* Now apply the translation and rotation to the ED structure */
3075 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3077 /* Find out how well we fit to the reference (just for output steps) */
3078 if (do_per_step(step, edi->outfrq) && MASTER(cr))
3082 /* Indices of reference and average structures are identical,
3083 * thus we can calculate the rmsd to SREF using xcoll */
3084 rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3088 /* We have to translate & rotate the reference atoms first */
3089 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3090 rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3094 /* update radsam references, when required */
3095 if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3097 project(buf->xcoll, edi);
3098 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3099 rad_project(edi, buf->xcoll, &edi->vecs.radfix);
3100 buf->oldrad = -1.e5;
3103 /* update radacc references, when required */
3104 if (do_per_step(step, iupdate) && step >= edi->presteps)
3106 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
3107 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3109 project(buf->xcoll, edi);
3110 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3115 buf->oldrad = edi->vecs.radacc.radius;
3119 /* apply the constraints */
3120 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3122 /* ED constraints should be applied already in the first MD step
3123 * (which is step 0), therefore we pass step+1 to the routine */
3124 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3127 /* write to edo, when required */
3128 if (do_per_step(step, edi->outfrq))
3130 project(buf->xcoll, edi);
3131 if (MASTER(cr) && !bSuppress)
3133 write_edo(edi, ed->edo, rmsdev);
3137 /* Copy back the positions unless monitoring only */
3138 if (ed_constraints(ed->eEDtype, edi))
3140 /* remove fitting */
3141 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3143 /* Copy the ED corrected positions into the coordinate array */
3144 /* Each node copies its local part. In the serial case, nat_loc is the
3145 * total number of ED atoms */
3146 for (i = 0; i < edi->sav.nr_loc; i++)
3148 /* Unshift local ED coordinate and store in x_unsh */
3149 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3150 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3152 /* dx is the ED correction to the positions: */
3153 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3157 /* dv is the ED correction to the velocity: */
3158 svmul(dt_1, dx, dv);
3159 /* apply the velocity correction: */
3160 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3162 /* Finally apply the position correction due to ED: */
3163 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3166 } /* END of if (edi->bNeedDoEdsam) */
3168 /* Prepare for the next ED group */
3169 edi = edi->next_edi;
3171 } /* END of loop over ED groups */