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     gmx_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     gmx_bool       fitmas;         /* true if trans fit with cm            */
 149     gmx_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     gmx_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     gmx_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     gmx_bool      bFirst;
 183     gmx_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     gmx_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     gmx_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     gmx_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, gmx_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,gmx_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,gmx_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                 /* Zero out values which were not scanned */
1349                 switch(nscan)
1350                 {
1351                     case 4:
1352                         /* Every 4 values read */
1353                         break;
1354                     case 3:
1355                         /* No value for slope, set to 0 */
1356                         refprojslope_dum = 0.0;
1357                         break;
1358                     case 2:
1359                         /* No values for reference projection and slope, set to 0 */
1360                         refproj_dum      = 0.0;
1361                         refprojslope_dum = 0.0;
1362                         break;
1363                     default:
1364                         gmx_fatal(FARGS,"Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1365                         break;
1366                 }
1367                 tvec->refproj[i]=refproj_dum;
1368                 tvec->refproj0[i]=refproj_dum;
1369                 tvec->refprojslope[i]=refprojslope_dum;
1370             }
1371             else /* Normal flooding */
1372             {
1373                 nscan = sscanf(line,"%d%lf",&idum,&rdum);
1374                 if (nscan != 2)
1375                     gmx_fatal(FARGS,"Expected 2 values for flooding vec: <nr> <stpsz>\n");
1376             }
1377             tvec->ieig[i]=idum;
1378             tvec->stpsz[i]=rdum;
1379         } /* end of loop over eigenvectors */
1380
1381         for(i=0; (i < tvec->neig); i++)
1382         {
1383             snew(tvec->vec[i],nr);
1384             scan_edvec(in,nr,tvec->vec[i]);
1385         }
1386     }
1387 }
1388
1389
1390 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1391 static void read_edvecs(FILE *in,int nr,t_edvecs *vecs)
1392 {
1393     read_edvec(in,nr,&vecs->mon   ,FALSE);
1394     read_edvec(in,nr,&vecs->linfix,FALSE);
1395     read_edvec(in,nr,&vecs->linacc,FALSE);
1396     read_edvec(in,nr,&vecs->radfix,FALSE);
1397     read_edvec(in,nr,&vecs->radacc,FALSE);
1398     read_edvec(in,nr,&vecs->radcon,FALSE);
1399 }
1400
1401
1402 /* Check if the same atom indices are used for reference and average positions */
1403 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1404 {
1405     int i;
1406
1407
1408     /* If the number of atoms differs between the two structures,
1409      * they cannot be identical */
1410     if (sref.nr != sav.nr)
1411         return FALSE;
1412
1413     /* Now that we know that both stuctures have the same number of atoms,
1414      * check if also the indices are identical */
1415     for (i=0; i < sav.nr; i++)
1416     {
1417         if (sref.anrs[i] != sav.anrs[i])
1418             return FALSE;
1419     }
1420     fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1421
1422     return TRUE;
1423 }
1424
1425
1426 static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int edi_nr, t_commrec *cr)
1427 {
1428     int readmagic;
1429     const int magic=669;
1430     gmx_bool bEOF;
1431
1432
1433     /* the edi file is not free format, so expect problems if the input is corrupt. */
1434
1435     /* check the magic number */
1436     readmagic=read_edint(in,&bEOF);
1437     /* Check whether we have reached the end of the input file */
1438     if (bEOF)
1439         return 0;
1440
1441     if (readmagic != magic)
1442     {
1443         if (readmagic==666 || readmagic==667 || readmagic==668)
1444             gmx_fatal(FARGS,"Wrong magic number: Use newest version of make_edi to produce edi file");
1445         else
1446             gmx_fatal(FARGS,"Wrong magic number %d in %s",readmagic,ed->edinam);
1447     }
1448
1449     /* check the number of atoms */
1450     edi->nini=read_edint(in,&bEOF);
1451     if (edi->nini != nr_mdatoms)
1452         gmx_fatal(FARGS,"Nr of atoms in %s (%d) does not match nr of md atoms (%d)",
1453                 ed->edinam,edi->nini,nr_mdatoms);
1454
1455     /* Done checking. For the rest we blindly trust the input */
1456     edi->fitmas          = read_checked_edint(in,"FITMAS");
1457     edi->pcamas          = read_checked_edint(in,"ANALYSIS_MAS");
1458     edi->outfrq          = read_checked_edint(in,"OUTFRQ");
1459     edi->maxedsteps      = read_checked_edint(in,"MAXLEN");
1460     edi->slope           = read_checked_edreal(in,"SLOPECRIT");
1461
1462     edi->presteps        = read_checked_edint(in,"PRESTEPS");
1463     edi->flood.deltaF0   = read_checked_edreal(in,"DELTA_F0");
1464     edi->flood.deltaF    = read_checked_edreal(in,"INIT_DELTA_F");
1465     edi->flood.tau       = read_checked_edreal(in,"TAU");
1466     edi->flood.constEfl  = read_checked_edreal(in,"EFL_NULL");
1467     edi->flood.alpha2    = read_checked_edreal(in,"ALPHA2");
1468     edi->flood.kT        = read_checked_edreal(in,"KT");
1469     edi->flood.bHarmonic = read_checked_edint(in,"HARMONIC");
1470     edi->flood.flood_id  = edi_nr;
1471     edi->sref.nr         = read_checked_edint(in,"NREF");
1472
1473     /* allocate space for reference positions and read them */
1474     snew(edi->sref.anrs,edi->sref.nr);
1475     snew(edi->sref.x   ,edi->sref.nr);
1476     if (PAR(cr))
1477         snew(edi->sref.x_old,edi->sref.nr);
1478     edi->sref.sqrtm    =NULL;
1479     read_edx(in,edi->sref.nr,edi->sref.anrs,edi->sref.x);
1480
1481     /* average positions. they define which atoms will be used for ED sampling */
1482     edi->sav.nr=read_checked_edint(in,"NAV");
1483     snew(edi->sav.anrs,edi->sav.nr);
1484     snew(edi->sav.x   ,edi->sav.nr);
1485     if (PAR(cr))
1486         snew(edi->sav.x_old,edi->sav.nr);
1487     read_edx(in,edi->sav.nr,edi->sav.anrs,edi->sav.x);
1488
1489     /* Check if the same atom indices are used for reference and average positions */
1490     edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1491
1492     /* eigenvectors */
1493     read_edvecs(in,edi->sav.nr,&edi->vecs);
1494     read_edvec(in,edi->sav.nr,&edi->flood.vecs,edi->flood.bHarmonic);
1495
1496     /* target positions */
1497     edi->star.nr=read_edint(in,&bEOF);
1498     if (edi->star.nr > 0)
1499     {
1500         snew(edi->star.anrs,edi->star.nr);
1501         snew(edi->star.x   ,edi->star.nr);
1502         edi->star.sqrtm    =NULL;
1503         read_edx(in,edi->star.nr,edi->star.anrs,edi->star.x);
1504     }
1505
1506     /* positions defining origin of expansion circle */
1507     edi->sori.nr=read_edint(in,&bEOF);
1508     if (edi->sori.nr > 0)
1509     {
1510         snew(edi->sori.anrs,edi->sori.nr);
1511         snew(edi->sori.x   ,edi->sori.nr);
1512         edi->sori.sqrtm    =NULL;
1513         read_edx(in,edi->sori.nr,edi->sori.anrs,edi->sori.x);
1514     }
1515
1516     /* all done */
1517     return 1;
1518 }
1519
1520
1521
1522 /* Read in the edi input file. Note that it may contain several ED data sets which were
1523  * achieved by concatenating multiple edi files. The standard case would be a single ED
1524  * data set, though. */
1525 static void read_edi_file(gmx_edsam_t ed, t_edpar *edi, int nr_mdatoms, t_commrec *cr)
1526 {
1527     FILE    *in;
1528     t_edpar *curr_edi,*last_edi;
1529     t_edpar *edi_read;
1530     int     edi_nr = 0;
1531
1532
1533     /* This routine is executed on the master only */
1534
1535     /* Open the .edi parameter input file */
1536     in = gmx_fio_fopen(ed->edinam,"r");
1537     fprintf(stderr, "ED: Reading edi file %s\n", ed->edinam);
1538
1539     /* Now read a sequence of ED input parameter sets from the edi file */
1540     curr_edi=edi;
1541     last_edi=edi;
1542     while( read_edi(in, ed, curr_edi, nr_mdatoms, edi_nr, cr) )
1543     {
1544         edi_nr++;
1545         /* Make shure that the number of atoms in each dataset is the same as in the tpr file */
1546         if (edi->nini != nr_mdatoms)
1547             gmx_fatal(FARGS,"edi file %s (dataset #%d) was made for %d atoms, but the simulation contains %d atoms.",
1548                     ed->edinam, edi_nr, edi->nini, nr_mdatoms);
1549         /* Since we arrived within this while loop we know that there is still another data set to be read in */
1550         /* We need to allocate space for the data: */
1551         snew(edi_read,1);
1552         /* Point the 'next_edi' entry to the next edi: */
1553         curr_edi->next_edi=edi_read;
1554         /* Keep the curr_edi pointer for the case that the next dataset is empty: */
1555         last_edi = curr_edi;
1556         /* Let's prepare to read in the next edi data set: */
1557         curr_edi = edi_read;
1558     }
1559     if (edi_nr == 0)
1560         gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", ed->edinam);
1561
1562     /* Terminate the edi dataset list with a NULL pointer: */
1563     last_edi->next_edi = NULL;
1564
1565     fprintf(stderr, "ED: Found %d ED dataset%s.\n", edi_nr, edi_nr>1? "s" : "");
1566
1567     /* Close the .edi file again */
1568     gmx_fio_fclose(in);
1569 }
1570
1571
1572 struct t_fit_to_ref {
1573     rvec *xcopy;       /* Working copy of the positions in fit_to_reference */
1574 };
1575
1576 /* Fit the current positions to the reference positions
1577  * Do not actually do the fit, just return rotation and translation.
1578  * Note that the COM of the reference structure was already put into
1579  * the origin by init_edi. */
1580 static void fit_to_reference(rvec      *xcoll,    /* The positions to be fitted */
1581                              rvec      transvec,  /* The translation vector */
1582                              matrix    rotmat,    /* The rotation matrix */
1583                              t_edpar   *edi)      /* Just needed for do_edfit */
1584 {
1585     rvec com;          /* center of mass */
1586     int  i;
1587     struct t_fit_to_ref *loc;
1588
1589
1590     GMX_MPE_LOG(ev_fit_to_reference_start);
1591
1592     /* Allocate memory the first time this routine is called for each edi dataset */
1593     if (NULL == edi->buf->fit_to_ref)
1594     {
1595         snew(edi->buf->fit_to_ref, 1);
1596         snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1597     }
1598     loc = edi->buf->fit_to_ref;
1599
1600     /* We do not touch the original positions but work on a copy. */
1601     for (i=0; i<edi->sref.nr; i++)
1602         copy_rvec(xcoll[i], loc->xcopy[i]);
1603
1604     /* Calculate the center of mass */
1605     get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1606
1607     transvec[XX] = -com[XX];
1608     transvec[YY] = -com[YY];
1609     transvec[ZZ] = -com[ZZ];
1610
1611     /* Subtract the center of mass from the copy */
1612     translate_x(loc->xcopy, edi->sref.nr, transvec);
1613
1614     /* Determine the rotation matrix */
1615     do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1616
1617     GMX_MPE_LOG(ev_fit_to_reference_finish);
1618 }
1619
1620
1621 static void translate_and_rotate(rvec *x,         /* The positions to be translated and rotated */
1622                                  int nat,         /* How many positions are there? */
1623                                  rvec transvec,   /* The translation vector */
1624                                  matrix rotmat)   /* The rotation matrix */
1625 {
1626     /* Translation */
1627     translate_x(x, nat, transvec);
1628
1629     /* Rotation */
1630     rotate_x(x, nat, rotmat);
1631 }
1632
1633
1634 /* Gets the rms deviation of the positions to the structure s */
1635 /* fit_to_structure has to be called before calling this routine! */
1636 static real rmsd_from_structure(rvec           *x,  /* The positions under consideration */
1637                                 struct gmx_edx *s)  /* The structure from which the rmsd shall be computed */
1638 {
1639     real  rmsd=0.0;
1640     int   i;
1641
1642
1643     for (i=0; i < s->nr; i++)
1644         rmsd += distance2(s->x[i], x[i]);
1645
1646     rmsd /= (real) s->nr;
1647     rmsd = sqrt(rmsd);
1648
1649     return rmsd;
1650 }
1651
1652
1653 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1654 {
1655     t_edpar *edi;
1656
1657
1658     if (ed->eEDtype != eEDnone)
1659     {
1660         /* Loop over ED datasets (usually there is just one dataset, though) */
1661         edi=ed->edpar;
1662         while (edi)
1663         {
1664             /* Local atoms of the reference structure (for fitting), need only be assembled
1665              * if their indices differ from the average ones */
1666             if (!edi->bRefEqAv)
1667                 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1668                         &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1669
1670             /* Local atoms of the average structure (on these ED will be performed) */
1671             dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1672                     &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1673
1674             /* Indicate that the ED shift vectors for this structure need to be updated
1675              * at the next call to communicate_group_positions, since obviously we are in a NS step */
1676             edi->buf->do_edsam->bUpdateShifts = TRUE;
1677
1678             /* Set the pointer to the next ED dataset (if any) */
1679             edi=edi->next_edi;
1680         }
1681     }
1682 }
1683
1684
1685 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1686 {
1687     int tx,ty,tz;
1688
1689
1690     GMX_MPE_LOG(ev_unshift_start);
1691
1692     tx=is[XX];
1693     ty=is[YY];
1694     tz=is[ZZ];
1695
1696     if(TRICLINIC(box))
1697     {
1698         xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1699         xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1700         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1701     } else
1702     {
1703         xu[XX] = x[XX]-tx*box[XX][XX];
1704         xu[YY] = x[YY]-ty*box[YY][YY];
1705         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1706     }
1707
1708     GMX_MPE_LOG(ev_unshift_finish);
1709 }
1710
1711
1712 static void do_linfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
1713 {
1714     int  i, j;
1715     real proj, add;
1716     rvec vec_dum;
1717
1718
1719     /* loop over linfix vectors */
1720     for (i=0; i<edi->vecs.linfix.neig; i++)
1721     {
1722         /* calculate the projection */
1723         proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1724
1725         /* calculate the correction */
1726         add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1727
1728         /* apply the correction */
1729         add /= edi->sav.sqrtm[i];
1730         for (j=0; j<edi->sav.nr; j++)
1731         {
1732             svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1733             rvec_inc(xcoll[j], vec_dum);
1734         }
1735     }
1736 }
1737
1738
1739 static void do_linacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1740 {
1741     int  i, j;
1742     real proj, add;
1743     rvec vec_dum;
1744
1745
1746     /* loop over linacc vectors */
1747     for (i=0; i<edi->vecs.linacc.neig; i++)
1748     {
1749         /* calculate the projection */
1750         proj=projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1751
1752         /* calculate the correction */
1753         add = 0.0;
1754         if (edi->vecs.linacc.stpsz[i] > 0.0)
1755         {
1756             if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1757                 add = edi->vecs.linacc.refproj[i] - proj;
1758         }
1759         if (edi->vecs.linacc.stpsz[i] < 0.0)
1760         {
1761             if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1762                 add = edi->vecs.linacc.refproj[i] - proj;
1763         }
1764
1765         /* apply the correction */
1766         add /= edi->sav.sqrtm[i];
1767         for (j=0; j<edi->sav.nr; j++)
1768         {
1769             svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
1770             rvec_inc(xcoll[j], vec_dum);
1771         }
1772
1773         /* new positions will act as reference */
1774         edi->vecs.linacc.refproj[i] = proj + add;
1775     }
1776 }
1777
1778
1779 static void do_radfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
1780 {
1781     int  i,j;
1782     real *proj, rad=0.0, ratio;
1783     rvec vec_dum;
1784
1785
1786     if (edi->vecs.radfix.neig == 0)
1787         return;
1788
1789     snew(proj, edi->vecs.radfix.neig);
1790
1791     /* loop over radfix vectors */
1792     for (i=0; i<edi->vecs.radfix.neig; i++)
1793     {
1794         /* calculate the projections, radius */
1795         proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
1796         rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
1797     }
1798
1799     rad   = sqrt(rad);
1800     ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
1801     edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
1802
1803     /* loop over radfix vectors */
1804     for (i=0; i<edi->vecs.radfix.neig; i++)
1805     {
1806         proj[i] -= edi->vecs.radfix.refproj[i];
1807
1808         /* apply the correction */
1809         proj[i] /= edi->sav.sqrtm[i];
1810         proj[i] *= ratio;
1811         for (j=0; j<edi->sav.nr; j++) {
1812             svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
1813             rvec_inc(xcoll[j], vec_dum);
1814         }
1815     }
1816
1817     sfree(proj);
1818 }
1819
1820
1821 static void do_radacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1822 {
1823     int  i,j;
1824     real *proj, rad=0.0, ratio=0.0;
1825     rvec vec_dum;
1826
1827
1828     if (edi->vecs.radacc.neig == 0)
1829         return;
1830
1831     snew(proj,edi->vecs.radacc.neig);
1832
1833     /* loop over radacc vectors */
1834     for (i=0; i<edi->vecs.radacc.neig; i++)
1835     {
1836         /* calculate the projections, radius */
1837         proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
1838         rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
1839     }
1840     rad = sqrt(rad);
1841
1842     /* only correct when radius decreased */
1843     if (rad < edi->vecs.radacc.radius)
1844     {
1845         ratio = edi->vecs.radacc.radius/rad - 1.0;
1846         rad   = edi->vecs.radacc.radius;
1847     }
1848     else
1849         edi->vecs.radacc.radius = rad;
1850
1851     /* loop over radacc vectors */
1852     for (i=0; i<edi->vecs.radacc.neig; i++)
1853     {
1854         proj[i] -= edi->vecs.radacc.refproj[i];
1855
1856         /* apply the correction */
1857         proj[i] /= edi->sav.sqrtm[i];
1858         proj[i] *= ratio;
1859         for (j=0; j<edi->sav.nr; j++)
1860         {
1861             svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
1862             rvec_inc(xcoll[j], vec_dum);
1863         }
1864     }
1865     sfree(proj);
1866 }
1867
1868
1869 struct t_do_radcon {
1870     real *proj;
1871 };
1872
1873 static void do_radcon(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1874 {
1875     int  i,j;
1876     real rad=0.0, ratio=0.0;
1877     struct t_do_radcon *loc;
1878     gmx_bool bFirst;
1879     rvec vec_dum;
1880
1881
1882     if(edi->buf->do_radcon != NULL)
1883     {
1884         bFirst = FALSE;
1885         loc    = edi->buf->do_radcon;
1886     }
1887     else
1888     {
1889         bFirst = TRUE;
1890         snew(edi->buf->do_radcon, 1);
1891     }
1892     loc = edi->buf->do_radcon;
1893
1894     if (edi->vecs.radcon.neig == 0)
1895         return;
1896
1897     if (bFirst)
1898         snew(loc->proj, edi->vecs.radcon.neig);
1899
1900     /* loop over radcon vectors */
1901     for (i=0; i<edi->vecs.radcon.neig; i++)
1902     {
1903         /* calculate the projections, radius */
1904         loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
1905         rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
1906     }
1907     rad = sqrt(rad);
1908     /* only correct when radius increased */
1909     if (rad > edi->vecs.radcon.radius)
1910     {
1911         ratio = edi->vecs.radcon.radius/rad - 1.0;
1912
1913         /* loop over radcon vectors */
1914         for (i=0; i<edi->vecs.radcon.neig; i++)
1915         {
1916             /* apply the correction */
1917             loc->proj[i] -= edi->vecs.radcon.refproj[i];
1918             loc->proj[i] /= edi->sav.sqrtm[i];
1919             loc->proj[i] *= ratio;
1920
1921             for (j=0; j<edi->sav.nr; j++)
1922             {
1923                 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
1924                 rvec_inc(xcoll[j], vec_dum);
1925             }
1926         }
1927     }
1928     else
1929         edi->vecs.radcon.radius = rad;
1930
1931     if (rad != edi->vecs.radcon.radius)
1932     {
1933         rad = 0.0;
1934         for (i=0; i<edi->vecs.radcon.neig; i++)
1935         {
1936             /* calculate the projections, radius */
1937             loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
1938             rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
1939         }
1940         rad = sqrt(rad);
1941     }
1942 }
1943
1944
1945 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step, t_commrec *cr)
1946 {
1947     int i;
1948
1949
1950     GMX_MPE_LOG(ev_ed_apply_cons_start);
1951
1952     /* subtract the average positions */
1953     for (i=0; i<edi->sav.nr; i++)
1954         rvec_dec(xcoll[i], edi->sav.x[i]);
1955
1956     /* apply the constraints */
1957     if (step >= 0)
1958         do_linfix(xcoll, edi, step, cr);
1959     do_linacc(xcoll, edi, cr);
1960     if (step >= 0)
1961         do_radfix(xcoll, edi, step, cr);
1962     do_radacc(xcoll, edi, cr);
1963     do_radcon(xcoll, edi, cr);
1964
1965     /* add back the average positions */
1966     for (i=0; i<edi->sav.nr; i++)
1967         rvec_inc(xcoll[i], edi->sav.x[i]);
1968
1969     GMX_MPE_LOG(ev_ed_apply_cons_finish);
1970 }
1971
1972
1973 /* Write out the projections onto the eigenvectors */
1974 static void write_edo(int nr_edi, t_edpar *edi, gmx_edsam_t ed, gmx_large_int_t step,real rmsd)
1975 {
1976     int i;
1977     char buf[22];
1978
1979
1980     if (edi->bNeedDoEdsam)
1981     {
1982         if (step == -1)
1983             fprintf(ed->edo, "Initial projections:\n");
1984         else
1985         {
1986             fprintf(ed->edo,"Step %s, ED #%d  ", gmx_step_str(step, buf), nr_edi);
1987             fprintf(ed->edo,"  RMSD %f nm\n",rmsd);
1988             if (ed->eEDtype == eEDflood)
1989                 fprintf(ed->edo, "  Efl=%f  deltaF=%f  Vfl=%f\n",edi->flood.Efl,edi->flood.deltaF,edi->flood.Vfl);
1990         }
1991
1992         if (edi->vecs.mon.neig)
1993         {
1994             fprintf(ed->edo,"  Monitor eigenvectors");
1995             for (i=0; i<edi->vecs.mon.neig; i++)
1996                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.mon.ieig[i],edi->vecs.mon.xproj[i]);
1997             fprintf(ed->edo,"\n");
1998         }
1999         if (edi->vecs.linfix.neig)
2000         {
2001             fprintf(ed->edo,"  Linfix  eigenvectors");
2002             for (i=0; i<edi->vecs.linfix.neig; i++)
2003                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.linfix.ieig[i],edi->vecs.linfix.xproj[i]);
2004             fprintf(ed->edo,"\n");
2005         }
2006         if (edi->vecs.linacc.neig)
2007         {
2008             fprintf(ed->edo,"  Linacc  eigenvectors");
2009             for (i=0; i<edi->vecs.linacc.neig; i++)
2010                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.linacc.ieig[i],edi->vecs.linacc.xproj[i]);
2011             fprintf(ed->edo,"\n");
2012         }
2013         if (edi->vecs.radfix.neig)
2014         {
2015             fprintf(ed->edo,"  Radfix  eigenvectors");
2016             for (i=0; i<edi->vecs.radfix.neig; i++)
2017                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radfix.ieig[i],edi->vecs.radfix.xproj[i]);
2018             fprintf(ed->edo,"\n");
2019             fprintf(ed->edo,"  fixed increment radius = %f\n", calc_radius(&edi->vecs.radfix));
2020         }
2021         if (edi->vecs.radacc.neig)
2022         {
2023             fprintf(ed->edo,"  Radacc  eigenvectors");
2024             for (i=0; i<edi->vecs.radacc.neig; i++)
2025                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radacc.ieig[i],edi->vecs.radacc.xproj[i]);
2026             fprintf(ed->edo,"\n");
2027             fprintf(ed->edo,"  acceptance radius      = %f\n", calc_radius(&edi->vecs.radacc));
2028         }
2029         if (edi->vecs.radcon.neig)
2030         {
2031             fprintf(ed->edo,"  Radcon  eigenvectors");
2032             for (i=0; i<edi->vecs.radcon.neig; i++)
2033                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radcon.ieig[i],edi->vecs.radcon.xproj[i]);
2034             fprintf(ed->edo,"\n");
2035             fprintf(ed->edo,"  contracting radius     = %f\n", calc_radius(&edi->vecs.radcon));
2036         }
2037     }
2038 }
2039
2040 /* Returns if any constraints are switched on */
2041 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2042 {
2043     if (edtype == eEDedsam || edtype == eEDflood)
2044     {
2045         return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2046                 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2047                 edi->vecs.radcon.neig);
2048     }
2049     return 0;
2050 }
2051
2052
2053 void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
2054                 t_inputrec  *ir,     /* input record                       */
2055                 t_commrec   *cr,     /* communication record               */
2056                 gmx_edsam_t ed,      /* contains all ED data               */
2057                 rvec        x[],     /* positions of the whole MD system   */
2058                 matrix      box)     /* the box                            */
2059 {
2060     t_edpar *edi = NULL;    /* points to a single edi data set */
2061     int     numedis=0;      /* keep track of the number of ED data sets in edi file */
2062     int     i,nr_edi;
2063     rvec    *x_pbc  = NULL; /* positions of the whole MD system with pbc removed  */
2064     rvec    *xfit   = NULL; /* the positions which will be fitted to the reference structure  */
2065     rvec    *xstart = NULL; /* the positions which are subject to ED sampling */
2066     rvec    fit_transvec;   /* translation ... */
2067     matrix  fit_rotmat;     /* ... and rotation from fit to reference structure */
2068
2069
2070     if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
2071         gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
2072
2073     GMX_MPE_LOG(ev_edsam_start);
2074
2075     if (MASTER(cr))
2076         fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2077
2078     /* Needed for initializing radacc radius in do_edsam */
2079     ed->bFirst = 1;
2080
2081     /* The input file is read by the master and the edi structures are
2082      * initialized here. Input is stored in ed->edpar. Then the edi
2083      * structures are transferred to the other nodes */
2084     if (MASTER(cr))
2085     {
2086         snew(ed->edpar,1);
2087         /* Read the whole edi file at once: */
2088         read_edi_file(ed,ed->edpar,mtop->natoms,cr);
2089
2090         /* Initialization for every ED/flooding dataset. Flooding uses one edi dataset per
2091          * flooding vector, Essential dynamics can be applied to more than one structure
2092          * as well, but will be done in the order given in the edi file, so
2093          * expect different results for different order of edi file concatenation! */
2094         edi=ed->edpar;
2095         while(edi != NULL)
2096         {
2097             init_edi(mtop,ir,cr,ed,edi);
2098
2099             /* Init flooding parameters if needed */
2100             init_flood(edi,ed,ir->delta_t,cr);
2101
2102             edi=edi->next_edi;
2103             numedis++;
2104         }
2105     }
2106
2107     /* The master does the work here. The other nodes get the positions
2108      * not before dd_partition_system which is called after init_edsam */
2109     if (MASTER(cr))
2110     {
2111         /* Remove pbc, make molecule whole.
2112          * When ir->bContinuation=TRUE this has already been done, but ok.
2113          */
2114         snew(x_pbc,mtop->natoms);
2115         m_rveccopy(mtop->natoms,x,x_pbc);
2116         do_pbc_first_mtop(NULL,ir->ePBC,box,mtop,x_pbc);
2117
2118         /* Reset pointer to first ED data set which contains the actual ED data */
2119         edi=ed->edpar;
2120
2121         /* Loop over all ED/flooding data sets (usually only one, though) */
2122         for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2123         {
2124             /* We use srenew to allocate memory since the size of the buffers
2125              * is likely to change with every ED dataset */
2126             srenew(xfit  , edi->sref.nr );
2127             srenew(xstart, edi->sav.nr  );
2128
2129             /* Extract the positions of the atoms to which will be fitted */
2130             for (i=0; i < edi->sref.nr; i++)
2131             {
2132                 copy_rvec(x_pbc[edi->sref.anrs[i]], xfit[i]);
2133
2134                 /* Save the sref positions such that in the next time step the molecule can
2135                  * be made whole again (in the parallel case) */
2136                 if (PAR(cr))
2137                     copy_rvec(xfit[i], edi->sref.x_old[i]);
2138             }
2139
2140             /* Extract the positions of the atoms subject to ED sampling */
2141             for (i=0; i < edi->sav.nr; i++)
2142             {
2143                 copy_rvec(x_pbc[edi->sav.anrs[i]], xstart[i]);
2144
2145                 /* Save the sav positions such that in the next time step the molecule can
2146                  * be made whole again (in the parallel case) */
2147                 if (PAR(cr))
2148                     copy_rvec(xstart[i], edi->sav.x_old[i]);
2149             }
2150
2151             /* Make the fit to the REFERENCE structure, get translation and rotation */
2152             fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2153
2154             /* Output how well we fit to the reference at the start */
2155             translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2156             fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm (dataset #%d)\n",
2157                     rmsd_from_structure(xfit, &edi->sref), nr_edi);
2158
2159             /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2160             translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2161
2162             /* calculate initial projections */
2163             project(xstart, edi);
2164
2165             /* process target structure, if required */
2166             if (edi->star.nr > 0)
2167             {
2168                 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
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 (eEDflood == ed->eEDtype)
2179                 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2180             if (edi->sori.nr > 0)
2181             {
2182                 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2183                 /* fit this structure to reference structure */
2184                 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2185                 /* do the fit */
2186                 translate_and_rotate(edi->sori.x, edi->sav.nr, fit_transvec, fit_rotmat);
2187                 rad_project(edi, edi->sori.x, &edi->vecs.radacc, cr);
2188                 rad_project(edi, edi->sori.x, &edi->vecs.radfix, cr);
2189                 if (ed->eEDtype == eEDflood)
2190                 {
2191                     fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2192                     /* Set center of flooding potential to the ORIGIN structure */
2193                     rad_project(edi, edi->sori.x, &edi->flood.vecs, cr);
2194                 }
2195             }
2196             else /* No origin structure given */
2197             {
2198                 rad_project(edi, xstart, &edi->vecs.radacc, cr);
2199                 rad_project(edi, xstart, &edi->vecs.radfix, cr);
2200                 if (ed->eEDtype == eEDflood)
2201                 {
2202                     if (edi->flood.bHarmonic)
2203                     {
2204                         fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2205                         for (i=0; i<edi->flood.vecs.neig; i++)
2206                             edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2207                     }
2208                     else
2209                     {
2210                         fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2211                         /* Set center of flooding potential to the center of the covariance matrix,
2212                          * i.e. the average structure, i.e. zero in the projected system */
2213                         for (i=0; i<edi->flood.vecs.neig; i++)
2214                             edi->flood.vecs.refproj[i] = 0.0;
2215                     }
2216                 }
2217             }
2218             /* For convenience, output the center of the flooding potential for the eigenvectors */
2219             if (eEDflood == ed->eEDtype)
2220             {
2221                 for (i=0; i<edi->flood.vecs.neig; i++)
2222                 {
2223                     fprintf(stdout, "ED: EV %d flooding potential center: %11.4e (adding %11.4e/timestep)\n",
2224                             i, edi->flood.vecs.refproj[i], edi->flood.vecs.refprojslope[i]);
2225                 }
2226             }
2227
2228             /* set starting projections for linsam */
2229             rad_project(edi, xstart, &edi->vecs.linacc, cr);
2230             rad_project(edi, xstart, &edi->vecs.linfix, cr);
2231
2232             /* Output to file, set the step to -1 so that write_edo knows it was called from init_edsam */
2233             if (ed->edo && !(ed->bStartFromCpt))
2234                 write_edo(nr_edi, edi, ed, -1, 0);
2235
2236             /* Prepare for the next edi data set: */
2237             edi=edi->next_edi;
2238         }
2239         /* Cleaning up on the master node: */
2240         sfree(x_pbc);
2241         sfree(xfit);
2242         sfree(xstart);
2243
2244     } /* end of MASTER only section */
2245
2246     if (PAR(cr))
2247     {
2248         /* First let everybody know how many ED data sets to expect */
2249         gmx_bcast(sizeof(numedis), &numedis, cr);
2250         /* Broadcast the essential dynamics / flooding data to all nodes */
2251         broadcast_ed_data(cr, ed, numedis);
2252     }
2253     else
2254     {
2255         /* In the single-CPU case, point the local atom numbers pointers to the global
2256          * one, so that we can use the same notation in serial and parallel case: */
2257
2258         /* Loop over all ED data sets (usually only one, though) */
2259         edi=ed->edpar;
2260         for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2261         {
2262             edi->sref.anrs_loc = edi->sref.anrs;
2263             edi->sav.anrs_loc  = edi->sav.anrs;
2264             edi->star.anrs_loc = edi->star.anrs;
2265             edi->sori.anrs_loc = edi->sori.anrs;
2266             /* For the same reason as above, make a dummy c_ind array: */
2267             snew(edi->sav.c_ind, edi->sav.nr);
2268             /* Initialize the array */
2269             for (i=0; i<edi->sav.nr; i++)
2270                 edi->sav.c_ind[i] = i;
2271             /* In the general case we will need a different-sized array for the reference indices: */
2272             if (!edi->bRefEqAv)
2273             {
2274                 snew(edi->sref.c_ind, edi->sref.nr);
2275                 for (i=0; i<edi->sref.nr; i++)
2276                     edi->sref.c_ind[i] = i;
2277             }
2278             /* Point to the very same array in case of other structures: */
2279             edi->star.c_ind = edi->sav.c_ind;
2280             edi->sori.c_ind = edi->sav.c_ind;
2281             /* In the serial case, the local number of atoms is the global one: */
2282             edi->sref.nr_loc = edi->sref.nr;
2283             edi->sav.nr_loc  = edi->sav.nr;
2284             edi->star.nr_loc = edi->star.nr;
2285             edi->sori.nr_loc = edi->sori.nr;
2286
2287             /* An on we go to the next edi dataset */
2288             edi=edi->next_edi;
2289         }
2290     }
2291
2292     /* Allocate space for ED buffer variables */
2293     /* Again, loop over ED data sets */
2294     edi=ed->edpar;
2295     for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2296     {
2297         /* Allocate space for ED buffer */
2298         snew(edi->buf, 1);
2299         snew(edi->buf->do_edsam, 1);
2300
2301         /* Space for collective ED buffer variables */
2302
2303         /* Collective positions of atoms with the average indices */
2304         snew(edi->buf->do_edsam->xcoll                  , edi->sav.nr);
2305         snew(edi->buf->do_edsam->shifts_xcoll           , edi->sav.nr); /* buffer for xcoll shifts */
2306         snew(edi->buf->do_edsam->extra_shifts_xcoll     , edi->sav.nr);
2307         /* Collective positions of atoms with the reference indices */
2308         if (!edi->bRefEqAv)
2309         {
2310             snew(edi->buf->do_edsam->xc_ref             , edi->sref.nr);
2311             snew(edi->buf->do_edsam->shifts_xc_ref      , edi->sref.nr); /* To store the shifts in */
2312             snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2313         }
2314
2315         /* Get memory for flooding forces */
2316         snew(edi->flood.forces_cartesian                , edi->sav.nr);
2317
2318 #ifdef DUMPEDI
2319         /* Dump it all into one file per process */
2320         dump_edi(edi, cr, nr_edi);
2321 #endif
2322
2323         /* An on we go to the next edi dataset */
2324         edi=edi->next_edi;
2325     }
2326
2327     /* Flush the edo file so that the user can check some things
2328      * when the simulation has started */
2329     if (ed->edo)
2330         fflush(ed->edo);
2331
2332     GMX_MPE_LOG(ev_edsam_finish);
2333 }
2334
2335
2336 void do_edsam(t_inputrec  *ir,
2337               gmx_large_int_t step,
2338               t_mdatoms   *md,
2339               t_commrec   *cr,
2340               rvec        xs[],   /* The local current positions on this processor */
2341               rvec        v[],    /* The velocities */
2342               matrix      box,
2343               gmx_edsam_t ed)
2344 {
2345     int     i,edinr,iupdate=500;
2346     matrix  rotmat;         /* rotation matrix */
2347     rvec    transvec;       /* translation vector */
2348     rvec    dv,dx,x_unsh;   /* tmp vectors for velocity, distance, unshifted x coordinate */
2349     real    dt_1;           /* 1/dt */
2350     struct t_do_edsam *buf;
2351     t_edpar *edi;
2352     real    rmsdev=-1;      /* RMSD from reference structure prior to applying the constraints */
2353     gmx_bool bSuppress=FALSE; /* Write .edo file on master? */
2354
2355
2356     /* Check if ED sampling has to be performed */
2357     if ( ed->eEDtype==eEDnone )
2358         return;
2359
2360     /* Suppress output on first call of do_edsam if
2361      * two-step sd2 integrator is used */
2362     if ( (ir->eI==eiSD2) && (v != NULL) )
2363         bSuppress = TRUE;
2364
2365     dt_1 = 1.0/ir->delta_t;
2366
2367     /* Loop over all ED datasets (usually one) */
2368     edi  = ed->edpar;
2369     edinr = 0;
2370     while (edi != NULL)
2371     {
2372         edinr++;
2373         if (edi->bNeedDoEdsam)
2374         {
2375
2376             buf=edi->buf->do_edsam;
2377
2378             if (ed->bFirst)
2379                 /* initialise radacc radius for slope criterion */
2380                 buf->oldrad=calc_radius(&edi->vecs.radacc);
2381
2382             /* Copy the positions into buf->xc* arrays and after ED
2383              * feed back corrections to the official positions */
2384
2385             /* Broadcast the ED positions such that every node has all of them
2386              * Every node contributes its local positions xs and stores it in
2387              * the collective buf->xcoll array. Note that for edinr > 1
2388              * xs could already have been modified by an earlier ED */
2389
2390             communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, buf->bUpdateShifts, xs,
2391                     edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old,  box);
2392
2393 #ifdef DEBUG_ED
2394             dump_xcoll(edi, buf, cr, step);
2395 #endif
2396             /* Only assembly reference positions if their indices differ from the average ones */
2397             if (!edi->bRefEqAv)
2398                 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, buf->bUpdateShifts, xs,
2399                         edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
2400
2401             /* If bUpdateShifts was TRUE then the shifts have just been updated in get_positions.
2402              * We do not need to uptdate the shifts until the next NS step */
2403             buf->bUpdateShifts = FALSE;
2404
2405             /* Now all nodes have all of the ED positions in edi->sav->xcoll,
2406              * as well as the indices in edi->sav.anrs */
2407
2408             /* Fit the reference indices to the reference structure */
2409             if (edi->bRefEqAv)
2410                 fit_to_reference(buf->xcoll , transvec, rotmat, edi);
2411             else
2412                 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
2413
2414             /* Now apply the translation and rotation to the ED structure */
2415             translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
2416
2417             /* Find out how well we fit to the reference (just for output steps) */
2418             if (do_per_step(step,edi->outfrq) && MASTER(cr))
2419             {
2420                 if (edi->bRefEqAv)
2421                 {
2422                     /* Indices of reference and average structures are identical,
2423                      * thus we can calculate the rmsd to SREF using xcoll */
2424                     rmsdev = rmsd_from_structure(buf->xcoll,&edi->sref);
2425                 }
2426                 else
2427                 {
2428                     /* We have to translate & rotate the reference atoms first */
2429                     translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
2430                     rmsdev = rmsd_from_structure(buf->xc_ref,&edi->sref);
2431                 }
2432             }
2433
2434             /* update radsam references, when required */
2435             if (do_per_step(step,edi->maxedsteps) && step >= edi->presteps)
2436             {
2437                 project(buf->xcoll, edi);
2438                 rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
2439                 rad_project(edi, buf->xcoll, &edi->vecs.radfix, cr);
2440                 buf->oldrad=-1.e5;
2441             }
2442
2443             /* update radacc references, when required */
2444             if (do_per_step(step,iupdate) && step >= edi->presteps)
2445             {
2446                 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
2447                 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
2448                 {
2449                     project(buf->xcoll, edi);
2450                     rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
2451                     buf->oldrad = 0.0;
2452                 } else
2453                     buf->oldrad = edi->vecs.radacc.radius;
2454             }
2455
2456             /* apply the constraints */
2457             if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
2458             {
2459                 /* ED constraints should be applied already in the first MD step
2460                  * (which is step 0), therefore we pass step+1 to the routine */
2461                 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step, cr);
2462             }
2463
2464             /* write to edo, when required */
2465             if (do_per_step(step,edi->outfrq))
2466             {
2467                 project(buf->xcoll, edi);
2468                 if (MASTER(cr) && !bSuppress)
2469                     write_edo(edinr, edi, ed, step, rmsdev);
2470             }
2471
2472             /* Copy back the positions unless monitoring only */
2473             if (ed_constraints(ed->eEDtype, edi))
2474             {
2475                 /* remove fitting */
2476                 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
2477
2478                 /* Copy the ED corrected positions into the coordinate array */
2479                 /* Each node copies its local part. In the serial case, nat_loc is the
2480                  * total number of ED atoms */
2481                 for (i=0; i<edi->sav.nr_loc; i++)
2482                 {
2483                     /* Unshift local ED coordinate and store in x_unsh */
2484                     ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
2485                                             buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
2486
2487                     /* dx is the ED correction to the positions: */
2488                     rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
2489
2490                     if (v != NULL)
2491                     {
2492                         /* dv is the ED correction to the velocity: */
2493                         svmul(dt_1, dx, dv);
2494                         /* apply the velocity correction: */
2495                         rvec_inc(v[edi->sav.anrs_loc[i]], dv);
2496                     }
2497                     /* Finally apply the position correction due to ED: */
2498                     copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
2499                 }
2500             }
2501         } /* END of if (edi->bNeedDoEdsam) */
2502
2503         /* Prepare for the next ED dataset */
2504         edi = edi->next_edi;
2505
2506     } /* END of loop over ED datasets */
2507
2508     ed->bFirst = FALSE;
2509
2510     GMX_MPE_LOG(ev_edsam_finish);
2511 }