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 <time.h>
  43 #include "typedefs.h"
  44 #include "string2.h"
  45 #include "smalloc.h"
  46 #include "names.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "txtdump.h"
  49 #include "vec.h"
  50 #include <time.h>
  51 #include "nrnb.h"
  52 #include "mshift.h"
  53 #include "mdrun.h"
  54 #include "update.h"
  55 #include "physics.h"
  56 #include "nrjac.h"
  57 #include "mtop_util.h"
  58 #include "gromacs/essentialdynamics/edsam.h"
  59 #include "gromacs/fileio/gmxfio.h"
  60 #include "xvgr.h"
  61 #include "gromacs/mdlib/groupcoord.h"
  62
  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         curr_edi = edi_read;
1798     }
1799     if (edi_nr == 0)
1800     {
1801         gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1802     }
1803
1804     /* Terminate the edi group list with a NULL pointer: */
1805     last_edi->next_edi = NULL;
1806
1807     fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1808
1809     /* Close the .edi file again */
1810     gmx_fio_fclose(in);
1811
1812     return edi_nr;
1813 }
1814
1815
1816 struct t_fit_to_ref {
1817     rvec *xcopy;       /* Working copy of the positions in fit_to_reference */
1818 };
1819
1820 /* Fit the current positions to the reference positions
1821  * Do not actually do the fit, just return rotation and translation.
1822  * Note that the COM of the reference structure was already put into
1823  * the origin by init_edi. */
1824 static void fit_to_reference(rvec      *xcoll,    /* The positions to be fitted */
1825                              rvec       transvec, /* The translation vector */
1826                              matrix     rotmat,   /* The rotation matrix */
1827                              t_edpar   *edi)      /* Just needed for do_edfit */
1828 {
1829     rvec                 com;                     /* center of mass */
1830     int                  i;
1831     struct t_fit_to_ref *loc;
1832
1833
1834     /* Allocate memory the first time this routine is called for each edi group */
1835     if (NULL == edi->buf->fit_to_ref)
1836     {
1837         snew(edi->buf->fit_to_ref, 1);
1838         snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1839     }
1840     loc = edi->buf->fit_to_ref;
1841
1842     /* We do not touch the original positions but work on a copy. */
1843     for (i = 0; i < edi->sref.nr; i++)
1844     {
1845         copy_rvec(xcoll[i], loc->xcopy[i]);
1846     }
1847
1848     /* Calculate the center of mass */
1849     get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1850
1851     transvec[XX] = -com[XX];
1852     transvec[YY] = -com[YY];
1853     transvec[ZZ] = -com[ZZ];
1854
1855     /* Subtract the center of mass from the copy */
1856     translate_x(loc->xcopy, edi->sref.nr, transvec);
1857
1858     /* Determine the rotation matrix */
1859     do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1860 }
1861
1862
1863 static void translate_and_rotate(rvec  *x,        /* The positions to be translated and rotated */
1864                                  int    nat,      /* How many positions are there? */
1865                                  rvec   transvec, /* The translation vector */
1866                                  matrix rotmat)   /* The rotation matrix */
1867 {
1868     /* Translation */
1869     translate_x(x, nat, transvec);
1870
1871     /* Rotation */
1872     rotate_x(x, nat, rotmat);
1873 }
1874
1875
1876 /* Gets the rms deviation of the positions to the structure s */
1877 /* fit_to_structure has to be called before calling this routine! */
1878 static real rmsd_from_structure(rvec           *x,  /* The positions under consideration */
1879                                 struct gmx_edx *s)  /* The structure from which the rmsd shall be computed */
1880 {
1881     real  rmsd = 0.0;
1882     int   i;
1883
1884
1885     for (i = 0; i < s->nr; i++)
1886     {
1887         rmsd += distance2(s->x[i], x[i]);
1888     }
1889
1890     rmsd /= (real) s->nr;
1891     rmsd  = sqrt(rmsd);
1892
1893     return rmsd;
1894 }
1895
1896
1897 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1898 {
1899     t_edpar *edi;
1900
1901
1902     if (ed->eEDtype != eEDnone)
1903     {
1904         /* Loop over ED groups */
1905         edi = ed->edpar;
1906         while (edi)
1907         {
1908             /* Local atoms of the reference structure (for fitting), need only be assembled
1909              * if their indices differ from the average ones */
1910             if (!edi->bRefEqAv)
1911             {
1912                 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1913                                             &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1914             }
1915
1916             /* Local atoms of the average structure (on these ED will be performed) */
1917             dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1918                                         &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1919
1920             /* Indicate that the ED shift vectors for this structure need to be updated
1921              * at the next call to communicate_group_positions, since obviously we are in a NS step */
1922             edi->buf->do_edsam->bUpdateShifts = TRUE;
1923
1924             /* Set the pointer to the next ED group (if any) */
1925             edi = edi->next_edi;
1926         }
1927     }
1928 }
1929
1930
1931 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1932 {
1933     int tx, ty, tz;
1934
1935
1936     tx = is[XX];
1937     ty = is[YY];
1938     tz = is[ZZ];
1939
1940     if (TRICLINIC(box))
1941     {
1942         xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1943         xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1944         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1945     }
1946     else
1947     {
1948         xu[XX] = x[XX]-tx*box[XX][XX];
1949         xu[YY] = x[YY]-ty*box[YY][YY];
1950         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1951     }
1952 }
1953
1954
1955 static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
1956 {
1957     int  i, j;
1958     real proj, add;
1959     rvec vec_dum;
1960
1961
1962     /* loop over linfix vectors */
1963     for (i = 0; i < edi->vecs.linfix.neig; i++)
1964     {
1965         /* calculate the projection */
1966         proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1967
1968         /* calculate the correction */
1969         add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1970
1971         /* apply the correction */
1972         add /= edi->sav.sqrtm[i];
1973         for (j = 0; j < edi->sav.nr; j++)
1974         {
1975             svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1976             rvec_inc(xcoll[j], vec_dum);
1977         }
1978     }
1979 }
1980
1981
1982 static void do_linacc(rvec *xcoll, t_edpar *edi)
1983 {
1984     int  i, j;
1985     real proj, add;
1986     rvec vec_dum;
1987
1988
1989     /* loop over linacc vectors */
1990     for (i = 0; i < edi->vecs.linacc.neig; i++)
1991     {
1992         /* calculate the projection */
1993         proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1994
1995         /* calculate the correction */
1996         add = 0.0;
1997         if (edi->vecs.linacc.stpsz[i] > 0.0)
1998         {
1999             if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
2000             {
2001                 add = edi->vecs.linacc.refproj[i] - proj;
2002             }
2003         }
2004         if (edi->vecs.linacc.stpsz[i] < 0.0)
2005         {
2006             if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
2007             {
2008                 add = edi->vecs.linacc.refproj[i] - proj;
2009             }
2010         }
2011
2012         /* apply the correction */
2013         add /= edi->sav.sqrtm[i];
2014         for (j = 0; j < edi->sav.nr; j++)
2015         {
2016             svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
2017             rvec_inc(xcoll[j], vec_dum);
2018         }
2019
2020         /* new positions will act as reference */
2021         edi->vecs.linacc.refproj[i] = proj + add;
2022     }
2023 }
2024
2025
2026 static void do_radfix(rvec *xcoll, t_edpar *edi)
2027 {
2028     int   i, j;
2029     real *proj, rad = 0.0, ratio;
2030     rvec  vec_dum;
2031
2032
2033     if (edi->vecs.radfix.neig == 0)
2034     {
2035         return;
2036     }
2037
2038     snew(proj, edi->vecs.radfix.neig);
2039
2040     /* loop over radfix vectors */
2041     for (i = 0; i < edi->vecs.radfix.neig; i++)
2042     {
2043         /* calculate the projections, radius */
2044         proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
2045         rad    += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
2046     }
2047
2048     rad                      = sqrt(rad);
2049     ratio                    = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2050     edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2051
2052     /* loop over radfix vectors */
2053     for (i = 0; i < edi->vecs.radfix.neig; i++)
2054     {
2055         proj[i] -= edi->vecs.radfix.refproj[i];
2056
2057         /* apply the correction */
2058         proj[i] /= edi->sav.sqrtm[i];
2059         proj[i] *= ratio;
2060         for (j = 0; j < edi->sav.nr; j++)
2061         {
2062             svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2063             rvec_inc(xcoll[j], vec_dum);
2064         }
2065     }
2066
2067     sfree(proj);
2068 }
2069
2070
2071 static void do_radacc(rvec *xcoll, t_edpar *edi)
2072 {
2073     int   i, j;
2074     real *proj, rad = 0.0, ratio = 0.0;
2075     rvec  vec_dum;
2076
2077
2078     if (edi->vecs.radacc.neig == 0)
2079     {
2080         return;
2081     }
2082
2083     snew(proj, edi->vecs.radacc.neig);
2084
2085     /* loop over radacc vectors */
2086     for (i = 0; i < edi->vecs.radacc.neig; i++)
2087     {
2088         /* calculate the projections, radius */
2089         proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2090         rad    += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
2091     }
2092     rad = sqrt(rad);
2093
2094     /* only correct when radius decreased */
2095     if (rad < edi->vecs.radacc.radius)
2096     {
2097         ratio = edi->vecs.radacc.radius/rad - 1.0;
2098         rad   = edi->vecs.radacc.radius;
2099     }
2100     else
2101     {
2102         edi->vecs.radacc.radius = rad;
2103     }
2104
2105     /* loop over radacc vectors */
2106     for (i = 0; i < edi->vecs.radacc.neig; i++)
2107     {
2108         proj[i] -= edi->vecs.radacc.refproj[i];
2109
2110         /* apply the correction */
2111         proj[i] /= edi->sav.sqrtm[i];
2112         proj[i] *= ratio;
2113         for (j = 0; j < edi->sav.nr; j++)
2114         {
2115             svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2116             rvec_inc(xcoll[j], vec_dum);
2117         }
2118     }
2119     sfree(proj);
2120 }
2121
2122
2123 struct t_do_radcon {
2124     real *proj;
2125 };
2126
2127 static void do_radcon(rvec *xcoll, t_edpar *edi)
2128 {
2129     int                 i, j;
2130     real                rad = 0.0, ratio = 0.0;
2131     struct t_do_radcon *loc;
2132     gmx_bool            bFirst;
2133     rvec                vec_dum;
2134
2135
2136     if (edi->buf->do_radcon != NULL)
2137     {
2138         bFirst = FALSE;
2139         loc    = edi->buf->do_radcon;
2140     }
2141     else
2142     {
2143         bFirst = TRUE;
2144         snew(edi->buf->do_radcon, 1);
2145     }
2146     loc = edi->buf->do_radcon;
2147
2148     if (edi->vecs.radcon.neig == 0)
2149     {
2150         return;
2151     }
2152
2153     if (bFirst)
2154     {
2155         snew(loc->proj, edi->vecs.radcon.neig);
2156     }
2157
2158     /* loop over radcon vectors */
2159     for (i = 0; i < edi->vecs.radcon.neig; i++)
2160     {
2161         /* calculate the projections, radius */
2162         loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2163         rad         += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2164     }
2165     rad = sqrt(rad);
2166     /* only correct when radius increased */
2167     if (rad > edi->vecs.radcon.radius)
2168     {
2169         ratio = edi->vecs.radcon.radius/rad - 1.0;
2170
2171         /* loop over radcon vectors */
2172         for (i = 0; i < edi->vecs.radcon.neig; i++)
2173         {
2174             /* apply the correction */
2175             loc->proj[i] -= edi->vecs.radcon.refproj[i];
2176             loc->proj[i] /= edi->sav.sqrtm[i];
2177             loc->proj[i] *= ratio;
2178
2179             for (j = 0; j < edi->sav.nr; j++)
2180             {
2181                 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2182                 rvec_inc(xcoll[j], vec_dum);
2183             }
2184         }
2185     }
2186     else
2187     {
2188         edi->vecs.radcon.radius = rad;
2189     }
2190
2191     if (rad != edi->vecs.radcon.radius)
2192     {
2193         rad = 0.0;
2194         for (i = 0; i < edi->vecs.radcon.neig; i++)
2195         {
2196             /* calculate the projections, radius */
2197             loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2198             rad         += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2199         }
2200         rad = sqrt(rad);
2201     }
2202 }
2203
2204
2205 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_int64_t step)
2206 {
2207     int i;
2208
2209
2210     /* subtract the average positions */
2211     for (i = 0; i < edi->sav.nr; i++)
2212     {
2213         rvec_dec(xcoll[i], edi->sav.x[i]);
2214     }
2215
2216     /* apply the constraints */
2217     if (step >= 0)
2218     {
2219         do_linfix(xcoll, edi, step);
2220     }
2221     do_linacc(xcoll, edi);
2222     if (step >= 0)
2223     {
2224         do_radfix(xcoll, edi);
2225     }
2226     do_radacc(xcoll, edi);
2227     do_radcon(xcoll, edi);
2228
2229     /* add back the average positions */
2230     for (i = 0; i < edi->sav.nr; i++)
2231     {
2232         rvec_inc(xcoll[i], edi->sav.x[i]);
2233     }
2234 }
2235
2236
2237 /* Write out the projections onto the eigenvectors. The order of output
2238  * corresponds to ed_output_legend() */
2239 static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
2240 {
2241     int i;
2242
2243
2244     /* Output how well we fit to the reference structure */
2245     fprintf(fp, EDcol_ffmt, rmsd);
2246
2247     for (i = 0; i < edi->vecs.mon.neig; i++)
2248     {
2249         fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
2250     }
2251
2252     for (i = 0; i < edi->vecs.linfix.neig; i++)
2253     {
2254         fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
2255     }
2256
2257     for (i = 0; i < edi->vecs.linacc.neig; i++)
2258     {
2259         fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
2260     }
2261
2262     for (i = 0; i < edi->vecs.radfix.neig; i++)
2263     {
2264         fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
2265     }
2266     if (edi->vecs.radfix.neig)
2267     {
2268         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
2269     }
2270
2271     for (i = 0; i < edi->vecs.radacc.neig; i++)
2272     {
2273         fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
2274     }
2275     if (edi->vecs.radacc.neig)
2276     {
2277         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
2278     }
2279
2280     for (i = 0; i < edi->vecs.radcon.neig; i++)
2281     {
2282         fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
2283     }
2284     if (edi->vecs.radcon.neig)
2285     {
2286         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
2287     }
2288 }
2289
2290 /* Returns if any constraints are switched on */
2291 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2292 {
2293     if (edtype == eEDedsam || edtype == eEDflood)
2294     {
2295         return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2296                 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2297                 edi->vecs.radcon.neig);
2298     }
2299     return 0;
2300 }
2301
2302
2303 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2304  * umbrella sampling simulations. */
2305 static void copyEvecReference(t_eigvec* floodvecs)
2306 {
2307     int i;
2308
2309
2310     if (NULL == floodvecs->refproj0)
2311     {
2312         snew(floodvecs->refproj0, floodvecs->neig);
2313     }
2314
2315     for (i = 0; i < floodvecs->neig; i++)
2316     {
2317         floodvecs->refproj0[i] = floodvecs->refproj[i];
2318     }
2319 }
2320
2321
2322 /* Call on MASTER only. Check whether the essential dynamics / flooding
2323  * groups of the checkpoint file are consistent with the provided .edi file. */
2324 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate)
2325 {
2326     t_edpar *edi = NULL;    /* points to a single edi data set */
2327     int      edinum;
2328
2329
2330     if (NULL == EDstate->nref || NULL == EDstate->nav)
2331     {
2332         gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2333                   "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2334                   "it must also continue with/without ED constraints when checkpointing.\n"
2335                   "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2336                   "from without a checkpoint.\n");
2337     }
2338
2339     edi    = ed->edpar;
2340     edinum = 0;
2341     while (edi != NULL)
2342     {
2343         /* Check number of atoms in the reference and average structures */
2344         if (EDstate->nref[edinum] != edi->sref.nr)
2345         {
2346             gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2347                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2348                       get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2349         }
2350         if (EDstate->nav[edinum] != edi->sav.nr)
2351         {
2352             gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2353                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2354                       get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2355         }
2356         edi = edi->next_edi;
2357         edinum++;
2358     }
2359
2360     if (edinum != EDstate->nED)
2361     {
2362         gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2363                   "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2364                   "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2365     }
2366 }
2367
2368
2369 /* The edsamstate struct stores the information we need to make the ED group
2370  * whole again after restarts from a checkpoint file. Here we do the following:
2371  * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2372  * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2373  * c) in any case, for subsequent checkpoint writing, we set the pointers in
2374  * edsamstate to the x_old arrays, which contain the correct PBC representation of
2375  * all ED structures at the last time step. */
2376 static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate)
2377 {
2378     int      i, nr_edi;
2379     t_edpar *edi;
2380
2381
2382     snew(EDstate->old_sref_p, EDstate->nED);
2383     snew(EDstate->old_sav_p, EDstate->nED);
2384
2385     /* If we did not read in a .cpt file, these arrays are not yet allocated */
2386     if (!EDstate->bFromCpt)
2387     {
2388         snew(EDstate->nref, EDstate->nED);
2389         snew(EDstate->nav, EDstate->nED);
2390     }
2391
2392     /* Loop over all ED/flooding data sets (usually only one, though) */
2393     edi = ed->edpar;
2394     for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2395     {
2396         /* We always need the last reference and average positions such that
2397          * in the next time step we can make the ED group whole again
2398          * if the atoms do not have the correct PBC representation */
2399         if (EDstate->bFromCpt)
2400         {
2401             /* Copy the last whole positions of reference and average group from .cpt */
2402             for (i = 0; i < edi->sref.nr; i++)
2403             {
2404                 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2405             }
2406             for (i = 0; i < edi->sav.nr; i++)
2407             {
2408                 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2409             }
2410         }
2411         else
2412         {
2413             EDstate->nref[nr_edi-1] = edi->sref.nr;
2414             EDstate->nav [nr_edi-1] = edi->sav.nr;
2415         }
2416
2417         /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2418         EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2419         EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2420
2421         edi = edi->next_edi;
2422     }
2423 }
2424
2425
2426 /* Adds 'buf' to 'str' */
2427 static void add_to_string(char **str, char *buf)
2428 {
2429     int len;
2430
2431
2432     len = strlen(*str) + strlen(buf) + 1;
2433     srenew(*str, len);
2434     strcat(*str, buf);
2435 }
2436
2437
2438 static void add_to_string_aligned(char **str, char *buf)
2439 {
2440     char buf_aligned[STRLEN];
2441
2442     sprintf(buf_aligned, EDcol_sfmt, buf);
2443     add_to_string(str, buf_aligned);
2444 }
2445
2446
2447 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, char *value, char *unit, char EDgroupchar)
2448 {
2449     char tmp[STRLEN], tmp2[STRLEN];
2450
2451
2452     sprintf(tmp, "%c %s", EDgroupchar, value);
2453     add_to_string_aligned(LegendStr, tmp);
2454     sprintf(tmp2, "%s (%s)", tmp, unit);
2455     (*setname)[*nsets] = strdup(tmp2);
2456     (*nsets)++;
2457 }
2458
2459
2460 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2461 {
2462     int  i;
2463     char tmp[STRLEN];
2464
2465
2466     for (i = 0; i < evec->neig; i++)
2467     {
2468         sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2469         nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2470     }
2471 }
2472
2473
2474 /* Makes a legend for the xvg output file. Call on MASTER only! */
2475 static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv)
2476 {
2477     t_edpar     *edi = NULL;
2478     int          i;
2479     int          nr_edi, nsets, n_flood, n_edsam;
2480     const char **setname;
2481     char         buf[STRLEN];
2482     char        *LegendStr = NULL;
2483
2484
2485     edi         = ed->edpar;
2486
2487     fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2488
2489     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2490     {
2491         fprintf(ed->edo, "#\n");
2492         fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2493         fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2494         fprintf(ed->edo, "#    monitor  : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig    != 1 ? "s" : "");
2495         fprintf(ed->edo, "#    LINFIX   : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2496         fprintf(ed->edo, "#    LINACC   : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2497         fprintf(ed->edo, "#    RADFIX   : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2498         fprintf(ed->edo, "#    RADACC   : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2499         fprintf(ed->edo, "#    RADCON   : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2500         fprintf(ed->edo, "#    FLOODING : %d vec%s  ", edi->flood.vecs.neig, edi->flood.vecs.neig  != 1 ? "s" : "");
2501
2502         if (edi->flood.vecs.neig)
2503         {
2504             /* If in any of the groups we find a flooding vector, flooding is turned on */
2505             ed->eEDtype = eEDflood;
2506
2507             /* Print what flavor of flooding we will do */
2508             if (0 == edi->flood.tau) /* constant flooding strength */
2509             {
2510                 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2511                 if (edi->flood.bHarmonic)
2512                 {
2513                     fprintf(ed->edo, ", harmonic");
2514                 }
2515             }
2516             else /* adaptive flooding */
2517             {
2518                 fprintf(ed->edo, ", adaptive");
2519             }
2520         }
2521         fprintf(ed->edo, "\n");
2522
2523         edi = edi->next_edi;
2524     }
2525
2526     /* Print a nice legend */
2527     snew(LegendStr, 1);
2528     LegendStr[0] = '\0';
2529     sprintf(buf, "#     %6s", "time");
2530     add_to_string(&LegendStr, buf);
2531
2532     /* Calculate the maximum number of columns we could end up with */
2533     edi     = ed->edpar;
2534     nsets   = 0;
2535     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2536     {
2537         nsets += 5 +edi->vecs.mon.neig
2538             +edi->vecs.linfix.neig
2539             +edi->vecs.linacc.neig
2540             +edi->vecs.radfix.neig
2541             +edi->vecs.radacc.neig
2542             +edi->vecs.radcon.neig
2543             + 6*edi->flood.vecs.neig;
2544         edi = edi->next_edi;
2545     }
2546     snew(setname, nsets);
2547
2548     /* In the mdrun time step in a first function call (do_flood()) the flooding
2549      * forces are calculated and in a second function call (do_edsam()) the
2550      * ED constraints. To get a corresponding legend, we need to loop twice
2551      * over the edi groups and output first the flooding, then the ED part */
2552
2553     /* The flooding-related legend entries, if flooding is done */
2554     nsets = 0;
2555     if (eEDflood == ed->eEDtype)
2556     {
2557         edi   = ed->edpar;
2558         for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2559         {
2560             /* Always write out the projection on the flooding EVs. Of course, this can also
2561              * be achieved with the monitoring option in do_edsam() (if switched on by the
2562              * user), but in that case the positions need to be communicated in do_edsam(),
2563              * which is not necessary when doing flooding only. */
2564             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2565
2566             for (i = 0; i < edi->flood.vecs.neig; i++)
2567             {
2568                 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2569                 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2570
2571                 /* Output the current reference projection if it changes with time;
2572                  * this can happen when flooding is used as harmonic restraint */
2573                 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
2574                 {
2575                     sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2576                     nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2577                 }
2578
2579                 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2580                 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2581                 {
2582                     sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2583                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2584                 }
2585
2586                 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2587                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2588
2589                 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2590                 {
2591                     sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2592                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2593                 }
2594
2595                 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2596                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2597             }
2598
2599             edi = edi->next_edi;
2600         } /* End of flooding-related legend entries */
2601     }
2602     n_flood = nsets;
2603
2604     /* Now the ED-related entries, if essential dynamics is done */
2605     edi         = ed->edpar;
2606     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2607     {
2608         if (bNeedDoEdsam(edi))  /* Only print ED legend if at least one ED option is on */
2609         {
2610             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2611
2612             /* Essential dynamics, projections on eigenvectors */
2613             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON"   );
2614             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2615             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2616             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2617             if (edi->vecs.radfix.neig)
2618             {
2619                 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2620             }
2621             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2622             if (edi->vecs.radacc.neig)
2623             {
2624                 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2625             }
2626             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2627             if (edi->vecs.radcon.neig)
2628             {
2629                 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2630             }
2631         }
2632         edi = edi->next_edi;
2633     } /* end of 'pure' essential dynamics legend entries */
2634     n_edsam = nsets - n_flood;
2635
2636     xvgr_legend(ed->edo, nsets, setname, oenv);
2637     sfree(setname);
2638
2639     fprintf(ed->edo, "#\n"
2640             "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2641             n_flood, 1 == n_flood ? "" : "s",
2642             n_edsam, 1 == n_edsam ? "" : "s");
2643     fprintf(ed->edo, "%s", LegendStr);
2644     sfree(LegendStr);
2645
2646     fflush(ed->edo);
2647 }
2648
2649
2650 /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle
2651  * contained in the input file, creates a NULL terminated list of t_edpar structures */
2652 void init_edsam(gmx_mtop_t   *mtop,
2653                 t_inputrec   *ir,
2654                 t_commrec    *cr,
2655                 gmx_edsam_t   ed,
2656                 rvec          x[],
2657                 matrix        box,
2658                 edsamstate_t *EDstate)
2659 {
2660     t_edpar *edi = NULL;                    /* points to a single edi data set */
2661     int      i, nr_edi, avindex;
2662     rvec    *x_pbc  = NULL;                 /* positions of the whole MD system with pbc removed  */
2663     rvec    *xfit   = NULL, *xstart = NULL; /* dummy arrays to determine initial RMSDs  */
2664     rvec     fit_transvec;                  /* translation ... */
2665     matrix   fit_rotmat;                    /* ... and rotation from fit to reference structure */
2666     rvec    *ref_x_old = NULL;              /* helper pointer */
2667
2668     if (MASTER(cr))
2669     {
2670         fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2671
2672         if (NULL == ed)
2673         {
2674             gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n"
2675                       "flooding simulation. Please also provide the correct .edi file with -ei.\n");
2676         }
2677     }
2678
2679     /* Needed for initializing radacc radius in do_edsam */
2680     ed->bFirst = TRUE;
2681
2682     /* The input file is read by the master and the edi structures are
2683      * initialized here. Input is stored in ed->edpar. Then the edi
2684      * structures are transferred to the other nodes */
2685     if (MASTER(cr))
2686     {
2687         /* Initialization for every ED/flooding group. Flooding uses one edi group per
2688          * flooding vector, Essential dynamics can be applied to more than one structure
2689          * as well, but will be done in the order given in the edi file, so
2690          * expect different results for different order of edi file concatenation! */
2691         edi = ed->edpar;
2692         while (edi != NULL)
2693         {
2694             init_edi(mtop, edi);
2695             init_flood(edi, ed, ir->delta_t);
2696             edi = edi->next_edi;
2697         }
2698     }
2699
2700     /* The master does the work here. The other nodes get the positions
2701      * not before dd_partition_system which is called after init_edsam */
2702     if (MASTER(cr))
2703     {
2704         if (!EDstate->bFromCpt)
2705         {
2706             /* Remove PBC, make molecule(s) subject to ED whole. */
2707             snew(x_pbc, mtop->natoms);
2708             m_rveccopy(mtop->natoms, x, x_pbc);
2709             do_pbc_first_mtop(NULL, ir->ePBC, box, mtop, x_pbc);
2710         }
2711         /* Reset pointer to first ED data set which contains the actual ED data */
2712         edi = ed->edpar;
2713         /* Loop over all ED/flooding data sets (usually only one, though) */
2714         for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2715         {
2716             /* For multiple ED groups we use the output frequency that was specified
2717              * in the first set */
2718             if (nr_edi > 1)
2719             {
2720                 edi->outfrq = ed->edpar->outfrq;
2721             }
2722
2723             /* Extract the initial reference and average positions. When starting
2724              * from .cpt, these have already been read into sref.x_old
2725              * in init_edsamstate() */
2726             if (!EDstate->bFromCpt)
2727             {
2728                 /* If this is the first run (i.e. no checkpoint present) we assume
2729                  * that the starting positions give us the correct PBC representation */
2730                 for (i = 0; i < edi->sref.nr; i++)
2731                 {
2732                     copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2733                 }
2734
2735                 for (i = 0; i < edi->sav.nr; i++)
2736                 {
2737                     copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2738                 }
2739             }
2740
2741             /* Now we have the PBC-correct start positions of the reference and
2742                average structure. We copy that over to dummy arrays on which we
2743                can apply fitting to print out the RMSD. We srenew the memory since
2744                the size of the buffers is likely different for every ED group */
2745             srenew(xfit, edi->sref.nr );
2746             srenew(xstart, edi->sav.nr  );
2747             if (edi->bRefEqAv)
2748             {
2749                 /* Reference indices are the same as average indices */
2750                 ref_x_old = edi->sav.x_old;
2751             }
2752             else
2753             {
2754                 ref_x_old = edi->sref.x_old;
2755             }
2756             copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr);
2757             copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2758
2759             /* Make the fit to the REFERENCE structure, get translation and rotation */
2760             fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2761
2762             /* Output how well we fit to the reference at the start */
2763             translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2764             fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2765                     rmsd_from_structure(xfit, &edi->sref));
2766             if (EDstate->nED > 1)
2767             {
2768                 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2769             }
2770             fprintf(stderr, "\n");
2771
2772             /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2773             translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2774
2775             /* calculate initial projections */
2776             project(xstart, edi);
2777
2778             /* For the target and origin structure both a reference (fit) and an
2779              * average structure can be provided in make_edi. If both structures
2780              * are the same, make_edi only stores one of them in the .edi file.
2781              * If they differ, first the fit and then the average structure is stored
2782              * in star (or sor), thus the number of entries in star/sor is
2783              * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2784              * the size of the average group. */
2785
2786             /* process target structure, if required */
2787             if (edi->star.nr > 0)
2788             {
2789                 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2790
2791                 /* get translation & rotation for fit of target structure to reference structure */
2792                 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2793                 /* do the fit */
2794                 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2795                 if (edi->star.nr == edi->sav.nr)
2796                 {
2797                     avindex = 0;
2798                 }
2799                 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2800                 {
2801                     /* The last sav.nr indices of the target structure correspond to
2802                      * the average structure, which must be projected */
2803                     avindex = edi->star.nr - edi->sav.nr;
2804                 }
2805                 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
2806             }
2807             else
2808             {
2809                 rad_project(edi, xstart, &edi->vecs.radcon);
2810             }
2811
2812             /* process structure that will serve as origin of expansion circle */
2813             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2814             {
2815                 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2816             }
2817
2818             if (edi->sori.nr > 0)
2819             {
2820                 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2821
2822                 /* fit this structure to reference structure */
2823                 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2824                 /* do the fit */
2825                 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2826                 if (edi->sori.nr == edi->sav.nr)
2827                 {
2828                     avindex = 0;
2829                 }
2830                 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2831                 {
2832                     /* For the projection, we need the last sav.nr indices of sori */
2833                     avindex = edi->sori.nr - edi->sav.nr;
2834                 }
2835
2836                 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2837                 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2838                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2839                 {
2840                     fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2841                     /* Set center of flooding potential to the ORIGIN structure */
2842                     rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
2843                     /* We already know that no (moving) reference position was provided,
2844                      * therefore we can overwrite refproj[0]*/
2845                     copyEvecReference(&edi->flood.vecs);
2846                 }
2847             }
2848             else /* No origin structure given */
2849             {
2850                 rad_project(edi, xstart, &edi->vecs.radacc);
2851                 rad_project(edi, xstart, &edi->vecs.radfix);
2852                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2853                 {
2854                     if (edi->flood.bHarmonic)
2855                     {
2856                         fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2857                         for (i = 0; i < edi->flood.vecs.neig; i++)
2858                         {
2859                             edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2860                         }
2861                     }
2862                     else
2863                     {
2864                         fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2865                         /* Set center of flooding potential to the center of the covariance matrix,
2866                          * i.e. the average structure, i.e. zero in the projected system */
2867                         for (i = 0; i < edi->flood.vecs.neig; i++)
2868                         {
2869                             edi->flood.vecs.refproj[i] = 0.0;
2870                         }
2871                     }
2872                 }
2873             }
2874             /* For convenience, output the center of the flooding potential for the eigenvectors */
2875             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2876             {
2877                 for (i = 0; i < edi->flood.vecs.neig; i++)
2878                 {
2879                     fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2880                     if (edi->flood.bHarmonic)
2881                     {
2882                         fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2883                     }
2884                     fprintf(stdout, "\n");
2885                 }
2886             }
2887
2888             /* set starting projections for linsam */
2889             rad_project(edi, xstart, &edi->vecs.linacc);
2890             rad_project(edi, xstart, &edi->vecs.linfix);
2891
2892             /* Prepare for the next edi data set: */
2893             edi = edi->next_edi;
2894         }
2895         /* Cleaning up on the master node: */
2896         if (!EDstate->bFromCpt)
2897         {
2898             sfree(x_pbc);
2899         }
2900         sfree(xfit);
2901         sfree(xstart);
2902
2903     } /* end of MASTER only section */
2904
2905     if (PAR(cr))
2906     {
2907         /* First let everybody know how many ED data sets to expect */
2908         gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr);
2909         /* Broadcast the essential dynamics / flooding data to all nodes */
2910         broadcast_ed_data(cr, ed, EDstate->nED);
2911     }
2912     else
2913     {
2914         /* In the single-CPU case, point the local atom numbers pointers to the global
2915          * one, so that we can use the same notation in serial and parallel case: */
2916
2917         /* Loop over all ED data sets (usually only one, though) */
2918         edi = ed->edpar;
2919         for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2920         {
2921             edi->sref.anrs_loc = edi->sref.anrs;
2922             edi->sav.anrs_loc  = edi->sav.anrs;
2923             edi->star.anrs_loc = edi->star.anrs;
2924             edi->sori.anrs_loc = edi->sori.anrs;
2925             /* For the same reason as above, make a dummy c_ind array: */
2926             snew(edi->sav.c_ind, edi->sav.nr);
2927             /* Initialize the array */
2928             for (i = 0; i < edi->sav.nr; i++)
2929             {
2930                 edi->sav.c_ind[i] = i;
2931             }
2932             /* In the general case we will need a different-sized array for the reference indices: */
2933             if (!edi->bRefEqAv)
2934             {
2935                 snew(edi->sref.c_ind, edi->sref.nr);
2936                 for (i = 0; i < edi->sref.nr; i++)
2937                 {
2938                     edi->sref.c_ind[i] = i;
2939                 }
2940             }
2941             /* Point to the very same array in case of other structures: */
2942             edi->star.c_ind = edi->sav.c_ind;
2943             edi->sori.c_ind = edi->sav.c_ind;
2944             /* In the serial case, the local number of atoms is the global one: */
2945             edi->sref.nr_loc = edi->sref.nr;
2946             edi->sav.nr_loc  = edi->sav.nr;
2947             edi->star.nr_loc = edi->star.nr;
2948             edi->sori.nr_loc = edi->sori.nr;
2949
2950             /* An on we go to the next ED group */
2951             edi = edi->next_edi;
2952         }
2953     }
2954
2955     /* Allocate space for ED buffer variables */
2956     /* Again, loop over ED data sets */
2957     edi = ed->edpar;
2958     for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2959     {
2960         /* Allocate space for ED buffer variables */
2961         snew_bc(cr, edi->buf, 1); /* MASTER has already allocated edi->buf in init_edi() */
2962         snew(edi->buf->do_edsam, 1);
2963
2964         /* Space for collective ED buffer variables */
2965
2966         /* Collective positions of atoms with the average indices */
2967         snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2968         snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr);            /* buffer for xcoll shifts */
2969         snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2970         /* Collective positions of atoms with the reference indices */
2971         if (!edi->bRefEqAv)
2972         {
2973             snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2974             snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr);       /* To store the shifts in */
2975             snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2976         }
2977
2978         /* Get memory for flooding forces */
2979         snew(edi->flood.forces_cartesian, edi->sav.nr);
2980
2981 #ifdef DUMPEDI
2982         /* Dump it all into one file per process */
2983         dump_edi(edi, cr, nr_edi);
2984 #endif
2985
2986         /* Next ED group */
2987         edi = edi->next_edi;
2988     }
2989
2990     /* Flush the edo file so that the user can check some things
2991      * when the simulation has started */
2992     if (ed->edo)
2993     {
2994         fflush(ed->edo);
2995     }
2996 }
2997
2998
2999 void do_edsam(t_inputrec     *ir,
3000               gmx_int64_t     step,
3001               t_commrec      *cr,
3002               rvec            xs[],
3003               rvec            v[],
3004               matrix          box,
3005               gmx_edsam_t     ed)
3006 {
3007     int                i, edinr, iupdate = 500;
3008     matrix             rotmat;         /* rotation matrix */
3009     rvec               transvec;       /* translation vector */
3010     rvec               dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
3011     real               dt_1;           /* 1/dt */
3012     struct t_do_edsam *buf;
3013     t_edpar           *edi;
3014     real               rmsdev    = -1;    /* RMSD from reference structure prior to applying the constraints */
3015     gmx_bool           bSuppress = FALSE; /* Write .xvg output file on master? */
3016
3017
3018     /* Check if ED sampling has to be performed */
3019     if (ed->eEDtype == eEDnone)
3020     {
3021         return;
3022     }
3023
3024     /* Suppress output on first call of do_edsam if
3025      * two-step sd2 integrator is used */
3026     if ( (ir->eI == eiSD2) && (v != NULL) )
3027     {
3028         bSuppress = TRUE;
3029     }
3030
3031     dt_1 = 1.0/ir->delta_t;
3032
3033     /* Loop over all ED groups (usually one) */
3034     edi   = ed->edpar;
3035     edinr = 0;
3036     while (edi != NULL)
3037     {
3038         edinr++;
3039         if (bNeedDoEdsam(edi))
3040         {
3041
3042             buf = edi->buf->do_edsam;
3043
3044             if (ed->bFirst)
3045             {
3046                 /* initialize radacc radius for slope criterion */
3047                 buf->oldrad = calc_radius(&edi->vecs.radacc);
3048             }
3049
3050             /* Copy the positions into buf->xc* arrays and after ED
3051              * feed back corrections to the official positions */
3052
3053             /* Broadcast the ED positions such that every node has all of them
3054              * Every node contributes its local positions xs and stores it in
3055              * the collective buf->xcoll array. Note that for edinr > 1
3056              * xs could already have been modified by an earlier ED */
3057
3058             communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3059                                         edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old,  box);
3060
3061             /* Only assembly reference positions if their indices differ from the average ones */
3062             if (!edi->bRefEqAv)
3063             {
3064                 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3065                                             edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3066             }
3067
3068             /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3069              * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3070              * set bUpdateShifts=TRUE in the parallel case. */
3071             buf->bUpdateShifts = FALSE;
3072
3073             /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3074              * as well as the indices in edi->sav.anrs */
3075
3076             /* Fit the reference indices to the reference structure */
3077             if (edi->bRefEqAv)
3078             {
3079                 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3080             }
3081             else
3082             {
3083                 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3084             }
3085
3086             /* Now apply the translation and rotation to the ED structure */
3087             translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3088
3089             /* Find out how well we fit to the reference (just for output steps) */
3090             if (do_per_step(step, edi->outfrq) && MASTER(cr))
3091             {
3092                 if (edi->bRefEqAv)
3093                 {
3094                     /* Indices of reference and average structures are identical,
3095                      * thus we can calculate the rmsd to SREF using xcoll */
3096                     rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3097                 }
3098                 else
3099                 {
3100                     /* We have to translate & rotate the reference atoms first */
3101                     translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3102                     rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3103                 }
3104             }
3105
3106             /* update radsam references, when required */
3107             if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3108             {
3109                 project(buf->xcoll, edi);
3110                 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3111                 rad_project(edi, buf->xcoll, &edi->vecs.radfix);
3112                 buf->oldrad = -1.e5;
3113             }
3114
3115             /* update radacc references, when required */
3116             if (do_per_step(step, iupdate) && step >= edi->presteps)
3117             {
3118                 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
3119                 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3120                 {
3121                     project(buf->xcoll, edi);
3122                     rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3123                     buf->oldrad = 0.0;
3124                 }
3125                 else
3126                 {
3127                     buf->oldrad = edi->vecs.radacc.radius;
3128                 }
3129             }
3130
3131             /* apply the constraints */
3132             if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3133             {
3134                 /* ED constraints should be applied already in the first MD step
3135                  * (which is step 0), therefore we pass step+1 to the routine */
3136                 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3137             }
3138
3139             /* write to edo, when required */
3140             if (do_per_step(step, edi->outfrq))
3141             {
3142                 project(buf->xcoll, edi);
3143                 if (MASTER(cr) && !bSuppress)
3144                 {
3145                     write_edo(edi, ed->edo, rmsdev);
3146                 }
3147             }
3148
3149             /* Copy back the positions unless monitoring only */
3150             if (ed_constraints(ed->eEDtype, edi))
3151             {
3152                 /* remove fitting */
3153                 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3154
3155                 /* Copy the ED corrected positions into the coordinate array */
3156                 /* Each node copies its local part. In the serial case, nat_loc is the
3157                  * total number of ED atoms */
3158                 for (i = 0; i < edi->sav.nr_loc; i++)
3159                 {
3160                     /* Unshift local ED coordinate and store in x_unsh */
3161                     ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3162                                             buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3163
3164                     /* dx is the ED correction to the positions: */
3165                     rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3166
3167                     if (v != NULL)
3168                     {
3169                         /* dv is the ED correction to the velocity: */
3170                         svmul(dt_1, dx, dv);
3171                         /* apply the velocity correction: */
3172                         rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3173                     }
3174                     /* Finally apply the position correction due to ED: */
3175                     copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3176                 }
3177             }
3178         } /* END of if ( bNeedDoEdsam(edi) ) */
3179
3180         /* Prepare for the next ED group */
3181         edi = edi->next_edi;
3182
3183     } /* END of loop over ED groups */
3184
3185     ed->bFirst = FALSE;
3186 }