src/gromacs/essentialdynamics/edsam.c

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