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