src/gromacs/mdlib/edsam.c

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