src/mdlib/edsam.c

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