src/gromacs/essentialdynamics/edsam.c

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