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,2021, 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 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,
 930                                 buf->xcoll,
 931                                 buf->shifts_xcoll,
 932                                 buf->extra_shifts_xcoll,
 933                                 bNS,
 934                                 x,
 935                                 edi->sav.nr,
 936                                 edi->sav.nr_loc,
 937                                 edi->sav.anrs_loc,
 938                                 edi->sav.c_ind,
 939                                 edi->sav.x_old,
 940                                 box);
 941
 942     /* Only assembly REFERENCE positions if their indices differ from the average ones */
 943     if (!edi->bRefEqAv)
 944     {
 945         communicate_group_positions(cr,
 946                                     buf->xc_ref,
 947                                     buf->shifts_xc_ref,
 948                                     buf->extra_shifts_xc_ref,
 949                                     bNS,
 950                                     x,
 951                                     edi->sref.nr,
 952                                     edi->sref.nr_loc,
 953                                     edi->sref.anrs_loc,
 954                                     edi->sref.c_ind,
 955                                     edi->sref.x_old,
 956                                     box);
 957     }
 958
 959     /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
 960      * We do not need to update the shifts until the next NS step */
 961     buf->bUpdateShifts = FALSE;
 962
 963     /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
 964      * as well as the indices in edi->sav.anrs */
 965
 966     /* Fit the reference indices to the reference structure */
 967     if (edi->bRefEqAv)
 968     {
 969         fit_to_reference(buf->xcoll, transvec, rotmat, edi);
 970     }
 971     else
 972     {
 973         fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
 974     }
 975
 976     /* Now apply the translation and rotation to the ED structure */
 977     translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
 978
 979     /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
 980     project_to_eigvectors(buf->xcoll, &edi->flood.vecs, *edi);
 981
 982     if (!edi->flood.bConstForce)
 983     {
 984         /* Compute Vfl(x) from flood.xproj */
 985         edi->flood.Vfl = flood_energy(edi, step);
 986
 987         update_adaption(edi);
 988
 989         /* Compute the flooding forces */
 990         flood_forces(edi);
 991     }
 992
 993     /* Translate them into cartesian positions */
 994     flood_blowup(*edi, edi->flood.forces_cartesian);
 995
 996     /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
 997     /* Each node rotates back its local forces */
 998     transpose(rotmat, tmat);
 999     rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
1000
1001     /* Finally add forces to the main force variable */
1002     for (i = 0; i < edi->sav.nr_loc; i++)
1003     {
1004         rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]);
1005     }
1006
1007     /* Output is written by the master process */
1008     if (do_per_step(step, edi->outfrq) && MASTER(cr))
1009     {
1010         /* Output how well we fit to the reference */
1011         if (edi->bRefEqAv)
1012         {
1013             /* Indices of reference and average structures are identical,
1014              * thus we can calculate the rmsd to SREF using xcoll */
1015             rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
1016         }
1017         else
1018         {
1019             /* We have to translate & rotate the reference atoms first */
1020             translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
1021             rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
1022         }
1023
1024         write_edo_flood(*edi, edo, rmsdev);
1025     }
1026 }
1027
1028
1029 /* Main flooding routine, called from do_force */
1030 void do_flood(const t_commrec*  cr,
1031               const t_inputrec& ir,
1032               const rvec        x[],
1033               rvec              force[],
1034               gmx_edsam*        ed,
1035               const matrix      box,
1036               int64_t           step,
1037               bool              bNS)
1038 {
1039     /* Write time to edo, when required. Output the time anyhow since we need
1040      * it in the output file for ED constraints. */
1041     if (MASTER(cr) && do_per_step(step, ed->edpar.begin()->outfrq))
1042     {
1043         fprintf(ed->edo, "\n%12f", ir.init_t + step * ir.delta_t);
1044     }
1045
1046     if (ed->eEDtype != EssentialDynamicsType::Flooding)
1047     {
1048         return;
1049     }
1050
1051     for (auto& edi : ed->edpar)
1052     {
1053         /* Call flooding for one matrix */
1054         if (edi.flood.vecs.neig)
1055         {
1056             do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS);
1057         }
1058     }
1059 }
1060
1061
1062 /* Called by init_edi, configure some flooding related variables and structures,
1063  * print headers to output files */
1064 static void init_flood(t_edpar* edi, gmx_edsam* ed, real dt)
1065 {
1066     int i;
1067
1068
1069     edi->flood.Efl = edi->flood.constEfl;
1070     edi->flood.Vfl = 0;
1071     edi->flood.dt  = dt;
1072
1073     if (edi->flood.vecs.neig)
1074     {
1075         /* If in any of the ED groups we find a flooding vector, flooding is turned on */
1076         ed->eEDtype = EssentialDynamicsType::Flooding;
1077
1078         fprintf(stderr,
1079                 "ED: Flooding %d eigenvector%s.\n",
1080                 edi->flood.vecs.neig,
1081                 edi->flood.vecs.neig > 1 ? "s" : "");
1082
1083         if (edi->flood.bConstForce)
1084         {
1085             /* We have used stpsz as a vehicle to carry the fproj values for constant
1086              * force flooding. Now we copy that to flood.vecs.fproj. Note that
1087              * in const force flooding, fproj is never changed. */
1088             for (i = 0; i < edi->flood.vecs.neig; i++)
1089             {
1090                 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1091
1092                 fprintf(stderr,
1093                         "ED: applying on eigenvector %d a constant force of %g\n",
1094                         edi->flood.vecs.ieig[i],
1095                         edi->flood.vecs.fproj[i]);
1096             }
1097         }
1098     }
1099 }
1100
1101
1102 #ifdef DEBUGHELPERS
1103 /*********** Energy book keeping ******/
1104 static void get_flood_enx_names(t_edpar* edi, char** names, int* nnames) /* get header of energies */
1105 {
1106     t_edpar* actual;
1107     int      count;
1108     char     buf[STRLEN];
1109     actual = edi;
1110     count  = 1;
1111     while (actual)
1112     {
1113         srenew(names, count);
1114         sprintf(buf, "Vfl_%d", count);
1115         names[count - 1] = gmx_strdup(buf);
1116         actual           = actual->next_edi;
1117         count++;
1118     }
1119     *nnames = count - 1;
1120 }
1121
1122
1123 static void get_flood_energies(t_edpar* edi, real Vfl[], int nnames)
1124 {
1125     /*fl has to be big enough to capture nnames-many entries*/
1126     t_edpar* actual;
1127     int      count;
1128
1129
1130     actual = edi;
1131     count  = 1;
1132     while (actual)
1133     {
1134         Vfl[count - 1] = actual->flood.Vfl;
1135         actual         = actual->next_edi;
1136         count++;
1137     }
1138     if (nnames != count - 1)
1139     {
1140         gmx_fatal(FARGS, "Number of energies is not consistent with t_edi structure");
1141     }
1142 }
1143 /************* END of FLOODING IMPLEMENTATION ****************************/
1144 #endif
1145
1146
1147 /* This function opens the ED input and output files, reads in all datasets it finds
1148  * in the input file, and cross-checks whether the .edi file information is consistent
1149  * with the essential dynamics data found in the checkpoint file (if present).
1150  * gmx make_edi can be used to create an .edi input file.
1151  */
1152 static std::unique_ptr<gmx::EssentialDynamics> ed_open(int                         natoms,
1153                                                        ObservablesHistory*         oh,
1154                                                        const char*                 ediFileName,
1155                                                        const char*                 edoFileName,
1156                                                        const gmx::StartingBehavior startingBehavior,
1157                                                        const gmx_output_env_t*     oenv,
1158                                                        const t_commrec*            cr)
1159 {
1160     auto edHandle = std::make_unique<gmx::EssentialDynamics>();
1161     auto ed       = edHandle->getLegacyED();
1162     /* We want to perform ED (this switch might later be upgraded to EssentialDynamicsType::Flooding) */
1163     ed->eEDtype = EssentialDynamicsType::EDSampling;
1164
1165     if (MASTER(cr))
1166     {
1167         // If we start from a checkpoint file, we already have an edsamHistory struct
1168         if (oh->edsamHistory == nullptr)
1169         {
1170             oh->edsamHistory = std::make_unique<edsamhistory_t>(edsamhistory_t{});
1171         }
1172         edsamhistory_t* EDstate = oh->edsamHistory.get();
1173
1174         /* Read the edi input file: */
1175         ed->edpar = read_edi_file(ediFileName, natoms);
1176
1177         /* Make sure the checkpoint was produced in a run using this .edi file */
1178         if (EDstate->bFromCpt)
1179         {
1180             crosscheck_edi_file_vs_checkpoint(*ed, EDstate);
1181         }
1182         else
1183         {
1184             EDstate->nED = ed->edpar.size();
1185         }
1186         init_edsamstate(*ed, EDstate);
1187
1188         /* The master opens the ED output file */
1189         if (startingBehavior == gmx::StartingBehavior::RestartWithAppending)
1190         {
1191             ed->edo = gmx_fio_fopen(edoFileName, "a+");
1192         }
1193         else
1194         {
1195             ed->edo = xvgropen(edoFileName,
1196                                "Essential dynamics / flooding output",
1197                                "Time (ps)",
1198                                "RMSDs (nm), projections on EVs (nm), ...",
1199                                oenv);
1200
1201             /* Make a descriptive legend */
1202             write_edo_legend(ed, EDstate->nED, oenv);
1203         }
1204     }
1205     return edHandle;
1206 }
1207
1208
1209 /* Broadcasts the structure data */
1210 static void bc_ed_positions(const t_commrec* cr, struct gmx_edx* s, EssentialDynamicsStructure stype)
1211 {
1212     snew_bc(MASTER(cr), s->anrs, s->nr); /* Index numbers     */
1213     snew_bc(MASTER(cr), s->x, s->nr);    /* Positions         */
1214     nblock_bc(cr->mpi_comm_mygroup, s->nr, s->anrs);
1215     nblock_bc(cr->mpi_comm_mygroup, s->nr, s->x);
1216
1217     /* For the average & reference structures we need an array for the collective indices,
1218      * and we need to broadcast the masses as well */
1219     if (stype == EssentialDynamicsStructure::Average || stype == EssentialDynamicsStructure::Reference)
1220     {
1221         /* We need these additional variables in the parallel case: */
1222         snew(s->c_ind, s->nr); /* Collective indices */
1223         /* Local atom indices get assigned in dd_make_local_group_indices.
1224          * There, also memory is allocated */
1225         s->nalloc_loc = 0;                    /* allocation size of s->anrs_loc */
1226         snew_bc(MASTER(cr), s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1227         nblock_bc(cr->mpi_comm_mygroup, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1228     }
1229
1230     /* broadcast masses for the reference structure (for mass-weighted fitting) */
1231     if (stype == EssentialDynamicsStructure::Reference)
1232     {
1233         snew_bc(MASTER(cr), s->m, s->nr);
1234         nblock_bc(cr->mpi_comm_mygroup, s->nr, s->m);
1235     }
1236
1237     /* For the average structure we might need the masses for mass-weighting */
1238     if (stype == EssentialDynamicsStructure::Average)
1239     {
1240         snew_bc(MASTER(cr), s->sqrtm, s->nr);
1241         nblock_bc(cr->mpi_comm_mygroup, s->nr, s->sqrtm);
1242         snew_bc(MASTER(cr), s->m, s->nr);
1243         nblock_bc(cr->mpi_comm_mygroup, s->nr, s->m);
1244     }
1245 }
1246
1247
1248 /* Broadcasts the eigenvector data */
1249 static void bc_ed_vecs(const t_commrec* cr, t_eigvec* ev, int length)
1250 {
1251     int i;
1252
1253     snew_bc(MASTER(cr), ev->ieig, ev->neig);    /* index numbers of eigenvector  */
1254     snew_bc(MASTER(cr), ev->stpsz, ev->neig);   /* stepsizes per eigenvector     */
1255     snew_bc(MASTER(cr), ev->xproj, ev->neig);   /* instantaneous x projection    */
1256     snew_bc(MASTER(cr), ev->fproj, ev->neig);   /* instantaneous f projection    */
1257     snew_bc(MASTER(cr), ev->refproj, ev->neig); /* starting or target projection */
1258
1259     nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->ieig);
1260     nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->stpsz);
1261     nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->xproj);
1262     nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->fproj);
1263     nblock_bc(cr->mpi_comm_mygroup, ev->neig, ev->refproj);
1264
1265     snew_bc(MASTER(cr), ev->vec, ev->neig); /* Eigenvector components        */
1266     for (i = 0; i < ev->neig; i++)
1267     {
1268         snew_bc(MASTER(cr), ev->vec[i], length);
1269         nblock_bc(cr->mpi_comm_mygroup, length, ev->vec[i]);
1270     }
1271 }
1272
1273
1274 /* Broadcasts the ED / flooding data to other nodes
1275  * and allocates memory where needed */
1276 static void broadcast_ed_data(const t_commrec* cr, gmx_edsam* ed)
1277 {
1278     /* Master lets the other nodes know if its ED only or also flooding */
1279     gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr->mpi_comm_mygroup);
1280
1281     int numedis = ed->edpar.size();
1282     /* First let everybody know how many ED data sets to expect */
1283     gmx_bcast(sizeof(numedis), &numedis, cr->mpi_comm_mygroup);
1284     nblock_abc(MASTER(cr), cr->mpi_comm_mygroup, numedis, &(ed->edpar));
1285     /* Now transfer the ED data set(s) */
1286     for (auto& edi : ed->edpar)
1287     {
1288         /* Broadcast a single ED data set */
1289         block_bc(cr->mpi_comm_mygroup, edi);
1290
1291         /* Broadcast positions */
1292         bc_ed_positions(cr, &(edi.sref), EssentialDynamicsStructure::Reference); /* reference positions (don't broadcast masses)    */
1293         bc_ed_positions(cr, &(edi.sav), EssentialDynamicsStructure::Average); /* average positions (do broadcast masses as well) */
1294         bc_ed_positions(cr, &(edi.star), EssentialDynamicsStructure::Target); /* target positions */
1295         bc_ed_positions(cr, &(edi.sori), EssentialDynamicsStructure::Origin); /* origin positions */
1296
1297         /* Broadcast eigenvectors */
1298         bc_ed_vecs(cr, &edi.vecs.mon, edi.sav.nr);
1299         bc_ed_vecs(cr, &edi.vecs.linfix, edi.sav.nr);
1300         bc_ed_vecs(cr, &edi.vecs.linacc, edi.sav.nr);
1301         bc_ed_vecs(cr, &edi.vecs.radfix, edi.sav.nr);
1302         bc_ed_vecs(cr, &edi.vecs.radacc, edi.sav.nr);
1303         bc_ed_vecs(cr, &edi.vecs.radcon, edi.sav.nr);
1304         /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1305         bc_ed_vecs(cr, &edi.flood.vecs, edi.sav.nr);
1306
1307         /* For harmonic restraints the reference projections can change with time */
1308         if (edi.flood.bHarmonic)
1309         {
1310             snew_bc(MASTER(cr), edi.flood.initialReferenceProjection, edi.flood.vecs.neig);
1311             snew_bc(MASTER(cr), edi.flood.referenceProjectionSlope, edi.flood.vecs.neig);
1312             nblock_bc(cr->mpi_comm_mygroup, edi.flood.vecs.neig, edi.flood.initialReferenceProjection);
1313             nblock_bc(cr->mpi_comm_mygroup, edi.flood.vecs.neig, edi.flood.referenceProjectionSlope);
1314         }
1315     }
1316 }
1317
1318
1319 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1320 static void init_edi(const gmx_mtop_t& mtop, t_edpar* edi)
1321 {
1322     int  i;
1323     real totalmass = 0.0;
1324     rvec com;
1325
1326     /* NOTE Init_edi is executed on the master process only
1327      * The initialized data sets are then transmitted to the
1328      * other nodes in broadcast_ed_data */
1329
1330     /* evaluate masses (reference structure) */
1331     snew(edi->sref.m, edi->sref.nr);
1332     int molb = 0;
1333     for (i = 0; i < edi->sref.nr; i++)
1334     {
1335         if (edi->fitmas)
1336         {
1337             edi->sref.m[i] = mtopGetAtomMass(&mtop, edi->sref.anrs[i], &molb);
1338         }
1339         else
1340         {
1341             edi->sref.m[i] = 1.0;
1342         }
1343
1344         /* Check that every m > 0. Bad things will happen otherwise. */
1345         if (edi->sref.m[i] <= 0.0)
1346         {
1347             gmx_fatal(FARGS,
1348                       "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1349                       "For a mass-weighted fit, all reference structure atoms need to have a mass "
1350                       ">0.\n"
1351                       "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1352                       "atoms from the reference structure by creating a proper index group.\n",
1353                       i,
1354                       edi->sref.anrs[i] + 1,
1355                       edi->sref.m[i]);
1356         }
1357
1358         totalmass += edi->sref.m[i];
1359     }
1360     edi->sref.mtot = totalmass;
1361
1362     /* Masses m and sqrt(m) for the average structure. Note that m
1363      * is needed if forces have to be evaluated in do_edsam */
1364     snew(edi->sav.sqrtm, edi->sav.nr);
1365     snew(edi->sav.m, edi->sav.nr);
1366     for (i = 0; i < edi->sav.nr; i++)
1367     {
1368         edi->sav.m[i] = mtopGetAtomMass(&mtop, edi->sav.anrs[i], &molb);
1369         if (edi->pcamas)
1370         {
1371             edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
1372         }
1373         else
1374         {
1375             edi->sav.sqrtm[i] = 1.0;
1376         }
1377
1378         /* Check that every m > 0. Bad things will happen otherwise. */
1379         if (edi->sav.sqrtm[i] <= 0.0)
1380         {
1381             gmx_fatal(FARGS,
1382                       "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1383                       "For ED with mass-weighting, all average structure atoms need to have a mass "
1384                       ">0.\n"
1385                       "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1386                       "atoms from the average structure by creating a proper index group.\n",
1387                       i,
1388                       edi->sav.anrs[i] + 1,
1389                       edi->sav.m[i]);
1390         }
1391     }
1392
1393     /* put reference structure in origin */
1394     get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1395     com[XX] = -com[XX];
1396     com[YY] = -com[YY];
1397     com[ZZ] = -com[ZZ];
1398     translate_x(edi->sref.x, edi->sref.nr, com);
1399
1400     /* Init ED buffer */
1401     snew(edi->buf, 1);
1402 }
1403
1404
1405 static void check(const char* line, const char* label)
1406 {
1407     if (!strstr(line, label))
1408     {
1409         gmx_fatal(FARGS,
1410                   "Could not find input parameter %s at expected position in edsam input-file "
1411                   "(.edi)\nline read instead is %s",
1412                   label,
1413                   line);
1414     }
1415 }
1416
1417
1418 static int read_checked_edint(FILE* file, const char* label)
1419 {
1420     char line[STRLEN + 1];
1421     int  idum;
1422
1423     fgets2(line, STRLEN, file);
1424     check(line, label);
1425     fgets2(line, STRLEN, file);
1426     sscanf(line, max_ev_fmt_d, &idum);
1427     return idum;
1428 }
1429
1430 static bool read_checked_edbool(FILE* file, const char* label)
1431 {
1432     return read_checked_edint(file, label) != 0;
1433 }
1434
1435 static int read_edint(FILE* file, bool* bEOF)
1436 {
1437     char  line[STRLEN + 1];
1438     int   idum;
1439     char* eof;
1440
1441
1442     eof = fgets2(line, STRLEN, file);
1443     if (eof == nullptr)
1444     {
1445         *bEOF = TRUE;
1446         return -1;
1447     }
1448     eof = fgets2(line, STRLEN, file);
1449     if (eof == nullptr)
1450     {
1451         *bEOF = TRUE;
1452         return -1;
1453     }
1454     sscanf(line, max_ev_fmt_d, &idum);
1455     *bEOF = FALSE;
1456     return idum;
1457 }
1458
1459
1460 static real read_checked_edreal(FILE* file, const char* label)
1461 {
1462     char   line[STRLEN + 1];
1463     double rdum;
1464
1465
1466     fgets2(line, STRLEN, file);
1467     check(line, label);
1468     fgets2(line, STRLEN, file);
1469     sscanf(line, max_ev_fmt_lf, &rdum);
1470     return static_cast<real>(rdum); /* always read as double and convert to single */
1471 }
1472
1473
1474 static void read_edx(FILE* file, int number, int* anrs, rvec* x)
1475 {
1476     int    i, j;
1477     char   line[STRLEN + 1];
1478     double d[3];
1479
1480
1481     for (i = 0; i < number; i++)
1482     {
1483         fgets2(line, STRLEN, file);
1484         sscanf(line, max_ev_fmt_dlflflf, &anrs[i], &d[0], &d[1], &d[2]);
1485         anrs[i]--; /* we are reading FORTRAN indices */
1486         for (j = 0; j < 3; j++)
1487         {
1488             x[i][j] = d[j]; /* always read as double and convert to single */
1489         }
1490     }
1491 }
1492
1493 namespace
1494 {
1495 /*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
1496  * \param[in] in the file to read from
1497  * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
1498  * \param[out] vec the eigen-vectors
1499  * \param[in] nEig the number of eigenvectors
1500  */
1501 void scan_edvec(FILE* in, int numAtoms, rvec*** vec, int nEig)
1502 {
1503     snew(*vec, nEig);
1504     for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
1505     {
1506         snew((*vec)[iEigenvector], numAtoms);
1507         for (int iAtom = 0; iAtom < numAtoms; iAtom++)
1508         {
1509             char   line[STRLEN + 1];
1510             double x, y, z;
1511             fgets2(line, STRLEN, in);
1512             sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
1513             (*vec)[iEigenvector][iAtom][XX] = x;
1514             (*vec)[iEigenvector][iAtom][YY] = y;
1515             (*vec)[iEigenvector][iAtom][ZZ] = z;
1516         }
1517     }
1518 }
1519 /*!\brief Read an essential dynamics eigenvector and corresponding step size.
1520  * \param[in] in input file
1521  * \param[in] nrAtoms number of atoms/coordinates
1522  * \param[out] tvec the eigenvector
1523  */
1524 void read_edvec(FILE* in, int nrAtoms, t_eigvec* tvec)
1525 {
1526     tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1527     if (tvec->neig <= 0)
1528     {
1529         return;
1530     }
1531
1532     snew(tvec->ieig, tvec->neig);
1533     snew(tvec->stpsz, tvec->neig);
1534     for (int i = 0; i < tvec->neig; i++)
1535     {
1536         char line[STRLEN + 1];
1537         fgets2(line, STRLEN, in);
1538         int       numEigenVectors;
1539         double    stepSize;
1540         const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
1541         if (nscan != 2)
1542         {
1543             gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1544         }
1545         tvec->ieig[i]  = numEigenVectors;
1546         tvec->stpsz[i] = stepSize;
1547     } /* end of loop over eigenvectors */
1548
1549     scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1550 }
1551 /*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
1552  *
1553  * Eigenvector and its intial reference and slope are stored on the
1554  * same line in the .edi format. To avoid re-winding the .edi file,
1555  * the reference eigen-vector and reference data are read in one go.
1556  *
1557  * \param[in] in input file
1558  * \param[in] nrAtoms number of atoms/coordinates
1559  * \param[out] tvec the eigenvector
1560  * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
1561  * \param[out] referenceProjectionSlope The slope of the reference projections.
1562  */
1563 bool readEdVecWithReferenceProjection(FILE*     in,
1564                                       int       nrAtoms,
1565                                       t_eigvec* tvec,
1566                                       real**    initialReferenceProjection,
1567                                       real**    referenceProjectionSlope)
1568 {
1569     bool providesReference = false;
1570     tvec->neig             = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1571     if (tvec->neig <= 0)
1572     {
1573         return false;
1574     }
1575
1576     snew(tvec->ieig, tvec->neig);
1577     snew(tvec->stpsz, tvec->neig);
1578     snew(*initialReferenceProjection, tvec->neig);
1579     snew(*referenceProjectionSlope, tvec->neig);
1580     for (int i = 0; i < tvec->neig; i++)
1581     {
1582         char line[STRLEN + 1];
1583         fgets2(line, STRLEN, in);
1584         int       numEigenVectors;
1585         double    stepSize            = 0.;
1586         double    referenceProjection = 0.;
1587         double    referenceSlope      = 0.;
1588         const int nscan               = sscanf(
1589                 line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize, &referenceProjection, &referenceSlope);
1590         /* Zero out values which were not scanned */
1591         switch (nscan)
1592         {
1593             case 4:
1594                 /* Every 4 values read, including reference position */
1595                 providesReference = true;
1596                 break;
1597             case 3:
1598                 /* A reference position is provided */
1599                 providesReference = true;
1600                 /* No value for slope, set to 0 */
1601                 referenceSlope = 0.0;
1602                 break;
1603             case 2:
1604                 /* No values for reference projection and slope, set to 0 */
1605                 referenceProjection = 0.0;
1606                 referenceSlope      = 0.0;
1607                 break;
1608             default:
1609                 gmx_fatal(FARGS,
1610                           "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> "
1611                           "<refproj> <refproj-slope>\n",
1612                           nscan);
1613         }
1614         (*initialReferenceProjection)[i] = referenceProjection;
1615         (*referenceProjectionSlope)[i]   = referenceSlope;
1616
1617         tvec->ieig[i]  = numEigenVectors;
1618         tvec->stpsz[i] = stepSize;
1619     } /* end of loop over eigenvectors */
1620
1621     scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
1622     return providesReference;
1623 }
1624
1625 /*!\brief Allocate working memory for the eigenvectors.
1626  * \param[in,out] tvec eigenvector for which memory will be allocated
1627  */
1628 void setup_edvec(t_eigvec* tvec)
1629 {
1630     snew(tvec->xproj, tvec->neig);
1631     snew(tvec->fproj, tvec->neig);
1632     snew(tvec->refproj, tvec->neig);
1633 }
1634 } // namespace
1635
1636
1637 /* Check if the same atom indices are used for reference and average positions */
1638 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1639 {
1640     int i;
1641
1642
1643     /* If the number of atoms differs between the two structures,
1644      * they cannot be identical */
1645     if (sref.nr != sav.nr)
1646     {
1647         return FALSE;
1648     }
1649
1650     /* Now that we know that both stuctures have the same number of atoms,
1651      * check if also the indices are identical */
1652     for (i = 0; i < sav.nr; i++)
1653     {
1654         if (sref.anrs[i] != sav.anrs[i])
1655         {
1656             return FALSE;
1657         }
1658     }
1659     fprintf(stderr,
1660             "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1661
1662     return TRUE;
1663 }
1664
1665 namespace
1666 {
1667
1668 //! Oldest (lowest) magic number indicating unsupported essential dynamics input
1669 constexpr int c_oldestUnsupportedMagicNumber = 666;
1670 //! Oldest (lowest) magic number indicating supported essential dynamics input
1671 constexpr int c_oldestSupportedMagicNumber = 669;
1672 //! Latest (highest) magic number indicating supported essential dynamics input
1673 constexpr int c_latestSupportedMagicNumber = 670;
1674
1675 /*!\brief Set up essential dynamics work parameters.
1676  * \param[in] edi Essential dynamics working structure.
1677  */
1678 void setup_edi(t_edpar* edi)
1679 {
1680     snew(edi->sref.x_old, edi->sref.nr);
1681     edi->sref.sqrtm = nullptr;
1682     snew(edi->sav.x_old, edi->sav.nr);
1683     if (edi->star.nr > 0)
1684     {
1685         edi->star.sqrtm = nullptr;
1686     }
1687     if (edi->sori.nr > 0)
1688     {
1689         edi->sori.sqrtm = nullptr;
1690     }
1691     setup_edvec(&edi->vecs.linacc);
1692     setup_edvec(&edi->vecs.mon);
1693     setup_edvec(&edi->vecs.linfix);
1694     setup_edvec(&edi->vecs.linacc);
1695     setup_edvec(&edi->vecs.radfix);
1696     setup_edvec(&edi->vecs.radacc);
1697     setup_edvec(&edi->vecs.radcon);
1698     setup_edvec(&edi->flood.vecs);
1699 }
1700
1701 /*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
1702  * \param[in] magicNumber indicates the essential dynamics input file version
1703  * \returns true if CONST_FORCE_FLOODING is supported
1704  */
1705 bool ediFileHasConstForceFlooding(int magicNumber)
1706 {
1707     return magicNumber > c_oldestSupportedMagicNumber;
1708 };
1709
1710 /*!\brief Reports reading success of the essential dynamics input file magic number.
1711  * \param[in] in pointer to input file
1712  * \param[out] magicNumber determines the edi file version
1713  * \returns true if setting file version from input file was successful.
1714  */
1715 bool ediFileHasValidDataSet(FILE* in, int* magicNumber)
1716 {
1717     /* the edi file is not free format, so expect problems if the input is corrupt. */
1718     bool reachedEndOfFile = true;
1719     /* check the magic number */
1720     *magicNumber = read_edint(in, &reachedEndOfFile);
1721     /* Check whether we have reached the end of the input file */
1722     return !reachedEndOfFile;
1723 }
1724
1725 /*!\brief Trigger fatal error if magic number is unsupported.
1726  * \param[in] magicNumber A magic number read from the edi file
1727  * \param[in] fn name of input file for error reporting
1728  */
1729 void exitWhenMagicNumberIsInvalid(int magicNumber, const char* fn)
1730 {
1731
1732     if (magicNumber < c_oldestSupportedMagicNumber || magicNumber > c_latestSupportedMagicNumber)
1733     {
1734         if (magicNumber >= c_oldestUnsupportedMagicNumber && magicNumber < c_oldestSupportedMagicNumber)
1735         {
1736             gmx_fatal(FARGS,
1737                       "Wrong magic number: Use newest version of make_edi to produce edi file");
1738         }
1739         else
1740         {
1741             gmx_fatal(FARGS, "Wrong magic number %d in %s", magicNumber, fn);
1742         }
1743     }
1744 }
1745
1746 /*!\brief reads an essential dynamics input file into a essential dynamics data structure.
1747  *
1748  * \param[in] in input file
1749  * \param[in] nr_mdatoms the number of md atoms, used for consistency checking
1750  * \param[in] hasConstForceFlooding determines if field "CONST_FORCE_FLOODING" is available in input file
1751  * \param[in] fn the file name of the input file for error reporting
1752  * \returns edi essential dynamics parameters
1753  */
1754 t_edpar read_edi(FILE* in, int nr_mdatoms, bool hasConstForceFlooding, const char* fn)
1755 {
1756     t_edpar edi;
1757     bool    bEOF;
1758     /* check the number of atoms */
1759     edi.nini = read_edint(in, &bEOF);
1760     if (edi.nini != nr_mdatoms)
1761     {
1762         gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi.nini, nr_mdatoms);
1763     }
1764
1765     /* Done checking. For the rest we blindly trust the input */
1766     edi.fitmas     = read_checked_edbool(in, "FITMAS");
1767     edi.pcamas     = read_checked_edbool(in, "ANALYSIS_MAS");
1768     edi.outfrq     = read_checked_edint(in, "OUTFRQ");
1769     edi.maxedsteps = read_checked_edint(in, "MAXLEN");
1770     edi.slope      = read_checked_edreal(in, "SLOPECRIT");
1771
1772     edi.presteps        = read_checked_edint(in, "PRESTEPS");
1773     edi.flood.deltaF0   = read_checked_edreal(in, "DELTA_F0");
1774     edi.flood.deltaF    = read_checked_edreal(in, "INIT_DELTA_F");
1775     edi.flood.tau       = read_checked_edreal(in, "TAU");
1776     edi.flood.constEfl  = read_checked_edreal(in, "EFL_NULL");
1777     edi.flood.alpha2    = read_checked_edreal(in, "ALPHA2");
1778     edi.flood.kT        = read_checked_edreal(in, "KT");
1779     edi.flood.bHarmonic = read_checked_edbool(in, "HARMONIC");
1780     if (hasConstForceFlooding)
1781     {
1782         edi.flood.bConstForce = read_checked_edbool(in, "CONST_FORCE_FLOODING");
1783     }
1784     else
1785     {
1786         edi.flood.bConstForce = FALSE;
1787     }
1788     edi.sref.nr = read_checked_edint(in, "NREF");
1789
1790     /* allocate space for reference positions and read them */
1791     snew(edi.sref.anrs, edi.sref.nr);
1792     snew(edi.sref.x, edi.sref.nr);
1793     read_edx(in, edi.sref.nr, edi.sref.anrs, edi.sref.x);
1794
1795     /* average positions. they define which atoms will be used for ED sampling */
1796     edi.sav.nr = read_checked_edint(in, "NAV");
1797     snew(edi.sav.anrs, edi.sav.nr);
1798     snew(edi.sav.x, edi.sav.nr);
1799     read_edx(in, edi.sav.nr, edi.sav.anrs, edi.sav.x);
1800
1801     /* Check if the same atom indices are used for reference and average positions */
1802     edi.bRefEqAv = check_if_same(edi.sref, edi.sav);
1803
1804     /* eigenvectors */
1805
1806     read_edvec(in, edi.sav.nr, &edi.vecs.mon);
1807     read_edvec(in, edi.sav.nr, &edi.vecs.linfix);
1808     read_edvec(in, edi.sav.nr, &edi.vecs.linacc);
1809     read_edvec(in, edi.sav.nr, &edi.vecs.radfix);
1810     read_edvec(in, edi.sav.nr, &edi.vecs.radacc);
1811     read_edvec(in, edi.sav.nr, &edi.vecs.radcon);
1812
1813     /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1814     bool bHaveReference = false;
1815     if (edi.flood.bHarmonic)
1816     {
1817         bHaveReference = readEdVecWithReferenceProjection(in,
1818                                                           edi.sav.nr,
1819                                                           &edi.flood.vecs,
1820                                                           &edi.flood.initialReferenceProjection,
1821                                                           &edi.flood.referenceProjectionSlope);
1822     }
1823     else
1824     {
1825         read_edvec(in, edi.sav.nr, &edi.flood.vecs);
1826     }
1827
1828     /* target positions */
1829     edi.star.nr = read_edint(in, &bEOF);
1830     if (edi.star.nr > 0)
1831     {
1832         snew(edi.star.anrs, edi.star.nr);
1833         snew(edi.star.x, edi.star.nr);
1834         read_edx(in, edi.star.nr, edi.star.anrs, edi.star.x);
1835     }
1836
1837     /* positions defining origin of expansion circle */
1838     edi.sori.nr = read_edint(in, &bEOF);
1839     if (edi.sori.nr > 0)
1840     {
1841         if (bHaveReference)
1842         {
1843             /* Both an -ori structure and a at least one manual reference point have been
1844              * specified. That's ambiguous and probably not intentional. */
1845             gmx_fatal(FARGS,
1846                       "ED: An origin structure has been provided and a at least one (moving) "
1847                       "reference\n"
1848                       "    point was manually specified in the edi file. That is ambiguous. "
1849                       "Aborting.\n");
1850         }
1851         snew(edi.sori.anrs, edi.sori.nr);
1852         snew(edi.sori.x, edi.sori.nr);
1853         read_edx(in, edi.sori.nr, edi.sori.anrs, edi.sori.x);
1854     }
1855     return edi;
1856 }
1857
1858 std::vector<t_edpar> read_edi_file(const char* fn, int nr_mdatoms)
1859 {
1860     std::vector<t_edpar> essentialDynamicsParameters;
1861     FILE*                in;
1862     /* This routine is executed on the master only */
1863
1864     /* Open the .edi parameter input file */
1865     in = gmx_fio_fopen(fn, "r");
1866     fprintf(stderr, "ED: Reading edi file %s\n", fn);
1867
1868     /* Now read a sequence of ED input parameter sets from the edi file */
1869     int ediFileMagicNumber;
1870     while (ediFileHasValidDataSet(in, &ediFileMagicNumber))
1871     {
1872         exitWhenMagicNumberIsInvalid(ediFileMagicNumber, fn);
1873         bool hasConstForceFlooding = ediFileHasConstForceFlooding(ediFileMagicNumber);
1874         auto edi                   = read_edi(in, nr_mdatoms, hasConstForceFlooding, fn);
1875         setup_edi(&edi);
1876         essentialDynamicsParameters.emplace_back(edi);
1877
1878         /* Since we arrived within this while loop we know that there is still another data set to be read in */
1879         /* We need to allocate space for the data: */
1880     }
1881     if (essentialDynamicsParameters.empty())
1882     {
1883         gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1884     }
1885
1886     fprintf(stderr,
1887             "ED: Found %zu ED group%s.\n",
1888             essentialDynamicsParameters.size(),
1889             essentialDynamicsParameters.size() > 1 ? "s" : "");
1890
1891     /* Close the .edi file again */
1892     gmx_fio_fclose(in);
1893
1894     return essentialDynamicsParameters;
1895 }
1896 } // anonymous namespace
1897
1898
1899 struct t_fit_to_ref
1900 {
1901     rvec* xcopy; /* Working copy of the positions in fit_to_reference */
1902 };
1903
1904 /* Fit the current positions to the reference positions
1905  * Do not actually do the fit, just return rotation and translation.
1906  * Note that the COM of the reference structure was already put into
1907  * the origin by init_edi. */
1908 static void fit_to_reference(rvec*    xcoll,    /* The positions to be fitted */
1909                              rvec     transvec, /* The translation vector */
1910                              matrix   rotmat,   /* The rotation matrix */
1911                              t_edpar* edi)      /* Just needed for do_edfit */
1912 {
1913     rvec                 com; /* center of mass */
1914     int                  i;
1915     struct t_fit_to_ref* loc;
1916
1917
1918     /* Allocate memory the first time this routine is called for each edi group */
1919     if (nullptr == edi->buf->fit_to_ref)
1920     {
1921         snew(edi->buf->fit_to_ref, 1);
1922         snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1923     }
1924     loc = edi->buf->fit_to_ref;
1925
1926     /* We do not touch the original positions but work on a copy. */
1927     for (i = 0; i < edi->sref.nr; i++)
1928     {
1929         copy_rvec(xcoll[i], loc->xcopy[i]);
1930     }
1931
1932     /* Calculate the center of mass */
1933     get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1934
1935     transvec[XX] = -com[XX];
1936     transvec[YY] = -com[YY];
1937     transvec[ZZ] = -com[ZZ];
1938
1939     /* Subtract the center of mass from the copy */
1940     translate_x(loc->xcopy, edi->sref.nr, transvec);
1941
1942     /* Determine the rotation matrix */
1943     do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1944 }
1945
1946
1947 static void translate_and_rotate(rvec*  x,        /* The positions to be translated and rotated */
1948                                  int    nat,      /* How many positions are there? */
1949                                  rvec   transvec, /* The translation vector */
1950                                  matrix rotmat)   /* The rotation matrix */
1951 {
1952     /* Translation */
1953     translate_x(x, nat, transvec);
1954
1955     /* Rotation */
1956     rotate_x(x, nat, rotmat);
1957 }
1958
1959
1960 /* Gets the rms deviation of the positions to the structure s */
1961 /* fit_to_structure has to be called before calling this routine! */
1962 static real rmsd_from_structure(rvec* x,           /* The positions under consideration */
1963                                 struct gmx_edx* s) /* The structure from which the rmsd shall be computed */
1964 {
1965     real rmsd = 0.0;
1966     int  i;
1967
1968
1969     for (i = 0; i < s->nr; i++)
1970     {
1971         rmsd += distance2(s->x[i], x[i]);
1972     }
1973
1974     rmsd /= static_cast<real>(s->nr);
1975     rmsd = sqrt(rmsd);
1976
1977     return rmsd;
1978 }
1979
1980
1981 void dd_make_local_ed_indices(gmx_domdec_t* dd, struct gmx_edsam* ed)
1982 {
1983     if (ed->eEDtype != EssentialDynamicsType::None)
1984     {
1985         /* Loop over ED groups */
1986
1987         for (auto& edi : ed->edpar)
1988         {
1989             /* Local atoms of the reference structure (for fitting), need only be assembled
1990              * if their indices differ from the average ones */
1991             if (!edi.bRefEqAv)
1992             {
1993                 dd_make_local_group_indices(dd->ga2la,
1994                                             edi.sref.nr,
1995                                             edi.sref.anrs,
1996                                             &edi.sref.nr_loc,
1997                                             &edi.sref.anrs_loc,
1998                                             &edi.sref.nalloc_loc,
1999                                             edi.sref.c_ind);
2000             }
2001
2002             /* Local atoms of the average structure (on these ED will be performed) */
2003             dd_make_local_group_indices(dd->ga2la,
2004                                         edi.sav.nr,
2005                                         edi.sav.anrs,
2006                                         &edi.sav.nr_loc,
2007                                         &edi.sav.anrs_loc,
2008                                         &edi.sav.nalloc_loc,
2009                                         edi.sav.c_ind);
2010
2011             /* Indicate that the ED shift vectors for this structure need to be updated
2012              * at the next call to communicate_group_positions, since obviously we are in a NS step */
2013             edi.buf->do_edsam->bUpdateShifts = TRUE;
2014         }
2015     }
2016 }
2017
2018
2019 static inline void ed_unshift_single_coord(const matrix box, const rvec x, const ivec is, rvec xu)
2020 {
2021     int tx, ty, tz;
2022
2023
2024     tx = is[XX];
2025     ty = is[YY];
2026     tz = is[ZZ];
2027
2028     if (TRICLINIC(box))
2029     {
2030         xu[XX] = x[XX] - tx * box[XX][XX] - ty * box[YY][XX] - tz * box[ZZ][XX];
2031         xu[YY] = x[YY] - ty * box[YY][YY] - tz * box[ZZ][YY];
2032         xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
2033     }
2034     else
2035     {
2036         xu[XX] = x[XX] - tx * box[XX][XX];
2037         xu[YY] = x[YY] - ty * box[YY][YY];
2038         xu[ZZ] = x[ZZ] - tz * box[ZZ][ZZ];
2039     }
2040 }
2041
2042 namespace
2043 {
2044 /*!\brief Apply fixed linear constraints to essential dynamics variable.
2045  * \param[in,out] xcoll The collected coordinates.
2046  * \param[in] edi the essential dynamics parameters
2047  * \param[in] step the current simulation step
2048  */
2049 void do_linfix(rvec* xcoll, const t_edpar& edi, int64_t step)
2050 {
2051     /* loop over linfix vectors */
2052     for (int i = 0; i < edi.vecs.linfix.neig; i++)
2053     {
2054         /* calculate the projection */
2055         real proj = projectx(edi, xcoll, edi.vecs.linfix.vec[i]);
2056
2057         /* calculate the correction */
2058         real preFactor = edi.vecs.linfix.refproj[i] + step * edi.vecs.linfix.stpsz[i] - proj;
2059
2060         /* apply the correction */
2061         preFactor /= edi.sav.sqrtm[i];
2062         for (int j = 0; j < edi.sav.nr; j++)
2063         {
2064             rvec differenceVector;
2065             svmul(preFactor, edi.vecs.linfix.vec[i][j], differenceVector);
2066             rvec_inc(xcoll[j], differenceVector);
2067         }
2068     }
2069 }
2070
2071 /*!\brief Apply acceptance linear constraints to essential dynamics variable.
2072  * \param[in,out] xcoll The collected coordinates.
2073  * \param[in] edi the essential dynamics parameters
2074  */
2075 void do_linacc(rvec* xcoll, t_edpar* edi)
2076 {
2077     /* loop over linacc vectors */
2078     for (int i = 0; i < edi->vecs.linacc.neig; i++)
2079     {
2080         /* calculate the projection */
2081         real proj = projectx(*edi, xcoll, edi->vecs.linacc.vec[i]);
2082
2083         /* calculate the correction */
2084         real preFactor = 0.0;
2085         if (edi->vecs.linacc.stpsz[i] > 0.0)
2086         {
2087             if ((proj - edi->vecs.linacc.refproj[i]) < 0.0)
2088             {
2089                 preFactor = edi->vecs.linacc.refproj[i] - proj;
2090             }
2091         }
2092         if (edi->vecs.linacc.stpsz[i] < 0.0)
2093         {
2094             if ((proj - edi->vecs.linacc.refproj[i]) > 0.0)
2095             {
2096                 preFactor = edi->vecs.linacc.refproj[i] - proj;
2097             }
2098         }
2099
2100         /* apply the correction */
2101         preFactor /= edi->sav.sqrtm[i];
2102         for (int j = 0; j < edi->sav.nr; j++)
2103         {
2104             rvec differenceVector;
2105             svmul(preFactor, edi->vecs.linacc.vec[i][j], differenceVector);
2106             rvec_inc(xcoll[j], differenceVector);
2107         }
2108
2109         /* new positions will act as reference */
2110         edi->vecs.linacc.refproj[i] = proj + preFactor;
2111     }
2112 }
2113 } // namespace
2114
2115
2116 static void do_radfix(rvec* xcoll, t_edpar* edi)
2117 {
2118     int   i, j;
2119     real *proj, rad = 0.0, ratio;
2120     rvec  vec_dum;
2121
2122
2123     if (edi->vecs.radfix.neig == 0)
2124     {
2125         return;
2126     }
2127
2128     snew(proj, edi->vecs.radfix.neig);
2129
2130     /* loop over radfix vectors */
2131     for (i = 0; i < edi->vecs.radfix.neig; i++)
2132     {
2133         /* calculate the projections, radius */
2134         proj[i] = projectx(*edi, xcoll, edi->vecs.radfix.vec[i]);
2135         rad += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
2136     }
2137
2138     rad   = sqrt(rad);
2139     ratio = (edi->vecs.radfix.stpsz[0] + edi->vecs.radfix.radius) / rad - 1.0;
2140     edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2141
2142     /* loop over radfix vectors */
2143     for (i = 0; i < edi->vecs.radfix.neig; i++)
2144     {
2145         proj[i] -= edi->vecs.radfix.refproj[i];
2146
2147         /* apply the correction */
2148         proj[i] /= edi->sav.sqrtm[i];
2149         proj[i] *= ratio;
2150         for (j = 0; j < edi->sav.nr; j++)
2151         {
2152             svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2153             rvec_inc(xcoll[j], vec_dum);
2154         }
2155     }
2156
2157     sfree(proj);
2158 }
2159
2160
2161 static void do_radacc(rvec* xcoll, t_edpar* edi)
2162 {
2163     int   i, j;
2164     real *proj, rad = 0.0, ratio = 0.0;
2165     rvec  vec_dum;
2166
2167
2168     if (edi->vecs.radacc.neig == 0)
2169     {
2170         return;
2171     }
2172
2173     snew(proj, edi->vecs.radacc.neig);
2174
2175     /* loop over radacc vectors */
2176     for (i = 0; i < edi->vecs.radacc.neig; i++)
2177     {
2178         /* calculate the projections, radius */
2179         proj[i] = projectx(*edi, xcoll, edi->vecs.radacc.vec[i]);
2180         rad += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2181     }
2182     rad = sqrt(rad);
2183
2184     /* only correct when radius decreased */
2185     if (rad < edi->vecs.radacc.radius)
2186     {
2187         ratio = edi->vecs.radacc.radius / rad - 1.0;
2188     }
2189     else
2190     {
2191         edi->vecs.radacc.radius = rad;
2192     }
2193
2194     /* loop over radacc vectors */
2195     for (i = 0; i < edi->vecs.radacc.neig; i++)
2196     {
2197         proj[i] -= edi->vecs.radacc.refproj[i];
2198
2199         /* apply the correction */
2200         proj[i] /= edi->sav.sqrtm[i];
2201         proj[i] *= ratio;
2202         for (j = 0; j < edi->sav.nr; j++)
2203         {
2204             svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2205             rvec_inc(xcoll[j], vec_dum);
2206         }
2207     }
2208     sfree(proj);
2209 }
2210
2211
2212 struct t_do_radcon
2213 {
2214     real* proj;
2215 };
2216
2217 static void do_radcon(rvec* xcoll, t_edpar* edi)
2218 {
2219     int                 i, j;
2220     real                rad = 0.0, ratio = 0.0;
2221     struct t_do_radcon* loc;
2222     gmx_bool            bFirst;
2223     rvec                vec_dum;
2224
2225
2226     if (edi->buf->do_radcon != nullptr)
2227     {
2228         bFirst = FALSE;
2229     }
2230     else
2231     {
2232         bFirst = TRUE;
2233         snew(edi->buf->do_radcon, 1);
2234     }
2235     loc = edi->buf->do_radcon;
2236
2237     if (edi->vecs.radcon.neig == 0)
2238     {
2239         return;
2240     }
2241
2242     if (bFirst)
2243     {
2244         snew(loc->proj, edi->vecs.radcon.neig);
2245     }
2246
2247     /* loop over radcon vectors */
2248     for (i = 0; i < edi->vecs.radcon.neig; i++)
2249     {
2250         /* calculate the projections, radius */
2251         loc->proj[i] = projectx(*edi, xcoll, edi->vecs.radcon.vec[i]);
2252         rad += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2253     }
2254     rad = sqrt(rad);
2255     /* only correct when radius increased */
2256     if (rad > edi->vecs.radcon.radius)
2257     {
2258         ratio = edi->vecs.radcon.radius / rad - 1.0;
2259
2260         /* loop over radcon vectors */
2261         for (i = 0; i < edi->vecs.radcon.neig; i++)
2262         {
2263             /* apply the correction */
2264             loc->proj[i] -= edi->vecs.radcon.refproj[i];
2265             loc->proj[i] /= edi->sav.sqrtm[i];
2266             loc->proj[i] *= ratio;
2267
2268             for (j = 0; j < edi->sav.nr; j++)
2269             {
2270                 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2271                 rvec_inc(xcoll[j], vec_dum);
2272             }
2273         }
2274     }
2275     else
2276     {
2277         edi->vecs.radcon.radius = rad;
2278     }
2279 }
2280
2281
2282 static void ed_apply_constraints(rvec* xcoll, t_edpar* edi, int64_t step)
2283 {
2284     int i;
2285
2286
2287     /* subtract the average positions */
2288     for (i = 0; i < edi->sav.nr; i++)
2289     {
2290         rvec_dec(xcoll[i], edi->sav.x[i]);
2291     }
2292
2293     /* apply the constraints */
2294     if (step >= 0)
2295     {
2296         do_linfix(xcoll, *edi, step);
2297     }
2298     do_linacc(xcoll, edi);
2299     if (step >= 0)
2300     {
2301         do_radfix(xcoll, edi);
2302     }
2303     do_radacc(xcoll, edi);
2304     do_radcon(xcoll, edi);
2305
2306     /* add back the average positions */
2307     for (i = 0; i < edi->sav.nr; i++)
2308     {
2309         rvec_inc(xcoll[i], edi->sav.x[i]);
2310     }
2311 }
2312
2313
2314 namespace
2315 {
2316 /*!\brief Write out the projections onto the eigenvectors.
2317  * The order of output corresponds to ed_output_legend().
2318  * \param[in] edi The essential dyanmics parameters
2319  * \param[in] fp The output file
2320  * \param[in] rmsd the rmsd to the reference structure
2321  */
2322 void write_edo(const t_edpar& edi, FILE* fp, real rmsd)
2323 {
2324     /* Output how well we fit to the reference structure */
2325     fprintf(fp, EDcol_ffmt, rmsd);
2326
2327     for (int i = 0; i < edi.vecs.mon.neig; i++)
2328     {
2329         fprintf(fp, EDcol_efmt, edi.vecs.mon.xproj[i]);
2330     }
2331
2332     for (int i = 0; i < edi.vecs.linfix.neig; i++)
2333     {
2334         fprintf(fp, EDcol_efmt, edi.vecs.linfix.xproj[i]);
2335     }
2336
2337     for (int i = 0; i < edi.vecs.linacc.neig; i++)
2338     {
2339         fprintf(fp, EDcol_efmt, edi.vecs.linacc.xproj[i]);
2340     }
2341
2342     for (int i = 0; i < edi.vecs.radfix.neig; i++)
2343     {
2344         fprintf(fp, EDcol_efmt, edi.vecs.radfix.xproj[i]);
2345     }
2346     if (edi.vecs.radfix.neig)
2347     {
2348         fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radfix)); /* fixed increment radius */
2349     }
2350
2351     for (int i = 0; i < edi.vecs.radacc.neig; i++)
2352     {
2353         fprintf(fp, EDcol_efmt, edi.vecs.radacc.xproj[i]);
2354     }
2355     if (edi.vecs.radacc.neig)
2356     {
2357         fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radacc)); /* acceptance radius */
2358     }
2359
2360     for (int i = 0; i < edi.vecs.radcon.neig; i++)
2361     {
2362         fprintf(fp, EDcol_efmt, edi.vecs.radcon.xproj[i]);
2363     }
2364     if (edi.vecs.radcon.neig)
2365     {
2366         fprintf(fp, EDcol_ffmt, calc_radius(edi.vecs.radcon)); /* contracting radius */
2367     }
2368 }
2369 } // namespace
2370
2371
2372 /* Returns if any constraints are switched on */
2373 static bool ed_constraints(EssentialDynamicsType edtype, const t_edpar& edi)
2374 {
2375     if (edtype == EssentialDynamicsType::EDSampling || edtype == EssentialDynamicsType::Flooding)
2376     {
2377         return ((edi.vecs.linfix.neig != 0) || (edi.vecs.linacc.neig != 0) || (edi.vecs.radfix.neig != 0)
2378                 || (edi.vecs.radacc.neig != 0) || (edi.vecs.radcon.neig != 0));
2379     }
2380     return false;
2381 }
2382
2383
2384 /* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for
2385  * flooding/ umbrella sampling simulations. */
2386 static void copyEvecReference(t_eigvec* floodvecs, real* initialReferenceProjection)
2387 {
2388     int i;
2389
2390
2391     if (nullptr == initialReferenceProjection)
2392     {
2393         snew(initialReferenceProjection, floodvecs->neig);
2394     }
2395
2396     for (i = 0; i < floodvecs->neig; i++)
2397     {
2398         initialReferenceProjection[i] = floodvecs->refproj[i];
2399     }
2400 }
2401
2402
2403 /* Call on MASTER only. Check whether the essential dynamics / flooding
2404  * groups of the checkpoint file are consistent with the provided .edi file. */
2405 static void crosscheck_edi_file_vs_checkpoint(const gmx_edsam& ed, edsamhistory_t* EDstate)
2406 {
2407     if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2408     {
2409         gmx_fatal(FARGS,
2410                   "Essential dynamics and flooding can only be switched on (or off) at the\n"
2411                   "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2412                   "it must also continue with/without ED constraints when checkpointing.\n"
2413                   "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2414                   "from without a checkpoint.\n");
2415     }
2416
2417     for (size_t edinum = 0; edinum < ed.edpar.size(); ++edinum)
2418     {
2419         /* Check number of atoms in the reference and average structures */
2420         if (EDstate->nref[edinum] != ed.edpar[edinum].sref.nr)
2421         {
2422             gmx_fatal(FARGS,
2423                       "The number of reference structure atoms in ED group %c is\n"
2424                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2425                       get_EDgroupChar(edinum + 1, 0),
2426                       EDstate->nref[edinum],
2427                       ed.edpar[edinum].sref.nr);
2428         }
2429         if (EDstate->nav[edinum] != ed.edpar[edinum].sav.nr)
2430         {
2431             gmx_fatal(FARGS,
2432                       "The number of average structure atoms in ED group %c is\n"
2433                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2434                       get_EDgroupChar(edinum + 1, 0),
2435                       EDstate->nav[edinum],
2436                       ed.edpar[edinum].sav.nr);
2437         }
2438     }
2439
2440     if (gmx::ssize(ed.edpar) != EDstate->nED)
2441     {
2442         gmx_fatal(FARGS,
2443                   "The number of essential dynamics / flooding groups is not consistent.\n"
2444                   "There are %d ED groups in the .cpt file, but %zu in the .edi file!\n"
2445                   "Are you sure this is the correct .edi file?\n",
2446                   EDstate->nED,
2447                   ed.edpar.size());
2448     }
2449 }
2450
2451
2452 /* The edsamstate struct stores the information we need to make the ED group
2453  * whole again after restarts from a checkpoint file. Here we do the following:
2454  * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2455  * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2456  * c) in any case, for subsequent checkpoint writing, we set the pointers in
2457  * edsamstate to the x_old arrays, which contain the correct PBC representation of
2458  * all ED structures at the last time step. */
2459 static void init_edsamstate(const gmx_edsam& ed, edsamhistory_t* EDstate)
2460 {
2461     snew(EDstate->old_sref_p, EDstate->nED);
2462     snew(EDstate->old_sav_p, EDstate->nED);
2463
2464     /* If we did not read in a .cpt file, these arrays are not yet allocated */
2465     if (!EDstate->bFromCpt)
2466     {
2467         snew(EDstate->nref, EDstate->nED);
2468         snew(EDstate->nav, EDstate->nED);
2469     }
2470
2471     /* Loop over all ED/flooding data sets (usually only one, though) */
2472     for (int nr_edi = 0; nr_edi < EDstate->nED; nr_edi++)
2473     {
2474         const auto& edi = ed.edpar[nr_edi];
2475         /* We always need the last reference and average positions such that
2476          * in the next time step we can make the ED group whole again
2477          * if the atoms do not have the correct PBC representation */
2478         if (EDstate->bFromCpt)
2479         {
2480             /* Copy the last whole positions of reference and average group from .cpt */
2481             for (int i = 0; i < edi.sref.nr; i++)
2482             {
2483                 copy_rvec(EDstate->old_sref[nr_edi][i], edi.sref.x_old[i]);
2484             }
2485             for (int i = 0; i < edi.sav.nr; i++)
2486             {
2487                 copy_rvec(EDstate->old_sav[nr_edi][i], edi.sav.x_old[i]);
2488             }
2489         }
2490         else
2491         {
2492             EDstate->nref[nr_edi] = edi.sref.nr;
2493             EDstate->nav[nr_edi]  = edi.sav.nr;
2494         }
2495
2496         /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2497         EDstate->old_sref_p[nr_edi] = edi.sref.x_old;
2498         EDstate->old_sav_p[nr_edi]  = edi.sav.x_old;
2499     }
2500 }
2501
2502
2503 /* Adds 'buf' to 'str' */
2504 static void add_to_string(char** str, const char* buf)
2505 {
2506     int len;
2507
2508
2509     len = strlen(*str) + strlen(buf) + 1;
2510     srenew(*str, len);
2511     strcat(*str, buf);
2512 }
2513
2514
2515 static void add_to_string_aligned(char** str, const char* buf)
2516 {
2517     char buf_aligned[STRLEN];
2518
2519     sprintf(buf_aligned, EDcol_sfmt, buf);
2520     add_to_string(str, buf_aligned);
2521 }
2522
2523
2524 static void nice_legend(const char*** setname,
2525                         int*          nsets,
2526                         char**        LegendStr,
2527                         const char*   value,
2528                         const char*   unit,
2529                         char          EDgroupchar)
2530 {
2531     auto tmp = gmx::formatString("%c %s", EDgroupchar, value);
2532     add_to_string_aligned(LegendStr, tmp.c_str());
2533     tmp += gmx::formatString(" (%s)", unit);
2534     (*setname)[*nsets] = gmx_strdup(tmp.c_str());
2535     (*nsets)++;
2536 }
2537
2538
2539 static void nice_legend_evec(const char*** setname,
2540                              int*          nsets,
2541                              char**        LegendStr,
2542                              t_eigvec*     evec,
2543                              char          EDgroupChar,
2544                              const char*   EDtype)
2545 {
2546     int  i;
2547     char tmp[STRLEN];
2548
2549
2550     for (i = 0; i < evec->neig; i++)
2551     {
2552         sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2553         nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2554     }
2555 }
2556
2557
2558 /* Makes a legend for the xvg output file. Call on MASTER only! */
2559 static void write_edo_legend(gmx_edsam* ed, int nED, const gmx_output_env_t* oenv)
2560 {
2561     int          i;
2562     int          nr_edi, nsets, n_flood, n_edsam;
2563     const char** setname;
2564     char         buf[STRLEN];
2565     char*        LegendStr = nullptr;
2566
2567
2568     auto edi = ed->edpar.begin();
2569
2570     fprintf(ed->edo, "# Output will be written every %d step%s\n", edi->outfrq, edi->outfrq != 1 ? "s" : "");
2571
2572     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2573     {
2574         fprintf(ed->edo, "#\n");
2575         fprintf(ed->edo,
2576                 "# Summary of applied con/restraints for the ED group %c\n",
2577                 get_EDgroupChar(nr_edi, nED));
2578         fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2579         fprintf(ed->edo, "#    monitor  : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : "");
2580         fprintf(ed->edo,
2581                 "#    LINFIX   : %d vec%s\n",
2582                 edi->vecs.linfix.neig,
2583                 edi->vecs.linfix.neig != 1 ? "s" : "");
2584         fprintf(ed->edo,
2585                 "#    LINACC   : %d vec%s\n",
2586                 edi->vecs.linacc.neig,
2587                 edi->vecs.linacc.neig != 1 ? "s" : "");
2588         fprintf(ed->edo,
2589                 "#    RADFIX   : %d vec%s\n",
2590                 edi->vecs.radfix.neig,
2591                 edi->vecs.radfix.neig != 1 ? "s" : "");
2592         fprintf(ed->edo,
2593                 "#    RADACC   : %d vec%s\n",
2594                 edi->vecs.radacc.neig,
2595                 edi->vecs.radacc.neig != 1 ? "s" : "");
2596         fprintf(ed->edo,
2597                 "#    RADCON   : %d vec%s\n",
2598                 edi->vecs.radcon.neig,
2599                 edi->vecs.radcon.neig != 1 ? "s" : "");
2600         fprintf(ed->edo, "#    FLOODING : %d vec%s  ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : "");
2601
2602         if (edi->flood.vecs.neig)
2603         {
2604             /* If in any of the groups we find a flooding vector, flooding is turned on */
2605             ed->eEDtype = EssentialDynamicsType::Flooding;
2606
2607             /* Print what flavor of flooding we will do */
2608             if (0 == edi->flood.tau) /* constant flooding strength */
2609             {
2610                 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2611                 if (edi->flood.bHarmonic)
2612                 {
2613                     fprintf(ed->edo, ", harmonic");
2614                 }
2615             }
2616             else /* adaptive flooding */
2617             {
2618                 fprintf(ed->edo, ", adaptive");
2619             }
2620         }
2621         fprintf(ed->edo, "\n");
2622
2623         ++edi;
2624     }
2625
2626     /* Print a nice legend */
2627     snew(LegendStr, 1);
2628     LegendStr[0] = '\0';
2629     sprintf(buf, "#     %6s", "time");
2630     add_to_string(&LegendStr, buf);
2631
2632     /* Calculate the maximum number of columns we could end up with */
2633     edi   = ed->edpar.begin();
2634     nsets = 0;
2635     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2636     {
2637         nsets += 5 + edi->vecs.mon.neig + edi->vecs.linfix.neig + edi->vecs.linacc.neig
2638                  + edi->vecs.radfix.neig + edi->vecs.radacc.neig + edi->vecs.radcon.neig
2639                  + 6 * edi->flood.vecs.neig;
2640         ++edi;
2641     }
2642     snew(setname, nsets);
2643
2644     /* In the mdrun time step in a first function call (do_flood()) the flooding
2645      * forces are calculated and in a second function call (do_edsam()) the
2646      * ED constraints. To get a corresponding legend, we need to loop twice
2647      * over the edi groups and output first the flooding, then the ED part */
2648
2649     /* The flooding-related legend entries, if flooding is done */
2650     nsets = 0;
2651     if (EssentialDynamicsType::Flooding == ed->eEDtype)
2652     {
2653         edi = ed->edpar.begin();
2654         for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2655         {
2656             /* Always write out the projection on the flooding EVs. Of course, this can also
2657              * be achieved with the monitoring option in do_edsam() (if switched on by the
2658              * user), but in that case the positions need to be communicated in do_edsam(),
2659              * which is not necessary when doing flooding only. */
2660             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2661
2662             for (i = 0; i < edi->flood.vecs.neig; i++)
2663             {
2664                 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2665                 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2666
2667                 /* Output the current reference projection if it changes with time;
2668                  * this can happen when flooding is used as harmonic restraint */
2669                 if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
2670                 {
2671                     sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2672                     nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2673                 }
2674
2675                 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2676                 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2677                 {
2678                     sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2679                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2680                 }
2681
2682                 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2683                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2684
2685                 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2686                 {
2687                     sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2688                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2689                 }
2690
2691                 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2692                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2693             }
2694
2695             ++edi;
2696         } /* End of flooding-related legend entries */
2697     }
2698     n_flood = nsets;
2699
2700     /* Now the ED-related entries, if essential dynamics is done */
2701     edi = ed->edpar.begin();
2702     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2703     {
2704         if (bNeedDoEdsam(*edi)) /* Only print ED legend if at least one ED option is on */
2705         {
2706             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED));
2707
2708             /* Essential dynamics, projections on eigenvectors */
2709             nice_legend_evec(&setname,
2710                              &nsets,
2711                              &LegendStr,
2712                              &edi->vecs.mon,
2713                              get_EDgroupChar(nr_edi, nED),
2714                              "MON");
2715             nice_legend_evec(&setname,
2716                              &nsets,
2717                              &LegendStr,
2718                              &edi->vecs.linfix,
2719                              get_EDgroupChar(nr_edi, nED),
2720                              "LINFIX");
2721             nice_legend_evec(&setname,
2722                              &nsets,
2723                              &LegendStr,
2724                              &edi->vecs.linacc,
2725                              get_EDgroupChar(nr_edi, nED),
2726                              "LINACC");
2727             nice_legend_evec(&setname,
2728                              &nsets,
2729                              &LegendStr,
2730                              &edi->vecs.radfix,
2731                              get_EDgroupChar(nr_edi, nED),
2732                              "RADFIX");
2733             if (edi->vecs.radfix.neig)
2734             {
2735                 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2736             }
2737             nice_legend_evec(&setname,
2738                              &nsets,
2739                              &LegendStr,
2740                              &edi->vecs.radacc,
2741                              get_EDgroupChar(nr_edi, nED),
2742                              "RADACC");
2743             if (edi->vecs.radacc.neig)
2744             {
2745                 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2746             }
2747             nice_legend_evec(&setname,
2748                              &nsets,
2749                              &LegendStr,
2750                              &edi->vecs.radcon,
2751                              get_EDgroupChar(nr_edi, nED),
2752                              "RADCON");
2753             if (edi->vecs.radcon.neig)
2754             {
2755                 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2756             }
2757         }
2758         ++edi;
2759     } /* end of 'pure' essential dynamics legend entries */
2760     n_edsam = nsets - n_flood;
2761
2762     xvgr_legend(ed->edo, nsets, setname, oenv);
2763     sfree(setname);
2764
2765     fprintf(ed->edo,
2766             "#\n"
2767             "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2768             n_flood,
2769             1 == n_flood ? "" : "s",
2770             n_edsam,
2771             1 == n_edsam ? "" : "s");
2772     fprintf(ed->edo, "%s", LegendStr);
2773     sfree(LegendStr);
2774
2775     fflush(ed->edo);
2776 }
2777
2778
2779 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2780  * contained in the input file, creates a NULL terminated list of t_edpar structures */
2781 std::unique_ptr<gmx::EssentialDynamics> init_edsam(const gmx::MDLogger&        mdlog,
2782                                                    const char*                 ediFileName,
2783                                                    const char*                 edoFileName,
2784                                                    const gmx_mtop_t&           mtop,
2785                                                    const t_inputrec&           ir,
2786                                                    const t_commrec*            cr,
2787                                                    gmx::Constraints*           constr,
2788                                                    const t_state*              globalState,
2789                                                    ObservablesHistory*         oh,
2790                                                    const gmx_output_env_t*     oenv,
2791                                                    const gmx::StartingBehavior startingBehavior)
2792 {
2793     int    i, avindex;
2794     rvec*  x_pbc = nullptr; /* positions of the whole MD system with pbc removed  */
2795     rvec * xfit = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs  */
2796     rvec   fit_transvec;                      /* translation ... */
2797     matrix fit_rotmat;                        /* ... and rotation from fit to reference structure */
2798     rvec*  ref_x_old = nullptr;               /* helper pointer */
2799
2800
2801     if (MASTER(cr))
2802     {
2803         fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2804
2805         if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName))
2806         {
2807             gmx_fatal(FARGS,
2808                       "The checkpoint file you provided is from an essential dynamics or flooding\n"
2809                       "simulation. Please also set the .edi file on the command line with -ei.\n");
2810         }
2811     }
2812     GMX_LOG(mdlog.info)
2813             .asParagraph()
2814             .appendText(
2815                     "Activating essential dynamics via files passed to "
2816                     "gmx mdrun -edi will change in future to be files passed to "
2817                     "gmx grompp and the related .mdp options may change also.");
2818
2819     /* Open input and output files, allocate space for ED data structure */
2820     auto edHandle = ed_open(mtop.natoms, oh, ediFileName, edoFileName, startingBehavior, oenv, cr);
2821     auto ed       = edHandle->getLegacyED();
2822     GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2823     constr->saveEdsamPointer(ed);
2824
2825     /* Needed for initializing radacc radius in do_edsam */
2826     ed->bFirst = TRUE;
2827
2828     /* The input file is read by the master and the edi structures are
2829      * initialized here. Input is stored in ed->edpar. Then the edi
2830      * structures are transferred to the other nodes */
2831     if (MASTER(cr))
2832     {
2833         /* Initialization for every ED/flooding group. Flooding uses one edi group per
2834          * flooding vector, Essential dynamics can be applied to more than one structure
2835          * as well, but will be done in the order given in the edi file, so
2836          * expect different results for different order of edi file concatenation! */
2837         for (auto& edi : ed->edpar)
2838         {
2839             init_edi(mtop, &edi);
2840             init_flood(&edi, ed, ir.delta_t);
2841         }
2842     }
2843
2844     /* The master does the work here. The other nodes get the positions
2845      * not before dd_partition_system which is called after init_edsam */
2846     if (MASTER(cr))
2847     {
2848         edsamhistory_t* EDstate = oh->edsamHistory.get();
2849
2850         if (!EDstate->bFromCpt)
2851         {
2852             /* Remove PBC, make molecule(s) subject to ED whole. */
2853             snew(x_pbc, mtop.natoms);
2854             copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop.natoms);
2855             do_pbc_first_mtop(nullptr, ir.pbcType, globalState->box, &mtop, x_pbc);
2856         }
2857         /* Reset pointer to first ED data set which contains the actual ED data */
2858         auto edi = ed->edpar.begin();
2859         GMX_ASSERT(!ed->edpar.empty(), "Essential dynamics parameters is undefined");
2860
2861         /* Loop over all ED/flooding data sets (usually only one, though) */
2862         for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2863         {
2864             /* For multiple ED groups we use the output frequency that was specified
2865              * in the first set */
2866             if (nr_edi > 1)
2867             {
2868                 edi->outfrq = ed->edpar.begin()->outfrq;
2869             }
2870
2871             /* Extract the initial reference and average positions. When starting
2872              * from .cpt, these have already been read into sref.x_old
2873              * in init_edsamstate() */
2874             if (!EDstate->bFromCpt)
2875             {
2876                 /* If this is the first run (i.e. no checkpoint present) we assume
2877                  * that the starting positions give us the correct PBC representation */
2878                 for (i = 0; i < edi->sref.nr; i++)
2879                 {
2880                     copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2881                 }
2882
2883                 for (i = 0; i < edi->sav.nr; i++)
2884                 {
2885                     copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2886                 }
2887             }
2888
2889             /* Now we have the PBC-correct start positions of the reference and
2890                average structure. We copy that over to dummy arrays on which we
2891                can apply fitting to print out the RMSD. We srenew the memory since
2892                the size of the buffers is likely different for every ED group */
2893             srenew(xfit, edi->sref.nr);
2894             srenew(xstart, edi->sav.nr);
2895             if (edi->bRefEqAv)
2896             {
2897                 /* Reference indices are the same as average indices */
2898                 ref_x_old = edi->sav.x_old;
2899             }
2900             else
2901             {
2902                 ref_x_old = edi->sref.x_old;
2903             }
2904             copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2905             copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2906
2907             /* Make the fit to the REFERENCE structure, get translation and rotation */
2908             fit_to_reference(xfit, fit_transvec, fit_rotmat, &(*edi));
2909
2910             /* Output how well we fit to the reference at the start */
2911             translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2912             fprintf(stderr,
2913                     "ED: Initial RMSD from reference after fit = %f nm",
2914                     rmsd_from_structure(xfit, &edi->sref));
2915             if (EDstate->nED > 1)
2916             {
2917                 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2918             }
2919             fprintf(stderr, "\n");
2920
2921             /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2922             translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2923
2924             /* calculate initial projections */
2925             project(xstart, &(*edi));
2926
2927             /* For the target and origin structure both a reference (fit) and an
2928              * average structure can be provided in make_edi. If both structures
2929              * are the same, make_edi only stores one of them in the .edi file.
2930              * If they differ, first the fit and then the average structure is stored
2931              * in star (or sor), thus the number of entries in star/sor is
2932              * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2933              * the size of the average group. */
2934
2935             /* process target structure, if required */
2936             if (edi->star.nr > 0)
2937             {
2938                 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2939
2940                 /* get translation & rotation for fit of target structure to reference structure */
2941                 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, &(*edi));
2942                 /* do the fit */
2943                 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2944                 if (edi->star.nr == edi->sav.nr)
2945                 {
2946                     avindex = 0;
2947                 }
2948                 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2949                 {
2950                     /* The last sav.nr indices of the target structure correspond to
2951                      * the average structure, which must be projected */
2952                     avindex = edi->star.nr - edi->sav.nr;
2953                 }
2954                 rad_project(*edi, &edi->star.x[avindex], &edi->vecs.radcon);
2955             }
2956             else
2957             {
2958                 rad_project(*edi, xstart, &edi->vecs.radcon);
2959             }
2960
2961             /* process structure that will serve as origin of expansion circle */
2962             if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2963             {
2964                 fprintf(stderr,
2965                         "ED: Setting center of flooding potential (0 = average structure)\n");
2966             }
2967
2968             if (edi->sori.nr > 0)
2969             {
2970                 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2971
2972                 /* fit this structure to reference structure */
2973                 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, &(*edi));
2974                 /* do the fit */
2975                 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2976                 if (edi->sori.nr == edi->sav.nr)
2977                 {
2978                     avindex = 0;
2979                 }
2980                 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2981                 {
2982                     /* For the projection, we need the last sav.nr indices of sori */
2983                     avindex = edi->sori.nr - edi->sav.nr;
2984                 }
2985
2986                 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2987                 rad_project(*edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2988                 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
2989                 {
2990                     fprintf(stderr,
2991                             "ED: The ORIGIN structure will define the flooding potential "
2992                             "center.\n");
2993                     /* Set center of flooding potential to the ORIGIN structure */
2994                     rad_project(*edi, &edi->sori.x[avindex], &edi->flood.vecs);
2995                     /* We already know that no (moving) reference position was provided,
2996                      * therefore we can overwrite refproj[0]*/
2997                     copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
2998                 }
2999             }
3000             else /* No origin structure given */
3001             {
3002                 rad_project(*edi, xstart, &edi->vecs.radacc);
3003                 rad_project(*edi, xstart, &edi->vecs.radfix);
3004                 if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
3005                 {
3006                     if (edi->flood.bHarmonic)
3007                     {
3008                         fprintf(stderr,
3009                                 "ED: A (possibly changing) ref. projection will define the "
3010                                 "flooding potential center.\n");
3011                         for (i = 0; i < edi->flood.vecs.neig; i++)
3012                         {
3013                             edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
3014                         }
3015                     }
3016                     else
3017                     {
3018                         fprintf(stderr,
3019                                 "ED: The AVERAGE structure will define the flooding potential "
3020                                 "center.\n");
3021                         /* Set center of flooding potential to the center of the covariance matrix,
3022                          * i.e. the average structure, i.e. zero in the projected system */
3023                         for (i = 0; i < edi->flood.vecs.neig; i++)
3024                         {
3025                             edi->flood.vecs.refproj[i] = 0.0;
3026                         }
3027                     }
3028                 }
3029             }
3030             /* For convenience, output the center of the flooding potential for the eigenvectors */
3031             if ((EssentialDynamicsType::Flooding == ed->eEDtype) && (!edi->flood.bConstForce))
3032             {
3033                 for (i = 0; i < edi->flood.vecs.neig; i++)
3034                 {
3035                     fprintf(stdout,
3036                             "ED: EV %d flooding potential center: %11.4e",
3037                             edi->flood.vecs.ieig[i],
3038                             edi->flood.vecs.refproj[i]);
3039                     if (edi->flood.bHarmonic)
3040                     {
3041                         fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.referenceProjectionSlope[i]);
3042                     }
3043                     fprintf(stdout, "\n");
3044                 }
3045             }
3046
3047             /* set starting projections for linsam */
3048             rad_project(*edi, xstart, &edi->vecs.linacc);
3049             rad_project(*edi, xstart, &edi->vecs.linfix);
3050
3051             /* Prepare for the next edi data set: */
3052             ++edi;
3053         }
3054         /* Cleaning up on the master node: */
3055         if (!EDstate->bFromCpt)
3056         {
3057             sfree(x_pbc);
3058         }
3059         sfree(xfit);
3060         sfree(xstart);
3061
3062     } /* end of MASTER only section */
3063
3064     if (PAR(cr))
3065     {
3066         /* Broadcast the essential dynamics / flooding data to all nodes */
3067         broadcast_ed_data(cr, ed);
3068     }
3069     else
3070     {
3071         /* In the single-CPU case, point the local atom numbers pointers to the global
3072          * one, so that we can use the same notation in serial and parallel case: */
3073         /* Loop over all ED data sets (usually only one, though) */
3074         for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
3075         {
3076             edi->sref.anrs_loc = edi->sref.anrs;
3077             edi->sav.anrs_loc  = edi->sav.anrs;
3078             edi->star.anrs_loc = edi->star.anrs;
3079             edi->sori.anrs_loc = edi->sori.anrs;
3080             /* For the same reason as above, make a dummy c_ind array: */
3081             snew(edi->sav.c_ind, edi->sav.nr);
3082             /* Initialize the array */
3083             for (i = 0; i < edi->sav.nr; i++)
3084             {
3085                 edi->sav.c_ind[i] = i;
3086             }
3087             /* In the general case we will need a different-sized array for the reference indices: */
3088             if (!edi->bRefEqAv)
3089             {
3090                 snew(edi->sref.c_ind, edi->sref.nr);
3091                 for (i = 0; i < edi->sref.nr; i++)
3092                 {
3093                     edi->sref.c_ind[i] = i;
3094                 }
3095             }
3096             /* Point to the very same array in case of other structures: */
3097             edi->star.c_ind = edi->sav.c_ind;
3098             edi->sori.c_ind = edi->sav.c_ind;
3099             /* In the serial case, the local number of atoms is the global one: */
3100             edi->sref.nr_loc = edi->sref.nr;
3101             edi->sav.nr_loc  = edi->sav.nr;
3102             edi->star.nr_loc = edi->star.nr;
3103             edi->sori.nr_loc = edi->sori.nr;
3104         }
3105     }
3106
3107     /* Allocate space for ED buffer variables */
3108     /* Again, loop over ED data sets */
3109     for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
3110     {
3111         /* Allocate space for ED buffer variables */
3112         snew_bc(MASTER(cr), edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
3113         snew(edi->buf->do_edsam, 1);
3114
3115         /* Space for collective ED buffer variables */
3116
3117         /* Collective positions of atoms with the average indices */
3118         snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
3119         snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr); /* buffer for xcoll shifts */
3120         snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
3121         /* Collective positions of atoms with the reference indices */
3122         if (!edi->bRefEqAv)
3123         {
3124             snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
3125             snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr); /* To store the shifts in */
3126             snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
3127         }
3128
3129         /* Get memory for flooding forces */
3130         snew(edi->flood.forces_cartesian, edi->sav.nr);
3131     }
3132
3133     /* Flush the edo file so that the user can check some things
3134      * when the simulation has started */
3135     if (ed->edo)
3136     {
3137         fflush(ed->edo);
3138     }
3139
3140     return edHandle;
3141 }
3142
3143
3144 void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[], rvec v[], const matrix box, gmx_edsam* ed)
3145 {
3146     int    i, edinr, iupdate = 500;
3147     matrix rotmat;         /* rotation matrix */
3148     rvec   transvec;       /* translation vector */
3149     rvec   dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3150     real   dt_1;           /* 1/dt */
3151     struct t_do_edsam* buf;
3152     real     rmsdev    = -1; /* RMSD from reference structure prior to applying the constraints */
3153     gmx_bool bSuppress = FALSE; /* Write .xvg output file on master? */
3154
3155
3156     /* Check if ED sampling has to be performed */
3157     if (ed->eEDtype == EssentialDynamicsType::None)
3158     {
3159         return;
3160     }
3161
3162     dt_1 = 1.0 / ir->delta_t;
3163
3164     /* Loop over all ED groups (usually one) */
3165     edinr = 0;
3166     for (auto& edi : ed->edpar)
3167     {
3168         edinr++;
3169         if (bNeedDoEdsam(edi))
3170         {
3171
3172             buf = edi.buf->do_edsam;
3173
3174             if (ed->bFirst)
3175             {
3176                 /* initialize radacc radius for slope criterion */
3177                 buf->oldrad = calc_radius(edi.vecs.radacc);
3178             }
3179
3180             /* Copy the positions into buf->xc* arrays and after ED
3181              * feed back corrections to the official positions */
3182
3183             /* Broadcast the ED positions such that every node has all of them
3184              * Every node contributes its local positions xs and stores it in
3185              * the collective buf->xcoll array. Note that for edinr > 1
3186              * xs could already have been modified by an earlier ED */
3187
3188             communicate_group_positions(cr,
3189                                         buf->xcoll,
3190                                         buf->shifts_xcoll,
3191                                         buf->extra_shifts_xcoll,
3192                                         PAR(cr) ? buf->bUpdateShifts : TRUE,
3193                                         xs,
3194                                         edi.sav.nr,
3195                                         edi.sav.nr_loc,
3196                                         edi.sav.anrs_loc,
3197                                         edi.sav.c_ind,
3198                                         edi.sav.x_old,
3199                                         box);
3200
3201             /* Only assembly reference positions if their indices differ from the average ones */
3202             if (!edi.bRefEqAv)
3203             {
3204                 communicate_group_positions(cr,
3205                                             buf->xc_ref,
3206                                             buf->shifts_xc_ref,
3207                                             buf->extra_shifts_xc_ref,
3208                                             PAR(cr) ? buf->bUpdateShifts : TRUE,
3209                                             xs,
3210                                             edi.sref.nr,
3211                                             edi.sref.nr_loc,
3212                                             edi.sref.anrs_loc,
3213                                             edi.sref.c_ind,
3214                                             edi.sref.x_old,
3215                                             box);
3216             }
3217
3218             /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3219              * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3220              * set bUpdateShifts=TRUE in the parallel case. */
3221             buf->bUpdateShifts = FALSE;
3222
3223             /* Now all nodes have all of the ED positions in edi.sav->xcoll,
3224              * as well as the indices in edi.sav.anrs */
3225
3226             /* Fit the reference indices to the reference structure */
3227             if (edi.bRefEqAv)
3228             {
3229                 fit_to_reference(buf->xcoll, transvec, rotmat, &edi);
3230             }
3231             else
3232             {
3233                 fit_to_reference(buf->xc_ref, transvec, rotmat, &edi);
3234             }
3235
3236             /* Now apply the translation and rotation to the ED structure */
3237             translate_and_rotate(buf->xcoll, edi.sav.nr, transvec, rotmat);
3238
3239             /* Find out how well we fit to the reference (just for output steps) */
3240             if (do_per_step(step, edi.outfrq) && MASTER(cr))
3241             {
3242                 if (edi.bRefEqAv)
3243                 {
3244                     /* Indices of reference and average structures are identical,
3245                      * thus we can calculate the rmsd to SREF using xcoll */
3246                     rmsdev = rmsd_from_structure(buf->xcoll, &edi.sref);
3247                 }
3248                 else
3249                 {
3250                     /* We have to translate & rotate the reference atoms first */
3251                     translate_and_rotate(buf->xc_ref, edi.sref.nr, transvec, rotmat);
3252                     rmsdev = rmsd_from_structure(buf->xc_ref, &edi.sref);
3253                 }
3254             }
3255
3256             /* update radsam references, when required */
3257             if (do_per_step(step, edi.maxedsteps) && step >= edi.presteps)
3258             {
3259                 project(buf->xcoll, &edi);
3260                 rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3261                 rad_project(edi, buf->xcoll, &edi.vecs.radfix);
3262                 buf->oldrad = -1.e5;
3263             }
3264
3265             /* update radacc references, when required */
3266             if (do_per_step(step, iupdate) && step >= edi.presteps)
3267             {
3268                 edi.vecs.radacc.radius = calc_radius(edi.vecs.radacc);
3269                 if (edi.vecs.radacc.radius - buf->oldrad < edi.slope)
3270                 {
3271                     project(buf->xcoll, &edi);
3272                     rad_project(edi, buf->xcoll, &edi.vecs.radacc);
3273                     buf->oldrad = 0.0;
3274                 }
3275                 else
3276                 {
3277                     buf->oldrad = edi.vecs.radacc.radius;
3278                 }
3279             }
3280
3281             /* apply the constraints */
3282             if (step >= edi.presteps && ed_constraints(ed->eEDtype, edi))
3283             {
3284                 /* ED constraints should be applied already in the first MD step
3285                  * (which is step 0), therefore we pass step+1 to the routine */
3286                 ed_apply_constraints(buf->xcoll, &edi, step + 1 - ir->init_step);
3287             }
3288
3289             /* write to edo, when required */
3290             if (do_per_step(step, edi.outfrq))
3291             {
3292                 project(buf->xcoll, &edi);
3293                 if (MASTER(cr) && !bSuppress)
3294                 {
3295                     write_edo(edi, ed->edo, rmsdev);
3296                 }
3297             }
3298
3299             /* Copy back the positions unless monitoring only */
3300             if (ed_constraints(ed->eEDtype, edi))
3301             {
3302                 /* remove fitting */
3303                 rmfit(edi.sav.nr, buf->xcoll, transvec, rotmat);
3304
3305                 /* Copy the ED corrected positions into the coordinate array */
3306                 /* Each node copies its local part. In the serial case, nat_loc is the
3307                  * total number of ED atoms */
3308                 for (i = 0; i < edi.sav.nr_loc; i++)
3309                 {
3310                     /* Unshift local ED coordinate and store in x_unsh */
3311                     ed_unshift_single_coord(
3312                             box, buf->xcoll[edi.sav.c_ind[i]], buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh);
3313
3314                     /* dx is the ED correction to the positions: */
3315                     rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx);
3316
3317                     if (v != nullptr)
3318                     {
3319                         /* dv is the ED correction to the velocity: */
3320                         svmul(dt_1, dx, dv);
3321                         /* apply the velocity correction: */
3322                         rvec_inc(v[edi.sav.anrs_loc[i]], dv);
3323                     }
3324                     /* Finally apply the position correction due to ED: */
3325                     copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]);
3326                 }
3327             }
3328         } /* END of if ( bNeedDoEdsam(edi) ) */
3329
3330         /* Prepare for the next ED group */
3331
3332     } /* END of loop over ED groups */
3333
3334     ed->bFirst = FALSE;
3335 }