src/gromacs/essentialdynamics/edsam.c

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