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