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