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