src/gromacs/essentialdynamics/edsam.c

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