src/gromacs/mdlib/edsam.c

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