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