src/gromacs/essentialdynamics/edsam.cpp

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