src/mdlib/edsam.c

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