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, 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                         break;
1514                 }
1515                 tvec->refproj[i]      = refproj_dum;
1516                 tvec->refproj0[i]     = refproj_dum;
1517                 tvec->refprojslope[i] = refprojslope_dum;
1518             }
1519             else /* Normal flooding */
1520             {
1521                 nscan = sscanf(line, max_ev_fmt_dlf, &idum, &rdum);
1522                 if (nscan != 2)
1523                 {
1524                     gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1525                 }
1526             }
1527             tvec->ieig[i]  = idum;
1528             tvec->stpsz[i] = rdum;
1529         } /* end of loop over eigenvectors */
1530
1531         for (i = 0; (i < tvec->neig); i++)
1532         {
1533             snew(tvec->vec[i], nr);
1534             scan_edvec(in, nr, tvec->vec[i]);
1535         }
1536     }
1537 }
1538
1539
1540 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1541 static void read_edvecs(FILE *in, int nr, t_edvecs *vecs)
1542 {
1543     gmx_bool bHaveReference = FALSE;
1544
1545
1546     read_edvec(in, nr, &vecs->mon, FALSE, &bHaveReference);
1547     read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1548     read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1549     read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1550     read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1551     read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1552 }
1553
1554
1555 /* Check if the same atom indices are used for reference and average positions */
1556 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1557 {
1558     int i;
1559
1560
1561     /* If the number of atoms differs between the two structures,
1562      * they cannot be identical */
1563     if (sref.nr != sav.nr)
1564     {
1565         return FALSE;
1566     }
1567
1568     /* Now that we know that both stuctures have the same number of atoms,
1569      * check if also the indices are identical */
1570     for (i = 0; i < sav.nr; i++)
1571     {
1572         if (sref.anrs[i] != sav.anrs[i])
1573         {
1574             return FALSE;
1575         }
1576     }
1577     fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1578
1579     return TRUE;
1580 }
1581
1582
1583 static int read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, const char *fn)
1584 {
1585     int       readmagic;
1586     const int magic = 670;
1587     gmx_bool  bEOF;
1588
1589     /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1590     gmx_bool bHaveReference = FALSE;
1591
1592
1593     /* the edi file is not free format, so expect problems if the input is corrupt. */
1594
1595     /* check the magic number */
1596     readmagic = read_edint(in, &bEOF);
1597     /* Check whether we have reached the end of the input file */
1598     if (bEOF)
1599     {
1600         return 0;
1601     }
1602
1603     if (readmagic != magic)
1604     {
1605         if (readmagic == 666 || readmagic == 667 || readmagic == 668)
1606         {
1607             gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1608         }
1609         else if (readmagic != 669)
1610         {
1611             gmx_fatal(FARGS, "Wrong magic number %d in %s", readmagic, fn);
1612         }
1613     }
1614
1615     /* check the number of atoms */
1616     edi->nini = read_edint(in, &bEOF);
1617     if (edi->nini != nr_mdatoms)
1618     {
1619         gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1620     }
1621
1622     /* Done checking. For the rest we blindly trust the input */
1623     edi->fitmas          = read_checked_edint(in, "FITMAS");
1624     edi->pcamas          = read_checked_edint(in, "ANALYSIS_MAS");
1625     edi->outfrq          = read_checked_edint(in, "OUTFRQ");
1626     edi->maxedsteps      = read_checked_edint(in, "MAXLEN");
1627     edi->slope           = read_checked_edreal(in, "SLOPECRIT");
1628
1629     edi->presteps        = read_checked_edint(in, "PRESTEPS");
1630     edi->flood.deltaF0   = read_checked_edreal(in, "DELTA_F0");
1631     edi->flood.deltaF    = read_checked_edreal(in, "INIT_DELTA_F");
1632     edi->flood.tau       = read_checked_edreal(in, "TAU");
1633     edi->flood.constEfl  = read_checked_edreal(in, "EFL_NULL");
1634     edi->flood.alpha2    = read_checked_edreal(in, "ALPHA2");
1635     edi->flood.kT        = read_checked_edreal(in, "KT");
1636     edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1637     if (readmagic > 669)
1638     {
1639         edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1640     }
1641     else
1642     {
1643         edi->flood.bConstForce = FALSE;
1644     }
1645     edi->sref.nr         = read_checked_edint(in, "NREF");
1646
1647     /* allocate space for reference positions and read them */
1648     snew(edi->sref.anrs, edi->sref.nr);
1649     snew(edi->sref.x, edi->sref.nr);
1650     snew(edi->sref.x_old, edi->sref.nr);
1651     edi->sref.sqrtm    = nullptr;
1652     read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1653
1654     /* average positions. they define which atoms will be used for ED sampling */
1655     edi->sav.nr = read_checked_edint(in, "NAV");
1656     snew(edi->sav.anrs, edi->sav.nr);
1657     snew(edi->sav.x, edi->sav.nr);
1658     snew(edi->sav.x_old, edi->sav.nr);
1659     read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1660
1661     /* Check if the same atom indices are used for reference and average positions */
1662     edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1663
1664     /* eigenvectors */
1665     read_edvecs(in, edi->sav.nr, &edi->vecs);
1666     read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
1667
1668     /* target positions */
1669     edi->star.nr = read_edint(in, &bEOF);
1670     if (edi->star.nr > 0)
1671     {
1672         snew(edi->star.anrs, edi->star.nr);
1673         snew(edi->star.x, edi->star.nr);
1674         edi->star.sqrtm    = nullptr;
1675         read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1676     }
1677
1678     /* positions defining origin of expansion circle */
1679     edi->sori.nr = read_edint(in, &bEOF);
1680     if (edi->sori.nr > 0)
1681     {
1682         if (bHaveReference)
1683         {
1684             /* Both an -ori structure and a at least one manual reference point have been
1685              * specified. That's ambiguous and probably not intentional. */
1686             gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1687                       "    point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1688         }
1689         snew(edi->sori.anrs, edi->sori.nr);
1690         snew(edi->sori.x, edi->sori.nr);
1691         edi->sori.sqrtm    = nullptr;
1692         read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1693     }
1694
1695     /* all done */
1696     return 1;
1697 }
1698
1699
1700
1701 /* Read in the edi input file. Note that it may contain several ED data sets which were
1702  * achieved by concatenating multiple edi files. The standard case would be a single ED
1703  * data set, though. */
1704 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1705 {
1706     FILE    *in;
1707     t_edpar *curr_edi, *last_edi;
1708     t_edpar *edi_read;
1709     int      edi_nr = 0;
1710
1711
1712     /* This routine is executed on the master only */
1713
1714     /* Open the .edi parameter input file */
1715     in = gmx_fio_fopen(fn, "r");
1716     fprintf(stderr, "ED: Reading edi file %s\n", fn);
1717
1718     /* Now read a sequence of ED input parameter sets from the edi file */
1719     curr_edi = edi;
1720     last_edi = edi;
1721     while (read_edi(in, curr_edi, nr_mdatoms, fn) )
1722     {
1723         edi_nr++;
1724
1725         /* Since we arrived within this while loop we know that there is still another data set to be read in */
1726         /* We need to allocate space for the data: */
1727         snew(edi_read, 1);
1728         /* Point the 'next_edi' entry to the next edi: */
1729         curr_edi->next_edi = edi_read;
1730         /* Keep the curr_edi pointer for the case that the next group is empty: */
1731         last_edi = curr_edi;
1732         /* Let's prepare to read in the next edi data set: */
1733         curr_edi = edi_read;
1734     }
1735     if (edi_nr == 0)
1736     {
1737         gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1738     }
1739
1740     /* Terminate the edi group list with a NULL pointer: */
1741     last_edi->next_edi = nullptr;
1742
1743     fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1744
1745     /* Close the .edi file again */
1746     gmx_fio_fclose(in);
1747
1748     return edi_nr;
1749 }
1750
1751
1752 struct t_fit_to_ref {
1753     rvec *xcopy;       /* Working copy of the positions in fit_to_reference */
1754 };
1755
1756 /* Fit the current positions to the reference positions
1757  * Do not actually do the fit, just return rotation and translation.
1758  * Note that the COM of the reference structure was already put into
1759  * the origin by init_edi. */
1760 static void fit_to_reference(rvec      *xcoll,    /* The positions to be fitted */
1761                              rvec       transvec, /* The translation vector */
1762                              matrix     rotmat,   /* The rotation matrix */
1763                              t_edpar   *edi)      /* Just needed for do_edfit */
1764 {
1765     rvec                 com;                     /* center of mass */
1766     int                  i;
1767     struct t_fit_to_ref *loc;
1768
1769
1770     /* Allocate memory the first time this routine is called for each edi group */
1771     if (nullptr == edi->buf->fit_to_ref)
1772     {
1773         snew(edi->buf->fit_to_ref, 1);
1774         snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1775     }
1776     loc = edi->buf->fit_to_ref;
1777
1778     /* We do not touch the original positions but work on a copy. */
1779     for (i = 0; i < edi->sref.nr; i++)
1780     {
1781         copy_rvec(xcoll[i], loc->xcopy[i]);
1782     }
1783
1784     /* Calculate the center of mass */
1785     get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1786
1787     transvec[XX] = -com[XX];
1788     transvec[YY] = -com[YY];
1789     transvec[ZZ] = -com[ZZ];
1790
1791     /* Subtract the center of mass from the copy */
1792     translate_x(loc->xcopy, edi->sref.nr, transvec);
1793
1794     /* Determine the rotation matrix */
1795     do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1796 }
1797
1798
1799 static void translate_and_rotate(rvec  *x,        /* The positions to be translated and rotated */
1800                                  int    nat,      /* How many positions are there? */
1801                                  rvec   transvec, /* The translation vector */
1802                                  matrix rotmat)   /* The rotation matrix */
1803 {
1804     /* Translation */
1805     translate_x(x, nat, transvec);
1806
1807     /* Rotation */
1808     rotate_x(x, nat, rotmat);
1809 }
1810
1811
1812 /* Gets the rms deviation of the positions to the structure s */
1813 /* fit_to_structure has to be called before calling this routine! */
1814 static real rmsd_from_structure(rvec           *x,  /* The positions under consideration */
1815                                 struct gmx_edx *s)  /* The structure from which the rmsd shall be computed */
1816 {
1817     real  rmsd = 0.0;
1818     int   i;
1819
1820
1821     for (i = 0; i < s->nr; i++)
1822     {
1823         rmsd += distance2(s->x[i], x[i]);
1824     }
1825
1826     rmsd /= (real) s->nr;
1827     rmsd  = sqrt(rmsd);
1828
1829     return rmsd;
1830 }
1831
1832
1833 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1834 {
1835     t_edpar *edi;
1836
1837
1838     if (ed->eEDtype != eEDnone)
1839     {
1840         /* Loop over ED groups */
1841         edi = ed->edpar;
1842         while (edi)
1843         {
1844             /* Local atoms of the reference structure (for fitting), need only be assembled
1845              * if their indices differ from the average ones */
1846             if (!edi->bRefEqAv)
1847             {
1848                 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1849                                             &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1850             }
1851
1852             /* Local atoms of the average structure (on these ED will be performed) */
1853             dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1854                                         &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1855
1856             /* Indicate that the ED shift vectors for this structure need to be updated
1857              * at the next call to communicate_group_positions, since obviously we are in a NS step */
1858             edi->buf->do_edsam->bUpdateShifts = TRUE;
1859
1860             /* Set the pointer to the next ED group (if any) */
1861             edi = edi->next_edi;
1862         }
1863     }
1864 }
1865
1866
1867 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1868 {
1869     int tx, ty, tz;
1870
1871
1872     tx = is[XX];
1873     ty = is[YY];
1874     tz = is[ZZ];
1875
1876     if (TRICLINIC(box))
1877     {
1878         xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1879         xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1880         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1881     }
1882     else
1883     {
1884         xu[XX] = x[XX]-tx*box[XX][XX];
1885         xu[YY] = x[YY]-ty*box[YY][YY];
1886         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1887     }
1888 }
1889
1890
1891 static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
1892 {
1893     int  i, j;
1894     real proj, add;
1895     rvec vec_dum;
1896
1897
1898     /* loop over linfix vectors */
1899     for (i = 0; i < edi->vecs.linfix.neig; i++)
1900     {
1901         /* calculate the projection */
1902         proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1903
1904         /* calculate the correction */
1905         add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1906
1907         /* apply the correction */
1908         add /= edi->sav.sqrtm[i];
1909         for (j = 0; j < edi->sav.nr; j++)
1910         {
1911             svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1912             rvec_inc(xcoll[j], vec_dum);
1913         }
1914     }
1915 }
1916
1917
1918 static void do_linacc(rvec *xcoll, t_edpar *edi)
1919 {
1920     int  i, j;
1921     real proj, add;
1922     rvec vec_dum;
1923
1924
1925     /* loop over linacc vectors */
1926     for (i = 0; i < edi->vecs.linacc.neig; i++)
1927     {
1928         /* calculate the projection */
1929         proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1930
1931         /* calculate the correction */
1932         add = 0.0;
1933         if (edi->vecs.linacc.stpsz[i] > 0.0)
1934         {
1935             if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1936             {
1937                 add = edi->vecs.linacc.refproj[i] - proj;
1938             }
1939         }
1940         if (edi->vecs.linacc.stpsz[i] < 0.0)
1941         {
1942             if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1943             {
1944                 add = edi->vecs.linacc.refproj[i] - proj;
1945             }
1946         }
1947
1948         /* apply the correction */
1949         add /= edi->sav.sqrtm[i];
1950         for (j = 0; j < edi->sav.nr; j++)
1951         {
1952             svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
1953             rvec_inc(xcoll[j], vec_dum);
1954         }
1955
1956         /* new positions will act as reference */
1957         edi->vecs.linacc.refproj[i] = proj + add;
1958     }
1959 }
1960
1961
1962 static void do_radfix(rvec *xcoll, t_edpar *edi)
1963 {
1964     int   i, j;
1965     real *proj, rad = 0.0, ratio;
1966     rvec  vec_dum;
1967
1968
1969     if (edi->vecs.radfix.neig == 0)
1970     {
1971         return;
1972     }
1973
1974     snew(proj, edi->vecs.radfix.neig);
1975
1976     /* loop over radfix vectors */
1977     for (i = 0; i < edi->vecs.radfix.neig; i++)
1978     {
1979         /* calculate the projections, radius */
1980         proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
1981         rad    += gmx::square(proj[i] - edi->vecs.radfix.refproj[i]);
1982     }
1983
1984     rad                      = sqrt(rad);
1985     ratio                    = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
1986     edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
1987
1988     /* loop over radfix vectors */
1989     for (i = 0; i < edi->vecs.radfix.neig; i++)
1990     {
1991         proj[i] -= edi->vecs.radfix.refproj[i];
1992
1993         /* apply the correction */
1994         proj[i] /= edi->sav.sqrtm[i];
1995         proj[i] *= ratio;
1996         for (j = 0; j < edi->sav.nr; j++)
1997         {
1998             svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
1999             rvec_inc(xcoll[j], vec_dum);
2000         }
2001     }
2002
2003     sfree(proj);
2004 }
2005
2006
2007 static void do_radacc(rvec *xcoll, t_edpar *edi)
2008 {
2009     int   i, j;
2010     real *proj, rad = 0.0, ratio = 0.0;
2011     rvec  vec_dum;
2012
2013
2014     if (edi->vecs.radacc.neig == 0)
2015     {
2016         return;
2017     }
2018
2019     snew(proj, edi->vecs.radacc.neig);
2020
2021     /* loop over radacc vectors */
2022     for (i = 0; i < edi->vecs.radacc.neig; i++)
2023     {
2024         /* calculate the projections, radius */
2025         proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2026         rad    += gmx::square(proj[i] - edi->vecs.radacc.refproj[i]);
2027     }
2028     rad = sqrt(rad);
2029
2030     /* only correct when radius decreased */
2031     if (rad < edi->vecs.radacc.radius)
2032     {
2033         ratio = edi->vecs.radacc.radius/rad - 1.0;
2034     }
2035     else
2036     {
2037         edi->vecs.radacc.radius = rad;
2038     }
2039
2040     /* loop over radacc vectors */
2041     for (i = 0; i < edi->vecs.radacc.neig; i++)
2042     {
2043         proj[i] -= edi->vecs.radacc.refproj[i];
2044
2045         /* apply the correction */
2046         proj[i] /= edi->sav.sqrtm[i];
2047         proj[i] *= ratio;
2048         for (j = 0; j < edi->sav.nr; j++)
2049         {
2050             svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2051             rvec_inc(xcoll[j], vec_dum);
2052         }
2053     }
2054     sfree(proj);
2055 }
2056
2057
2058 struct t_do_radcon {
2059     real *proj;
2060 };
2061
2062 static void do_radcon(rvec *xcoll, t_edpar *edi)
2063 {
2064     int                 i, j;
2065     real                rad = 0.0, ratio = 0.0;
2066     struct t_do_radcon *loc;
2067     gmx_bool            bFirst;
2068     rvec                vec_dum;
2069
2070
2071     if (edi->buf->do_radcon != nullptr)
2072     {
2073         bFirst = FALSE;
2074     }
2075     else
2076     {
2077         bFirst = TRUE;
2078         snew(edi->buf->do_radcon, 1);
2079     }
2080     loc = edi->buf->do_radcon;
2081
2082     if (edi->vecs.radcon.neig == 0)
2083     {
2084         return;
2085     }
2086
2087     if (bFirst)
2088     {
2089         snew(loc->proj, edi->vecs.radcon.neig);
2090     }
2091
2092     /* loop over radcon vectors */
2093     for (i = 0; i < edi->vecs.radcon.neig; i++)
2094     {
2095         /* calculate the projections, radius */
2096         loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2097         rad         += gmx::square(loc->proj[i] - edi->vecs.radcon.refproj[i]);
2098     }
2099     rad = sqrt(rad);
2100     /* only correct when radius increased */
2101     if (rad > edi->vecs.radcon.radius)
2102     {
2103         ratio = edi->vecs.radcon.radius/rad - 1.0;
2104
2105         /* loop over radcon vectors */
2106         for (i = 0; i < edi->vecs.radcon.neig; i++)
2107         {
2108             /* apply the correction */
2109             loc->proj[i] -= edi->vecs.radcon.refproj[i];
2110             loc->proj[i] /= edi->sav.sqrtm[i];
2111             loc->proj[i] *= ratio;
2112
2113             for (j = 0; j < edi->sav.nr; j++)
2114             {
2115                 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2116                 rvec_inc(xcoll[j], vec_dum);
2117             }
2118         }
2119
2120     }
2121     else
2122     {
2123         edi->vecs.radcon.radius = rad;
2124     }
2125
2126 }
2127
2128
2129 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
2130 {
2131     int i;
2132
2133
2134     /* subtract the average positions */
2135     for (i = 0; i < edi->sav.nr; i++)
2136     {
2137         rvec_dec(xcoll[i], edi->sav.x[i]);
2138     }
2139
2140     /* apply the constraints */
2141     if (step >= 0)
2142     {
2143         do_linfix(xcoll, edi, step);
2144     }
2145     do_linacc(xcoll, edi);
2146     if (step >= 0)
2147     {
2148         do_radfix(xcoll, edi);
2149     }
2150     do_radacc(xcoll, edi);
2151     do_radcon(xcoll, edi);
2152
2153     /* add back the average positions */
2154     for (i = 0; i < edi->sav.nr; i++)
2155     {
2156         rvec_inc(xcoll[i], edi->sav.x[i]);
2157     }
2158 }
2159
2160
2161 /* Write out the projections onto the eigenvectors. The order of output
2162  * corresponds to ed_output_legend() */
2163 static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
2164 {
2165     int i;
2166
2167
2168     /* Output how well we fit to the reference structure */
2169     fprintf(fp, EDcol_ffmt, rmsd);
2170
2171     for (i = 0; i < edi->vecs.mon.neig; i++)
2172     {
2173         fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
2174     }
2175
2176     for (i = 0; i < edi->vecs.linfix.neig; i++)
2177     {
2178         fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
2179     }
2180
2181     for (i = 0; i < edi->vecs.linacc.neig; i++)
2182     {
2183         fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
2184     }
2185
2186     for (i = 0; i < edi->vecs.radfix.neig; i++)
2187     {
2188         fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
2189     }
2190     if (edi->vecs.radfix.neig)
2191     {
2192         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
2193     }
2194
2195     for (i = 0; i < edi->vecs.radacc.neig; i++)
2196     {
2197         fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
2198     }
2199     if (edi->vecs.radacc.neig)
2200     {
2201         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
2202     }
2203
2204     for (i = 0; i < edi->vecs.radcon.neig; i++)
2205     {
2206         fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
2207     }
2208     if (edi->vecs.radcon.neig)
2209     {
2210         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
2211     }
2212 }
2213
2214 /* Returns if any constraints are switched on */
2215 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2216 {
2217     if (edtype == eEDedsam || edtype == eEDflood)
2218     {
2219         return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2220                 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2221                 edi->vecs.radcon.neig);
2222     }
2223     return 0;
2224 }
2225
2226
2227 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2228  * umbrella sampling simulations. */
2229 static void copyEvecReference(t_eigvec* floodvecs)
2230 {
2231     int i;
2232
2233
2234     if (nullptr == floodvecs->refproj0)
2235     {
2236         snew(floodvecs->refproj0, floodvecs->neig);
2237     }
2238
2239     for (i = 0; i < floodvecs->neig; i++)
2240     {
2241         floodvecs->refproj0[i] = floodvecs->refproj[i];
2242     }
2243 }
2244
2245
2246 /* Call on MASTER only. Check whether the essential dynamics / flooding
2247  * groups of the checkpoint file are consistent with the provided .edi file. */
2248 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamhistory_t *EDstate)
2249 {
2250     t_edpar *edi = nullptr;    /* points to a single edi data set */
2251     int      edinum;
2252
2253
2254     if (nullptr == EDstate->nref || nullptr == EDstate->nav)
2255     {
2256         gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2257                   "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2258                   "it must also continue with/without ED constraints when checkpointing.\n"
2259                   "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2260                   "from without a checkpoint.\n");
2261     }
2262
2263     edi    = ed->edpar;
2264     edinum = 0;
2265     while (edi != nullptr)
2266     {
2267         /* Check number of atoms in the reference and average structures */
2268         if (EDstate->nref[edinum] != edi->sref.nr)
2269         {
2270             gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2271                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2272                       get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2273         }
2274         if (EDstate->nav[edinum] != edi->sav.nr)
2275         {
2276             gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2277                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2278                       get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2279         }
2280         edi = edi->next_edi;
2281         edinum++;
2282     }
2283
2284     if (edinum != EDstate->nED)
2285     {
2286         gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2287                   "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2288                   "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2289     }
2290 }
2291
2292
2293 /* The edsamstate struct stores the information we need to make the ED group
2294  * whole again after restarts from a checkpoint file. Here we do the following:
2295  * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2296  * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2297  * c) in any case, for subsequent checkpoint writing, we set the pointers in
2298  * edsamstate to the x_old arrays, which contain the correct PBC representation of
2299  * all ED structures at the last time step. */
2300 static void init_edsamstate(gmx_edsam_t ed, edsamhistory_t *EDstate)
2301 {
2302     int      i, nr_edi;
2303     t_edpar *edi;
2304
2305
2306     snew(EDstate->old_sref_p, EDstate->nED);
2307     snew(EDstate->old_sav_p, EDstate->nED);
2308
2309     /* If we did not read in a .cpt file, these arrays are not yet allocated */
2310     if (!EDstate->bFromCpt)
2311     {
2312         snew(EDstate->nref, EDstate->nED);
2313         snew(EDstate->nav, EDstate->nED);
2314     }
2315
2316     /* Loop over all ED/flooding data sets (usually only one, though) */
2317     edi = ed->edpar;
2318     for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2319     {
2320         /* We always need the last reference and average positions such that
2321          * in the next time step we can make the ED group whole again
2322          * if the atoms do not have the correct PBC representation */
2323         if (EDstate->bFromCpt)
2324         {
2325             /* Copy the last whole positions of reference and average group from .cpt */
2326             for (i = 0; i < edi->sref.nr; i++)
2327             {
2328                 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2329             }
2330             for (i = 0; i < edi->sav.nr; i++)
2331             {
2332                 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2333             }
2334         }
2335         else
2336         {
2337             EDstate->nref[nr_edi-1] = edi->sref.nr;
2338             EDstate->nav [nr_edi-1] = edi->sav.nr;
2339         }
2340
2341         /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2342         EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2343         EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2344
2345         edi = edi->next_edi;
2346     }
2347 }
2348
2349
2350 /* Adds 'buf' to 'str' */
2351 static void add_to_string(char **str, char *buf)
2352 {
2353     int len;
2354
2355
2356     len = strlen(*str) + strlen(buf) + 1;
2357     srenew(*str, len);
2358     strcat(*str, buf);
2359 }
2360
2361
2362 static void add_to_string_aligned(char **str, char *buf)
2363 {
2364     char buf_aligned[STRLEN];
2365
2366     sprintf(buf_aligned, EDcol_sfmt, buf);
2367     add_to_string(str, buf_aligned);
2368 }
2369
2370
2371 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, const char *value, const char *unit, char EDgroupchar)
2372 {
2373     char tmp[STRLEN], tmp2[STRLEN];
2374
2375
2376     sprintf(tmp, "%c %s", EDgroupchar, value);
2377     add_to_string_aligned(LegendStr, tmp);
2378     sprintf(tmp2, "%s (%s)", tmp, unit);
2379     (*setname)[*nsets] = gmx_strdup(tmp2);
2380     (*nsets)++;
2381 }
2382
2383
2384 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2385 {
2386     int  i;
2387     char tmp[STRLEN];
2388
2389
2390     for (i = 0; i < evec->neig; i++)
2391     {
2392         sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2393         nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2394     }
2395 }
2396
2397
2398 /* Makes a legend for the xvg output file. Call on MASTER only! */
2399 static void write_edo_legend(gmx_edsam_t ed, int nED, const gmx_output_env_t *oenv)
2400 {
2401     t_edpar     *edi = nullptr;
2402     int          i;
2403     int          nr_edi, nsets, n_flood, n_edsam;
2404     const char **setname;
2405     char         buf[STRLEN];
2406     char        *LegendStr = nullptr;
2407
2408
2409     edi         = ed->edpar;
2410
2411     fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2412
2413     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2414     {
2415         fprintf(ed->edo, "#\n");
2416         fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2417         fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2418         fprintf(ed->edo, "#    monitor  : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig    != 1 ? "s" : "");
2419         fprintf(ed->edo, "#    LINFIX   : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2420         fprintf(ed->edo, "#    LINACC   : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2421         fprintf(ed->edo, "#    RADFIX   : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2422         fprintf(ed->edo, "#    RADACC   : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2423         fprintf(ed->edo, "#    RADCON   : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2424         fprintf(ed->edo, "#    FLOODING : %d vec%s  ", edi->flood.vecs.neig, edi->flood.vecs.neig  != 1 ? "s" : "");
2425
2426         if (edi->flood.vecs.neig)
2427         {
2428             /* If in any of the groups we find a flooding vector, flooding is turned on */
2429             ed->eEDtype = eEDflood;
2430
2431             /* Print what flavor of flooding we will do */
2432             if (0 == edi->flood.tau) /* constant flooding strength */
2433             {
2434                 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2435                 if (edi->flood.bHarmonic)
2436                 {
2437                     fprintf(ed->edo, ", harmonic");
2438                 }
2439             }
2440             else /* adaptive flooding */
2441             {
2442                 fprintf(ed->edo, ", adaptive");
2443             }
2444         }
2445         fprintf(ed->edo, "\n");
2446
2447         edi = edi->next_edi;
2448     }
2449
2450     /* Print a nice legend */
2451     snew(LegendStr, 1);
2452     LegendStr[0] = '\0';
2453     sprintf(buf, "#     %6s", "time");
2454     add_to_string(&LegendStr, buf);
2455
2456     /* Calculate the maximum number of columns we could end up with */
2457     edi     = ed->edpar;
2458     nsets   = 0;
2459     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2460     {
2461         nsets += 5 +edi->vecs.mon.neig
2462             +edi->vecs.linfix.neig
2463             +edi->vecs.linacc.neig
2464             +edi->vecs.radfix.neig
2465             +edi->vecs.radacc.neig
2466             +edi->vecs.radcon.neig
2467             + 6*edi->flood.vecs.neig;
2468         edi = edi->next_edi;
2469     }
2470     snew(setname, nsets);
2471
2472     /* In the mdrun time step in a first function call (do_flood()) the flooding
2473      * forces are calculated and in a second function call (do_edsam()) the
2474      * ED constraints. To get a corresponding legend, we need to loop twice
2475      * over the edi groups and output first the flooding, then the ED part */
2476
2477     /* The flooding-related legend entries, if flooding is done */
2478     nsets = 0;
2479     if (eEDflood == ed->eEDtype)
2480     {
2481         edi   = ed->edpar;
2482         for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2483         {
2484             /* Always write out the projection on the flooding EVs. Of course, this can also
2485              * be achieved with the monitoring option in do_edsam() (if switched on by the
2486              * user), but in that case the positions need to be communicated in do_edsam(),
2487              * which is not necessary when doing flooding only. */
2488             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2489
2490             for (i = 0; i < edi->flood.vecs.neig; i++)
2491             {
2492                 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2493                 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2494
2495                 /* Output the current reference projection if it changes with time;
2496                  * this can happen when flooding is used as harmonic restraint */
2497                 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
2498                 {
2499                     sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2500                     nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2501                 }
2502
2503                 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2504                 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2505                 {
2506                     sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2507                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2508                 }
2509
2510                 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2511                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2512
2513                 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2514                 {
2515                     sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2516                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2517                 }
2518
2519                 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2520                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2521             }
2522
2523             edi = edi->next_edi;
2524         } /* End of flooding-related legend entries */
2525     }
2526     n_flood = nsets;
2527
2528     /* Now the ED-related entries, if essential dynamics is done */
2529     edi         = ed->edpar;
2530     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2531     {
2532         if (bNeedDoEdsam(edi))  /* Only print ED legend if at least one ED option is on */
2533         {
2534             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2535
2536             /* Essential dynamics, projections on eigenvectors */
2537             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON"   );
2538             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2539             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2540             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2541             if (edi->vecs.radfix.neig)
2542             {
2543                 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2544             }
2545             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2546             if (edi->vecs.radacc.neig)
2547             {
2548                 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2549             }
2550             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2551             if (edi->vecs.radcon.neig)
2552             {
2553                 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2554             }
2555         }
2556         edi = edi->next_edi;
2557     } /* end of 'pure' essential dynamics legend entries */
2558     n_edsam = nsets - n_flood;
2559
2560     xvgr_legend(ed->edo, nsets, setname, oenv);
2561     sfree(setname);
2562
2563     fprintf(ed->edo, "#\n"
2564             "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2565             n_flood, 1 == n_flood ? "" : "s",
2566             n_edsam, 1 == n_edsam ? "" : "s");
2567     fprintf(ed->edo, "%s", LegendStr);
2568     sfree(LegendStr);
2569
2570     fflush(ed->edo);
2571 }
2572
2573
2574 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2575  * contained in the input file, creates a NULL terminated list of t_edpar structures */
2576 gmx_edsam_t init_edsam(
2577         const char             *ediFileName,
2578         const char             *edoFileName,
2579         const gmx_mtop_t       *mtop,
2580         const t_inputrec       *ir,
2581         const t_commrec        *cr,
2582         gmx::Constraints       *constr,
2583         const t_state          *globalState,
2584         ObservablesHistory     *oh,
2585         const gmx_output_env_t *oenv,
2586         gmx_bool                bAppend)
2587 {
2588     t_edpar *edi = nullptr;                       /* points to a single edi data set */
2589     int      i, avindex;
2590     int      nED    = 0;                          /* Number of ED data sets */
2591     rvec    *x_pbc  = nullptr;                    /* positions of the whole MD system with pbc removed  */
2592     rvec    *xfit   = nullptr, *xstart = nullptr; /* dummy arrays to determine initial RMSDs  */
2593     rvec     fit_transvec;                        /* translation ... */
2594     matrix   fit_rotmat;                          /* ... and rotation from fit to reference structure */
2595     rvec    *ref_x_old = nullptr;                 /* helper pointer */
2596
2597
2598     if (MASTER(cr))
2599     {
2600         fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2601
2602         if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName) )
2603         {
2604             gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or flooding\n"
2605                       "simulation. Please also set the .edi file on the command line with -ei.\n");
2606         }
2607     }
2608
2609     /* Open input and output files, allocate space for ED data structure */
2610     gmx_edsam_t ed = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr);
2611     GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
2612     constr->saveEdsamPointer(ed);
2613
2614     /* Needed for initializing radacc radius in do_edsam */
2615     ed->bFirst = TRUE;
2616
2617     /* The input file is read by the master and the edi structures are
2618      * initialized here. Input is stored in ed->edpar. Then the edi
2619      * structures are transferred to the other nodes */
2620     if (MASTER(cr))
2621     {
2622         /* Initialization for every ED/flooding group. Flooding uses one edi group per
2623          * flooding vector, Essential dynamics can be applied to more than one structure
2624          * as well, but will be done in the order given in the edi file, so
2625          * expect different results for different order of edi file concatenation! */
2626         edi = ed->edpar;
2627         while (edi != nullptr)
2628         {
2629             init_edi(mtop, edi);
2630             init_flood(edi, ed, ir->delta_t);
2631             edi = edi->next_edi;
2632         }
2633     }
2634
2635     /* The master does the work here. The other nodes get the positions
2636      * not before dd_partition_system which is called after init_edsam */
2637     if (MASTER(cr))
2638     {
2639         edsamhistory_t *EDstate = oh->edsamHistory.get();
2640
2641         if (!EDstate->bFromCpt)
2642         {
2643             /* Remove PBC, make molecule(s) subject to ED whole. */
2644             snew(x_pbc, mtop->natoms);
2645             copy_rvecn(as_rvec_array(globalState->x.data()), x_pbc, 0, mtop->natoms);
2646             do_pbc_first_mtop(nullptr, ir->ePBC, globalState->box, mtop, x_pbc);
2647         }
2648         /* Reset pointer to first ED data set which contains the actual ED data */
2649         edi = ed->edpar;
2650         GMX_ASSERT(nullptr != edi,
2651                    "The pointer to the essential dynamics parameters is undefined");
2652
2653         nED = EDstate->nED;
2654         /* Loop over all ED/flooding data sets (usually only one, though) */
2655         for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2656         {
2657             /* For multiple ED groups we use the output frequency that was specified
2658              * in the first set */
2659             if (nr_edi > 1)
2660             {
2661                 edi->outfrq = ed->edpar->outfrq;
2662             }
2663
2664             /* Extract the initial reference and average positions. When starting
2665              * from .cpt, these have already been read into sref.x_old
2666              * in init_edsamstate() */
2667             if (!EDstate->bFromCpt)
2668             {
2669                 /* If this is the first run (i.e. no checkpoint present) we assume
2670                  * that the starting positions give us the correct PBC representation */
2671                 for (i = 0; i < edi->sref.nr; i++)
2672                 {
2673                     copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2674                 }
2675
2676                 for (i = 0; i < edi->sav.nr; i++)
2677                 {
2678                     copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2679                 }
2680             }
2681
2682             /* Now we have the PBC-correct start positions of the reference and
2683                average structure. We copy that over to dummy arrays on which we
2684                can apply fitting to print out the RMSD. We srenew the memory since
2685                the size of the buffers is likely different for every ED group */
2686             srenew(xfit, edi->sref.nr );
2687             srenew(xstart, edi->sav.nr  );
2688             if (edi->bRefEqAv)
2689             {
2690                 /* Reference indices are the same as average indices */
2691                 ref_x_old = edi->sav.x_old;
2692             }
2693             else
2694             {
2695                 ref_x_old = edi->sref.x_old;
2696             }
2697             copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2698             copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2699
2700             /* Make the fit to the REFERENCE structure, get translation and rotation */
2701             fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2702
2703             /* Output how well we fit to the reference at the start */
2704             translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2705             fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2706                     rmsd_from_structure(xfit, &edi->sref));
2707             if (EDstate->nED > 1)
2708             {
2709                 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2710             }
2711             fprintf(stderr, "\n");
2712
2713             /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2714             translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2715
2716             /* calculate initial projections */
2717             project(xstart, edi);
2718
2719             /* For the target and origin structure both a reference (fit) and an
2720              * average structure can be provided in make_edi. If both structures
2721              * are the same, make_edi only stores one of them in the .edi file.
2722              * If they differ, first the fit and then the average structure is stored
2723              * in star (or sor), thus the number of entries in star/sor is
2724              * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2725              * the size of the average group. */
2726
2727             /* process target structure, if required */
2728             if (edi->star.nr > 0)
2729             {
2730                 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2731
2732                 /* get translation & rotation for fit of target structure to reference structure */
2733                 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2734                 /* do the fit */
2735                 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2736                 if (edi->star.nr == edi->sav.nr)
2737                 {
2738                     avindex = 0;
2739                 }
2740                 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2741                 {
2742                     /* The last sav.nr indices of the target structure correspond to
2743                      * the average structure, which must be projected */
2744                     avindex = edi->star.nr - edi->sav.nr;
2745                 }
2746                 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
2747             }
2748             else
2749             {
2750                 rad_project(edi, xstart, &edi->vecs.radcon);
2751             }
2752
2753             /* process structure that will serve as origin of expansion circle */
2754             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2755             {
2756                 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2757             }
2758
2759             if (edi->sori.nr > 0)
2760             {
2761                 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2762
2763                 /* fit this structure to reference structure */
2764                 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2765                 /* do the fit */
2766                 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2767                 if (edi->sori.nr == edi->sav.nr)
2768                 {
2769                     avindex = 0;
2770                 }
2771                 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2772                 {
2773                     /* For the projection, we need the last sav.nr indices of sori */
2774                     avindex = edi->sori.nr - edi->sav.nr;
2775                 }
2776
2777                 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2778                 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2779                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2780                 {
2781                     fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2782                     /* Set center of flooding potential to the ORIGIN structure */
2783                     rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
2784                     /* We already know that no (moving) reference position was provided,
2785                      * therefore we can overwrite refproj[0]*/
2786                     copyEvecReference(&edi->flood.vecs);
2787                 }
2788             }
2789             else /* No origin structure given */
2790             {
2791                 rad_project(edi, xstart, &edi->vecs.radacc);
2792                 rad_project(edi, xstart, &edi->vecs.radfix);
2793                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2794                 {
2795                     if (edi->flood.bHarmonic)
2796                     {
2797                         fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2798                         for (i = 0; i < edi->flood.vecs.neig; i++)
2799                         {
2800                             edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2801                         }
2802                     }
2803                     else
2804                     {
2805                         fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2806                         /* Set center of flooding potential to the center of the covariance matrix,
2807                          * i.e. the average structure, i.e. zero in the projected system */
2808                         for (i = 0; i < edi->flood.vecs.neig; i++)
2809                         {
2810                             edi->flood.vecs.refproj[i] = 0.0;
2811                         }
2812                     }
2813                 }
2814             }
2815             /* For convenience, output the center of the flooding potential for the eigenvectors */
2816             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2817             {
2818                 for (i = 0; i < edi->flood.vecs.neig; i++)
2819                 {
2820                     fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2821                     if (edi->flood.bHarmonic)
2822                     {
2823                         fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2824                     }
2825                     fprintf(stdout, "\n");
2826                 }
2827             }
2828
2829             /* set starting projections for linsam */
2830             rad_project(edi, xstart, &edi->vecs.linacc);
2831             rad_project(edi, xstart, &edi->vecs.linfix);
2832
2833             /* Prepare for the next edi data set: */
2834             edi = edi->next_edi;
2835         }
2836         /* Cleaning up on the master node: */
2837         if (!EDstate->bFromCpt)
2838         {
2839             sfree(x_pbc);
2840         }
2841         sfree(xfit);
2842         sfree(xstart);
2843
2844     } /* end of MASTER only section */
2845
2846     if (PAR(cr))
2847     {
2848         /* First let everybody know how many ED data sets to expect */
2849         gmx_bcast(sizeof(nED), &nED, cr);
2850         /* Broadcast the essential dynamics / flooding data to all nodes */
2851         broadcast_ed_data(cr, ed, nED);
2852     }
2853     else
2854     {
2855         /* In the single-CPU case, point the local atom numbers pointers to the global
2856          * one, so that we can use the same notation in serial and parallel case: */
2857         /* Loop over all ED data sets (usually only one, though) */
2858         edi = ed->edpar;
2859         for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
2860         {
2861             edi->sref.anrs_loc = edi->sref.anrs;
2862             edi->sav.anrs_loc  = edi->sav.anrs;
2863             edi->star.anrs_loc = edi->star.anrs;
2864             edi->sori.anrs_loc = edi->sori.anrs;
2865             /* For the same reason as above, make a dummy c_ind array: */
2866             snew(edi->sav.c_ind, edi->sav.nr);
2867             /* Initialize the array */
2868             for (i = 0; i < edi->sav.nr; i++)
2869             {
2870                 edi->sav.c_ind[i] = i;
2871             }
2872             /* In the general case we will need a different-sized array for the reference indices: */
2873             if (!edi->bRefEqAv)
2874             {
2875                 snew(edi->sref.c_ind, edi->sref.nr);
2876                 for (i = 0; i < edi->sref.nr; i++)
2877                 {
2878                     edi->sref.c_ind[i] = i;
2879                 }
2880             }
2881             /* Point to the very same array in case of other structures: */
2882             edi->star.c_ind = edi->sav.c_ind;
2883             edi->sori.c_ind = edi->sav.c_ind;
2884             /* In the serial case, the local number of atoms is the global one: */
2885             edi->sref.nr_loc = edi->sref.nr;
2886             edi->sav.nr_loc  = edi->sav.nr;
2887             edi->star.nr_loc = edi->star.nr;
2888             edi->sori.nr_loc = edi->sori.nr;
2889
2890             /* An on we go to the next ED group */
2891             edi = edi->next_edi;
2892         }
2893     }
2894
2895     /* Allocate space for ED buffer variables */
2896     /* Again, loop over ED data sets */
2897     edi = ed->edpar;
2898     for (int nr_edi = 1; nr_edi <= nED; nr_edi++)
2899     {
2900         /* Allocate space for ED buffer variables */
2901         snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
2902         snew(edi->buf->do_edsam, 1);
2903
2904         /* Space for collective ED buffer variables */
2905
2906         /* Collective positions of atoms with the average indices */
2907         snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2908         snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr);            /* buffer for xcoll shifts */
2909         snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2910         /* Collective positions of atoms with the reference indices */
2911         if (!edi->bRefEqAv)
2912         {
2913             snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2914             snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr);       /* To store the shifts in */
2915             snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2916         }
2917
2918         /* Get memory for flooding forces */
2919         snew(edi->flood.forces_cartesian, edi->sav.nr);
2920
2921         if (gmx_debug_at)
2922         {
2923             /* Dump it all into one file per process */
2924             dump_edi(edi, cr, nr_edi);
2925         }
2926
2927         /* Next ED group */
2928         edi = edi->next_edi;
2929     }
2930
2931     /* Flush the edo file so that the user can check some things
2932      * when the simulation has started */
2933     if (ed->edo)
2934     {
2935         fflush(ed->edo);
2936     }
2937
2938     return ed;
2939 }
2940
2941
2942 void do_edsam(const t_inputrec *ir,
2943               gmx_int64_t       step,
2944               const t_commrec  *cr,
2945               rvec              xs[],
2946               rvec              v[],
2947               matrix            box,
2948               gmx_edsam_t       ed)
2949 {
2950     int                i, edinr, iupdate = 500;
2951     matrix             rotmat;         /* rotation matrix */
2952     rvec               transvec;       /* translation vector */
2953     rvec               dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
2954     real               dt_1;           /* 1/dt */
2955     struct t_do_edsam *buf;
2956     t_edpar           *edi;
2957     real               rmsdev    = -1;    /* RMSD from reference structure prior to applying the constraints */
2958     gmx_bool           bSuppress = FALSE; /* Write .xvg output file on master? */
2959
2960
2961     /* Check if ED sampling has to be performed */
2962     if (ed->eEDtype == eEDnone)
2963     {
2964         return;
2965     }
2966
2967     dt_1 = 1.0/ir->delta_t;
2968
2969     /* Loop over all ED groups (usually one) */
2970     edi   = ed->edpar;
2971     edinr = 0;
2972     while (edi != nullptr)
2973     {
2974         edinr++;
2975         if (bNeedDoEdsam(edi))
2976         {
2977
2978             buf = edi->buf->do_edsam;
2979
2980             if (ed->bFirst)
2981             {
2982                 /* initialize radacc radius for slope criterion */
2983                 buf->oldrad = calc_radius(&edi->vecs.radacc);
2984             }
2985
2986             /* Copy the positions into buf->xc* arrays and after ED
2987              * feed back corrections to the official positions */
2988
2989             /* Broadcast the ED positions such that every node has all of them
2990              * Every node contributes its local positions xs and stores it in
2991              * the collective buf->xcoll array. Note that for edinr > 1
2992              * xs could already have been modified by an earlier ED */
2993
2994             communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
2995                                         edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old,  box);
2996
2997             /* Only assembly reference positions if their indices differ from the average ones */
2998             if (!edi->bRefEqAv)
2999             {
3000                 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3001                                             edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3002             }
3003
3004             /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3005              * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3006              * set bUpdateShifts=TRUE in the parallel case. */
3007             buf->bUpdateShifts = FALSE;
3008
3009             /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3010              * as well as the indices in edi->sav.anrs */
3011
3012             /* Fit the reference indices to the reference structure */
3013             if (edi->bRefEqAv)
3014             {
3015                 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3016             }
3017             else
3018             {
3019                 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3020             }
3021
3022             /* Now apply the translation and rotation to the ED structure */
3023             translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3024
3025             /* Find out how well we fit to the reference (just for output steps) */
3026             if (do_per_step(step, edi->outfrq) && MASTER(cr))
3027             {
3028                 if (edi->bRefEqAv)
3029                 {
3030                     /* Indices of reference and average structures are identical,
3031                      * thus we can calculate the rmsd to SREF using xcoll */
3032                     rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3033                 }
3034                 else
3035                 {
3036                     /* We have to translate & rotate the reference atoms first */
3037                     translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3038                     rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3039                 }
3040             }
3041
3042             /* update radsam references, when required */
3043             if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3044             {
3045                 project(buf->xcoll, edi);
3046                 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3047                 rad_project(edi, buf->xcoll, &edi->vecs.radfix);
3048                 buf->oldrad = -1.e5;
3049             }
3050
3051             /* update radacc references, when required */
3052             if (do_per_step(step, iupdate) && step >= edi->presteps)
3053             {
3054                 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
3055                 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3056                 {
3057                     project(buf->xcoll, edi);
3058                     rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3059                     buf->oldrad = 0.0;
3060                 }
3061                 else
3062                 {
3063                     buf->oldrad = edi->vecs.radacc.radius;
3064                 }
3065             }
3066
3067             /* apply the constraints */
3068             if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3069             {
3070                 /* ED constraints should be applied already in the first MD step
3071                  * (which is step 0), therefore we pass step+1 to the routine */
3072                 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3073             }
3074
3075             /* write to edo, when required */
3076             if (do_per_step(step, edi->outfrq))
3077             {
3078                 project(buf->xcoll, edi);
3079                 if (MASTER(cr) && !bSuppress)
3080                 {
3081                     write_edo(edi, ed->edo, rmsdev);
3082                 }
3083             }
3084
3085             /* Copy back the positions unless monitoring only */
3086             if (ed_constraints(ed->eEDtype, edi))
3087             {
3088                 /* remove fitting */
3089                 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3090
3091                 /* Copy the ED corrected positions into the coordinate array */
3092                 /* Each node copies its local part. In the serial case, nat_loc is the
3093                  * total number of ED atoms */
3094                 for (i = 0; i < edi->sav.nr_loc; i++)
3095                 {
3096                     /* Unshift local ED coordinate and store in x_unsh */
3097                     ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3098                                             buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3099
3100                     /* dx is the ED correction to the positions: */
3101                     rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3102
3103                     if (v != nullptr)
3104                     {
3105                         /* dv is the ED correction to the velocity: */
3106                         svmul(dt_1, dx, dv);
3107                         /* apply the velocity correction: */
3108                         rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3109                     }
3110                     /* Finally apply the position correction due to ED: */
3111                     copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3112                 }
3113             }
3114         } /* END of if ( bNeedDoEdsam(edi) ) */
3115
3116         /* Prepare for the next ED group */
3117         edi = edi->next_edi;
3118
3119     } /* END of loop over ED groups */
3120
3121     ed->bFirst = FALSE;
3122 }
3123
3124 void done_ed(gmx_edsam_t *ed)
3125 {
3126     if (*ed)
3127     {
3128         /* ed->edo is opened sometimes with xvgropen, sometimes with
3129          * gmx_fio_fopen, so we use the least common denominator for
3130          * closing. */
3131         gmx_fio_fclose((*ed)->edo);
3132     }
3133
3134     /* TODO deallocate ed and set pointer to NULL */
3135 }