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