src/gromacs/essentialdynamics/edsam.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "edsam.h"
  41
  42 #include <cmath>
  43 #include <cstdio>
  44 #include <cstring>
  45 #include <ctime>
  46
  47 #include <memory>
  48
  49 #include "gromacs/commandline/filenm.h"
  50 #include "gromacs/domdec/domdec_struct.h"
  51 #include "gromacs/fileio/gmxfio.h"
  52 #include "gromacs/fileio/xvgr.h"
  53 #include "gromacs/gmxlib/network.h"
  54 #include "gromacs/gmxlib/nrnb.h"
  55 #include "gromacs/linearalgebra/nrjac.h"
  56 #include "gromacs/math/functions.h"
  57 #include "gromacs/math/utilities.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/math/vectypes.h"
  60 #include "gromacs/mdlib/broadcaststructs.h"
  61 #include "gromacs/mdlib/constr.h"
  62 #include "gromacs/mdlib/groupcoord.h"
  63 #include "gromacs/mdlib/stat.h"
  64 #include "gromacs/mdlib/update.h"
  65 #include "gromacs/mdrunutility/handlerestart.h"
  66 #include "gromacs/mdtypes/commrec.h"
  67 #include "gromacs/mdtypes/edsamhistory.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/md_enums.h"
  70 #include "gromacs/mdtypes/observableshistory.h"
  71 #include "gromacs/mdtypes/state.h"
  72 #include "gromacs/pbcutil/pbc.h"
  73 #include "gromacs/topology/mtop_lookup.h"
  74 #include "gromacs/topology/topology.h"
  75 #include "gromacs/utility/cstringutil.h"
  76 #include "gromacs/utility/fatalerror.h"
  77 #include "gromacs/utility/gmxassert.h"
  78 #include "gromacs/utility/logger.h"
  79 #include "gromacs/utility/smalloc.h"
  80 #include "gromacs/utility/strconvert.h"
  81
  82 namespace
  83 {
  84 /*! \brief Identifies the type of ED: none, normal ED, flooding. */
  85 enum class EssentialDynamicsType
  86 {
  87     //! No essential dynamics
  88     None,
  89     //! Essential dynamics sampling
  90     EDSampling,
  91     //! Essential dynamics flooding
  92     Flooding
  93 };
  94
  95 /*! \brief Identify on which structure to operate. */
  96 enum class EssentialDynamicsStructure
  97 {
  98     //! Reference structure
  99     Reference,
 100     //! Average Structure
 101     Average,
 102     //! Origin Structure
 103     Origin,
 104     //! Target Structure
 105     Target
 106 };
 107
 108 /*! \brief Essential dynamics vector.
 109  * TODO: split into working data and input data
 110  * NOTE: called eigvec, because traditionally eigenvectors from PCA analysis
 111  * were used as essential dynamics vectors, however, vectors used for ED need
 112  * not have any special properties
 113  */
 114 struct t_eigvec
 115 {
 116     //! nr of eigenvectors
 117     int neig = 0;
 118     //! index nrs of eigenvectors
 119     int* ieig = nullptr;
 120     //! stepsizes (per eigenvector)
 121     real* stpsz = nullptr;
 122     //! eigenvector components
 123     rvec** vec = nullptr;
 124     //! instantaneous x projections
 125     real* xproj = nullptr;
 126     //! instantaneous f projections
 127     real* fproj = nullptr;
 128     //! instantaneous radius
 129     real radius = 0.;
 130     //! starting or target projections
 131     real* refproj = nullptr;
 132 };
 133
 134 /*! \brief Essential dynamics vectors per method implementation.
 135  */
 136 struct t_edvecs
 137 {
 138     //! only monitored, no constraints
 139     t_eigvec mon = {};
 140     //! fixed linear constraints
 141     t_eigvec linfix = {};
 142     //! acceptance linear constraints
 143     t_eigvec linacc = {};
 144     //! fixed radial constraints (exp)
 145     t_eigvec radfix = {};
 146     //! acceptance radial constraints (exp)
 147     t_eigvec radacc = {};
 148     //! acceptance rad. contraction constr.
 149     t_eigvec radcon = {};
 150 };
 151
 152 /*! \brief Essential dynamics flooding parameters and work data.
 153  * TODO: split into working data and input parameters
 154  * NOTE: The implementation here follows:
 155  * O.E. Lange, L.V. Schafer, and H. Grubmuller,
 156  * “Flooding in GROMACS: Accelerated barrier crossings in molecular dynamics,”
 157  *  J. Comp. Chem., 27 1693–1702 (2006)
 158  */
 159 struct t_edflood
 160 {
 161     /*! \brief Target destabilisation free energy.
 162      * Controls the flooding potential strength.
 163      * Linked to the expected speed-up of mean first passage time out of flooded minimum */
 164     real deltaF0 = 0;
 165     //! Do not calculate a flooding potential, instead flood with a constant force
 166     bool bConstForce = false;
 167     //! Relaxation time scale for slowest flooded degree of freedom
 168     real tau = 0;
 169     //! Current estimated destabilisation free energy
 170     real deltaF = 0;
 171     //! Flooding energy, acting as a proportionality factor for the flooding potential
 172     real Efl = 0;
 173     //! Boltzmann constant times temperature, provided by user
 174     real kT = 0;
 175     //! The flooding potential
 176     real Vfl = 0;
 177     //! Integration time step
 178     real dt = 0;
 179     //! Inital flooding strenth
 180     real constEfl = 0;
 181     //! Empirical global scaling parameter of the essential dynamics vectors.
 182     real alpha2 = 0;
 183     //! The forces from flooding in atom coordinate space (in contrast to projections onto essential dynamics vectors)
 184     rvec* forces_cartesian = nullptr;
 185     //! The vectors along which to flood
 186     t_eigvec vecs = {};
 187     //! Use flooding for harmonic restraint on the eigenvector
 188     bool bHarmonic = false;
 189     //! The initial reference projection of the flooding vectors. Only with harmonic restraint.
 190     real* initialReferenceProjection = nullptr;
 191     //! The current reference projection is the initialReferenceProjection + step * slope. Only with harmonic restraint.
 192     real* referenceProjectionSlope = nullptr;
 193 };
 194 } // namespace
 195
 196
 197 /* This type is for the average, reference, target, and origin structure    */
 198 struct gmx_edx
 199 {
 200     int  nr         = 0;       /* number of atoms this structure contains  */
 201     int  nr_loc     = 0;       /* number of atoms on local node            */
 202     int* anrs       = nullptr; /* atom index numbers                       */
 203     int* anrs_loc   = nullptr; /* local atom index numbers                 */
 204     int  nalloc_loc = 0;       /* allocation size of anrs_loc              */
 205     int* c_ind      = nullptr; /* at which position of the whole anrs
 206                                 * array is a local atom?, i.e.
 207                                 * c_ind[0...nr_loc-1] gives the atom index
 208                                 * with respect to the collective
 209                                 * anrs[0...nr-1] array                     */
 210     rvec* x     = nullptr;     /* positions for this structure             */
 211     rvec* x_old = nullptr;     /* Last positions which have the correct PBC
 212                                   representation of the ED group. In
 213                                   combination with keeping track of the
 214                                   shift vectors, the ED group can always
 215                                   be made whole                            */
 216     real* m     = nullptr;     /* masses                                   */
 217     real  mtot  = 0.;          /* total mass (only used in sref)           */
 218     real* sqrtm = nullptr;     /* sqrt of the masses used for mass-
 219                                 * weighting of analysis (only used in sav) */
 220 };
 221
 222
 223 typedef struct edpar
 224 {
 225     int      nini     = 0;     /* total Nr of atoms                    */
 226     gmx_bool fitmas   = false; /* true if trans fit with cm            */
 227     gmx_bool pcamas   = false; /* true if mass-weighted PCA            */
 228     int      presteps = 0;     /* number of steps to run without any
 229                                 *    perturbations ... just monitoring */
 230     int outfrq     = 0;        /* freq (in steps) of writing to edo    */
 231     int maxedsteps = 0;        /* max nr of steps per cycle            */
 232
 233     /* all gmx_edx datasets are copied to all nodes in the parallel case   */
 234     gmx_edx sref = {};         /* reference positions, to these fitting
 235                                 * will be done                         */
 236     gmx_bool bRefEqAv = false; /* If true, reference & average indices
 237                                 * are the same. Used for optimization  */
 238     gmx_edx sav  = {};         /* average positions                    */
 239     gmx_edx star = {};         /* target positions                     */
 240     gmx_edx sori = {};         /* origin positions                     */
 241
 242     t_edvecs vecs  = {}; /* eigenvectors                         */
 243     real     slope = 0;  /* minimal slope in acceptance radexp   */
 244
 245     t_edflood           flood = {};      /* parameters especially for flooding   */
 246     struct t_ed_buffer* buf   = nullptr; /* handle to local buffers              */
 247 } t_edpar;
 248
 249
 250 struct gmx_edsam
 251 {
 252     ~gmx_edsam();
 253     //! The type of ED
 254     EssentialDynamicsType eEDtype = EssentialDynamicsType::None;
 255     //! output file pointer
 256     FILE*                edo = nullptr;
 257     std::vector<t_edpar> edpar;
 258     gmx_bool             bFirst = false;
 259 };
 260 gmx_edsam::~gmx_edsam()
 261 {
 262     if (edo)
 263     {
 264         /* edo is opened sometimes with xvgropen, sometimes with
 265          * gmx_fio_fopen, so we use the least common denominator for
 266          * closing. */
 267         gmx_fio_fclose(edo);
 268     }
 269 }
 270
 271 struct t_do_edsam
 272 {
 273     matrix old_rotmat;
 274     real   oldrad;
 275     rvec   old_transvec, older_transvec, transvec_compact;
 276     rvec*  xcoll;                 /* Positions from all nodes, this is the
 277                                      collective set we work on.
 278                                      These are the positions of atoms with
 279                                      average structure indices */
 280     rvec*    xc_ref;              /* same but with reference structure indices */
 281     ivec*    shifts_xcoll;        /* Shifts for xcoll  */
 282     ivec*    extra_shifts_xcoll;  /* xcoll shift changes since last NS step */
 283     ivec*    shifts_xc_ref;       /* Shifts for xc_ref */
 284     ivec*    extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
 285     gmx_bool bUpdateShifts;       /* TRUE in NS steps to indicate that the
 286                                      ED shifts for this ED group need to
 287                                      be updated */
 288 };
 289
 290
 291 /* definition of ED buffer structure */
 292 struct t_ed_buffer
 293 {
 294     struct t_fit_to_ref* fit_to_ref;
 295     struct t_do_edfit*   do_edfit;
 296     struct t_do_edsam*   do_edsam;
 297     struct t_do_radcon*  do_radcon;
 298 };
 299
 300 namespace gmx
 301 {
 302 class EssentialDynamics::Impl
 303 {
 304 public:
 305     // TODO: move all fields from gmx_edsam here and remove gmx_edsam
 306     gmx_edsam essentialDynamics_;
 307 };
 308 EssentialDynamics::EssentialDynamics() : impl_(new Impl) {}
 309 EssentialDynamics::~EssentialDynamics() = default;
 310
 311 gmx_edsam* EssentialDynamics::getLegacyED()
 312 {
 313     return &impl_->essentialDynamics_;
 314 }
 315 } // namespace gmx
 316
 317 /* Function declarations */
 318 static void fit_to_reference(rvec* xcoll, rvec transvec, matrix rotmat, t_edpar* edi);
 319 static void translate_and_rotate(rvec* x, int nat, rvec transvec, matrix rotmat);
 320 static real rmsd_from_structure(rvec* x, struct gmx_edx* s);
 321 namespace
 322 {
 323 /*! \brief Read in the essential dynamics input file and return its contents.
 324  * \param[in] fn the file name of the edi file to be read
 325  * \param[in] nr_mdatoms the number of atoms in the simulation
 326  * \returns A vector of containing the essentiail dyanmics parameters.
 327  * NOTE: edi files may that it may contain several ED data sets from concatenated edi files.
 328  * The standard case would be a single ED data set, though. */
 329 std::vector<t_edpar> read_edi_file(const char* fn, int nr_mdatoms);
 330 } // namespace
 331 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam& ed, edsamhistory_t* EDstate);
 332 static void init_edsamstate(const gmx_edsam& ed, edsamhistory_t* EDstate);
 333 static void write_edo_legend(gmx_edsam* ed, int nED, const gmx_output_env_t* oenv);
 334 /* End function declarations */
 335
 336 /* Define formats for the column width in the output file */
 337 const char EDcol_sfmt[] = "%17s";
 338 const char EDcol_efmt[] = "%17.5e";
 339 const char EDcol_ffmt[] = "%17f";
 340
 341 /* Define a formats for reading, otherwise cppcheck complains for scanf without width limits */
 342 const char max_ev_fmt_d[]       = "%7d"; /* Number of eigenvectors. 9,999,999 should be enough */
 343 const char max_ev_fmt_lf[]      = "%12lf";
 344 const char max_ev_fmt_dlf[]     = "%7d%12lf";
 345 const char max_ev_fmt_dlflflf[] = "%7d%12lf%12lf%12lf";
 346 const char max_ev_fmt_lelele[]  = "%12le%12le%12le";
 347
 348 namespace
 349 {
 350 /*!\brief Determines whether to perform essential dynamics constraints other than flooding.
 351  * \param[in] edi the essential dynamics parameters
 352  * \returns true if essential dyanmics constraints need to be performed
 353  */
 354 bool bNeedDoEdsam(const t_edpar& edi)
 355 {
 356     return (edi.vecs.mon.neig != 0) || (edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0)
 357            || (edi.vecs.radfix.neig != 0) || (edi.vecs.radacc.neig != 0) || (edi.vecs.radcon.neig != 0);
 358 }
 359 } // namespace
 360
 361
 362 /* Multiple ED groups will be labeled with letters instead of numbers
 363  * to avoid confusion with eigenvector indices */
 364 static char get_EDgroupChar(int nr_edi, int nED)
 365 {
 366     if (nED == 1)
 367     {
 368         return ' ';
 369     }
 370
 371     /* nr_edi = 1 -> A
 372      * nr_edi = 2 -> B ...
 373      */
 374     return 'A' + nr_edi - 1;
 375 }
 376
 377 namespace
 378 {
 379 /*! \brief The mass-weighted inner product of two coordinate vectors.
 380  * Does not subtract average positions, projection on single eigenvector is returned
 381  * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
 382  * Average position is subtracted in ed_apply_constraints prior to calling projectx
 383  * \param[in] edi Essential dynamics parameters
 384  * \param[in] xcoll vector of atom coordinates
 385  * \param[in] vec vector of coordinates to project onto
 386  * \return mass-weighted projection.
 387  */
 388 real projectx(const t_edpar& edi, rvec* xcoll, rvec* vec)
 389 {
 390     int  i;
 391     real proj = 0.0;
 392
 393
 394     for (i = 0; i < edi.sav.nr; i++)
 395     {
 396         proj += edi.sav.sqrtm[i] * iprod(vec[i], xcoll[i]);
 397     }
 398
 399     return proj;
 400 }
 401 /*!\brief Project coordinates onto vector after substracting average position.
 402  * projection is stored in vec->refproj which is used for radacc, radfix,
 403  * radcon and center of flooding potential.
 404  * Average positions are first substracted from x, then added back again.
 405  * \param[in] edi essential dynamics parameters with average position
 406  * \param[in] x Coordinates to be projected
 407  * \param[out] vec eigenvector, radius and refproj are overwritten here
 408  */
 409 void rad_project(const t_edpar& edi, rvec* x, t_eigvec* vec)
 410 {
 411     int  i;
 412     real rad = 0.0;
 413
 414     /* Subtract average positions */
 415     for (i = 0; i < edi.sav.nr; i++)
 416     {
 417         rvec_dec(x[i], edi.sav.x[i]);
 418     }
 419
 420     for (i = 0; i < vec->neig; i++)
 421     {
 422         vec->refproj[i] = projectx(edi, x, vec->vec[i]);
 423         rad += gmx::square((vec->refproj[i] - vec->xproj[i]));
 424     }
 425     vec->radius = sqrt(rad);
 426
 427     /* Add average positions */
 428     for (i = 0; i < edi.sav.nr; i++)
 429     {
 430         rvec_inc(x[i], edi.sav.x[i]);
 431     }
 432 }
 433
 434 /*!\brief Projects coordinates onto eigenvectors and stores result in vec->xproj.
 435  * Mass-weighting is applied. Subtracts average positions prior to projection and add
 436  * them afterwards to retain the unchanged vector.
 437  * \param[in] x The coordinates to project to an eigenvector
 438  * \param[in,out] vec The eigenvectors
 439  * \param[in] edi essential dynamics parameters holding average structure and masses
 440  */
 441 void project_to_eigvectors(rvec* x, t_eigvec* vec, const t_edpar& edi)
 442 {
 443     if (!vec->neig)
 444     {
 445         return;
 446     }
 447
 448     /* Subtract average positions */
 449     for (int i = 0; i < edi.sav.nr; i++)
 450     {
 451         rvec_dec(x[i], edi.sav.x[i]);
 452     }
 453
 454     for (int i = 0; i < vec->neig; i++)
 455     {
 456         vec->xproj[i] = projectx(edi, x, vec->vec[i]);
 457     }
 458
 459     /* Add average positions */
 460     for (int i = 0; i < edi.sav.nr; i++)
 461     {
 462         rvec_inc(x[i], edi.sav.x[i]);
 463     }
 464 }
 465 } // namespace
 466
 467 /* Project vector x onto all edi->vecs (mon, linfix,...) */
 468 static void project(rvec*    x,   /* positions to project */
 469                     t_edpar* edi) /* edi data set */
 470 {
 471     /* It is not more work to subtract the average position in every
 472      * subroutine again, because these routines are rarely used simultaneously */
 473     project_to_eigvectors(x, &edi->vecs.mon, *edi);
 474     project_to_eigvectors(x, &edi->vecs.linfix, *edi);
 475     project_to_eigvectors(x, &edi->vecs.linacc, *edi);
 476     project_to_eigvectors(x, &edi->vecs.radfix, *edi);
 477     project_to_eigvectors(x, &edi->vecs.radacc, *edi);
 478     project_to_eigvectors(x, &edi->vecs.radcon, *edi);
 479 }
 480
 481 namespace
 482 {
 483 /*!\brief Evaluates the distance from reference to current eigenvector projection.
 484  * \param[in] vec eigenvector
 485  * \returns distance
 486  */
 487 real calc_radius(const t_eigvec& vec)
 488 {
 489     real rad = 0.0;
 490
 491     for (int i = 0; i < vec.neig; i++)
 492     {
 493         rad += gmx::square((vec.refproj[i] - vec.xproj[i]));
 494     }
 495
 496     return rad = sqrt(rad);
 497 }
 498 } // namespace
 499
 500 struct t_do_edfit
 501 {
 502     double** omega;
 503     double** om;
 504 };
 505
 506 static void do_edfit(int natoms, rvec* xp, rvec* x, matrix R, t_edpar* edi)
 507 {
 508     /* this is a copy of do_fit with some modifications */
 509     int    c, r, n, j, i, irot;
 510     double d[6], xnr, xpc;
 511     matrix vh, vk, u;
 512     int    index;
 513     real   max_d;
 514
 515     struct t_do_edfit* loc;
 516     gmx_bool           bFirst;
 517
 518     if (edi->buf->do_edfit != nullptr)
 519     {
 520         bFirst = FALSE;
 521     }
 522     else
 523     {
 524         bFirst = TRUE;
 525         snew(edi->buf->do_edfit, 1);
 526     }
 527     loc = edi->buf->do_edfit;
 528
 529     if (bFirst)
 530     {
 531         snew(loc->omega, 2 * DIM);
 532         snew(loc->om, 2 * DIM);
 533         for (i = 0; i < 2 * DIM; i++)
 534         {
 535             snew(loc->omega[i], 2 * DIM);
 536             snew(loc->om[i], 2 * DIM);
 537         }
 538     }
 539
 540     for (i = 0; (i < 6); i++)
 541     {
 542         d[i] = 0;
 543         for (j = 0; (j < 6); j++)
 544         {
 545             loc->omega[i][j] = 0;
 546             loc->om[i][j]    = 0;
 547         }
 548     }
 549
 550     /* calculate the matrix U */
 551     clear_mat(u);
 552     for (n = 0; (n < natoms); n++)
 553     {
 554         for (c = 0; (c < DIM); c++)
 555         {
 556             xpc = xp[n][c];
 557             for (r = 0; (r < DIM); r++)
 558             {
 559                 xnr = x[n][r];
 560                 u[c][r] += xnr * xpc;
 561             }
 562         }
 563     }
 564
 565     /* construct loc->omega */
 566     /* loc->omega is symmetric -> loc->omega==loc->omega' */
 567     for (r = 0; (r < 6); r++)
 568     {
 569         for (c = 0; (c <= r); c++)
 570         {
 571             if ((r >= 3) && (c < 3))
 572             {
 573                 loc->omega[r][c] = u[r - 3][c];
 574                 loc->omega[c][r] = u[r - 3][c];
 575             }
 576             else
 577             {
 578                 loc->omega[r][c] = 0;
 579                 loc->omega[c][r] = 0;
 580             }
 581         }
 582     }
 583
 584     /* determine h and k */
 585     jacobi(loc->omega, 6, d, loc->om, &irot);
 586
 587     if (irot == 0)
 588     {
 589         fprintf(stderr, "IROT=0\n");
 590     }
 591
 592     index = 0; /* For the compiler only */
 593
 594     for (j = 0; (j < 3); j++)
 595     {
 596         max_d = -1000;
 597         for (i = 0; (i < 6); i++)
 598         {
 599             if (d[i] > max_d)
 600             {
 601                 max_d = d[i];
 602                 index = i;
 603             }
 604         }
 605         d[index] = -10000;
 606         for (i = 0; (i < 3); i++)
 607         {
 608             vh[j][i] = M_SQRT2 * loc->om[i][index];
 609             vk[j][i] = M_SQRT2 * loc->om[i + DIM][index];
 610         }
 611     }
 612
 613     /* determine R */
 614     for (c = 0; (c < 3); c++)
 615     {
 616         for (r = 0; (r < 3); r++)
 617         {
 618             R[c][r] = vk[0][r] * vh[0][c] + vk[1][r] * vh[1][c] + vk[2][r] * vh[2][c];
 619         }
 620     }
 621     if (det(R) < 0)
 622     {
 623         for (c = 0; (c < 3); c++)
 624         {
 625             for (r = 0; (r < 3); r++)
 626             {
 627                 R[c][r] = vk[0][r] * vh[0][c] + vk[1][r] * vh[1][c] - vk[2][r] * vh[2][c];
 628             }
 629         }
 630     }
 631 }
 632
 633
 634 static void rmfit(int nat, rvec* xcoll, const rvec transvec, matrix rotmat)
 635 {
 636     rvec   vec;
 637     matrix tmat;
 638
 639
 640     /* Remove rotation.
 641      * The inverse rotation is described by the transposed rotation matrix */
 642     transpose(rotmat, tmat);
 643     rotate_x(xcoll, nat, tmat);
 644
 645     /* Remove translation */
 646     vec[XX] = -transvec[XX];
 647     vec[YY] = -transvec[YY];
 648     vec[ZZ] = -transvec[ZZ];
 649     translate_x(xcoll, nat, vec);
 650 }
 651
 652
 653 /**********************************************************************************
 654  ******************** FLOODING ****************************************************
 655  **********************************************************************************
 656
 657    The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
 658    The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
 659    read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
 660
 661    do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
 662    the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
 663
 664    since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
 665    edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
 666
 667    do_flood makes a copy of the positions,
 668    fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
 669    space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
 670    fit. Then do_flood adds these forces to the forcefield-forces
 671    (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
 672
 673    To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
 674    structure is projected to the system of eigenvectors and then this position in the subspace is used as
 675    center of the flooding potential.   If the option is not used, the center will be zero in the subspace,
 676    i.e. the average structure as given in the make_edi file.
 677
 678    To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
 679    signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
 680    Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
 681    so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
 682    further adaption is applied, Efl will stay constant at zero.
 683
 684    To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
 685    used as spring constants for the harmonic potential.
 686    Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
 687    as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
 688
 689    To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
 690    the routine read_edi_file reads all of theses flooding files.
 691    The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
 692    calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
 693    edi there is no interdependence whatsoever. The forces are added together.
 694
 695    To write energies into the .edr file, call the function
 696         get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
 697    and call
 698         get_flood_energies(real Vfl[],int nnames);
 699
 700    TODO:
 701    - 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.
 702
 703    Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
 704    two edsam files from two peptide chains
 705  */
 706
 707 // TODO split this into multiple files
 708 namespace
 709 {
 710 /*!\brief Output flooding simulation settings and results to file.
 711  * \param[in] edi Essential dynamics input parameters
 712  * \param[in] fp output file
 713  * \param[in] rmsd rmsd to reference structure
 714  */
 715 void write_edo_flood(const t_edpar& edi, FILE* fp, real rmsd)
 716 {
 717     /* Output how well we fit to the reference structure */
 718     fprintf(fp, EDcol_ffmt, rmsd);
 719
 720     for (int i = 0; i < edi.flood.vecs.neig; i++)
 721     {
 722         fprintf(fp, EDcol_efmt, edi.flood.vecs.xproj[i]);
 723
 724         /* Check whether the reference projection changes with time (this can happen
 725          * in case flooding is used as harmonic restraint). If so, output the
 726          * current reference projection */
 727         if (edi.flood.bHarmonic && edi.flood.referenceProjectionSlope[i] != 0.0)
 728         {
 729             fprintf(fp, EDcol_efmt, edi.flood.vecs.refproj[i]);
 730         }
 731
 732         /* Output Efl if we are doing adaptive flooding */
 733         if (0 != edi.flood.tau)
 734         {
 735             fprintf(fp, EDcol_efmt, edi.flood.Efl);
 736         }
 737         fprintf(fp, EDcol_efmt, edi.flood.Vfl);
 738
 739         /* Output deltaF if we are doing adaptive flooding */
 740         if (0 != edi.flood.tau)
 741         {
 742             fprintf(fp, EDcol_efmt, edi.flood.deltaF);
 743         }
 744         fprintf(fp, EDcol_efmt, edi.flood.vecs.fproj[i]);
 745     }
 746 }
 747 } // namespace
 748
 749
 750 /* From flood.xproj compute the Vfl(x) at this point */
 751 static real flood_energy(t_edpar* edi, int64_t step)
 752 {
 753     /* compute flooding energy Vfl
 754        Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
 755        \lambda_i is the reciprocal eigenvalue 1/\sigma_i
 756          it is already computed by make_edi and stored in stpsz[i]
 757        bHarmonic:
 758        Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
 759      */
 760     real sum;
 761     real Vfl;
 762     int  i;
 763
 764
 765     /* Each time this routine is called (i.e. each time step), we add a small
 766      * value to the reference projection. This way a harmonic restraint towards
 767      * a moving reference is realized. If no value for the additive constant
 768      * is provided in the edi file, the reference will not change. */
 769     if (edi->flood.bHarmonic)
 770     {
 771         for (i = 0; i < edi->flood.vecs.neig; i++)
 772         {
 773             edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i]
 774                                          + step * edi->flood.referenceProjectionSlope[i];
 775         }
 776     }
 777
 778     sum = 0.0;
 779     /* Compute sum which will be the exponent of the exponential */
 780     for (i = 0; i < edi->flood.vecs.neig; i++)
 781     {
 782         /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
 783         sum += edi->flood.vecs.stpsz[i] * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i])
 784                * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i]);
 785     }
 786
 787     /* Compute the Gauss function*/
 788     if (edi->flood.bHarmonic)
 789     {
 790         Vfl = -0.5 * edi->flood.Efl * sum; /* minus sign because Efl is negative, if restrain is on. */
 791     }
 792     else
 793     {
 794         Vfl = edi->flood.Efl != 0
 795                       ? edi->flood.Efl
 796                                 * std::exp(-edi->flood.kT / 2 / edi->flood.Efl / edi->flood.alpha2 * sum)
 797                       : 0;
 798     }
 799
 800     return Vfl;
 801 }
 802
 803
 804 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
 805 static void flood_forces(t_edpar* edi)
 806 {
 807     /* compute the forces in the subspace of the flooding eigenvectors
 808      * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
 809
 810     int  i;
 811     real energy = edi->flood.Vfl;
 812
 813
 814     if (edi->flood.bHarmonic)
 815     {
 816         for (i = 0; i < edi->flood.vecs.neig; i++)
 817         {
 818             edi->flood.vecs.fproj[i] = edi->flood.Efl * edi->flood.vecs.stpsz[i]
 819                                        * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i]);
 820         }
 821     }
 822     else
 823     {
 824         for (i = 0; i < edi->flood.vecs.neig; i++)
 825         {
 826             /* if Efl is zero the forces are zero if not use the formula */
 827             edi->flood.vecs.fproj[i] =
 828                     edi->flood.Efl != 0
 829                             ? edi->flood.kT / edi->flood.Efl / edi->flood.alpha2 * energy
 830                                       * edi->flood.vecs.stpsz[i]
 831                                       * (edi->flood.vecs.xproj[i] - edi->flood.vecs.refproj[i])
 832                             : 0;
 833         }
 834     }
 835 }
 836
 837 namespace
 838 {
 839 /*!\brief Raise forces from subspace into cartesian space.
 840  * This function lifts the forces from the subspace to the cartesian space
 841  * all the values not contained in the subspace are assumed to be zero and then
 842  * a coordinate transformation from eigenvector to cartesian vectors is performed
 843  * The nonexistent values don't have to be set to zero explicitly, they would occur
 844  * as zero valued summands, hence we just stop to compute this part of the sum.
 845  * For every atom we add all the contributions to this atom from all the different eigenvectors.
 846  * NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
 847  * field forces_cart prior the computation, but we compute the forces separately
 848  * to have them accessible for diagnostics
 849  *
 850  * \param[in] edi Essential dynamics input parameters
 851  * \param[out] forces_cart The cartesian forces
 852  */
 853
 854 void flood_blowup(const t_edpar& edi, rvec* forces_cart)
 855 {
 856     const real* forces_sub = edi.flood.vecs.fproj;
 857     /* Calculate the cartesian forces for the local atoms */
 858
 859     /* Clear forces first */
 860     for (int j = 0; j < edi.sav.nr_loc; j++)
 861     {
 862         clear_rvec(forces_cart[j]);
 863     }
 864
 865     /* Now compute atomwise */
 866     for (int j = 0; j < edi.sav.nr_loc; j++)
 867     {
 868         /* Compute forces_cart[edi.sav.anrs[j]] */
 869         for (int eig = 0; eig < edi.flood.vecs.neig; eig++)
 870         {
 871             rvec addedForce;
 872             /* Force vector is force * eigenvector (compute only atom j) */
 873             svmul(forces_sub[eig], edi.flood.vecs.vec[eig][edi.sav.c_ind[j]], addedForce);
 874             /* Add this vector to the cartesian forces */
 875             rvec_inc(forces_cart[j], addedForce);
 876         }
 877     }
 878 }
 879
 880 } // namespace
 881
 882
 883 /* Update the values of Efl, deltaF depending on tau and Vfl */
 884 static void update_adaption(t_edpar* edi)
 885 {
 886     /* this function updates the parameter Efl and deltaF according to the rules given in
 887      * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
 888      * J. chem Phys. */
 889
 890     if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau) > 0.00000001)
 891     {
 892         edi->flood.Efl = edi->flood.Efl
 893                          + edi->flood.dt / edi->flood.tau * (edi->flood.deltaF0 - edi->flood.deltaF);
 894         /* check if restrain (inverted flooding) -> don't let EFL become positive */
 895         if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001)
 896         {
 897             edi->flood.Efl = 0;
 898         }
 899
 900         edi->flood.deltaF = (1 - edi->flood.dt / edi->flood.tau) * edi->flood.deltaF
 901                             + edi->flood.dt / edi->flood.tau * edi->flood.Vfl;
 902     }
 903 }
 904
 905
 906 static void do_single_flood(FILE*            edo,
 907                             const rvec       x[],
 908                             rvec             force[],
 909                             t_edpar*         edi,
 910                             int64_t          step,
 911                             const matrix     box,
 912                             const t_commrec* cr,
 913                             gmx_bool         bNS) /* Are we in a neighbor searching step? */
 914 {
 915     int                i;
 916     matrix             rotmat;   /* rotation matrix */
 917     matrix             tmat;     /* inverse rotation */
 918     rvec               transvec; /* translation vector */
 919     real               rmsdev;
 920     struct t_do_edsam* buf;
 921
 922
 923     buf = edi->buf->do_edsam;
 924
 925
 926     /* Broadcast the positions of the AVERAGE structure such that they are known on
 927      * every processor. Each node contributes its local positions x and stores them in
 928      * the collective ED array buf->xcoll */
 929     communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
 930                                 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind,
 931                                 edi->sav.x_old, box);
 932
 933     /* Only assembly REFERENCE positions if their indices differ from the average ones */
 934     if (!edi->bRefEqAv)
 935     {
 936         communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref,
 937                                     bNS, x, edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc,
 938                                     edi->sref.c_ind, edi->sref.x_old, box);
 939     }
 940
 941     /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
 942      * We do not need to update the shifts until the next NS step */
 943     buf->bUpdateShifts = FALSE;
 944
 945     /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
 946      * as well as the indices in edi->sav.anrs */
 947
 948     /* Fit the reference indices to the reference structure */
 949     if (edi->bRefEqAv)
 950     {
 951         fit_to_reference(buf->xcoll, transvec, rotmat, edi);
 952     }
 953     else
 954     {
 955         fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
 956     }
 957
 958     /* Now apply the translation and rotation to the ED structure */
 959     translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
 960
 961     /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
 962     project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
 963
 964     if (!edi->flood.bConstForce)
 965     {
 966         /* Compute Vfl(x) from flood.xproj */
 967         edi->flood.Vfl = flood_energy(edi, step);
 968
 969         update_adaption(edi);
 970
 971         /* Compute the flooding forces */
 972         flood_forces(edi);
 973     }
 974
 975     /* Translate them into cartesian positions */
 976     flood_blowup(*edi, edi->flood.forces_cartesian);
 977
 978     /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
 979     /* Each node rotates back its local forces */
 980     transpose(rotmat, tmat);
 981     rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
 982
 983     /* Finally add forces to the main force variable */
 984     for (i = 0; i < edi->sav.nr_loc; i++)
 985     {
 986         rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
 987     }
 988
 989     /* Output is written by the master process */
 990     if (do_per_step(step, edi->outfrq) && MASTER(cr))
 991     {
 992         /* Output how well we fit to the reference */
 993         if (edi->bRefEqAv)
 994         {
 995             /* Indices of reference and average structures are identical,
 996              * thus we can calculate the rmsd to SREF using xcoll */
 997             rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
 998         }
 999         else
1000         {
1001             /* We have to translate & rotate the reference atoms first */
1002             translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1003             rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1004         }
1005
1006         write_edo_flood(*edi, edo, rmsdev);
1007     }
1008 }
1009
1010
1011 /* Main flooding routine, called from do_force */
1012 extern void do_flood(const t_commrec*  cr,
1013                      const t_inputrec* ir,
1014                      const rvec        x[],
1015                      rvec              force[],
1016                      gmx_edsam*        ed,
1017                      const matrix      box,
1018                      int64_t           step,
1019                      gmx_bool          bNS)
1020 {
1021     /* Write time to edo, when required. Output the time anyhow since we need
1022      * it in the output file for ED constraints. */
1023     if (MASTER(cr) && do_per_step(step, ed->edpar.begin()->outfrq))
1024     {
1025         fprintf(ed->edo, "\n%12f", ir->init_t + step * ir->delta_t);
1026     }
1027
1028     if (ed->eEDtype != EssentialDynamicsType::Flooding)
1029     {
1030         return;
1031     }
1032
1033     for (auto& edi : ed->edpar)
1034     {
1035         /* Call flooding for one matrix */
1036         if (edi.flood.vecs.neig)
1037         {
1038             do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS);
1039         }
1040     }
1041 }
1042
1043
1044 /* Called by init_edi, configure some flooding related variables and structures,
1045  * print headers to output files */
1046 static void init_flood(t_edpar* edi, gmx_edsam* ed, real dt)
1047 {
1048     int i;
1049
1050
1051     edi->flood.Efl = edi->flood.constEfl;
1052     edi->flood.Vfl = 0;
1053     edi->flood.dt  = dt;
1054
1055     if (edi->flood.vecs.neig)
1056     {
1057         /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1058         ed->eEDtype = EssentialDynamicsType::Flooding;
1059
1060         fprintf(stderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig,
1061                 edi->flood.vecs.neig > 1 ? "s" : "");
1062
1063         if (edi->flood.bConstForce)
1064         {
1065             /* We have used stpsz as a vehicle to carry the fproj values for constant
1066              * force flooding. Now we copy that to flood.vecs.fproj. Note that
1067              * in const force flooding, fproj is never changed. */
1068             for (i = 0; i < edi->flood.vecs.neig; i++)
1069             {
1070                 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1071
1072                 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1073                         edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1074             }
1075         }
1076     }
1077 }
1078
1079
1080 #ifdef DEBUGHELPERS
1081 /*********** Energy book keeping ******/
1082 static void get_flood_enx_names(t_edpar* edi, char** names, int* nnames) /* get header of energies */
1083 {
1084     t_edpar* actual;
1085     int      count;
1086     char     buf[STRLEN];
1087     actual = edi;
1088     count  = 1;
1089     while (actual)
1090     {
1091         srenew(names, count);
1092         sprintf(buf, "Vfl_%d", count);
1093         names[count - 1] = gmx_strdup(buf);
1094         actual           = actual->next_edi;
1095         count++;
1096     }
1097     *nnames = count - 1;
1098 }
1099
1100
1101 static void get_flood_energies(t_edpar* edi, real Vfl[], int nnames)
1102 {
1103     /*fl has to be big enough to capture nnames-many entries*/
1104     t_edpar* actual;
1105     int      count;
1106
1107
1108     actual = edi;
1109     count  = 1;
1110     while (actual)
1111     {
1112         Vfl[count - 1] = actual->flood.Vfl;
1113         actual         = actual->next_edi;
1114         count++;
1115     }
1116     if (nnames != count - 1)
1117     {
1118         gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1119     }
1120 }
1121 /************* END of FLOODING IMPLEMENTATION ****************************/
1122 #endif
1123
1124
1125 /* This function opens the ED input and output files, reads in all datasets it finds
1126  * in the input file, and cross-checks whether the .edi file information is consistent
1127  * with the essential dynamics data found in the checkpoint file (if present).
1128  * gmx make_edi can be used to create an .edi input file.
1129  */
1130 static std::unique_ptr<gmx::EssentialDynamics> ed_open(int                         natoms,
1131                                                        ObservablesHistory*         oh,
1132                                                        const char*                 ediFileName,
1133                                                        const char*                 edoFileName,
1134                                                        const gmx::StartingBehavior startingBehavior,
1135                                                        const gmx_output_env_t*     oenv,
1136                                                        const t_commrec*            cr)
1137 {
1138     auto edHandle = std::make_unique<gmx::EssentialDynamics>();
1139     auto ed       = edHandle->getLegacyED();
1140     /* We want to perform ED (this switch might later be upgraded to EssentialDynamicsType::Flooding) */
1141     ed->eEDtype = EssentialDynamicsType::EDSampling;
1142
1143     if (MASTER(cr))
1144     {
1145         // If we start from a checkpoint file, we already have an edsamHistory struct
1146         if (oh->edsamHistory == nullptr)
1147         {
1148             oh->edsamHistory = std::make_unique<edsamhistory_t>(edsamhistory_t{});
1149         }
1150         edsamhistory_t* EDstate = oh->edsamHistory.get();
1151
1152         /* Read the edi input file: */
1153         ed->edpar = read_edi_file(ediFileName, natoms);
1154
1155         /* Make sure the checkpoint was produced in a run using this .edi file */
1156         if (EDstate->bFromCpt)
1157         {
1158             crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1159         }
1160         else
1161         {
1162             EDstate->nED = ed->edpar.size();
1163         }
1164         init_edsamstate(*ed, EDstate);
1165
1166         /* The master opens the ED output file */
1167         if (startingBehavior == gmx::StartingBehavior::RestartWithAppending)
1168         {
1169             ed->edo = gmx_fio_fopen(edoFileName, "a+");
1170         }
1171         else
1172         {
1173             ed->edo = xvgropen(edoFileName, "Essential dynamics / flooding output", "Time (ps)",
1174                                "RMSDs (nm), projections on EVs (nm), ...", oenv);
1175
1176             /* Make a descriptive legend */
1177             write_edo_legend(ed, EDstate->nED, oenv);
1178         }
1179     }
1180     return edHandle;
1181 }
1182
1183
1184 /* Broadcasts the structure data */
1185 static void bc_ed_positions(const t_commrec* cr, struct gmx_edx* s, EssentialDynamicsStructure stype)
1186 {
1187     snew_bc(cr, s->anrs, s->nr); /* Index numbers     */
1188     snew_bc(cr, s->x, s->nr);    /* Positions         */
1189     nblock_bc(cr, s->nr, s->anrs);
1190     nblock_bc(cr, s->nr, s->x);
1191
1192     /* For the average & reference structures we need an array for the collective indices,
1193      * and we need to broadcast the masses as well */
1194     if (stype == EssentialDynamicsStructure::Average || stype == EssentialDynamicsStructure::Reference)
1195     {
1196         /* We need these additional variables in the parallel case: */
1197         snew(s->c_ind, s->nr); /* Collective indices */
1198         /* Local atom indices get assigned in dd_make_local_group_indices.
1199          * There, also memory is allocated */
1200         s->nalloc_loc = 0;            /* allocation size of s->anrs_loc */
1201         snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1202         nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1203     }
1204
1205     /* broadcast masses for the reference structure (for mass-weighted fitting) */
1206     if (stype == EssentialDynamicsStructure::Reference)
1207     {
1208         snew_bc(cr, s->m, s->nr);
1209         nblock_bc(cr, s->nr, s->m);
1210     }
1211
1212     /* For the average structure we might need the masses for mass-weighting */
1213     if (stype == EssentialDynamicsStructure::Average)
1214     {
1215         snew_bc(cr, s->sqrtm, s->nr);
1216         nblock_bc(cr, s->nr, s->sqrtm);
1217         snew_bc(cr, s->m, s->nr);
1218         nblock_bc(cr, s->nr, s->m);
1219     }
1220 }
1221
1222
1223 /* Broadcasts the eigenvector data */
1224 static void bc_ed_vecs(const t_commrec* cr, t_eigvec* ev, int length)
1225 {
1226     int i;
1227
1228     snew_bc(cr, ev->ieig, ev->neig);    /* index numbers of eigenvector  */
1229     snew_bc(cr, ev->stpsz, ev->neig);   /* stepsizes per eigenvector     */
1230     snew_bc(cr, ev->xproj, ev->neig);   /* instantaneous x projection    */
1231     snew_bc(cr, ev->fproj, ev->neig);   /* instantaneous f projection    */
1232     snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1233
1234     nblock_bc(cr, ev->neig, ev->ieig);
1235     nblock_bc(cr, ev->neig, ev->stpsz);
1236     nblock_bc(cr, ev->neig, ev->xproj);
1237     nblock_bc(cr, ev->neig, ev->fproj);
1238     nblock_bc(cr, ev->neig, ev->refproj);
1239
1240     snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components        */
1241     for (i = 0; i < ev->neig; i++)
1242     {
1243         snew_bc(cr, ev->vec[i], length);
1244         nblock_bc(cr, length, ev->vec[i]);
1245     }
1246 }
1247
1248
1249 /* Broadcasts the ED / flooding data to other nodes
1250  * and allocates memory where needed */
1251 static void broadcast_ed_data(const t_commrec* cr, gmx_edsam* ed)
1252 {
1253     /* Master lets the other nodes know if its ED only or also flooding */
1254     gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1255
1256     int numedis = ed->edpar.size();
1257     /* First let everybody know how many ED data sets to expect */
1258     gmx_bcast(sizeof(numedis), &numedis, cr);
1259     nblock_abc(cr, numedis, &(ed->edpar));
1260     /* Now transfer the ED data set(s) */
1261     for (auto& edi : ed->edpar)
1262     {
1263         /* Broadcast a single ED data set */
1264         block_bc(cr, edi);
1265
1266         /* Broadcast positions */
1267         bc_ed_positions(cr, &(edi.sref),
1268                         EssentialDynamicsStructure::Reference); /* reference positions (don't broadcast masses)    */
1269         bc_ed_positions(cr, &(edi.sav),
1270                         EssentialDynamicsStructure::Average); /* average positions (do broadcast masses as well) */
1271         bc_ed_positions(cr, &(edi.star), EssentialDynamicsStructure::Target); /* target positions */
1272         bc_ed_positions(cr, &(edi.sori), EssentialDynamicsStructure::Origin); /* origin positions */
1273
1274         /* Broadcast eigenvectors */
1275         bc_ed_vecs(cr, &edi.vecs.mon, edi.sav.nr);
1276         bc_ed_vecs(cr, &edi.vecs.linfix, edi.sav.nr);
1277         bc_ed_vecs(cr, &edi.vecs.linacc, edi.sav.nr);
1278         bc_ed_vecs(cr, &edi.vecs.radfix, edi.sav.nr);
1279         bc_ed_vecs(cr, &edi.vecs.radacc, edi.sav.nr);
1280         bc_ed_vecs(cr, &edi.vecs.radcon, edi.sav.nr);
1281         /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1282         bc_ed_vecs(cr, &edi.flood.vecs, edi.sav.nr);
1283
1284         /* For harmonic restraints the reference projections can change with time */
1285         if (edi.flood.bHarmonic)
1286         {
1287             snew_bc(cr, edi.flood.initialReferenceProjection, edi.flood.vecs.neig);
1288             snew_bc(cr, edi.flood.referenceProjectionSlope, edi.flood.vecs.neig);
1289             nblock_bc(cr, edi.flood.vecs.neig, edi.flood.initialReferenceProjection);
1290             nblock_bc(cr, edi.flood.vecs.neig, edi.flood.referenceProjectionSlope);
1291         }
1292     }
1293 }
1294
1295
1296 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1297 static void init_edi(const gmx_mtop_t* mtop, t_edpar* edi)
1298 {
1299     int  i;
1300     real totalmass = 0.0;
1301     rvec com;
1302
1303     /* NOTE Init_edi is executed on the master process only
1304      * The initialized data sets are then transmitted to the
1305      * other nodes in broadcast_ed_data */
1306
1307     /* evaluate masses (reference structure) */
1308     snew(edi->sref.m, edi->sref.nr);
1309     int molb = 0;
1310     for (i = 0; i < edi->sref.nr; i++)
1311     {
1312         if (edi->fitmas)
1313         {
1314             edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
1315         }
1316         else
1317         {
1318             edi->sref.m[i] = 1.0;
1319         }
1320
1321         /* Check that every m > 0. Bad things will happen otherwise. */
1322         if (edi->sref.m[i] <= 0.0)
1323         {
1324             gmx_fatal(FARGS,
1325                       "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1326                       "For a mass-weighted fit, all reference structure atoms need to have a mass "
1327                       ">0.\n"
1328                       "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1329                       "atoms from the reference structure by creating a proper index group.\n",
1330                       i, edi->sref.anrs[i] + 1, edi->sref.m[i]);
1331         }
1332
1333         totalmass += edi->sref.m[i];
1334     }
1335     edi->sref.mtot = totalmass;
1336
1337     /* Masses m and sqrt(m) for the average structure. Note that m
1338      * is needed if forces have to be evaluated in do_edsam */
1339     snew(edi->sav.sqrtm, edi->sav.nr);
1340     snew(edi->sav.m, edi->sav.nr);
1341     for (i = 0; i < edi->sav.nr; i++)
1342     {
1343         edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
1344         if (edi->pcamas)
1345         {
1346             edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1347         }
1348         else
1349         {
1350             edi->sav.sqrtm[i] = 1.0;
1351         }
1352
1353         /* Check that every m > 0. Bad things will happen otherwise. */
1354         if (edi->sav.sqrtm[i] <= 0.0)
1355         {
1356             gmx_fatal(FARGS,
1357                       "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1358                       "For ED with mass-weighting, all average structure atoms need to have a mass "
1359                       ">0.\n"
1360                       "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1361                       "atoms from the average structure by creating a proper index group.\n",
1362                       i, edi->sav.anrs[i] + 1, edi->sav.m[i]);
1363         }
1364     }
1365
1366     /* put reference structure in origin */
1367     get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1368     com[XX] = -com[XX];
1369     com[YY] = -com[YY];
1370     com[ZZ] = -com[ZZ];
1371     translate_x(edi->sref.x, edi->sref.nr, com);
1372
1373     /* Init ED buffer */
1374     snew(edi->buf, 1);
1375 }
1376
1377
1378 static void check(const char* line, const char* label)
1379 {
1380     if (!strstr(line, label))
1381     {
1382         gmx_fatal(FARGS,
1383                   "Could not find input parameter %s at expected position in edsam input-file "
1384                   "(.edi)\nline read instead is %s",
1385                   label, line);
1386     }
1387 }
1388
1389
1390 static int read_checked_edint(FILE* file, const char* label)
1391 {
1392     char line[STRLEN + 1];
1393     int  idum;
1394
1395     fgets2(line, STRLEN, file);
1396     check(line, label);
1397     fgets2(line, STRLEN, file);
1398     sscanf(line, max_ev_fmt_d, &idum);
1399     return idum;
1400 }
1401
1402 static bool read_checked_edbool(FILE* file, const char* label)
1403 {
1404     return read_checked_edint(file, label) != 0;
1405 }
1406
1407 static int read_edint(FILE* file, bool* bEOF)
1408 {
1409     char  line[STRLEN + 1];
1410     int   idum;
1411     char* eof;
1412
1413
1414     eof = fgets2(line, STRLEN, file);
1415     if (eof == nullptr)
1416     {
1417         *bEOF = TRUE;
1418         return -1;
1419     }
1420     eof = fgets2(line, STRLEN, file);
1421     if (eof == nullptr)
1422     {
1423         *bEOF = TRUE;
1424         return -1;
1425     }
1426     sscanf(line, max_ev_fmt_d, &idum);
1427     *bEOF = FALSE;
1428     return idum;
1429 }
1430
1431
1432 static real read_checked_edreal(FILE* file, const char* label)
1433 {
1434     char   line[STRLEN + 1];
1435     double rdum;
1436
1437
1438     fgets2(line, STRLEN, file);
1439     check(line, label);
1440     fgets2(line, STRLEN, file);
1441     sscanf(line, max_ev_fmt_lf, &rdum);
1442     return static_cast<real>(rdum); /* always read as double and convert to single */
1443 }
1444
1445
1446 static void read_edx(FILE* file, int number, int* anrs, rvec* x)
1447 {
1448     int    i, j;
1449     char   line[STRLEN + 1];
1450     double d[3];
1451
1452
1453     for (i = 0; i < number; i++)
1454     {
1455         fgets2(line, STRLEN, file);
1456         sscanf(line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1457         anrs[i]--; /* we are reading FORTRAN indices */
1458         for (j = 0; j < 3; j++)
1459         {
1460             x[i][j] = d[j]; /* always read as double and convert to single */
1461         }
1462     }
1463 }
1464
1465 namespace
1466 {
1467 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1468  * \param[in] in the file to read from
1469  * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1470  * \param[out] vec the eigen-vectors
1471  * \param[in] nEig the number of eigenvectors
1472  */
1473 void scan_edvec(FILE* in, int numAtoms, rvec*** vec, int nEig)
1474 {
1475     snew(*vec, nEig);
1476     for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1477     {
1478         snew((*vec)[iEigenvector], numAtoms);
1479         for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1480         {
1481             char   line[STRLEN + 1];
1482             double x, y, z;
1483             fgets2(line, STRLEN, in);
1484             sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1485             (*vec)[iEigenvector][iAtom][XX] = x;
1486             (*vec)[iEigenvector][iAtom][YY] = y;
1487             (*vec)[iEigenvector][iAtom][ZZ] = z;
1488         }
1489     }
1490 }
1491 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1492  * \param[in] in input file
1493  * \param[in] nrAtoms number of atoms/coordinates
1494  * \param[out] tvec the eigenvector
1495  */
1496 void read_edvec(FILE* in, int nrAtoms, t_eigvec* tvec)
1497 {
1498     tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1499     if (tvec->neig <= 0)
1500     {
1501         return;
1502     }
1503
1504     snew(tvec->ieig, tvec->neig);
1505     snew(tvec->stpsz, tvec->neig);
1506     for (int i = 0; i < tvec->neig; i++)
1507     {
1508         char line[STRLEN + 1];
1509         fgets2(line, STRLEN, in);
1510         int       numEigenVectors;
1511         double    stepSize;
1512         const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1513         if (nscan != 2)
1514         {
1515             gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1516         }
1517         tvec->ieig[i]  = numEigenVectors;
1518         tvec->stpsz[i] = stepSize;
1519     } /* end of loop over eigenvectors */
1520
1521     scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1522 }
1523 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1524  *
1525  * Eigenvector and its intial reference and slope are stored on the
1526  * same line in the .edi format. To avoid re-winding the .edi file,
1527  * the reference eigen-vector and reference data are read in one go.
1528  *
1529  * \param[in] in input file
1530  * \param[in] nrAtoms number of atoms/coordinates
1531  * \param[out] tvec the eigenvector
1532  * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1533  * \param[out] referenceProjectionSlope The slope of the reference projections.
1534  */
1535 bool readEdVecWithReferenceProjection(FILE*     in,
1536                                       int       nrAtoms,
1537                                       t_eigvec* tvec,
1538                                       real**    initialReferenceProjection,
1539                                       real**    referenceProjectionSlope)
1540 {
1541     bool providesReference = false;
1542     tvec->neig             = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1543     if (tvec->neig <= 0)
1544     {
1545         return false;
1546     }
1547
1548     snew(tvec->ieig, tvec->neig);
1549     snew(tvec->stpsz, tvec->neig);
1550     snew(*initialReferenceProjection, tvec->neig);
1551     snew(*referenceProjectionSlope, tvec->neig);
1552     for (int i = 0; i < tvec->neig; i++)
1553     {
1554         char line[STRLEN + 1];
1555         fgets2(line, STRLEN, in);
1556         int       numEigenVectors;
1557         double    stepSize            = 0.;
1558         double    referenceProjection = 0.;
1559         double    referenceSlope      = 0.;
1560         const int nscan = sscanf(line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize,
1561                                  &referenceProjection, &referenceSlope);
1562         /* Zero out values which were not scanned */
1563         switch (nscan)
1564         {
1565             case 4:
1566                 /* Every 4 values read, including reference position */
1567                 providesReference = true;
1568                 break;
1569             case 3:
1570                 /* A reference position is provided */
1571                 providesReference = true;
1572                 /* No value for slope, set to 0 */
1573                 referenceSlope = 0.0;
1574                 break;
1575             case 2:
1576                 /* No values for reference projection and slope, set to 0 */
1577                 referenceProjection = 0.0;
1578                 referenceSlope      = 0.0;
1579                 break;
1580             default:
1581                 gmx_fatal(FARGS,
1582                           "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> "
1583                           "<refproj> <refproj-slope>\n",
1584                           nscan);
1585         }
1586         (*initialReferenceProjection)[i] = referenceProjection;
1587         (*referenceProjectionSlope)[i]   = referenceSlope;
1588
1589         tvec->ieig[i]  = numEigenVectors;
1590         tvec->stpsz[i] = stepSize;
1591     } /* end of loop over eigenvectors */
1592
1593     scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1594     return providesReference;
1595 }
1596
1597 /*!\brief Allocate working memory for the eigenvectors.
1598  * \param[in,out] tvec eigenvector for which memory will be allocated
1599  */
1600 void setup_edvec(t_eigvec* tvec)
1601 {
1602     snew(tvec->xproj, tvec->neig);
1603     snew(tvec->fproj, tvec->neig);
1604     snew(tvec->refproj, tvec->neig);
1605 }
1606 } // namespace
1607
1608
1609 /* Check if the same atom indices are used for reference and average positions */
1610 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1611 {
1612     int i;
1613
1614
1615     /* If the number of atoms differs between the two structures,
1616      * they cannot be identical */
1617     if (sref.nr != sav.nr)
1618     {
1619         return FALSE;
1620     }
1621
1622     /* Now that we know that both stuctures have the same number of atoms,
1623      * check if also the indices are identical */
1624     for (i = 0; i < sav.nr; i++)
1625     {
1626         if (sref.anrs[i] != sav.anrs[i])
1627         {
1628             return FALSE;
1629         }
1630     }
1631     fprintf(stderr,
1632             "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1633
1634     return TRUE;
1635 }
1636
1637 namespace
1638 {
1639
1640 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1641 constexpr int c_oldestUnsupportedMagicNumber = 666;
1642 //! Oldest (lowest) magic number indicating supported essential dynamics input
1643 constexpr int c_oldestSupportedMagicNumber = 669;
1644 //! Latest (highest) magic number indicating supported essential dynamics input
1645 constexpr int c_latestSupportedMagicNumber = 670;
1646
1647 /*!\brief Set up essential dynamics work parameters.
1648  * \param[in] edi Essential dynamics working structure.
1649  */
1650 void setup_edi(t_edpar* edi)
1651 {
1652     snew(edi->sref.x_old, edi->sref.nr);
1653     edi->sref.sqrtm = nullptr;
1654     snew(edi->sav.x_old, edi->sav.nr);
1655     if (edi->star.nr > 0)
1656     {
1657         edi->star.sqrtm = nullptr;
1658     }
1659     if (edi->sori.nr > 0)
1660     {
1661         edi->sori.sqrtm = nullptr;
1662     }
1663     setup_edvec(&edi->vecs.linacc);
1664     setup_edvec(&edi->vecs.mon);
1665     setup_edvec(&edi->vecs.linfix);
1666     setup_edvec(&edi->vecs.linacc);
1667     setup_edvec(&edi->vecs.radfix);
1668     setup_edvec(&edi->vecs.radacc);
1669     setup_edvec(&edi->vecs.radcon);
1670     setup_edvec(&edi->flood.vecs);
1671 }
1672
1673 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1674  * \param[in] magicNumber indicates the essential dynamics input file version
1675  * \returns true if CONST_FORCE_FLOODING is supported
1676  */
1677 bool ediFileHasConstForceFlooding(int magicNumber)
1678 {
1679     return magicNumber > c_oldestSupportedMagicNumber;
1680 };
1681
1682 /*!\brief Reports reading success of the essential dynamics input file magic number.
1683  * \param[in] in pointer to input file
1684  * \param[out] magicNumber determines the edi file version
1685  * \returns true if setting file version from input file was successful.
1686  */
1687 bool ediFileHasValidDataSet(FILE* in, int* magicNumber)
1688 {
1689     /* the edi file is not free format, so expect problems if the input is corrupt. */
1690     bool reachedEndOfFile = true;
1691     /* check the magic number */
1692     *magicNumber = read_edint(in, &reachedEndOfFile);
1693     /* Check whether we have reached the end of the input file */
1694     return !reachedEndOfFile;
1695 }
1696
1697 /*!\brief Trigger fatal error if magic number is unsupported.
1698  * \param[in] magicNumber A magic number read from the edi file
1699  * \param[in] fn name of input file for error reporting
1700  */
1701 void exitWhenMagicNumberIsInvalid(int magicNumber, const char* fn)
1702 {
1703
1704     if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1705     {
1706         if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1707         {
1708             gmx_fatal(FARGS,
1709                       "Wrong magic number: Use newest version of make_edi to produce edi file");
1710         }
1711         else
1712         {
1713             gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1714         }
1715     }
1716 }
1717
1718 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1719  *
1720  * \param[in] in input file
1721  * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1722  * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1723  * \param[in] fn the file name of the input file for error reporting
1724  * \returns edi essential dynamics parameters
1725  */
1726 t_edpar read_edi(FILE* in, int nr_mdatoms, bool hasConstForceFlooding, const char* fn)
1727 {
1728     t_edpar edi;
1729     bool    bEOF;
1730     /* check the number of atoms */
1731     edi.nini = read_edint(in, &bEOF);
1732     if (edi.nini != nr_mdatoms)
1733     {
1734         gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi.nini,
1735                   nr_mdatoms);
1736     }
1737
1738     /* Done checking. For the rest we blindly trust the input */
1739     edi.fitmas     = read_checked_edbool(in, "FITMAS");
1740     edi.pcamas     = read_checked_edbool(in, "ANALYSIS_MAS");
1741     edi.outfrq     = read_checked_edint(in, "OUTFRQ");
1742     edi.maxedsteps = read_checked_edint(in, "MAXLEN");
1743     edi.slope      = read_checked_edreal(in, "SLOPECRIT");
1744
1745     edi.presteps        = read_checked_edint(in, "PRESTEPS");
1746     edi.flood.deltaF0   = read_checked_edreal(in, "DELTA_F0");
1747     edi.flood.deltaF    = read_checked_edreal(in, "INIT_DELTA_F");
1748     edi.flood.tau       = read_checked_edreal(in, "TAU");
1749     edi.flood.constEfl  = read_checked_edreal(in, "EFL_NULL");
1750     edi.flood.alpha2    = read_checked_edreal(in, "ALPHA2");
1751     edi.flood.kT        = read_checked_edreal(in, "KT");
1752     edi.flood.bHarmonic = read_checked_edbool(in, "HARMONIC");
1753     if (hasConstForceFlooding)
1754     {
1755         edi.flood.bConstForce = read_checked_edbool(in, "CONST_FORCE_FLOODING");
1756     }
1757     else
1758     {
1759         edi.flood.bConstForce = FALSE;
1760     }
1761     edi.sref.nr = read_checked_edint(in, "NREF");
1762
1763     /* allocate space for reference positions and read them */
1764     snew(edi.sref.anrs, edi.sref.nr);
1765     snew(edi.sref.x, edi.sref.nr);
1766     read_edx(in, edi.sref.nr, edi.sref.anrs, edi.sref.x);
1767
1768     /* average positions. they define which atoms will be used for ED sampling */
1769     edi.sav.nr = read_checked_edint(in, "NAV");
1770     snew(edi.sav.anrs, edi.sav.nr);
1771     snew(edi.sav.x, edi.sav.nr);
1772     read_edx(in, edi.sav.nr, edi.sav.anrs, edi.sav.x);
1773
1774     /* Check if the same atom indices are used for reference and average positions */
1775     edi.bRefEqAv = check_if_same(edi.sref, edi.sav);
1776
1777     /* eigenvectors */
1778
1779     read_edvec(in, edi.sav.nr, &edi.vecs.mon);
1780     read_edvec(in, edi.sav.nr, &edi.vecs.linfix);
1781     read_edvec(in, edi.sav.nr, &edi.vecs.linacc);
1782     read_edvec(in, edi.sav.nr, &edi.vecs.radfix);
1783     read_edvec(in, edi.sav.nr, &edi.vecs.radacc);
1784     read_edvec(in, edi.sav.nr, &edi.vecs.radcon);
1785
1786     /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1787     bool bHaveReference = false;
1788     if (edi.flood.bHarmonic)
1789     {
1790         bHaveReference = readEdVecWithReferenceProjection(in, edi.sav.nr, &edi.flood.vecs,
1791                                                           &edi.flood.initialReferenceProjection,
1792                                                           &edi.flood.referenceProjectionSlope);
1793     }
1794     else
1795     {
1796         read_edvec(in, edi.sav.nr, &edi.flood.vecs);
1797     }
1798
1799     /* target positions */
1800     edi.star.nr = read_edint(in, &bEOF);
1801     if (edi.star.nr > 0)
1802     {
1803         snew(edi.star.anrs, edi.star.nr);
1804         snew(edi.star.x, edi.star.nr);
1805         read_edx(in, edi.star.nr, edi.star.anrs, edi.star.x);
1806     }
1807
1808     /* positions defining origin of expansion circle */
1809     edi.sori.nr = read_edint(in, &bEOF);
1810     if (edi.sori.nr > 0)
1811     {
1812         if (bHaveReference)
1813         {
1814             /* Both an -ori structure and a at least one manual reference point have been
1815              * specified. That's ambiguous and probably not intentional. */
1816             gmx_fatal(FARGS,
1817                       "ED: An origin structure has been provided and a at least one (moving) "
1818                       "reference\n"
1819                       "    point was manually specified in the edi file. That is ambiguous. "
1820                       "Aborting.\n");
1821         }
1822         snew(edi.sori.anrs, edi.sori.nr);
1823         snew(edi.sori.x, edi.sori.nr);
1824         read_edx(in, edi.sori.nr, edi.sori.anrs, edi.sori.x);
1825     }
1826     return edi;
1827 }
1828
1829 std::vector<t_edpar> read_edi_file(const char* fn, int nr_mdatoms)
1830 {
1831     std::vector<t_edpar> essentialDynamicsParameters;
1832     FILE*                in;
1833     /* This routine is executed on the master only */
1834
1835     /* Open the .edi parameter input file */
1836     in = gmx_fio_fopen(fn, "r");
1837     fprintf(stderr, "ED: Reading edi file %s\n", fn);
1838
1839     /* Now read a sequence of ED input parameter sets from the edi file */
1840     int ediFileMagicNumber;
1841     while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1842     {
1843         exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1844         bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1845         auto edi                   = read_edi(in, nr_mdatoms, hasConstForceFlooding, fn);
1846         setup_edi(&edi);
1847         essentialDynamicsParameters.emplace_back(edi);
1848
1849         /* Since we arrived within this while loop we know that there is still another data set to be read in */
1850         /* We need to allocate space for the data: */
1851     }
1852     if (essentialDynamicsParameters.empty())
1853     {
1854         gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1855     }
1856
1857     fprintf(stderr, "ED: Found %zu ED group%s.\n", essentialDynamicsParameters.size(),
1858             essentialDynamicsParameters.size() > 1 ? "s" : "");
1859
1860     /* Close the .edi file again */
1861     gmx_fio_fclose(in);
1862
1863     return essentialDynamicsParameters;
1864 }
1865 } // anonymous namespace
1866
1867
1868 struct t_fit_to_ref
1869 {
1870     rvec* xcopy; /* Working copy of the positions in fit_to_reference */
1871 };
1872
1873 /* Fit the current positions to the reference positions
1874  * Do not actually do the fit, just return rotation and translation.
1875  * Note that the COM of the reference structure was already put into
1876  * the origin by init_edi. */
1877 static void fit_to_reference(rvec*    xcoll,    /* The positions to be fitted */
1878                              rvec     transvec, /* The translation vector */
1879                              matrix   rotmat,   /* The rotation matrix */
1880                              t_edpar* edi)      /* Just needed for do_edfit */
1881 {
1882     rvec                 com; /* center of mass */
1883     int                  i;
1884     struct t_fit_to_ref* loc;
1885
1886
1887     /* Allocate memory the first time this routine is called for each edi group */
1888     if (nullptr == edi->buf->fit_to_ref)
1889     {
1890         snew(edi->buf->fit_to_ref, 1);
1891         snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1892     }
1893     loc = edi->buf->fit_to_ref;
1894
1895     /* We do not touch the original positions but work on a copy. */
1896     for (i = 0; i < edi->sref.nr; i++)
1897     {
1898         copy_rvec(xcoll[i], loc->xcopy[i]);
1899     }
1900
1901     /* Calculate the center of mass */
1902     get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1903
1904     transvec[XX] = -com[XX];
1905     transvec[YY] = -com[YY];
1906     transvec[ZZ] = -com[ZZ];
1907
1908     /* Subtract the center of mass from the copy */
1909     translate_x(loc->xcopy, edi->sref.nr, transvec);
1910
1911     /* Determine the rotation matrix */
1912     do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1913 }
1914
1915
1916 static void translate_and_rotate(rvec*  x,        /* The positions to be translated and rotated */
1917                                  int    nat,      /* How many positions are there? */
1918                                  rvec   transvec, /* The translation vector */
1919                                  matrix rotmat)   /* The rotation matrix */
1920 {
1921     /* Translation */
1922     translate_x(x, nat, transvec);
1923
1924     /* Rotation */
1925     rotate_x(x, nat, rotmat);
1926 }
1927
1928
1929 /* Gets the rms deviation of the positions to the structure s */
1930 /* fit_to_structure has to be called before calling this routine! */
1931 static real rmsd_from_structure(rvec* x,           /* The positions under consideration */
1932                                 struct gmx_edx* s) /* The structure from which the rmsd shall be computed */
1933 {
1934     real rmsd = 0.0;
1935     int  i;
1936
1937
1938     for (i = 0; i < s->nr; i++)
1939     {
1940         rmsd += distance2(s->x[i], x[i]);
1941     }
1942
1943     rmsd /= static_cast<real>(s->nr);
1944     rmsd = sqrt(rmsd);
1945
1946     return rmsd;
1947 }
1948
1949
1950 void dd_make_local_ed_indices(gmx_domdec_t* dd, struct gmx_edsam* ed)
1951 {
1952     if (ed->eEDtype != EssentialDynamicsType::None)
1953     {
1954         /* Loop over ED groups */
1955
1956         for (auto& edi : ed->edpar)
1957         {
1958             /* Local atoms of the reference structure (for fitting), need only be assembled
1959              * if their indices differ from the average ones */
1960             if (!edi.bRefEqAv)
1961             {
1962                 dd_make_local_group_indices(dd->ga2la, edi.sref.nr, edi.sref.anrs, &edi.sref.nr_loc,
1963                                             &edi.sref.anrs_loc, &edi.sref.nalloc_loc, edi.sref.c_ind);
1964             }
1965
1966             /* Local atoms of the average structure (on these ED will be performed) */
1967             dd_make_local_group_indices(dd->ga2la, edi.sav.nr, edi.sav.anrs, &edi.sav.nr_loc,
1968                                         &edi.sav.anrs_loc, &edi.sav.nalloc_loc, edi.sav.c_ind);
1969
1970             /* Indicate that the ED shift vectors for this structure need to be updated
1971              * at the next call to communicate_group_positions, since obviously we are in a NS step */
1972             edi.buf->do_edsam->bUpdateShifts = TRUE;
1973         }
1974     }
1975 }
1976
1977
1978 static inline void ed_unshift_single_coord(const matrix box, const rvec x, const ivec is, rvec xu)
1979 {
1980     int tx, ty, tz;
1981
1982
1983     tx = is[XX];
1984     ty = is[YY];
1985     tz = is[ZZ];
1986
1987     if (TRICLINIC(box))
1988     {
1989         xu[XX] = x[XX] - tx * box[XX][XX] - ty * box[YY][XX] - tz * box[ZZ][XX];
1990         xu[YY] = x[YY] - ty * box[YY][YY] - tz * box[ZZ][YY];
1991         xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
1992     }
1993     else
1994     {
1995         xu[XX] = x[XX] - tx * box[XX][XX];
1996         xu[YY] = x[YY] - ty * box[YY][YY];
1997         xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
1998     }
1999 }
2000
2001 namespace
2002 {
2003 /*!\brief Apply fixed linear constraints to essential dynamics variable.
2004  * \param[in,out] xcoll The collected coordinates.
2005  * \param[in] edi the essential dynamics parameters
2006  * \param[in] step the current simulation step
2007  */
2008 void do_linfix(rvec* xcoll, const t_edpar& edi, int64_t step)
2009 {
2010     /* loop over linfix vectors */
2011     for (int i = 0; i < edi.vecs.linfix.neig; i++)
2012     {
2013         /* calculate the projection */
2014         real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
2015
2016         /* calculate the correction */
2017         real preFactor = edi.vecs.linfix.refproj[i] + step * edi.vecs.linfix.stpsz[i] - proj;
2018
2019         /* apply the correction */
2020         preFactor /= edi.sav.sqrtm[i];
2021         for (int j = 0; j < edi.sav.nr; j++)
2022         {
2023             rvec differenceVector;
2024             svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
2025             rvec_inc(xcoll[j], differenceVector);
2026         }
2027     }
2028 }
2029
2030 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
2031  * \param[in,out] xcoll The collected coordinates.
2032  * \param[in] edi the essential dynamics parameters
2033  */
2034 void do_linacc(rvec* xcoll, t_edpar* edi)
2035 {
2036     /* loop over linacc vectors */
2037     for (int i = 0; i < edi->vecs.linacc.neig; i++)
2038     {
2039         /* calculate the projection */
2040         real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
2041
2042         /* calculate the correction */
2043         real preFactor = 0.0;
2044         if (edi->vecs.linacc.stpsz[i] > 0.0)
2045         {
2046             if ((proj - edi->vecs.linacc.refproj[i]) < 0.0)
2047             {
2048                 preFactor = edi->vecs.linacc.refproj[i] - proj;
2049             }
2050         }
2051         if (edi->vecs.linacc.stpsz[i] < 0.0)
2052         {
2053             if ((proj - edi->vecs.linacc.refproj[i]) > 0.0)
2054             {
2055                 preFactor = edi->vecs.linacc.refproj[i] - proj;
2056             }
2057         }
2058
2059         /* apply the correction */
2060         preFactor /= edi->sav.sqrtm[i];
2061         for (int j = 0; j < edi->sav.nr; j++)
2062         {
2063             rvec differenceVector;
2064             svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2065             rvec_inc(xcoll[j], differenceVector);
2066         }
2067
2068         /* new positions will act as reference */
2069         edi->vecs.linacc.refproj[i] = proj + preFactor;
2070     }
2071 }
2072 } // namespace
2073
2074
2075 static void do_radfix(rvec* xcoll, t_edpar* edi)
2076 {
2077     int   i, j;
2078     real *proj, rad = 0.0, ratio;
2079     rvec  vec_dum;
2080
2081
2082     if (edi->vecs.radfix.neig == 0)
2083     {
2084         return;
2085     }
2086
2087     snew(proj, edi->vecs.radfix.neig);
2088
2089     /* loop over radfix vectors */
2090     for (i = 0; i < edi->vecs.radfix.neig; i++)
2091     {
2092         /* calculate the projections, radius */
2093         proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2094         rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2095     }
2096
2097     rad   = sqrt(rad);
2098     ratio = (edi->vecs.radfix.stpsz[0] + edi->vecs.radfix.radius) / rad - 1.0;
2099     edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2100
2101     /* loop over radfix vectors */
2102     for (i = 0; i < edi->vecs.radfix.neig; i++)
2103     {
2104         proj[i] -= edi->vecs.radfix.refproj[i];
2105
2106         /* apply the correction */
2107         proj[i] /= edi->sav.sqrtm[i];
2108         proj[i] *= ratio;
2109         for (j = 0; j < edi->sav.nr; j++)
2110         {
2111             svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2112             rvec_inc(xcoll[j], vec_dum);
2113         }
2114     }
2115
2116     sfree(proj);
2117 }
2118
2119
2120 static void do_radacc(rvec* xcoll, t_edpar* edi)
2121 {
2122     int   i, j;
2123     real *proj, rad = 0.0, ratio = 0.0;
2124     rvec  vec_dum;
2125
2126
2127     if (edi->vecs.radacc.neig == 0)
2128     {
2129         return;
2130     }
2131
2132     snew(proj, edi->vecs.radacc.neig);
2133
2134     /* loop over radacc vectors */
2135     for (i = 0; i < edi->vecs.radacc.neig; i++)
2136     {
2137         /* calculate the projections, radius */
2138         proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2139         rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2140     }
2141     rad = sqrt(rad);
2142
2143     /* only correct when radius decreased */
2144     if (rad < edi->vecs.radacc.radius)
2145     {
2146         ratio = edi->vecs.radacc.radius / rad - 1.0;
2147     }
2148     else
2149     {
2150         edi->vecs.radacc.radius = rad;
2151     }
2152
2153     /* loop over radacc vectors */
2154     for (i = 0; i < edi->vecs.radacc.neig; i++)
2155     {
2156         proj[i] -= edi->vecs.radacc.refproj[i];
2157
2158         /* apply the correction */
2159         proj[i] /= edi->sav.sqrtm[i];
2160         proj[i] *= ratio;
2161         for (j = 0; j < edi->sav.nr; j++)
2162         {
2163             svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2164             rvec_inc(xcoll[j], vec_dum);
2165         }
2166     }
2167     sfree(proj);
2168 }
2169
2170
2171 struct t_do_radcon
2172 {
2173     real* proj;
2174 };
2175
2176 static void do_radcon(rvec* xcoll, t_edpar* edi)
2177 {
2178     int                 i, j;
2179     real                rad = 0.0, ratio = 0.0;
2180     struct t_do_radcon* loc;
2181     gmx_bool            bFirst;
2182     rvec                vec_dum;
2183
2184
2185     if (edi->buf->do_radcon != nullptr)
2186     {
2187         bFirst = FALSE;
2188     }
2189     else
2190     {
2191         bFirst = TRUE;
2192         snew(edi->buf->do_radcon, 1);
2193     }
2194     loc = edi->buf->do_radcon;
2195
2196     if (edi->vecs.radcon.neig == 0)
2197     {
2198         return;
2199     }
2200
2201     if (bFirst)
2202     {
2203         snew(loc->proj, edi->vecs.radcon.neig);
2204     }
2205
2206     /* loop over radcon vectors */
2207     for (i = 0; i < edi->vecs.radcon.neig; i++)
2208     {
2209         /* calculate the projections, radius */
2210         loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2211         rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2212     }
2213     rad = sqrt(rad);
2214     /* only correct when radius increased */
2215     if (rad > edi->vecs.radcon.radius)
2216     {
2217         ratio = edi->vecs.radcon.radius / rad - 1.0;
2218
2219         /* loop over radcon vectors */
2220         for (i = 0; i < edi->vecs.radcon.neig; i++)
2221         {
2222             /* apply the correction */
2223             loc->proj[i] -= edi->vecs.radcon.refproj[i];
2224             loc->proj[i] /= edi->sav.sqrtm[i];
2225             loc->proj[i] *= ratio;
2226
2227             for (j = 0; j < edi->sav.nr; j++)
2228             {
2229                 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2230                 rvec_inc(xcoll[j], vec_dum);
2231             }
2232         }
2233     }
2234     else
2235     {
2236         edi->vecs.radcon.radius = rad;
2237     }
2238 }
2239
2240
2241 static void ed_apply_constraints(rvec* xcoll, t_edpar* edi, int64_t step)
2242 {
2243     int i;
2244
2245
2246     /* subtract the average positions */
2247     for (i = 0; i < edi->sav.nr; i++)
2248     {
2249         rvec_dec(xcoll[i], edi->sav.x[i]);
2250     }
2251
2252     /* apply the constraints */
2253     if (step >= 0)
2254     {
2255         do_linfix(xcoll, *edi, step);
2256     }
2257     do_linacc(xcoll, edi);
2258     if (step >= 0)
2259     {
2260         do_radfix(xcoll, edi);
2261     }
2262     do_radacc(xcoll, edi);
2263     do_radcon(xcoll, edi);
2264
2265     /* add back the average positions */
2266     for (i = 0; i < edi->sav.nr; i++)
2267     {
2268         rvec_inc(xcoll[i], edi->sav.x[i]);
2269     }
2270 }
2271
2272
2273 namespace
2274 {
2275 /*!\brief Write out the projections onto the eigenvectors.
2276  * The order of output corresponds to ed_output_legend().
2277  * \param[in] edi The essential dyanmics parameters
2278  * \param[in] fp The output file
2279  * \param[in] rmsd the rmsd to the reference structure
2280  */
2281 void write_edo(const t_edpar& edi, FILE* fp, real rmsd)
2282 {
2283     /* Output how well we fit to the reference structure */
2284     fprintf(fp, EDcol_ffmt, rmsd);
2285
2286     for (int i = 0; i < edi.vecs.mon.neig; i++)
2287     {
2288         fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2289     }
2290
2291     for (int i = 0; i < edi.vecs.linfix.neig; i++)
2292     {
2293         fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2294     }
2295
2296     for (int i = 0; i < edi.vecs.linacc.neig; i++)
2297     {
2298         fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2299     }
2300
2301     for (int i = 0; i < edi.vecs.radfix.neig; i++)
2302     {
2303         fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2304     }
2305     if (edi.vecs.radfix.neig)
2306     {
2307         fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2308     }
2309
2310     for (int i = 0; i < edi.vecs.radacc.neig; i++)
2311     {
2312         fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2313     }
2314     if (edi.vecs.radacc.neig)
2315     {
2316         fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2317     }
2318
2319     for (int i = 0; i < edi.vecs.radcon.neig; i++)
2320     {
2321         fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2322     }
2323     if (edi.vecs.radcon.neig)
2324     {
2325         fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2326     }
2327 }
2328 } // namespace
2329
2330
2331 /* Returns if any constraints are switched on */
2332 static bool ed_constraints(EssentialDynamicsType edtype, const t_edpar& edi)
2333 {
2334     if (edtype == EssentialDynamicsType::EDSampling || edtype == EssentialDynamicsType::Flooding)
2335     {
2336         return ((edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0) || (edi.vecs.radfix.neig != 0)
2337                 || (edi.vecs.radacc.neig != 0) || (edi.vecs.radcon.neig != 0));
2338     }
2339     return false;
2340 }
2341
2342
2343 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for
2344  * flooding/ umbrella sampling simulations. */
2345 static void copyEvecReference(t_eigvec* floodvecs, real* initialReferenceProjection)
2346 {
2347     int i;
2348
2349
2350     if (nullptr == initialReferenceProjection)
2351     {
2352         snew(initialReferenceProjection, floodvecs->neig);
2353     }
2354
2355     for (i = 0; i < floodvecs->neig; i++)
2356     {
2357         initialReferenceProjection[i] = floodvecs->refproj[i];
2358     }
2359 }
2360
2361
2362 /* Call on MASTER only. Check whether the essential dynamics / flooding
2363  * groups of the checkpoint file are consistent with the provided .edi file. */
2364 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam& ed, edsamhistory_t* EDstate)
2365 {
2366     if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2367     {
2368         gmx_fatal(FARGS,
2369                   "Essential dynamics and flooding can only be switched on (or off) at the\n"
2370                   "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2371                   "it must also continue with/without ED constraints when checkpointing.\n"
2372                   "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2373                   "from without a checkpoint.\n");
2374     }
2375
2376     for (size_t edinum = 0; edinum < ed.edpar.size(); ++edinum)
2377     {
2378         /* Check number of atoms in the reference and average structures */
2379         if (EDstate->nref[edinum] != ed.edpar[edinum].sref.nr)
2380         {
2381             gmx_fatal(FARGS,
2382                       "The number of reference structure atoms in ED group %c is\n"
2383                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2384                       get_EDgroupChar(edinum + 1, 0), EDstate->nref[edinum], ed.edpar[edinum].sref.nr);
2385         }
2386         if (EDstate->nav[edinum] != ed.edpar[edinum].sav.nr)
2387         {
2388             gmx_fatal(FARGS,
2389                       "The number of average structure atoms in ED group %c is\n"
2390                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2391                       get_EDgroupChar(edinum + 1, 0), EDstate->nav[edinum], ed.edpar[edinum].sav.nr);
2392         }
2393     }
2394
2395     if (gmx::ssize(ed.edpar) != EDstate->nED)
2396     {
2397         gmx_fatal(FARGS,
2398                   "The number of essential dynamics / flooding groups is not consistent.\n"
2399                   "There are %d ED groups in the .cpt file, but %zu in the .edi file!\n"
2400                   "Are you sure this is the correct .edi file?\n",
2401                   EDstate->nED, ed.edpar.size());
2402     }
2403 }
2404
2405
2406 /* The edsamstate struct stores the information we need to make the ED group
2407  * whole again after restarts from a checkpoint file. Here we do the following:
2408  * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2409  * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2410  * c) in any case, for subsequent checkpoint writing, we set the pointers in
2411  * edsamstate to the x_old arrays, which contain the correct PBC representation of
2412  * all ED structures at the last time step. */
2413 static void init_edsamstate(const gmx_edsam& ed, edsamhistory_t* EDstate)
2414 {
2415     snew(EDstate->old_sref_p, EDstate->nED);
2416     snew(EDstate->old_sav_p, EDstate->nED);
2417
2418     /* If we did not read in a .cpt file, these arrays are not yet allocated */
2419     if (!EDstate->bFromCpt)
2420     {
2421         snew(EDstate->nref, EDstate->nED);
2422         snew(EDstate->nav, EDstate->nED);
2423     }
2424
2425     /* Loop over all ED/flooding data sets (usually only one, though) */
2426     for (int nr_edi = 0; nr_edi < EDstate->nED; nr_edi++)
2427     {
2428         const auto& edi = ed.edpar[nr_edi];
2429         /* We always need the last reference and average positions such that
2430          * in the next time step we can make the ED group whole again
2431          * if the atoms do not have the correct PBC representation */
2432         if (EDstate->bFromCpt)
2433         {
2434             /* Copy the last whole positions of reference and average group from .cpt */
2435             for (int i = 0; i < edi.sref.nr; i++)
2436             {
2437                 copy_rvec(EDstate->old_sref[nr_edi][i], edi.sref.x_old[i]);
2438             }
2439             for (int i = 0; i < edi.sav.nr; i++)
2440             {
2441                 copy_rvec(EDstate->old_sav[nr_edi][i], edi.sav.x_old[i]);
2442             }
2443         }
2444         else
2445         {
2446             EDstate->nref[nr_edi] = edi.sref.nr;
2447             EDstate->nav[nr_edi]  = edi.sav.nr;
2448         }
2449
2450         /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2451         EDstate->old_sref_p[nr_edi] = edi.sref.x_old;
2452         EDstate->old_sav_p[nr_edi]  = edi.sav.x_old;
2453     }
2454 }
2455
2456
2457 /* Adds 'buf' to 'str' */
2458 static void add_to_string(char** str, const char* buf)
2459 {
2460     int len;
2461
2462
2463     len = strlen(*str) + strlen(buf) + 1;
2464     srenew(*str, len);
2465     strcat(*str, buf);
2466 }
2467
2468
2469 static void add_to_string_aligned(char** str, const char* buf)
2470 {
2471     char buf_aligned[STRLEN];
2472
2473     sprintf(buf_aligned, EDcol_sfmt, buf);
2474     add_to_string(str, buf_aligned);
2475 }
2476
2477
2478 static void nice_legend(const char*** setname,
2479                         int*          nsets,
2480                         char**        LegendStr,
2481                         const char*   value,
2482                         const char*   unit,
2483                         char          EDgroupchar)
2484 {
2485     auto tmp = gmx::formatString("%c %s", EDgroupchar, value);
2486     add_to_string_aligned(LegendStr, tmp.c_str());
2487     tmp += gmx::formatString(" (%s)", unit);
2488     (*setname)[*nsets] = gmx_strdup(tmp.c_str());
2489     (*nsets)++;
2490 }
2491
2492
2493 static void nice_legend_evec(const char*** setname,
2494                              int*          nsets,
2495                              char**        LegendStr,
2496                              t_eigvec*     evec,
2497                              char          EDgroupChar,
2498                              const char*   EDtype)
2499 {
2500     int  i;
2501     char tmp[STRLEN];
2502
2503
2504     for (i = 0; i < evec->neig; i++)
2505     {
2506         sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2507         nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2508     }
2509 }
2510
2511
2512 /* Makes a legend for the xvg output file. Call on MASTER only! */
2513 static void write_edo_legend(gmx_edsam* ed, int nED, const gmx_output_env_t* oenv)
2514 {
2515     int          i;
2516     int          nr_edi, nsets, n_flood, n_edsam;
2517     const char** setname;
2518     char         buf[STRLEN];
2519     char*        LegendStr = nullptr;
2520
2521
2522     auto edi = ed->edpar.begin();
2523
2524     fprintf(ed->edo, "# Output will be written every %d step%s\n", edi->outfrq,
2525             edi->outfrq != 1 ? "s" : "");
2526
2527     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2528     {
2529         fprintf(ed->edo, "#\n");
2530         fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n",
2531                 get_EDgroupChar(nr_edi, nED));
2532         fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2533         fprintf(ed->edo, "#    monitor  : %d vec%s\n", edi->vecs.mon.neig,
2534                 edi->vecs.mon.neig != 1 ? "s" : "");
2535         fprintf(ed->edo, "#    LINFIX   : %d vec%s\n", edi->vecs.linfix.neig,
2536                 edi->vecs.linfix.neig != 1 ? "s" : "");
2537         fprintf(ed->edo, "#    LINACC   : %d vec%s\n", edi->vecs.linacc.neig,
2538                 edi->vecs.linacc.neig != 1 ? "s" : "");
2539         fprintf(ed->edo, "#    RADFIX   : %d vec%s\n", edi->vecs.radfix.neig,
2540                 edi->vecs.radfix.neig != 1 ? "s" : "");
2541         fprintf(ed->edo, "#    RADACC   : %d vec%s\n", edi->vecs.radacc.neig,
2542                 edi->vecs.radacc.neig != 1 ? "s" : "");
2543         fprintf(ed->edo, "#    RADCON   : %d vec%s\n", edi->vecs.radcon.neig,
2544                 edi->vecs.radcon.neig != 1 ? "s" : "");
2545         fprintf(ed->edo, "#    FLOODING : %d vec%s  ", edi->flood.vecs.neig,
2546                 edi->flood.vecs.neig != 1 ? "s" : "");
2547
2548         if (edi->flood.vecs.neig)
2549         {
2550             /* If in any of the groups we find a flooding vector, flooding is turned on */
2551             ed->eEDtype = EssentialDynamicsType::Flooding;
2552
2553             /* Print what flavor of flooding we will do */
2554             if (0 == edi->flood.tau) /* constant flooding strength */
2555             {
2556                 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2557                 if (edi->flood.bHarmonic)
2558                 {
2559                     fprintf(ed->edo, ", harmonic");
2560                 }
2561             }
2562             else /* adaptive flooding */
2563             {
2564                 fprintf(ed->edo, ", adaptive");
2565             }
2566         }
2567         fprintf(ed->edo, "\n");
2568
2569         ++edi;
2570     }
2571
2572     /* Print a nice legend */
2573     snew(LegendStr, 1);
2574     LegendStr[0] = '\0';
2575     sprintf(buf, "#     %6s", "time");
2576     add_to_string(&LegendStr, buf);
2577
2578     /* Calculate the maximum number of columns we could end up with */
2579     edi   = ed->edpar.begin();
2580     nsets = 0;
2581     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2582     {
2583         nsets += 5 + edi->vecs.mon.neig + edi->vecs.linfix.neig + edi->vecs.linacc.neig
2584                  + edi->vecs.radfix.neig + edi->vecs.radacc.neig + edi->vecs.radcon.neig
2585                  + 6 * edi->flood.vecs.neig;
2586         ++edi;
2587     }
2588     snew(setname, nsets);
2589
2590     /* In the mdrun time step in a first function call (do_flood()) the flooding
2591      * forces are calculated and in a second function call (do_edsam()) the
2592      * ED constraints. To get a corresponding legend, we need to loop twice
2593      * over the edi groups and output first the flooding, then the ED part */
2594
2595     /* The flooding-related legend entries, if flooding is done */
2596     nsets = 0;
2597     if (EssentialDynamicsType::Flooding == ed->eEDtype)
2598     {
2599         edi = ed->edpar.begin();
2600         for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2601         {
2602             /* Always write out the projection on the flooding EVs. Of course, this can also
2603              * be achieved with the monitoring option in do_edsam() (if switched on by the
2604              * user), but in that case the positions need to be communicated in do_edsam(),
2605              * which is not necessary when doing flooding only. */
2606             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2607
2608             for (i = 0; i < edi->flood.vecs.neig; i++)
2609             {
2610                 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2611                 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2612
2613                 /* Output the current reference projection if it changes with time;
2614                  * this can happen when flooding is used as harmonic restraint */
2615                 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2616                 {
2617                     sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2618                     nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2619                 }
2620
2621                 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2622                 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2623                 {
2624                     sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2625                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2626                 }
2627
2628                 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2629                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2630
2631                 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2632                 {
2633                     sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2634                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2635                 }
2636
2637                 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2638                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2639             }
2640
2641             ++edi;
2642         } /* End of flooding-related legend entries */
2643     }
2644     n_flood = nsets;
2645
2646     /* Now the ED-related entries, if essential dynamics is done */
2647     edi = ed->edpar.begin();
2648     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2649     {
2650         if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2651         {
2652             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2653
2654             /* Essential dynamics, projections on eigenvectors */
2655             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon,
2656                              get_EDgroupChar(nr_edi, nED), "MON");
2657             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix,
2658                              get_EDgroupChar(nr_edi, nED), "LINFIX");
2659             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc,
2660                              get_EDgroupChar(nr_edi, nED), "LINACC");
2661             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix,
2662                              get_EDgroupChar(nr_edi, nED), "RADFIX");
2663             if (edi->vecs.radfix.neig)
2664             {
2665                 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm",
2666                             get_EDgroupChar(nr_edi, nED));
2667             }
2668             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc,
2669                              get_EDgroupChar(nr_edi, nED), "RADACC");
2670             if (edi->vecs.radacc.neig)
2671             {
2672                 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm",
2673                             get_EDgroupChar(nr_edi, nED));
2674             }
2675             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon,
2676                              get_EDgroupChar(nr_edi, nED), "RADCON");
2677             if (edi->vecs.radcon.neig)
2678             {
2679                 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm",
2680                             get_EDgroupChar(nr_edi, nED));
2681             }
2682         }
2683         ++edi;
2684     } /* end of 'pure' essential dynamics legend entries */
2685     n_edsam = nsets - n_flood;
2686
2687     xvgr_legend(ed->edo, nsets, setname, oenv);
2688     sfree(setname);
2689
2690     fprintf(ed->edo,
2691             "#\n"
2692             "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2693             n_flood, 1 == n_flood ? "" : "s", n_edsam, 1 == n_edsam ? "" : "s");
2694     fprintf(ed->edo, "%s", LegendStr);
2695     sfree(LegendStr);
2696
2697     fflush(ed->edo);
2698 }
2699
2700
2701 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2702  * contained in the input file, creates a NULL terminated list of t_edpar structures */
2703 std::unique_ptr<gmx::EssentialDynamics> init_edsam(const gmx::MDLogger&        mdlog,
2704                                                    const char*                 ediFileName,
2705                                                    const char*                 edoFileName,
2706                                                    const gmx_mtop_t*           mtop,
2707                                                    const t_inputrec*           ir,
2708                                                    const t_commrec*            cr,
2709                                                    gmx::Constraints*           constr,
2710                                                    const t_state*              globalState,
2711                                                    ObservablesHistory*         oh,
2712                                                    const gmx_output_env_t*     oenv,
2713                                                    const gmx::StartingBehavior startingBehavior)
2714 {
2715     int    i, avindex;
2716     rvec*  x_pbc = nullptr; /* positions of the whole MD system with pbc removed  */
2717     rvec * xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs  */
2718     rvec   fit_transvec;                      /* translation ... */
2719     matrix fit_rotmat;                        /* ... and rotation from fit to reference structure */
2720     rvec*  ref_x_old = nullptr;               /* helper pointer */
2721
2722
2723     if (MASTER(cr))
2724     {
2725         fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2726
2727         if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName))
2728         {
2729             gmx_fatal(FARGS,
2730                       "The checkpoint file you provided is from an essential dynamics or flooding\n"
2731                       "simulation. Please also set the .edi file on the command line with -ei.\n");
2732         }
2733     }
2734     GMX_LOG(mdlog.info)
2735             .asParagraph()
2736             .appendText(
2737                     "Activating essential dynamics via files passed to "
2738                     "gmx mdrun -edi will change in future to be files passed to "
2739                     "gmx grompp and the related .mdp options may change also.");
2740
2741     /* Open input and output files, allocate space for ED data structure */
2742     auto edHandle = ed_open(mtop->natoms, oh, ediFileName, edoFileName, startingBehavior, oenv, cr);
2743     auto ed       = edHandle->getLegacyED();
2744     GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2745     constr->saveEdsamPointer(ed);
2746
2747     /* Needed for initializing radacc radius in do_edsam */
2748     ed->bFirst = TRUE;
2749
2750     /* The input file is read by the master and the edi structures are
2751      * initialized here. Input is stored in ed->edpar. Then the edi
2752      * structures are transferred to the other nodes */
2753     if (MASTER(cr))
2754     {
2755         /* Initialization for every ED/flooding group. Flooding uses one edi group per
2756          * flooding vector, Essential dynamics can be applied to more than one structure
2757          * as well, but will be done in the order given in the edi file, so
2758          * expect different results for different order of edi file concatenation! */
2759         for (auto& edi : ed->edpar)
2760         {
2761             init_edi(mtop, &edi);
2762             init_flood(&edi, ed, ir->delta_t);
2763         }
2764     }
2765
2766     /* The master does the work here. The other nodes get the positions
2767      * not before dd_partition_system which is called after init_edsam */
2768     if (MASTER(cr))
2769     {
2770         edsamhistory_t* EDstate = oh->edsamHistory.get();
2771
2772         if (!EDstate->bFromCpt)
2773         {
2774             /* Remove PBC, make molecule(s) subject to ED whole. */
2775             snew(x_pbc, mtop->natoms);
2776             copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop->natoms);
2777             do_pbc_first_mtop(nullptr, ir->pbcType, globalState->box, mtop, x_pbc);
2778         }
2779         /* Reset pointer to first ED data set which contains the actual ED data */
2780         auto edi = ed->edpar.begin();
2781         GMX_ASSERT(!ed->edpar.empty(), "Essential dynamics parameters is undefined");
2782
2783         /* Loop over all ED/flooding data sets (usually only one, though) */
2784         for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2785         {
2786             /* For multiple ED groups we use the output frequency that was specified
2787              * in the first set */
2788             if (nr_edi > 1)
2789             {
2790                 edi->outfrq = ed->edpar.begin()->outfrq;
2791             }
2792
2793             /* Extract the initial reference and average positions. When starting
2794              * from .cpt, these have already been read into sref.x_old
2795              * in init_edsamstate() */
2796             if (!EDstate->bFromCpt)
2797             {
2798                 /* If this is the first run (i.e. no checkpoint present) we assume
2799                  * that the starting positions give us the correct PBC representation */
2800                 for (i = 0; i < edi->sref.nr; i++)
2801                 {
2802                     copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2803                 }
2804
2805                 for (i = 0; i < edi->sav.nr; i++)
2806                 {
2807                     copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2808                 }
2809             }
2810
2811             /* Now we have the PBC-correct start positions of the reference and
2812                average structure. We copy that over to dummy arrays on which we
2813                can apply fitting to print out the RMSD. We srenew the memory since
2814                the size of the buffers is likely different for every ED group */
2815             srenew(xfit, edi->sref.nr);
2816             srenew(xstart, edi->sav.nr);
2817             if (edi->bRefEqAv)
2818             {
2819                 /* Reference indices are the same as average indices */
2820                 ref_x_old = edi->sav.x_old;
2821             }
2822             else
2823             {
2824                 ref_x_old = edi->sref.x_old;
2825             }
2826             copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2827             copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2828
2829             /* Make the fit to the REFERENCE structure, get translation and rotation */
2830             fit_to_reference(xfit, fit_transvec, fit_rotmat, &(*edi));
2831
2832             /* Output how well we fit to the reference at the start */
2833             translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2834             fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2835                     rmsd_from_structure(xfit, &edi->sref));
2836             if (EDstate->nED > 1)
2837             {
2838                 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2839             }
2840             fprintf(stderr, "\n");
2841
2842             /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2843             translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2844
2845             /* calculate initial projections */
2846             project(xstart, &(*edi));
2847
2848             /* For the target and origin structure both a reference (fit) and an
2849              * average structure can be provided in make_edi. If both structures
2850              * are the same, make_edi only stores one of them in the .edi file.
2851              * If they differ, first the fit and then the average structure is stored
2852              * in star (or sor), thus the number of entries in star/sor is
2853              * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2854              * the size of the average group. */
2855
2856             /* process target structure, if required */
2857             if (edi->star.nr > 0)
2858             {
2859                 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2860
2861                 /* get translation & rotation for fit of target structure to reference structure */
2862                 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, &(*edi));
2863                 /* do the fit */
2864                 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2865                 if (edi->star.nr == edi->sav.nr)
2866                 {
2867                     avindex = 0;
2868                 }
2869                 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2870                 {
2871                     /* The last sav.nr indices of the target structure correspond to
2872                      * the average structure, which must be projected */
2873                     avindex = edi->star.nr - edi->sav.nr;
2874                 }
2875                 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2876             }
2877             else
2878             {
2879                 rad_project(*edi, xstart, &edi->vecs.radcon);
2880             }
2881
2882             /* process structure that will serve as origin of expansion circle */
2883             if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2884             {
2885                 fprintf(stderr,
2886                         "ED: Setting center of flooding potential (0 = average structure)\n");
2887             }
2888
2889             if (edi->sori.nr > 0)
2890             {
2891                 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2892
2893                 /* fit this structure to reference structure */
2894                 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, &(*edi));
2895                 /* do the fit */
2896                 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2897                 if (edi->sori.nr == edi->sav.nr)
2898                 {
2899                     avindex = 0;
2900                 }
2901                 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2902                 {
2903                     /* For the projection, we need the last sav.nr indices of sori */
2904                     avindex = edi->sori.nr - edi->sav.nr;
2905                 }
2906
2907                 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2908                 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2909                 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2910                 {
2911                     fprintf(stderr,
2912                             "ED: The ORIGIN structure will define the flooding potential "
2913                             "center.\n");
2914                     /* Set center of flooding potential to the ORIGIN structure */
2915                     rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2916                     /* We already know that no (moving) reference position was provided,
2917                      * therefore we can overwrite refproj[0]*/
2918                     copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
2919                 }
2920             }
2921             else /* No origin structure given */
2922             {
2923                 rad_project(*edi, xstart, &edi->vecs.radacc);
2924                 rad_project(*edi, xstart, &edi->vecs.radfix);
2925                 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2926                 {
2927                     if (edi->flood.bHarmonic)
2928                     {
2929                         fprintf(stderr,
2930                                 "ED: A (possibly changing) ref. projection will define the "
2931                                 "flooding potential center.\n");
2932                         for (i = 0; i < edi->flood.vecs.neig; i++)
2933                         {
2934                             edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
2935                         }
2936                     }
2937                     else
2938                     {
2939                         fprintf(stderr,
2940                                 "ED: The AVERAGE structure will define the flooding potential "
2941                                 "center.\n");
2942                         /* Set center of flooding potential to the center of the covariance matrix,
2943                          * i.e. the average structure, i.e. zero in the projected system */
2944                         for (i = 0; i < edi->flood.vecs.neig; i++)
2945                         {
2946                             edi->flood.vecs.refproj[i] = 0.0;
2947                         }
2948                     }
2949                 }
2950             }
2951             /* For convenience, output the center of the flooding potential for the eigenvectors */
2952             if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2953             {
2954                 for (i = 0; i < edi->flood.vecs.neig; i++)
2955                 {
2956                     fprintf(stdout, "ED: EV %d flooding potential center: %11.4e",
2957                             edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2958                     if (edi->flood.bHarmonic)
2959                     {
2960                         fprintf(stdout, " (adding %11.4e/timestep)",
2961                                 edi->flood.referenceProjectionSlope[i]);
2962                     }
2963                     fprintf(stdout, "\n");
2964                 }
2965             }
2966
2967             /* set starting projections for linsam */
2968             rad_project(*edi, xstart, &edi->vecs.linacc);
2969             rad_project(*edi, xstart, &edi->vecs.linfix);
2970
2971             /* Prepare for the next edi data set: */
2972             ++edi;
2973         }
2974         /* Cleaning up on the master node: */
2975         if (!EDstate->bFromCpt)
2976         {
2977             sfree(x_pbc);
2978         }
2979         sfree(xfit);
2980         sfree(xstart);
2981
2982     } /* end of MASTER only section */
2983
2984     if (PAR(cr))
2985     {
2986         /* Broadcast the essential dynamics / flooding data to all nodes */
2987         broadcast_ed_data(cr, ed);
2988     }
2989     else
2990     {
2991         /* In the single-CPU case, point the local atom numbers pointers to the global
2992          * one, so that we can use the same notation in serial and parallel case: */
2993         /* Loop over all ED data sets (usually only one, though) */
2994         for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
2995         {
2996             edi->sref.anrs_loc = edi->sref.anrs;
2997             edi->sav.anrs_loc  = edi->sav.anrs;
2998             edi->star.anrs_loc = edi->star.anrs;
2999             edi->sori.anrs_loc = edi->sori.anrs;
3000             /* For the same reason as above, make a dummy c_ind array: */
3001             snew(edi->sav.c_ind, edi->sav.nr);
3002             /* Initialize the array */
3003             for (i = 0; i < edi->sav.nr; i++)
3004             {
3005                 edi->sav.c_ind[i] = i;
3006             }
3007             /* In the general case we will need a different-sized array for the reference indices: */
3008             if (!edi->bRefEqAv)
3009             {
3010                 snew(edi->sref.c_ind, edi->sref.nr);
3011                 for (i = 0; i < edi->sref.nr; i++)
3012                 {
3013                     edi->sref.c_ind[i] = i;
3014                 }
3015             }
3016             /* Point to the very same array in case of other structures: */
3017             edi->star.c_ind = edi->sav.c_ind;
3018             edi->sori.c_ind = edi->sav.c_ind;
3019             /* In the serial case, the local number of atoms is the global one: */
3020             edi->sref.nr_loc = edi->sref.nr;
3021             edi->sav.nr_loc  = edi->sav.nr;
3022             edi->star.nr_loc = edi->star.nr;
3023             edi->sori.nr_loc = edi->sori.nr;
3024         }
3025     }
3026
3027     /* Allocate space for ED buffer variables */
3028     /* Again, loop over ED data sets */
3029     for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
3030     {
3031         /* Allocate space for ED buffer variables */
3032         snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
3033         snew(edi->buf->do_edsam, 1);
3034
3035         /* Space for collective ED buffer variables */
3036
3037         /* Collective positions of atoms with the average indices */
3038         snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
3039         snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
3040         snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
3041         /* Collective positions of atoms with the reference indices */
3042         if (!edi->bRefEqAv)
3043         {
3044             snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
3045             snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
3046             snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
3047         }
3048
3049         /* Get memory for flooding forces */
3050         snew(edi->flood.forces_cartesian, edi->sav.nr);
3051     }
3052
3053     /* Flush the edo file so that the user can check some things
3054      * when the simulation has started */
3055     if (ed->edo)
3056     {
3057         fflush(ed->edo);
3058     }
3059
3060     return edHandle;
3061 }
3062
3063
3064 void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[], rvec v[], const matrix box, gmx_edsam* ed)
3065 {
3066     int    i, edinr, iupdate = 500;
3067     matrix rotmat;         /* rotation matrix */
3068     rvec   transvec;       /* translation vector */
3069     rvec   dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3070     real   dt_1;           /* 1/dt */
3071     struct t_do_edsam* buf;
3072     real     rmsdev    = -1; /* RMSD from reference structure prior to applying the constraints */
3073     gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3074
3075
3076     /* Check if ED sampling has to be performed */
3077     if (ed->eEDtype == EssentialDynamicsType::None)
3078     {
3079         return;
3080     }
3081
3082     dt_1 = 1.0 / ir->delta_t;
3083
3084     /* Loop over all ED groups (usually one) */
3085     edinr = 0;
3086     for (auto& edi : ed->edpar)
3087     {
3088         edinr++;
3089         if (bNeedDoEdsam(edi))
3090         {
3091
3092             buf = edi.buf->do_edsam;
3093
3094             if (ed->bFirst)
3095             {
3096                 /* initialize radacc radius for slope criterion */
3097                 buf->oldrad = calc_radius(edi.vecs.radacc);
3098             }
3099
3100             /* Copy the positions into buf->xc* arrays and after ED
3101              * feed back corrections to the official positions */
3102
3103             /* Broadcast the ED positions such that every node has all of them
3104              * Every node contributes its local positions xs and stores it in
3105              * the collective buf->xcoll array. Note that for edinr > 1
3106              * xs could already have been modified by an earlier ED */
3107
3108             communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll,
3109                                         PAR(cr) ? buf->bUpdateShifts : TRUE, xs, edi.sav.nr, edi.sav.nr_loc,
3110                                         edi.sav.anrs_loc, edi.sav.c_ind, edi.sav.x_old, box);
3111
3112             /* Only assembly reference positions if their indices differ from the average ones */
3113             if (!edi.bRefEqAv)
3114             {
3115                 communicate_group_positions(
3116                         cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref,
3117                         PAR(cr) ? buf->bUpdateShifts : TRUE, xs, edi.sref.nr, edi.sref.nr_loc,
3118                         edi.sref.anrs_loc, edi.sref.c_ind, edi.sref.x_old, box);
3119             }
3120
3121             /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3122              * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3123              * set bUpdateShifts=TRUE in the parallel case. */
3124             buf->bUpdateShifts = FALSE;
3125
3126             /* Now all nodes have all of the ED positions in edi.sav->xcoll,
3127              * as well as the indices in edi.sav.anrs */
3128
3129             /* Fit the reference indices to the reference structure */
3130             if (edi.bRefEqAv)
3131             {
3132                 fit_to_reference(buf->xcoll, transvec, rotmat, &edi);
3133             }
3134             else
3135             {
3136                 fit_to_reference(buf->xc_ref, transvec, rotmat, &edi);
3137             }
3138
3139             /* Now apply the translation and rotation to the ED structure */
3140             translate_and_rotate(buf->xcoll, edi.sav.nr, transvec, rotmat);
3141
3142             /* Find out how well we fit to the reference (just for output steps) */
3143             if (do_per_step(step, edi.outfrq) && MASTER(cr))
3144             {
3145                 if (edi.bRefEqAv)
3146                 {
3147                     /* Indices of reference and average structures are identical,
3148                      * thus we can calculate the rmsd to SREF using xcoll */
3149                     rmsdev = rmsd_from_structure(buf->xcoll, &edi.sref);
3150                 }
3151                 else
3152                 {
3153                     /* We have to translate & rotate the reference atoms first */
3154                     translate_and_rotate(buf->xc_ref, edi.sref.nr, transvec, rotmat);
3155                     rmsdev = rmsd_from_structure(buf->xc_ref, &edi.sref);
3156                 }
3157             }
3158
3159             /* update radsam references, when required */
3160             if (do_per_step(step, edi.maxedsteps) && step >= edi.presteps)
3161             {
3162                 project(buf->xcoll, &edi);
3163                 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3164                 rad_project(edi, buf->xcoll, &edi.vecs.radfix);
3165                 buf->oldrad = -1.e5;
3166             }
3167
3168             /* update radacc references, when required */
3169             if (do_per_step(step, iupdate) && step >= edi.presteps)
3170             {
3171                 edi.vecs.radacc.radius = calc_radius(edi.vecs.radacc);
3172                 if (edi.vecs.radacc.radius - buf->oldrad < edi.slope)
3173                 {
3174                     project(buf->xcoll, &edi);
3175                     rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3176                     buf->oldrad = 0.0;
3177                 }
3178                 else
3179                 {
3180                     buf->oldrad = edi.vecs.radacc.radius;
3181                 }
3182             }
3183
3184             /* apply the constraints */
3185             if (step >= edi.presteps && ed_constraints(ed->eEDtype, edi))
3186             {
3187                 /* ED constraints should be applied already in the first MD step
3188                  * (which is step 0), therefore we pass step+1 to the routine */
3189                 ed_apply_constraints(buf->xcoll, &edi, step + 1 - ir->init_step);
3190             }
3191
3192             /* write to edo, when required */
3193             if (do_per_step(step, edi.outfrq))
3194             {
3195                 project(buf->xcoll, &edi);
3196                 if (MASTER(cr) && !bSuppress)
3197                 {
3198                     write_edo(edi, ed->edo, rmsdev);
3199                 }
3200             }
3201
3202             /* Copy back the positions unless monitoring only */
3203             if (ed_constraints(ed->eEDtype, edi))
3204             {
3205                 /* remove fitting */
3206                 rmfit(edi.sav.nr, buf->xcoll, transvec, rotmat);
3207
3208                 /* Copy the ED corrected positions into the coordinate array */
3209                 /* Each node copies its local part. In the serial case, nat_loc is the
3210                  * total number of ED atoms */
3211                 for (i = 0; i < edi.sav.nr_loc; i++)
3212                 {
3213                     /* Unshift local ED coordinate and store in x_unsh */
3214                     ed_unshift_single_coord(box, buf->xcoll[edi.sav.c_ind[i]],
3215                                             buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh);
3216
3217                     /* dx is the ED correction to the positions: */
3218                     rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx);
3219
3220                     if (v != nullptr)
3221                     {
3222                         /* dv is the ED correction to the velocity: */
3223                         svmul(dt_1, dx, dv);
3224                         /* apply the velocity correction: */
3225                         rvec_inc(v[edi.sav.anrs_loc[i]], dv);
3226                     }
3227                     /* Finally apply the position correction due to ED: */
3228                     copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]);
3229                 }
3230             }
3231         } /* END of if ( bNeedDoEdsam(edi) ) */
3232
3233         /* Prepare for the next ED group */
3234
3235     } /* END of loop over ED groups */
3236
3237     ed->bFirst = FALSE;
3238 }