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 "gromacs/fileio/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 "gromacs/fileio/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     alook = gmx_mtop_atomlookup_init(mtop);
1346
1347     /* evaluate masses (reference structure) */
1348     snew(edi->sref.m, edi->sref.nr);
1349     for (i = 0; i < edi->sref.nr; i++)
1350     {
1351         if (edi->fitmas)
1352         {
1353             gmx_mtop_atomnr_to_atom(alook, edi->sref.anrs[i], &atom);
1354             edi->sref.m[i] = atom->m;
1355         }
1356         else
1357         {
1358             edi->sref.m[i] = 1.0;
1359         }
1360
1361         /* Check that every m > 0. Bad things will happen otherwise. */
1362         if (edi->sref.m[i] <= 0.0)
1363         {
1364             gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1365                       "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1366                       "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1367                       "atoms from the reference structure by creating a proper index group.\n",
1368                       i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1369         }
1370
1371         totalmass += edi->sref.m[i];
1372     }
1373     edi->sref.mtot = totalmass;
1374
1375     /* Masses m and sqrt(m) for the average structure. Note that m
1376      * is needed if forces have to be evaluated in do_edsam */
1377     snew(edi->sav.sqrtm, edi->sav.nr );
1378     snew(edi->sav.m, edi->sav.nr );
1379     for (i = 0; i < edi->sav.nr; i++)
1380     {
1381         gmx_mtop_atomnr_to_atom(alook, edi->sav.anrs[i], &atom);
1382         edi->sav.m[i] = atom->m;
1383         if (edi->pcamas)
1384         {
1385             edi->sav.sqrtm[i] = sqrt(atom->m);
1386         }
1387         else
1388         {
1389             edi->sav.sqrtm[i] = 1.0;
1390         }
1391
1392         /* Check that every m > 0. Bad things will happen otherwise. */
1393         if (edi->sav.sqrtm[i] <= 0.0)
1394         {
1395             gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1396                       "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1397                       "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1398                       "atoms from the average structure by creating a proper index group.\n",
1399                       i, edi->sav.anrs[i]+1, atom->m);
1400         }
1401     }
1402
1403     gmx_mtop_atomlookup_destroy(alook);
1404
1405     /* put reference structure in origin */
1406     get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1407     com[XX] = -com[XX];
1408     com[YY] = -com[YY];
1409     com[ZZ] = -com[ZZ];
1410     translate_x(edi->sref.x, edi->sref.nr, com);
1411
1412     /* Init ED buffer */
1413     snew(edi->buf, 1);
1414 }
1415
1416
1417 static void check(const char *line, const char *label)
1418 {
1419     if (!strstr(line, label))
1420     {
1421         gmx_fatal(FARGS, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line);
1422     }
1423 }
1424
1425
1426 static int read_checked_edint(FILE *file, const char *label)
1427 {
1428     char line[STRLEN+1];
1429     int  idum;
1430
1431
1432     fgets2 (line, STRLEN, file);
1433     check(line, label);
1434     fgets2 (line, STRLEN, file);
1435     sscanf (line, "%d", &idum);
1436     return idum;
1437 }
1438
1439
1440 static int read_edint(FILE *file, gmx_bool *bEOF)
1441 {
1442     char  line[STRLEN+1];
1443     int   idum;
1444     char *eof;
1445
1446
1447     eof = fgets2 (line, STRLEN, file);
1448     if (eof == NULL)
1449     {
1450         *bEOF = TRUE;
1451         return -1;
1452     }
1453     eof = fgets2 (line, STRLEN, file);
1454     if (eof == NULL)
1455     {
1456         *bEOF = TRUE;
1457         return -1;
1458     }
1459     sscanf (line, "%d", &idum);
1460     *bEOF = FALSE;
1461     return idum;
1462 }
1463
1464
1465 static real read_checked_edreal(FILE *file, const char *label)
1466 {
1467     char   line[STRLEN+1];
1468     double rdum;
1469
1470
1471     fgets2 (line, STRLEN, file);
1472     check(line, label);
1473     fgets2 (line, STRLEN, file);
1474     sscanf (line, "%lf", &rdum);
1475     return (real) rdum; /* always read as double and convert to single */
1476 }
1477
1478
1479 static void read_edx(FILE *file, int number, int *anrs, rvec *x)
1480 {
1481     int    i, j;
1482     char   line[STRLEN+1];
1483     double d[3];
1484
1485
1486     for (i = 0; i < number; i++)
1487     {
1488         fgets2 (line, STRLEN, file);
1489         sscanf (line, "%d%lf%lf%lf", &anrs[i], &d[0], &d[1], &d[2]);
1490         anrs[i]--; /* we are reading FORTRAN indices */
1491         for (j = 0; j < 3; j++)
1492         {
1493             x[i][j] = d[j]; /* always read as double and convert to single */
1494         }
1495     }
1496 }
1497
1498
1499 static void scan_edvec(FILE *in, int nr, rvec *vec)
1500 {
1501     char   line[STRLEN+1];
1502     int    i;
1503     double x, y, z;
1504
1505
1506     for (i = 0; (i < nr); i++)
1507     {
1508         fgets2 (line, STRLEN, in);
1509         sscanf (line, "%le%le%le", &x, &y, &z);
1510         vec[i][XX] = x;
1511         vec[i][YY] = y;
1512         vec[i][ZZ] = z;
1513     }
1514 }
1515
1516
1517 static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1518 {
1519     int    i, idum, nscan;
1520     double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0;
1521     char   line[STRLEN+1];
1522
1523
1524     tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
1525     if (tvec->neig > 0)
1526     {
1527         snew(tvec->ieig, tvec->neig);
1528         snew(tvec->stpsz, tvec->neig);
1529         snew(tvec->vec, tvec->neig);
1530         snew(tvec->xproj, tvec->neig);
1531         snew(tvec->fproj, tvec->neig);
1532         snew(tvec->refproj, tvec->neig);
1533         if (bReadRefproj)
1534         {
1535             snew(tvec->refproj0, tvec->neig);
1536             snew(tvec->refprojslope, tvec->neig);
1537         }
1538
1539         for (i = 0; (i < tvec->neig); i++)
1540         {
1541             fgets2 (line, STRLEN, in);
1542             if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1543             {
1544                 nscan = sscanf(line, "%d%lf%lf%lf", &idum, &rdum, &refproj_dum, &refprojslope_dum);
1545                 /* Zero out values which were not scanned */
1546                 switch (nscan)
1547                 {
1548                     case 4:
1549                         /* Every 4 values read, including reference position */
1550                         *bHaveReference = TRUE;
1551                         break;
1552                     case 3:
1553                         /* A reference position is provided */
1554                         *bHaveReference = TRUE;
1555                         /* No value for slope, set to 0 */
1556                         refprojslope_dum = 0.0;
1557                         break;
1558                     case 2:
1559                         /* No values for reference projection and slope, set to 0 */
1560                         refproj_dum      = 0.0;
1561                         refprojslope_dum = 0.0;
1562                         break;
1563                     default:
1564                         gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1565                         break;
1566                 }
1567                 tvec->refproj[i]      = refproj_dum;
1568                 tvec->refproj0[i]     = refproj_dum;
1569                 tvec->refprojslope[i] = refprojslope_dum;
1570             }
1571             else /* Normal flooding */
1572             {
1573                 nscan = sscanf(line, "%d%lf", &idum, &rdum);
1574                 if (nscan != 2)
1575                 {
1576                     gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
1577                 }
1578             }
1579             tvec->ieig[i]  = idum;
1580             tvec->stpsz[i] = rdum;
1581         } /* end of loop over eigenvectors */
1582
1583         for (i = 0; (i < tvec->neig); i++)
1584         {
1585             snew(tvec->vec[i], nr);
1586             scan_edvec(in, nr, tvec->vec[i]);
1587         }
1588     }
1589 }
1590
1591
1592 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1593 static void read_edvecs(FILE *in, int nr, t_edvecs *vecs)
1594 {
1595     gmx_bool bHaveReference = FALSE;
1596
1597
1598     read_edvec(in, nr, &vecs->mon, FALSE, &bHaveReference);
1599     read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1600     read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1601     read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1602     read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1603     read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1604 }
1605
1606
1607 /* Check if the same atom indices are used for reference and average positions */
1608 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1609 {
1610     int i;
1611
1612
1613     /* If the number of atoms differs between the two structures,
1614      * they cannot be identical */
1615     if (sref.nr != sav.nr)
1616     {
1617         return FALSE;
1618     }
1619
1620     /* Now that we know that both stuctures have the same number of atoms,
1621      * check if also the indices are identical */
1622     for (i = 0; i < sav.nr; i++)
1623     {
1624         if (sref.anrs[i] != sav.anrs[i])
1625         {
1626             return FALSE;
1627         }
1628     }
1629     fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1630
1631     return TRUE;
1632 }
1633
1634
1635 static int read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, const char *fn)
1636 {
1637     int       readmagic;
1638     const int magic = 670;
1639     gmx_bool  bEOF;
1640
1641     /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1642     gmx_bool bHaveReference = FALSE;
1643
1644
1645     /* the edi file is not free format, so expect problems if the input is corrupt. */
1646
1647     /* check the magic number */
1648     readmagic = read_edint(in, &bEOF);
1649     /* Check whether we have reached the end of the input file */
1650     if (bEOF)
1651     {
1652         return 0;
1653     }
1654
1655     if (readmagic != magic)
1656     {
1657         if (readmagic == 666 || readmagic == 667 || readmagic == 668)
1658         {
1659             gmx_fatal(FARGS, "Wrong magic number: Use newest version of make_edi to produce edi file");
1660         }
1661         else if (readmagic != 669)
1662         {
1663             gmx_fatal(FARGS, "Wrong magic number %d in %s", readmagic, fn);
1664         }
1665     }
1666
1667     /* check the number of atoms */
1668     edi->nini = read_edint(in, &bEOF);
1669     if (edi->nini != nr_mdatoms)
1670     {
1671         gmx_fatal(FARGS, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms);
1672     }
1673
1674     /* Done checking. For the rest we blindly trust the input */
1675     edi->fitmas          = read_checked_edint(in, "FITMAS");
1676     edi->pcamas          = read_checked_edint(in, "ANALYSIS_MAS");
1677     edi->outfrq          = read_checked_edint(in, "OUTFRQ");
1678     edi->maxedsteps      = read_checked_edint(in, "MAXLEN");
1679     edi->slope           = read_checked_edreal(in, "SLOPECRIT");
1680
1681     edi->presteps        = read_checked_edint(in, "PRESTEPS");
1682     edi->flood.deltaF0   = read_checked_edreal(in, "DELTA_F0");
1683     edi->flood.deltaF    = read_checked_edreal(in, "INIT_DELTA_F");
1684     edi->flood.tau       = read_checked_edreal(in, "TAU");
1685     edi->flood.constEfl  = read_checked_edreal(in, "EFL_NULL");
1686     edi->flood.alpha2    = read_checked_edreal(in, "ALPHA2");
1687     edi->flood.kT        = read_checked_edreal(in, "KT");
1688     edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC");
1689     if (readmagic > 669)
1690     {
1691         edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING");
1692     }
1693     else
1694     {
1695         edi->flood.bConstForce = FALSE;
1696     }
1697     edi->sref.nr         = read_checked_edint(in, "NREF");
1698
1699     /* allocate space for reference positions and read them */
1700     snew(edi->sref.anrs, edi->sref.nr);
1701     snew(edi->sref.x, edi->sref.nr);
1702     snew(edi->sref.x_old, edi->sref.nr);
1703     edi->sref.sqrtm    = NULL;
1704     read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x);
1705
1706     /* average positions. they define which atoms will be used for ED sampling */
1707     edi->sav.nr = read_checked_edint(in, "NAV");
1708     snew(edi->sav.anrs, edi->sav.nr);
1709     snew(edi->sav.x, edi->sav.nr);
1710     snew(edi->sav.x_old, edi->sav.nr);
1711     read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x);
1712
1713     /* Check if the same atom indices are used for reference and average positions */
1714     edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1715
1716     /* eigenvectors */
1717     read_edvecs(in, edi->sav.nr, &edi->vecs);
1718     read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
1719
1720     /* target positions */
1721     edi->star.nr = read_edint(in, &bEOF);
1722     if (edi->star.nr > 0)
1723     {
1724         snew(edi->star.anrs, edi->star.nr);
1725         snew(edi->star.x, edi->star.nr);
1726         edi->star.sqrtm    = NULL;
1727         read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x);
1728     }
1729
1730     /* positions defining origin of expansion circle */
1731     edi->sori.nr = read_edint(in, &bEOF);
1732     if (edi->sori.nr > 0)
1733     {
1734         if (bHaveReference)
1735         {
1736             /* Both an -ori structure and a at least one manual reference point have been
1737              * specified. That's ambiguous and probably not intentional. */
1738             gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1739                       "    point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1740         }
1741         snew(edi->sori.anrs, edi->sori.nr);
1742         snew(edi->sori.x, edi->sori.nr);
1743         edi->sori.sqrtm    = NULL;
1744         read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x);
1745     }
1746
1747     /* all done */
1748     return 1;
1749 }
1750
1751
1752
1753 /* Read in the edi input file. Note that it may contain several ED data sets which were
1754  * achieved by concatenating multiple edi files. The standard case would be a single ED
1755  * data set, though. */
1756 static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
1757 {
1758     FILE    *in;
1759     t_edpar *curr_edi, *last_edi;
1760     t_edpar *edi_read;
1761     int      edi_nr = 0;
1762
1763
1764     /* This routine is executed on the master only */
1765
1766     /* Open the .edi parameter input file */
1767     in = gmx_fio_fopen(fn, "r");
1768     fprintf(stderr, "ED: Reading edi file %s\n", fn);
1769
1770     /* Now read a sequence of ED input parameter sets from the edi file */
1771     curr_edi = edi;
1772     last_edi = edi;
1773     while (read_edi(in, curr_edi, nr_mdatoms, fn) )
1774     {
1775         edi_nr++;
1776
1777         /* Since we arrived within this while loop we know that there is still another data set to be read in */
1778         /* We need to allocate space for the data: */
1779         snew(edi_read, 1);
1780         /* Point the 'next_edi' entry to the next edi: */
1781         curr_edi->next_edi = edi_read;
1782         /* Keep the curr_edi pointer for the case that the next group is empty: */
1783         last_edi = curr_edi;
1784         /* Let's prepare to read in the next edi data set: */
1785         curr_edi = edi_read;
1786     }
1787     if (edi_nr == 0)
1788     {
1789         gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
1790     }
1791
1792     /* Terminate the edi group list with a NULL pointer: */
1793     last_edi->next_edi = NULL;
1794
1795     fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : "");
1796
1797     /* Close the .edi file again */
1798     gmx_fio_fclose(in);
1799
1800     return edi_nr;
1801 }
1802
1803
1804 struct t_fit_to_ref {
1805     rvec *xcopy;       /* Working copy of the positions in fit_to_reference */
1806 };
1807
1808 /* Fit the current positions to the reference positions
1809  * Do not actually do the fit, just return rotation and translation.
1810  * Note that the COM of the reference structure was already put into
1811  * the origin by init_edi. */
1812 static void fit_to_reference(rvec      *xcoll,    /* The positions to be fitted */
1813                              rvec       transvec, /* The translation vector */
1814                              matrix     rotmat,   /* The rotation matrix */
1815                              t_edpar   *edi)      /* Just needed for do_edfit */
1816 {
1817     rvec                 com;                     /* center of mass */
1818     int                  i;
1819     struct t_fit_to_ref *loc;
1820
1821
1822     /* Allocate memory the first time this routine is called for each edi group */
1823     if (NULL == edi->buf->fit_to_ref)
1824     {
1825         snew(edi->buf->fit_to_ref, 1);
1826         snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1827     }
1828     loc = edi->buf->fit_to_ref;
1829
1830     /* We do not touch the original positions but work on a copy. */
1831     for (i = 0; i < edi->sref.nr; i++)
1832     {
1833         copy_rvec(xcoll[i], loc->xcopy[i]);
1834     }
1835
1836     /* Calculate the center of mass */
1837     get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1838
1839     transvec[XX] = -com[XX];
1840     transvec[YY] = -com[YY];
1841     transvec[ZZ] = -com[ZZ];
1842
1843     /* Subtract the center of mass from the copy */
1844     translate_x(loc->xcopy, edi->sref.nr, transvec);
1845
1846     /* Determine the rotation matrix */
1847     do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1848 }
1849
1850
1851 static void translate_and_rotate(rvec  *x,        /* The positions to be translated and rotated */
1852                                  int    nat,      /* How many positions are there? */
1853                                  rvec   transvec, /* The translation vector */
1854                                  matrix rotmat)   /* The rotation matrix */
1855 {
1856     /* Translation */
1857     translate_x(x, nat, transvec);
1858
1859     /* Rotation */
1860     rotate_x(x, nat, rotmat);
1861 }
1862
1863
1864 /* Gets the rms deviation of the positions to the structure s */
1865 /* fit_to_structure has to be called before calling this routine! */
1866 static real rmsd_from_structure(rvec           *x,  /* The positions under consideration */
1867                                 struct gmx_edx *s)  /* The structure from which the rmsd shall be computed */
1868 {
1869     real  rmsd = 0.0;
1870     int   i;
1871
1872
1873     for (i = 0; i < s->nr; i++)
1874     {
1875         rmsd += distance2(s->x[i], x[i]);
1876     }
1877
1878     rmsd /= (real) s->nr;
1879     rmsd  = sqrt(rmsd);
1880
1881     return rmsd;
1882 }
1883
1884
1885 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1886 {
1887     t_edpar *edi;
1888
1889
1890     if (ed->eEDtype != eEDnone)
1891     {
1892         /* Loop over ED groups */
1893         edi = ed->edpar;
1894         while (edi)
1895         {
1896             /* Local atoms of the reference structure (for fitting), need only be assembled
1897              * if their indices differ from the average ones */
1898             if (!edi->bRefEqAv)
1899             {
1900                 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1901                                             &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1902             }
1903
1904             /* Local atoms of the average structure (on these ED will be performed) */
1905             dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1906                                         &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1907
1908             /* Indicate that the ED shift vectors for this structure need to be updated
1909              * at the next call to communicate_group_positions, since obviously we are in a NS step */
1910             edi->buf->do_edsam->bUpdateShifts = TRUE;
1911
1912             /* Set the pointer to the next ED group (if any) */
1913             edi = edi->next_edi;
1914         }
1915     }
1916 }
1917
1918
1919 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1920 {
1921     int tx, ty, tz;
1922
1923
1924     tx = is[XX];
1925     ty = is[YY];
1926     tz = is[ZZ];
1927
1928     if (TRICLINIC(box))
1929     {
1930         xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1931         xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1932         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1933     }
1934     else
1935     {
1936         xu[XX] = x[XX]-tx*box[XX][XX];
1937         xu[YY] = x[YY]-ty*box[YY][YY];
1938         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1939     }
1940 }
1941
1942
1943 static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_large_int_t step)
1944 {
1945     int  i, j;
1946     real proj, add;
1947     rvec vec_dum;
1948
1949
1950     /* loop over linfix vectors */
1951     for (i = 0; i < edi->vecs.linfix.neig; i++)
1952     {
1953         /* calculate the projection */
1954         proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1955
1956         /* calculate the correction */
1957         add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1958
1959         /* apply the correction */
1960         add /= edi->sav.sqrtm[i];
1961         for (j = 0; j < edi->sav.nr; j++)
1962         {
1963             svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1964             rvec_inc(xcoll[j], vec_dum);
1965         }
1966     }
1967 }
1968
1969
1970 static void do_linacc(rvec *xcoll, t_edpar *edi)
1971 {
1972     int  i, j;
1973     real proj, add;
1974     rvec vec_dum;
1975
1976
1977     /* loop over linacc vectors */
1978     for (i = 0; i < edi->vecs.linacc.neig; i++)
1979     {
1980         /* calculate the projection */
1981         proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1982
1983         /* calculate the correction */
1984         add = 0.0;
1985         if (edi->vecs.linacc.stpsz[i] > 0.0)
1986         {
1987             if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1988             {
1989                 add = edi->vecs.linacc.refproj[i] - proj;
1990             }
1991         }
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
2000         /* apply the correction */
2001         add /= edi->sav.sqrtm[i];
2002         for (j = 0; j < edi->sav.nr; j++)
2003         {
2004             svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
2005             rvec_inc(xcoll[j], vec_dum);
2006         }
2007
2008         /* new positions will act as reference */
2009         edi->vecs.linacc.refproj[i] = proj + add;
2010     }
2011 }
2012
2013
2014 static void do_radfix(rvec *xcoll, t_edpar *edi)
2015 {
2016     int   i, j;
2017     real *proj, rad = 0.0, ratio;
2018     rvec  vec_dum;
2019
2020
2021     if (edi->vecs.radfix.neig == 0)
2022     {
2023         return;
2024     }
2025
2026     snew(proj, edi->vecs.radfix.neig);
2027
2028     /* loop over radfix vectors */
2029     for (i = 0; i < edi->vecs.radfix.neig; i++)
2030     {
2031         /* calculate the projections, radius */
2032         proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
2033         rad    += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
2034     }
2035
2036     rad                      = sqrt(rad);
2037     ratio                    = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
2038     edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
2039
2040     /* loop over radfix vectors */
2041     for (i = 0; i < edi->vecs.radfix.neig; i++)
2042     {
2043         proj[i] -= edi->vecs.radfix.refproj[i];
2044
2045         /* apply the correction */
2046         proj[i] /= edi->sav.sqrtm[i];
2047         proj[i] *= ratio;
2048         for (j = 0; j < edi->sav.nr; j++)
2049         {
2050             svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2051             rvec_inc(xcoll[j], vec_dum);
2052         }
2053     }
2054
2055     sfree(proj);
2056 }
2057
2058
2059 static void do_radacc(rvec *xcoll, t_edpar *edi)
2060 {
2061     int   i, j;
2062     real *proj, rad = 0.0, ratio = 0.0;
2063     rvec  vec_dum;
2064
2065
2066     if (edi->vecs.radacc.neig == 0)
2067     {
2068         return;
2069     }
2070
2071     snew(proj, edi->vecs.radacc.neig);
2072
2073     /* loop over radacc vectors */
2074     for (i = 0; i < edi->vecs.radacc.neig; i++)
2075     {
2076         /* calculate the projections, radius */
2077         proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2078         rad    += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
2079     }
2080     rad = sqrt(rad);
2081
2082     /* only correct when radius decreased */
2083     if (rad < edi->vecs.radacc.radius)
2084     {
2085         ratio = edi->vecs.radacc.radius/rad - 1.0;
2086         rad   = edi->vecs.radacc.radius;
2087     }
2088     else
2089     {
2090         edi->vecs.radacc.radius = rad;
2091     }
2092
2093     /* loop over radacc vectors */
2094     for (i = 0; i < edi->vecs.radacc.neig; i++)
2095     {
2096         proj[i] -= edi->vecs.radacc.refproj[i];
2097
2098         /* apply the correction */
2099         proj[i] /= edi->sav.sqrtm[i];
2100         proj[i] *= ratio;
2101         for (j = 0; j < edi->sav.nr; j++)
2102         {
2103             svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2104             rvec_inc(xcoll[j], vec_dum);
2105         }
2106     }
2107     sfree(proj);
2108 }
2109
2110
2111 struct t_do_radcon {
2112     real *proj;
2113 };
2114
2115 static void do_radcon(rvec *xcoll, t_edpar *edi)
2116 {
2117     int                 i, j;
2118     real                rad = 0.0, ratio = 0.0;
2119     struct t_do_radcon *loc;
2120     gmx_bool            bFirst;
2121     rvec                vec_dum;
2122
2123
2124     if (edi->buf->do_radcon != NULL)
2125     {
2126         bFirst = FALSE;
2127         loc    = edi->buf->do_radcon;
2128     }
2129     else
2130     {
2131         bFirst = TRUE;
2132         snew(edi->buf->do_radcon, 1);
2133     }
2134     loc = edi->buf->do_radcon;
2135
2136     if (edi->vecs.radcon.neig == 0)
2137     {
2138         return;
2139     }
2140
2141     if (bFirst)
2142     {
2143         snew(loc->proj, edi->vecs.radcon.neig);
2144     }
2145
2146     /* loop over radcon vectors */
2147     for (i = 0; i < edi->vecs.radcon.neig; i++)
2148     {
2149         /* calculate the projections, radius */
2150         loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2151         rad         += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2152     }
2153     rad = sqrt(rad);
2154     /* only correct when radius increased */
2155     if (rad > edi->vecs.radcon.radius)
2156     {
2157         ratio = edi->vecs.radcon.radius/rad - 1.0;
2158
2159         /* loop over radcon vectors */
2160         for (i = 0; i < edi->vecs.radcon.neig; i++)
2161         {
2162             /* apply the correction */
2163             loc->proj[i] -= edi->vecs.radcon.refproj[i];
2164             loc->proj[i] /= edi->sav.sqrtm[i];
2165             loc->proj[i] *= ratio;
2166
2167             for (j = 0; j < edi->sav.nr; j++)
2168             {
2169                 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2170                 rvec_inc(xcoll[j], vec_dum);
2171             }
2172         }
2173     }
2174     else
2175     {
2176         edi->vecs.radcon.radius = rad;
2177     }
2178
2179     if (rad != edi->vecs.radcon.radius)
2180     {
2181         rad = 0.0;
2182         for (i = 0; i < edi->vecs.radcon.neig; i++)
2183         {
2184             /* calculate the projections, radius */
2185             loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2186             rad         += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2187         }
2188         rad = sqrt(rad);
2189     }
2190 }
2191
2192
2193 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step)
2194 {
2195     int i;
2196
2197
2198     /* subtract the average positions */
2199     for (i = 0; i < edi->sav.nr; i++)
2200     {
2201         rvec_dec(xcoll[i], edi->sav.x[i]);
2202     }
2203
2204     /* apply the constraints */
2205     if (step >= 0)
2206     {
2207         do_linfix(xcoll, edi, step);
2208     }
2209     do_linacc(xcoll, edi);
2210     if (step >= 0)
2211     {
2212         do_radfix(xcoll, edi);
2213     }
2214     do_radacc(xcoll, edi);
2215     do_radcon(xcoll, edi);
2216
2217     /* add back the average positions */
2218     for (i = 0; i < edi->sav.nr; i++)
2219     {
2220         rvec_inc(xcoll[i], edi->sav.x[i]);
2221     }
2222 }
2223
2224
2225 /* Write out the projections onto the eigenvectors. The order of output
2226  * corresponds to ed_output_legend() */
2227 static void write_edo(t_edpar *edi, FILE *fp, real rmsd)
2228 {
2229     int i;
2230
2231
2232     /* Output how well we fit to the reference structure */
2233     fprintf(fp, EDcol_ffmt, rmsd);
2234
2235     for (i = 0; i < edi->vecs.mon.neig; i++)
2236     {
2237         fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
2238     }
2239
2240     for (i = 0; i < edi->vecs.linfix.neig; i++)
2241     {
2242         fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
2243     }
2244
2245     for (i = 0; i < edi->vecs.linacc.neig; i++)
2246     {
2247         fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
2248     }
2249
2250     for (i = 0; i < edi->vecs.radfix.neig; i++)
2251     {
2252         fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
2253     }
2254     if (edi->vecs.radfix.neig)
2255     {
2256         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
2257     }
2258
2259     for (i = 0; i < edi->vecs.radacc.neig; i++)
2260     {
2261         fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
2262     }
2263     if (edi->vecs.radacc.neig)
2264     {
2265         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
2266     }
2267
2268     for (i = 0; i < edi->vecs.radcon.neig; i++)
2269     {
2270         fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
2271     }
2272     if (edi->vecs.radcon.neig)
2273     {
2274         fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
2275     }
2276 }
2277
2278 /* Returns if any constraints are switched on */
2279 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2280 {
2281     if (edtype == eEDedsam || edtype == eEDflood)
2282     {
2283         return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2284                 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2285                 edi->vecs.radcon.neig);
2286     }
2287     return 0;
2288 }
2289
2290
2291 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2292  * umbrella sampling simulations. */
2293 static void copyEvecReference(t_eigvec* floodvecs)
2294 {
2295     int i;
2296
2297
2298     if (NULL == floodvecs->refproj0)
2299     {
2300         snew(floodvecs->refproj0, floodvecs->neig);
2301     }
2302
2303     for (i = 0; i < floodvecs->neig; i++)
2304     {
2305         floodvecs->refproj0[i] = floodvecs->refproj[i];
2306     }
2307 }
2308
2309
2310 /* Call on MASTER only. Check whether the essential dynamics / flooding
2311  * groups of the checkpoint file are consistent with the provided .edi file. */
2312 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate)
2313 {
2314     t_edpar *edi = NULL;    /* points to a single edi data set */
2315     int      edinum;
2316
2317
2318     if (NULL == EDstate->nref || NULL == EDstate->nav)
2319     {
2320         gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2321                   "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2322                   "it must also continue with/without ED constraints when checkpointing.\n"
2323                   "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2324                   "from without a checkpoint.\n");
2325     }
2326
2327     edi    = ed->edpar;
2328     edinum = 0;
2329     while (edi != NULL)
2330     {
2331         /* Check number of atoms in the reference and average structures */
2332         if (EDstate->nref[edinum] != edi->sref.nr)
2333         {
2334             gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
2335                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2336                       get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
2337         }
2338         if (EDstate->nav[edinum] != edi->sav.nr)
2339         {
2340             gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
2341                       "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2342                       get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
2343         }
2344         edi = edi->next_edi;
2345         edinum++;
2346     }
2347
2348     if (edinum != EDstate->nED)
2349     {
2350         gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
2351                   "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
2352                   "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
2353     }
2354 }
2355
2356
2357 /* The edsamstate struct stores the information we need to make the ED group
2358  * whole again after restarts from a checkpoint file. Here we do the following:
2359  * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2360  * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2361  * c) in any case, for subsequent checkpoint writing, we set the pointers in
2362  * edsamstate to the x_old arrays, which contain the correct PBC representation of
2363  * all ED structures at the last time step. */
2364 static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate)
2365 {
2366     int      i, nr_edi;
2367     t_edpar *edi;
2368
2369
2370     snew(EDstate->old_sref_p, EDstate->nED);
2371     snew(EDstate->old_sav_p, EDstate->nED);
2372
2373     /* If we did not read in a .cpt file, these arrays are not yet allocated */
2374     if (!EDstate->bFromCpt)
2375     {
2376         snew(EDstate->nref, EDstate->nED);
2377         snew(EDstate->nav, EDstate->nED);
2378     }
2379
2380     /* Loop over all ED/flooding data sets (usually only one, though) */
2381     edi = ed->edpar;
2382     for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2383     {
2384         /* We always need the last reference and average positions such that
2385          * in the next time step we can make the ED group whole again
2386          * if the atoms do not have the correct PBC representation */
2387         if (EDstate->bFromCpt)
2388         {
2389             /* Copy the last whole positions of reference and average group from .cpt */
2390             for (i = 0; i < edi->sref.nr; i++)
2391             {
2392                 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2393             }
2394             for (i = 0; i < edi->sav.nr; i++)
2395             {
2396                 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2397             }
2398         }
2399         else
2400         {
2401             EDstate->nref[nr_edi-1] = edi->sref.nr;
2402             EDstate->nav [nr_edi-1] = edi->sav.nr;
2403         }
2404
2405         /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2406         EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2407         EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old;
2408
2409         edi = edi->next_edi;
2410     }
2411 }
2412
2413
2414 /* Adds 'buf' to 'str' */
2415 static void add_to_string(char **str, char *buf)
2416 {
2417     int len;
2418
2419
2420     len = strlen(*str) + strlen(buf) + 1;
2421     srenew(*str, len);
2422     strcat(*str, buf);
2423 }
2424
2425
2426 static void add_to_string_aligned(char **str, char *buf)
2427 {
2428     char buf_aligned[STRLEN];
2429
2430     sprintf(buf_aligned, EDcol_sfmt, buf);
2431     add_to_string(str, buf_aligned);
2432 }
2433
2434
2435 static void nice_legend(const char ***setname, int *nsets, char **LegendStr, char *value, char *unit, char EDgroupchar)
2436 {
2437     char tmp[STRLEN], tmp2[STRLEN];
2438
2439
2440     sprintf(tmp, "%c %s", EDgroupchar, value);
2441     add_to_string_aligned(LegendStr, tmp);
2442     sprintf(tmp2, "%s (%s)", tmp, unit);
2443     (*setname)[*nsets] = strdup(tmp2);
2444     (*nsets)++;
2445 }
2446
2447
2448 static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
2449 {
2450     int  i;
2451     char tmp[STRLEN];
2452
2453
2454     for (i = 0; i < evec->neig; i++)
2455     {
2456         sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
2457         nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
2458     }
2459 }
2460
2461
2462 /* Makes a legend for the xvg output file. Call on MASTER only! */
2463 static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv)
2464 {
2465     t_edpar     *edi = NULL;
2466     int          i;
2467     int          nr_edi, nsets, n_flood, n_edsam;
2468     const char **setname;
2469     char         buf[STRLEN];
2470     char        *LegendStr = NULL;
2471
2472
2473     edi         = ed->edpar;
2474
2475     fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : "");
2476
2477     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2478     {
2479         /* Remember for each ED group whether we have to do essential dynamics
2480          * constraints or possibly only flooding */
2481         edi->bNeedDoEdsam = edi->vecs.mon.neig
2482             || edi->vecs.linfix.neig
2483             || edi->vecs.linacc.neig
2484             || edi->vecs.radfix.neig
2485             || edi->vecs.radacc.neig
2486             || edi->vecs.radcon.neig;
2487
2488         fprintf(ed->edo, "#\n");
2489         fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
2490         fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
2491         fprintf(ed->edo, "#    monitor  : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig    != 1 ? "s" : "");
2492         fprintf(ed->edo, "#    LINFIX   : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : "");
2493         fprintf(ed->edo, "#    LINACC   : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : "");
2494         fprintf(ed->edo, "#    RADFIX   : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : "");
2495         fprintf(ed->edo, "#    RADACC   : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : "");
2496         fprintf(ed->edo, "#    RADCON   : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : "");
2497         fprintf(ed->edo, "#    FLOODING : %d vec%s  ", edi->flood.vecs.neig, edi->flood.vecs.neig  != 1 ? "s" : "");
2498
2499         if (edi->flood.vecs.neig)
2500         {
2501             /* If in any of the groups we find a flooding vector, flooding is turned on */
2502             ed->eEDtype = eEDflood;
2503
2504             /* Print what flavor of flooding we will do */
2505             if (0 == edi->flood.tau) /* constant flooding strength */
2506             {
2507                 fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
2508                 if (edi->flood.bHarmonic)
2509                 {
2510                     fprintf(ed->edo, ", harmonic");
2511                 }
2512             }
2513             else /* adaptive flooding */
2514             {
2515                 fprintf(ed->edo, ", adaptive");
2516             }
2517         }
2518         fprintf(ed->edo, "\n");
2519
2520         edi = edi->next_edi;
2521     }
2522
2523     /* Print a nice legend */
2524     snew(LegendStr, 1);
2525     LegendStr[0] = '\0';
2526     sprintf(buf, "#     %6s", "time");
2527     add_to_string(&LegendStr, buf);
2528
2529     /* Calculate the maximum number of columns we could end up with */
2530     edi     = ed->edpar;
2531     nsets   = 0;
2532     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2533     {
2534         nsets += 5 +edi->vecs.mon.neig
2535             +edi->vecs.linfix.neig
2536             +edi->vecs.linacc.neig
2537             +edi->vecs.radfix.neig
2538             +edi->vecs.radacc.neig
2539             +edi->vecs.radcon.neig
2540             + 6*edi->flood.vecs.neig;
2541         edi = edi->next_edi;
2542     }
2543     snew(setname, nsets);
2544
2545     /* In the mdrun time step in a first function call (do_flood()) the flooding
2546      * forces are calculated and in a second function call (do_edsam()) the
2547      * ED constraints. To get a corresponding legend, we need to loop twice
2548      * over the edi groups and output first the flooding, then the ED part */
2549
2550     /* The flooding-related legend entries, if flooding is done */
2551     nsets = 0;
2552     if (eEDflood == ed->eEDtype)
2553     {
2554         edi   = ed->edpar;
2555         for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2556         {
2557             /* Always write out the projection on the flooding EVs. Of course, this can also
2558              * be achieved with the monitoring option in do_edsam() (if switched on by the
2559              * user), but in that case the positions need to be communicated in do_edsam(),
2560              * which is not necessary when doing flooding only. */
2561             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2562
2563             for (i = 0; i < edi->flood.vecs.neig; i++)
2564             {
2565                 sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
2566                 nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2567
2568                 /* Output the current reference projection if it changes with time;
2569                  * this can happen when flooding is used as harmonic restraint */
2570                 if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
2571                 {
2572                     sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
2573                     nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
2574                 }
2575
2576                 /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
2577                 if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
2578                 {
2579                     sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
2580                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2581                 }
2582
2583                 sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
2584                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2585
2586                 if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
2587                 {
2588                     sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
2589                     nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
2590                 }
2591
2592                 sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
2593                 nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
2594             }
2595
2596             edi = edi->next_edi;
2597         } /* End of flooding-related legend entries */
2598     }
2599     n_flood = nsets;
2600
2601     /* Now the ED-related entries, if essential dynamics is done */
2602     edi         = ed->edpar;
2603     for (nr_edi = 1; nr_edi <= nED; nr_edi++)
2604     {
2605         if (edi->bNeedDoEdsam) /* Only print ED legend if at least one ED option is on */
2606         {
2607             nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
2608
2609             /* Essential dynamics, projections on eigenvectors */
2610             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON"   );
2611             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
2612             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
2613             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
2614             if (edi->vecs.radfix.neig)
2615             {
2616                 nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
2617             }
2618             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
2619             if (edi->vecs.radacc.neig)
2620             {
2621                 nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
2622             }
2623             nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
2624             if (edi->vecs.radcon.neig)
2625             {
2626                 nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
2627             }
2628         }
2629         edi = edi->next_edi;
2630     } /* end of 'pure' essential dynamics legend entries */
2631     n_edsam = nsets - n_flood;
2632
2633     xvgr_legend(ed->edo, nsets, setname, oenv);
2634     sfree(setname);
2635
2636     fprintf(ed->edo, "#\n"
2637             "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
2638             n_flood, 1 == n_flood ? "" : "s",
2639             n_edsam, 1 == n_edsam ? "" : "s");
2640     fprintf(ed->edo, "%s", LegendStr);
2641     sfree(LegendStr);
2642
2643     fflush(ed->edo);
2644 }
2645
2646
2647 void init_edsam(gmx_mtop_t   *mtop,  /* global topology                    */
2648                 t_inputrec   *ir,    /* input record                       */
2649                 t_commrec    *cr,    /* communication record               */
2650                 gmx_edsam_t   ed,    /* contains all ED data               */
2651                 rvec          x[],   /* positions of the whole MD system   */
2652                 matrix        box,   /* the box                            */
2653                 edsamstate_t *EDstate)
2654 {
2655     t_edpar *edi = NULL;                    /* points to a single edi data set */
2656     int      i, nr_edi, avindex;
2657     rvec    *x_pbc  = NULL;                 /* positions of the whole MD system with pbc removed  */
2658     rvec    *xfit   = NULL, *xstart = NULL; /* dummy arrays to determine initial RMSDs  */
2659     rvec     fit_transvec;                  /* translation ... */
2660     matrix   fit_rotmat;                    /* ... and rotation from fit to reference structure */
2661
2662
2663     if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
2664     {
2665         gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
2666     }
2667
2668     if (MASTER(cr))
2669     {
2670         fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2671
2672         if (NULL == ed)
2673         {
2674             gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n"
2675                       "flooding simulation. Please also provide the correct .edi file with -ei.\n");
2676         }
2677     }
2678
2679     /* Needed for initializing radacc radius in do_edsam */
2680     ed->bFirst = TRUE;
2681
2682     /* The input file is read by the master and the edi structures are
2683      * initialized here. Input is stored in ed->edpar. Then the edi
2684      * structures are transferred to the other nodes */
2685     if (MASTER(cr))
2686     {
2687         /* Initialization for every ED/flooding group. Flooding uses one edi group per
2688          * flooding vector, Essential dynamics can be applied to more than one structure
2689          * as well, but will be done in the order given in the edi file, so
2690          * expect different results for different order of edi file concatenation! */
2691         edi = ed->edpar;
2692         while (edi != NULL)
2693         {
2694             init_edi(mtop, edi);
2695             init_flood(edi, ed, ir->delta_t);
2696             edi = edi->next_edi;
2697         }
2698     }
2699
2700     /* The master does the work here. The other nodes get the positions
2701      * not before dd_partition_system which is called after init_edsam */
2702     if (MASTER(cr))
2703     {
2704         /* Remove pbc, make molecule whole.
2705          * When ir->bContinuation=TRUE this has already been done, but ok.
2706          */
2707         snew(x_pbc, mtop->natoms);
2708         m_rveccopy(mtop->natoms, x, x_pbc);
2709         do_pbc_first_mtop(NULL, ir->ePBC, box, mtop, x_pbc);
2710
2711         /* Reset pointer to first ED data set which contains the actual ED data */
2712         edi = ed->edpar;
2713         /* Loop over all ED/flooding data sets (usually only one, though) */
2714         for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2715         {
2716             /* For multiple ED groups we use the output frequency that was specified
2717              * in the first set */
2718             if (nr_edi > 1)
2719             {
2720                 edi->outfrq = ed->edpar->outfrq;
2721             }
2722
2723             /* Extract the initial reference and average positions. When starting
2724              * from .cpt, these have already been read into sref.x_old
2725              * in init_edsamstate() */
2726             if (!EDstate->bFromCpt)
2727             {
2728                 /* If this is the first run (i.e. no checkpoint present) we assume
2729                  * that the starting positions give us the correct PBC representation */
2730                 for (i = 0; i < edi->sref.nr; i++)
2731                 {
2732                     copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2733                 }
2734
2735                 for (i = 0; i < edi->sav.nr; i++)
2736                 {
2737                     copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2738                 }
2739             }
2740
2741             /* Now we have the PBC-correct start positions of the reference and
2742                average structure. We copy that over to dummy arrays on which we
2743                can apply fitting to print out the RMSD. We srenew the memory since
2744                the size of the buffers is likely different for every ED group */
2745             srenew(xfit, edi->sref.nr );
2746             srenew(xstart, edi->sav.nr  );
2747             copy_rvecn(edi->sref.x_old, xfit, 0, edi->sref.nr);
2748             copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2749
2750             /* Make the fit to the REFERENCE structure, get translation and rotation */
2751             fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2752
2753             /* Output how well we fit to the reference at the start */
2754             translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2755             fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
2756                     rmsd_from_structure(xfit, &edi->sref));
2757             if (EDstate->nED > 1)
2758             {
2759                 fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
2760             }
2761             fprintf(stderr, "\n");
2762
2763             /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2764             translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2765
2766             /* calculate initial projections */
2767             project(xstart, edi);
2768
2769             /* For the target and origin structure both a reference (fit) and an
2770              * average structure can be provided in make_edi. If both structures
2771              * are the same, make_edi only stores one of them in the .edi file.
2772              * If they differ, first the fit and then the average structure is stored
2773              * in star (or sor), thus the number of entries in star/sor is
2774              * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2775              * the size of the average group. */
2776
2777             /* process target structure, if required */
2778             if (edi->star.nr > 0)
2779             {
2780                 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2781
2782                 /* get translation & rotation for fit of target structure to reference structure */
2783                 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2784                 /* do the fit */
2785                 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2786                 if (edi->star.nr == edi->sav.nr)
2787                 {
2788                     avindex = 0;
2789                 }
2790                 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2791                 {
2792                     /* The last sav.nr indices of the target structure correspond to
2793                      * the average structure, which must be projected */
2794                     avindex = edi->star.nr - edi->sav.nr;
2795                 }
2796                 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
2797             }
2798             else
2799             {
2800                 rad_project(edi, xstart, &edi->vecs.radcon);
2801             }
2802
2803             /* process structure that will serve as origin of expansion circle */
2804             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2805             {
2806                 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2807             }
2808
2809             if (edi->sori.nr > 0)
2810             {
2811                 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2812
2813                 /* fit this structure to reference structure */
2814                 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2815                 /* do the fit */
2816                 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2817                 if (edi->sori.nr == edi->sav.nr)
2818                 {
2819                     avindex = 0;
2820                 }
2821                 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2822                 {
2823                     /* For the projection, we need the last sav.nr indices of sori */
2824                     avindex = edi->sori.nr - edi->sav.nr;
2825                 }
2826
2827                 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
2828                 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
2829                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2830                 {
2831                     fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2832                     /* Set center of flooding potential to the ORIGIN structure */
2833                     rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
2834                     /* We already know that no (moving) reference position was provided,
2835                      * therefore we can overwrite refproj[0]*/
2836                     copyEvecReference(&edi->flood.vecs);
2837                 }
2838             }
2839             else /* No origin structure given */
2840             {
2841                 rad_project(edi, xstart, &edi->vecs.radacc);
2842                 rad_project(edi, xstart, &edi->vecs.radfix);
2843                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2844                 {
2845                     if (edi->flood.bHarmonic)
2846                     {
2847                         fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2848                         for (i = 0; i < edi->flood.vecs.neig; i++)
2849                         {
2850                             edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2851                         }
2852                     }
2853                     else
2854                     {
2855                         fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2856                         /* Set center of flooding potential to the center of the covariance matrix,
2857                          * i.e. the average structure, i.e. zero in the projected system */
2858                         for (i = 0; i < edi->flood.vecs.neig; i++)
2859                         {
2860                             edi->flood.vecs.refproj[i] = 0.0;
2861                         }
2862                     }
2863                 }
2864             }
2865             /* For convenience, output the center of the flooding potential for the eigenvectors */
2866             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2867             {
2868                 for (i = 0; i < edi->flood.vecs.neig; i++)
2869                 {
2870                     fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2871                     if (edi->flood.bHarmonic)
2872                     {
2873                         fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2874                     }
2875                     fprintf(stdout, "\n");
2876                 }
2877             }
2878
2879             /* set starting projections for linsam */
2880             rad_project(edi, xstart, &edi->vecs.linacc);
2881             rad_project(edi, xstart, &edi->vecs.linfix);
2882
2883             /* Prepare for the next edi data set: */
2884             edi = edi->next_edi;
2885         }
2886         /* Cleaning up on the master node: */
2887         sfree(x_pbc);
2888         sfree(xfit);
2889         sfree(xstart);
2890
2891     } /* end of MASTER only section */
2892
2893     if (PAR(cr))
2894     {
2895         /* First let everybody know how many ED data sets to expect */
2896         gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr);
2897         /* Broadcast the essential dynamics / flooding data to all nodes */
2898         broadcast_ed_data(cr, ed, EDstate->nED);
2899     }
2900     else
2901     {
2902         /* In the single-CPU case, point the local atom numbers pointers to the global
2903          * one, so that we can use the same notation in serial and parallel case: */
2904
2905         /* Loop over all ED data sets (usually only one, though) */
2906         edi = ed->edpar;
2907         for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2908         {
2909             edi->sref.anrs_loc = edi->sref.anrs;
2910             edi->sav.anrs_loc  = edi->sav.anrs;
2911             edi->star.anrs_loc = edi->star.anrs;
2912             edi->sori.anrs_loc = edi->sori.anrs;
2913             /* For the same reason as above, make a dummy c_ind array: */
2914             snew(edi->sav.c_ind, edi->sav.nr);
2915             /* Initialize the array */
2916             for (i = 0; i < edi->sav.nr; i++)
2917             {
2918                 edi->sav.c_ind[i] = i;
2919             }
2920             /* In the general case we will need a different-sized array for the reference indices: */
2921             if (!edi->bRefEqAv)
2922             {
2923                 snew(edi->sref.c_ind, edi->sref.nr);
2924                 for (i = 0; i < edi->sref.nr; i++)
2925                 {
2926                     edi->sref.c_ind[i] = i;
2927                 }
2928             }
2929             /* Point to the very same array in case of other structures: */
2930             edi->star.c_ind = edi->sav.c_ind;
2931             edi->sori.c_ind = edi->sav.c_ind;
2932             /* In the serial case, the local number of atoms is the global one: */
2933             edi->sref.nr_loc = edi->sref.nr;
2934             edi->sav.nr_loc  = edi->sav.nr;
2935             edi->star.nr_loc = edi->star.nr;
2936             edi->sori.nr_loc = edi->sori.nr;
2937
2938             /* An on we go to the next ED group */
2939             edi = edi->next_edi;
2940         }
2941     }
2942
2943     /* Allocate space for ED buffer variables */
2944     /* Again, loop over ED data sets */
2945     edi = ed->edpar;
2946     for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2947     {
2948         /* Allocate space for ED buffer */
2949         snew(edi->buf, 1);
2950         snew(edi->buf->do_edsam, 1);
2951
2952         /* Space for collective ED buffer variables */
2953
2954         /* Collective positions of atoms with the average indices */
2955         snew(edi->buf->do_edsam->xcoll, edi->sav.nr);
2956         snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr);            /* buffer for xcoll shifts */
2957         snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr);
2958         /* Collective positions of atoms with the reference indices */
2959         if (!edi->bRefEqAv)
2960         {
2961             snew(edi->buf->do_edsam->xc_ref, edi->sref.nr);
2962             snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr);       /* To store the shifts in */
2963             snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2964         }
2965
2966         /* Get memory for flooding forces */
2967         snew(edi->flood.forces_cartesian, edi->sav.nr);
2968
2969 #ifdef DUMPEDI
2970         /* Dump it all into one file per process */
2971         dump_edi(edi, cr, nr_edi);
2972 #endif
2973
2974         /* Next ED group */
2975         edi = edi->next_edi;
2976     }
2977
2978     /* Flush the edo file so that the user can check some things
2979      * when the simulation has started */
2980     if (ed->edo)
2981     {
2982         fflush(ed->edo);
2983     }
2984 }
2985
2986
2987 void do_edsam(t_inputrec     *ir,
2988               gmx_large_int_t step,
2989               t_commrec      *cr,
2990               rvec            xs[], /* The local current positions on this processor */
2991               rvec            v[],  /* The velocities */
2992               matrix          box,
2993               gmx_edsam_t     ed)
2994 {
2995     int                i, edinr, iupdate = 500;
2996     matrix             rotmat;         /* rotation matrix */
2997     rvec               transvec;       /* translation vector */
2998     rvec               dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
2999     real               dt_1;           /* 1/dt */
3000     struct t_do_edsam *buf;
3001     t_edpar           *edi;
3002     real               rmsdev    = -1;    /* RMSD from reference structure prior to applying the constraints */
3003     gmx_bool           bSuppress = FALSE; /* Write .xvg output file on master? */
3004
3005
3006     /* Check if ED sampling has to be performed */
3007     if (ed->eEDtype == eEDnone)
3008     {
3009         return;
3010     }
3011
3012     /* Suppress output on first call of do_edsam if
3013      * two-step sd2 integrator is used */
3014     if ( (ir->eI == eiSD2) && (v != NULL) )
3015     {
3016         bSuppress = TRUE;
3017     }
3018
3019     dt_1 = 1.0/ir->delta_t;
3020
3021     /* Loop over all ED groups (usually one) */
3022     edi   = ed->edpar;
3023     edinr = 0;
3024     while (edi != NULL)
3025     {
3026         edinr++;
3027         if (edi->bNeedDoEdsam)
3028         {
3029
3030             buf = edi->buf->do_edsam;
3031
3032             if (ed->bFirst)
3033             {
3034                 /* initialise radacc radius for slope criterion */
3035                 buf->oldrad = calc_radius(&edi->vecs.radacc);
3036             }
3037
3038             /* Copy the positions into buf->xc* arrays and after ED
3039              * feed back corrections to the official positions */
3040
3041             /* Broadcast the ED positions such that every node has all of them
3042              * Every node contributes its local positions xs and stores it in
3043              * the collective buf->xcoll array. Note that for edinr > 1
3044              * xs could already have been modified by an earlier ED */
3045
3046             communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3047                                         edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old,  box);
3048
3049             /* Only assembly reference positions if their indices differ from the average ones */
3050             if (!edi->bRefEqAv)
3051             {
3052                 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
3053                                             edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
3054             }
3055
3056             /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
3057              * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
3058              * set bUpdateShifts=TRUE in the parallel case. */
3059             buf->bUpdateShifts = FALSE;
3060
3061             /* Now all nodes have all of the ED positions in edi->sav->xcoll,
3062              * as well as the indices in edi->sav.anrs */
3063
3064             /* Fit the reference indices to the reference structure */
3065             if (edi->bRefEqAv)
3066             {
3067                 fit_to_reference(buf->xcoll, transvec, rotmat, edi);
3068             }
3069             else
3070             {
3071                 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
3072             }
3073
3074             /* Now apply the translation and rotation to the ED structure */
3075             translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
3076
3077             /* Find out how well we fit to the reference (just for output steps) */
3078             if (do_per_step(step, edi->outfrq) && MASTER(cr))
3079             {
3080                 if (edi->bRefEqAv)
3081                 {
3082                     /* Indices of reference and average structures are identical,
3083                      * thus we can calculate the rmsd to SREF using xcoll */
3084                     rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref);
3085                 }
3086                 else
3087                 {
3088                     /* We have to translate & rotate the reference atoms first */
3089                     translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
3090                     rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref);
3091                 }
3092             }
3093
3094             /* update radsam references, when required */
3095             if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps)
3096             {
3097                 project(buf->xcoll, edi);
3098                 rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3099                 rad_project(edi, buf->xcoll, &edi->vecs.radfix);
3100                 buf->oldrad = -1.e5;
3101             }
3102
3103             /* update radacc references, when required */
3104             if (do_per_step(step, iupdate) && step >= edi->presteps)
3105             {
3106                 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
3107                 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
3108                 {
3109                     project(buf->xcoll, edi);
3110                     rad_project(edi, buf->xcoll, &edi->vecs.radacc);
3111                     buf->oldrad = 0.0;
3112                 }
3113                 else
3114                 {
3115                     buf->oldrad = edi->vecs.radacc.radius;
3116                 }
3117             }
3118
3119             /* apply the constraints */
3120             if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
3121             {
3122                 /* ED constraints should be applied already in the first MD step
3123                  * (which is step 0), therefore we pass step+1 to the routine */
3124                 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
3125             }
3126
3127             /* write to edo, when required */
3128             if (do_per_step(step, edi->outfrq))
3129             {
3130                 project(buf->xcoll, edi);
3131                 if (MASTER(cr) && !bSuppress)
3132                 {
3133                     write_edo(edi, ed->edo, rmsdev);
3134                 }
3135             }
3136
3137             /* Copy back the positions unless monitoring only */
3138             if (ed_constraints(ed->eEDtype, edi))
3139             {
3140                 /* remove fitting */
3141                 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
3142
3143                 /* Copy the ED corrected positions into the coordinate array */
3144                 /* Each node copies its local part. In the serial case, nat_loc is the
3145                  * total number of ED atoms */
3146                 for (i = 0; i < edi->sav.nr_loc; i++)
3147                 {
3148                     /* Unshift local ED coordinate and store in x_unsh */
3149                     ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
3150                                             buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
3151
3152                     /* dx is the ED correction to the positions: */
3153                     rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
3154
3155                     if (v != NULL)
3156                     {
3157                         /* dv is the ED correction to the velocity: */
3158                         svmul(dt_1, dx, dv);
3159                         /* apply the velocity correction: */
3160                         rvec_inc(v[edi->sav.anrs_loc[i]], dv);
3161                     }
3162                     /* Finally apply the position correction due to ED: */
3163                     copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
3164                 }
3165             }
3166         } /* END of if (edi->bNeedDoEdsam) */
3167
3168         /* Prepare for the next ED group */
3169         edi = edi->next_edi;
3170
3171     } /* END of loop over ED groups */
3172
3173     ed->bFirst = FALSE;
3174 }