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