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