src/mdlib/edsam.c

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