src/gromacs/mdlib/edsam.c

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