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