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         gmx_bool bNS)       /* Are we in a neighbor searching step? */
 861 {
 862     int i;
 863     matrix  rotmat;         /* rotation matrix */
 864     matrix  tmat;           /* inverse rotation */
 865     rvec    transvec;       /* translation vector */
 866     struct t_do_edsam *buf;
 867
 868
 869     buf=edi->buf->do_edsam;
 870
 871
 872     /* Broadcast the positions of the AVERAGE structure such that they are known on
 873      * every processor. Each node contributes its local positions x and stores them in
 874      * the collective ED array buf->xcoll */
 875     communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
 876                     edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
 877
 878     /* Only assembly REFERENCE positions if their indices differ from the average ones */
 879     if (!edi->bRefEqAv)
 880         communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
 881                 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
 882
 883     /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
 884      * We do not need to update the shifts until the next NS step */
 885     buf->bUpdateShifts = FALSE;
 886
 887     /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
 888      * as well as the indices in edi->sav.anrs */
 889
 890     /* Fit the reference indices to the reference structure */
 891     if (edi->bRefEqAv)
 892         fit_to_reference(buf->xcoll , transvec, rotmat, edi);
 893     else
 894         fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
 895
 896     /* Now apply the translation and rotation to the ED structure */
 897     translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
 898
 899     /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
 900     project_to_eigvectors(buf->xcoll,&edi->flood.vecs,edi);
 901
 902     if (FALSE == edi->flood.bConstForce)
 903     {
 904         /* Compute Vfl(x) from flood.xproj */
 905         edi->flood.Vfl = flood_energy(edi, step);
 906
 907         update_adaption(edi);
 908
 909         /* Compute the flooding forces */
 910         flood_forces(edi);
 911     }
 912
 913     /* Translate them into cartesian positions */
 914     flood_blowup(edi, edi->flood.forces_cartesian);
 915
 916     /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
 917     /* Each node rotates back its local forces */
 918     transpose(rotmat,tmat);
 919     rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
 920
 921     /* Finally add forces to the main force variable */
 922     for (i=0; i<edi->sav.nr_loc; i++)
 923         rvec_inc(force[edi->sav.anrs_loc[i]],edi->flood.forces_cartesian[i]);
 924
 925     /* Output is written by the master process */
 926     if (do_per_step(step,edi->outfrq) && MASTER(cr))
 927         write_edo_flood(edi,edo,step);
 928 }
 929
 930
 931 /* Main flooding routine, called from do_force */
 932 extern void do_flood(
 933         FILE            *log,    /* md.log file */
 934         t_commrec       *cr,     /* Communication record */
 935         rvec            x[],     /* Positions on the local processor */
 936         rvec            force[], /* forcefield forces, to these the flooding forces are added */
 937         gmx_edsam_t     ed,      /* ed data structure contains all ED and flooding datasets */
 938         matrix          box,     /* the box */
 939         gmx_large_int_t step,    /* The relative time step since ir->init_step is already subtracted */
 940         gmx_bool        bNS)     /* Are we in a neighbor searching step? */
 941 {
 942     t_edpar *edi;
 943
 944
 945     if (ed->eEDtype != eEDflood)
 946         return;
 947
 948     edi = ed->edpar;
 949     while (edi)
 950     {
 951         /* Call flooding for one matrix */
 952         if (edi->flood.vecs.neig)
 953             do_single_flood(ed->edo,x,force,edi,step,box,cr,bNS);
 954         edi = edi->next_edi;
 955     }
 956 }
 957
 958
 959 /* Called by init_edi, configure some flooding related variables and structures,
 960  * print headers to output files */
 961 static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt, t_commrec *cr)
 962 {
 963     int i;
 964
 965
 966     edi->flood.Efl = edi->flood.constEfl;
 967     edi->flood.Vfl = 0;
 968     edi->flood.dt  = dt;
 969
 970     if (edi->flood.vecs.neig)
 971     {
 972         /* If in any of the datasets we find a flooding vector, flooding is turned on */
 973         ed->eEDtype = eEDflood;
 974
 975         fprintf(stderr,"ED: Flooding of matrix %d is switched on.\n", edi->flood.flood_id);
 976
 977         if (edi->flood.bConstForce)
 978         {
 979             /* We have used stpsz as a vehicle to carry the fproj values for constant
 980              * force flooding. Now we copy that to flood.vecs.fproj. Note that
 981              * in const force flooding, fproj is never changed. */
 982             for (i=0; i<edi->flood.vecs.neig; i++)
 983             {
 984                 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
 985
 986                 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
 987                         edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
 988             }
 989         }
 990         fprintf(ed->edo,"FL_HEADER: Flooding of matrix %d is switched on! The flooding output will have the following format:\n",
 991                 edi->flood.flood_id);
 992         fprintf(ed->edo,"FL_HEADER: Step     Efl          Vfl       deltaF\n");
 993     }
 994 }
 995
 996
 997 #ifdef DEBUGHELPERS
 998 /*********** Energy book keeping ******/
 999 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames)  /* get header of energies */
1000 {
1001     t_edpar *actual;
1002     int count;
1003     char buf[STRLEN];
1004     actual=edi;
1005     count = 1;
1006     while (actual)
1007     {
1008         srenew(names,count);
1009         sprintf(buf,"Vfl_%d",count);
1010         names[count-1]=strdup(buf);
1011         actual=actual->next_edi;
1012         count++;
1013     }
1014     *nnames=count-1;
1015 }
1016
1017
1018 static void get_flood_energies(t_edpar *edi, real Vfl[],int nnames)
1019 {
1020     /*fl has to be big enough to capture nnames-many entries*/
1021     t_edpar *actual;
1022     int count;
1023
1024
1025     actual=edi;
1026     count = 1;
1027     while (actual)
1028     {
1029         Vfl[count-1]=actual->flood.Vfl;
1030         actual=actual->next_edi;
1031         count++;
1032     }
1033     if (nnames!=count-1)
1034         gmx_fatal(FARGS,"Number of energies is not consistent with t_edi structure");
1035 }
1036 /************* END of FLOODING IMPLEMENTATION ****************************/
1037 #endif
1038
1039
1040 gmx_edsam_t ed_open(int nfile,const t_filenm fnm[],unsigned long Flags,t_commrec *cr)
1041 {
1042     gmx_edsam_t ed;
1043
1044
1045     /* Allocate space for the ED data structure */
1046     snew(ed, 1);
1047
1048     /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1049     ed->eEDtype = eEDedsam;
1050
1051     if (MASTER(cr))
1052     {
1053         /* Open .edi input file: */
1054         ed->edinam=ftp2fn(efEDI,nfile,fnm);
1055         /* The master opens the .edo output file */
1056         fprintf(stderr,"ED sampling will be performed!\n");
1057         ed->edonam = ftp2fn(efEDO,nfile,fnm);
1058         ed->edo    = gmx_fio_fopen(ed->edonam,(Flags & MD_APPENDFILES)? "a+" : "w+");
1059         ed->bStartFromCpt = Flags & MD_STARTFROMCPT;
1060     }
1061     return ed;
1062 }
1063
1064
1065 /* Broadcasts the structure data */
1066 static void bc_ed_positions(t_commrec *cr, struct gmx_edx *s, int stype)
1067 {
1068     snew_bc(cr, s->anrs, s->nr   );    /* Index numbers     */
1069     snew_bc(cr, s->x   , s->nr   );    /* Positions         */
1070     nblock_bc(cr, s->nr, s->anrs );
1071     nblock_bc(cr, s->nr, s->x    );
1072
1073     /* For the average & reference structures we need an array for the collective indices,
1074      * and we need to broadcast the masses as well */
1075     if (stype == eedAV || stype == eedREF)
1076     {
1077         /* We need these additional variables in the parallel case: */
1078         snew(s->c_ind    , s->nr   );   /* Collective indices */
1079         /* Local atom indices get assigned in dd_make_local_group_indices.
1080          * There, also memory is allocated */
1081         s->nalloc_loc = 0;              /* allocation size of s->anrs_loc */
1082         snew_bc(cr, s->x_old, s->nr);   /* To be able to always make the ED molecule whole, ...        */
1083         nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1084     }
1085
1086     /* broadcast masses for the reference structure (for mass-weighted fitting) */
1087     if (stype == eedREF)
1088     {
1089         snew_bc(cr, s->m, s->nr);
1090         nblock_bc(cr, s->nr, s->m);
1091     }
1092
1093     /* For the average structure we might need the masses for mass-weighting */
1094     if (stype == eedAV)
1095     {
1096         snew_bc(cr, s->sqrtm, s->nr);
1097         nblock_bc(cr, s->nr, s->sqrtm);
1098         snew_bc(cr, s->m, s->nr);
1099         nblock_bc(cr, s->nr, s->m);
1100     }
1101 }
1102
1103
1104 /* Broadcasts the eigenvector data */
1105 static void bc_ed_vecs(t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
1106 {
1107     int i;
1108
1109     snew_bc(cr, ev->ieig   , ev->neig);  /* index numbers of eigenvector  */
1110     snew_bc(cr, ev->stpsz  , ev->neig);  /* stepsizes per eigenvector     */
1111     snew_bc(cr, ev->xproj  , ev->neig);  /* instantaneous x projection    */
1112     snew_bc(cr, ev->fproj  , ev->neig);  /* instantaneous f projection    */
1113     snew_bc(cr, ev->refproj, ev->neig);  /* starting or target projection */
1114
1115     nblock_bc(cr, ev->neig, ev->ieig   );
1116     nblock_bc(cr, ev->neig, ev->stpsz  );
1117     nblock_bc(cr, ev->neig, ev->xproj  );
1118     nblock_bc(cr, ev->neig, ev->fproj  );
1119     nblock_bc(cr, ev->neig, ev->refproj);
1120
1121     snew_bc(cr, ev->vec, ev->neig);      /* Eigenvector components        */
1122     for (i=0; i<ev->neig; i++)
1123     {
1124         snew_bc(cr, ev->vec[i], length);
1125         nblock_bc(cr, length, ev->vec[i]);
1126     }
1127
1128     /* For harmonic restraints the reference projections can change with time */
1129     if (bHarmonic)
1130     {
1131         snew_bc(cr, ev->refproj0    , ev->neig);
1132         snew_bc(cr, ev->refprojslope, ev->neig);
1133         nblock_bc(cr, ev->neig, ev->refproj0    );
1134         nblock_bc(cr, ev->neig, ev->refprojslope);
1135     }
1136 }
1137
1138
1139 /* Broadcasts the ED / flooding data to other nodes
1140  * and allocates memory where needed */
1141 static void broadcast_ed_data(t_commrec *cr, gmx_edsam_t ed, int numedis)
1142 {
1143     int     nr;
1144     t_edpar *edi;
1145
1146
1147     /* Master lets the other nodes know if its ED only or also flooding */
1148     gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1149
1150     snew_bc(cr, ed->edpar,1);
1151     /* Now transfer the ED data set(s) */
1152     edi = ed->edpar;
1153     for (nr=0; nr<numedis; nr++)
1154     {
1155         /* Broadcast a single ED data set */
1156         block_bc(cr, *edi);
1157
1158         /* Broadcast positions */
1159         bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses)    */
1160         bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1161         bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions                                */
1162         bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions                                */
1163
1164         /* Broadcast eigenvectors */
1165         bc_ed_vecs(cr, &edi->vecs.mon   , edi->sav.nr, FALSE);
1166         bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
1167         bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
1168         bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
1169         bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
1170         bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
1171         /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1172         bc_ed_vecs(cr, &edi->flood.vecs,  edi->sav.nr, edi->flood.bHarmonic);
1173
1174         /* Set the pointer to the next ED dataset */
1175         if (edi->next_edi)
1176         {
1177           snew_bc(cr, edi->next_edi, 1);
1178           edi = edi->next_edi;
1179         }
1180     }
1181 }
1182
1183
1184 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1185 static void init_edi(gmx_mtop_t *mtop,t_inputrec *ir,
1186                      t_commrec *cr,gmx_edsam_t ed,t_edpar *edi)
1187 {
1188     int  i;
1189     real totalmass = 0.0;
1190     rvec com;
1191     t_atom *atom;
1192
1193     /* NOTE Init_edi is executed on the master process only
1194      * The initialized data sets are then transmitted to the
1195      * other nodes in broadcast_ed_data */
1196
1197     edi->bNeedDoEdsam = edi->vecs.mon.neig
1198                      || edi->vecs.linfix.neig
1199                      || edi->vecs.linacc.neig
1200                      || edi->vecs.radfix.neig
1201                      || edi->vecs.radacc.neig
1202                      || edi->vecs.radcon.neig;
1203
1204     /* evaluate masses (reference structure) */
1205     snew(edi->sref.m, edi->sref.nr);
1206     for (i = 0; i < edi->sref.nr; i++)
1207     {
1208         if (edi->fitmas)
1209         {
1210             gmx_mtop_atomnr_to_atom(mtop,edi->sref.anrs[i],&atom);
1211             edi->sref.m[i] = atom->m;
1212         }
1213         else
1214         {
1215             edi->sref.m[i] = 1.0;
1216         }
1217
1218         /* Check that every m > 0. Bad things will happen otherwise. */
1219         if (edi->sref.m[i] <= 0.0)
1220         {
1221             gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1222                              "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1223                              "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1224                              "atoms from the reference structure by creating a proper index group.\n",
1225                       i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1226         }
1227
1228         totalmass += edi->sref.m[i];
1229     }
1230     edi->sref.mtot = totalmass;
1231
1232     /* Masses m and sqrt(m) for the average structure. Note that m
1233      * is needed if forces have to be evaluated in do_edsam */
1234     snew(edi->sav.sqrtm, edi->sav.nr );
1235     snew(edi->sav.m    , edi->sav.nr );
1236     for (i = 0; i < edi->sav.nr; i++)
1237     {
1238         gmx_mtop_atomnr_to_atom(mtop,edi->sav.anrs[i],&atom);
1239         edi->sav.m[i] = atom->m;
1240         if (edi->pcamas)
1241         {
1242             edi->sav.sqrtm[i] = sqrt(atom->m);
1243         }
1244         else
1245         {
1246             edi->sav.sqrtm[i] = 1.0;
1247         }
1248
1249         /* Check that every m > 0. Bad things will happen otherwise. */
1250         if (edi->sav.sqrtm[i] <= 0.0)
1251         {
1252             gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1253                              "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1254                              "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1255                              "atoms from the average structure by creating a proper index group.\n",
1256                       i, edi->sav.anrs[i]+1, atom->m);
1257         }
1258     }
1259
1260     /* put reference structure in origin */
1261     get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1262     com[XX] = -com[XX];
1263     com[YY] = -com[YY];
1264     com[ZZ] = -com[ZZ];
1265     translate_x(edi->sref.x, edi->sref.nr, com);
1266
1267     /* Init ED buffer */
1268     snew(edi->buf, 1);
1269 }
1270
1271
1272 static void check(const char *line, const char *label)
1273 {
1274     if (!strstr(line,label))
1275         gmx_fatal(FARGS,"Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s",label,line);
1276 }
1277
1278
1279 static int read_checked_edint(FILE *file,const char *label)
1280 {
1281     char line[STRLEN+1];
1282     int idum;
1283
1284
1285     fgets2 (line,STRLEN,file);
1286     check(line,label);
1287     fgets2 (line,STRLEN,file);
1288     sscanf (line,"%d",&idum);
1289     return idum;
1290 }
1291
1292
1293 static int read_edint(FILE *file,gmx_bool *bEOF)
1294 {
1295     char line[STRLEN+1];
1296     int idum;
1297     char *eof;
1298
1299
1300     eof=fgets2 (line,STRLEN,file);
1301     if (eof==NULL)
1302     {
1303         *bEOF = TRUE;
1304         return -1;
1305     }
1306     eof=fgets2 (line,STRLEN,file);
1307     if (eof==NULL)
1308     {
1309         *bEOF = TRUE;
1310         return -1;
1311     }
1312     sscanf (line,"%d",&idum);
1313     *bEOF = FALSE;
1314     return idum;
1315 }
1316
1317
1318 static real read_checked_edreal(FILE *file,const char *label)
1319 {
1320     char line[STRLEN+1];
1321     double rdum;
1322
1323
1324     fgets2 (line,STRLEN,file);
1325     check(line,label);
1326     fgets2 (line,STRLEN,file);
1327     sscanf (line,"%lf",&rdum);
1328     return (real) rdum; /* always read as double and convert to single */
1329 }
1330
1331
1332 static void read_edx(FILE *file,int number,int *anrs,rvec *x)
1333 {
1334     int i,j;
1335     char line[STRLEN+1];
1336     double d[3];
1337
1338
1339     for(i=0; i<number; i++)
1340     {
1341         fgets2 (line,STRLEN,file);
1342         sscanf (line,"%d%lf%lf%lf",&anrs[i],&d[0],&d[1],&d[2]);
1343         anrs[i]--; /* we are reading FORTRAN indices */
1344         for(j=0; j<3; j++)
1345             x[i][j]=d[j]; /* always read as double and convert to single */
1346     }
1347 }
1348
1349
1350 static void scan_edvec(FILE *in,int nr,rvec *vec)
1351 {
1352     char line[STRLEN+1];
1353     int i;
1354     double x,y,z;
1355
1356
1357     for(i=0; (i < nr); i++)
1358     {
1359         fgets2 (line,STRLEN,in);
1360         sscanf (line,"%le%le%le",&x,&y,&z);
1361         vec[i][XX]=x;
1362         vec[i][YY]=y;
1363         vec[i][ZZ]=z;
1364     }
1365 }
1366
1367
1368 static void read_edvec(FILE *in,int nr,t_eigvec *tvec,gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1369 {
1370     int i,idum,nscan;
1371     double rdum,refproj_dum=0.0,refprojslope_dum=0.0;
1372     char line[STRLEN+1];
1373
1374
1375     tvec->neig=read_checked_edint(in,"NUMBER OF EIGENVECTORS");
1376     if (tvec->neig >0)
1377     {
1378         snew(tvec->ieig   ,tvec->neig);
1379         snew(tvec->stpsz  ,tvec->neig);
1380         snew(tvec->vec    ,tvec->neig);
1381         snew(tvec->xproj  ,tvec->neig);
1382         snew(tvec->fproj  ,tvec->neig);
1383         snew(tvec->refproj,tvec->neig);
1384         if (bReadRefproj)
1385         {
1386             snew(tvec->refproj0    ,tvec->neig);
1387             snew(tvec->refprojslope,tvec->neig);
1388         }
1389
1390         for(i=0; (i < tvec->neig); i++)
1391         {
1392             fgets2 (line,STRLEN,in);
1393             if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1394             {
1395                 nscan = sscanf(line,"%d%lf%lf%lf",&idum,&rdum,&refproj_dum,&refprojslope_dum);
1396                 /* Zero out values which were not scanned */
1397                 switch(nscan)
1398                 {
1399                     case 4:
1400                         /* Every 4 values read, including reference position */
1401                         *bHaveReference = TRUE;
1402                         break;
1403                     case 3:
1404                         /* A reference position is provided */
1405                         *bHaveReference = TRUE;
1406                         /* No value for slope, set to 0 */
1407                         refprojslope_dum = 0.0;
1408                         break;
1409                     case 2:
1410                         /* No values for reference projection and slope, set to 0 */
1411                         refproj_dum      = 0.0;
1412                         refprojslope_dum = 0.0;
1413                         break;
1414                     default:
1415                         gmx_fatal(FARGS,"Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1416                         break;
1417                 }
1418                 tvec->refproj[i]=refproj_dum;
1419                 tvec->refproj0[i]=refproj_dum;
1420                 tvec->refprojslope[i]=refprojslope_dum;
1421             }
1422             else /* Normal flooding */
1423             {
1424                 nscan = sscanf(line,"%d%lf",&idum,&rdum);
1425                 if (nscan != 2)
1426                     gmx_fatal(FARGS,"Expected 2 values for flooding vec: <nr> <stpsz>\n");
1427             }
1428             tvec->ieig[i]=idum;
1429             tvec->stpsz[i]=rdum;
1430         } /* end of loop over eigenvectors */
1431
1432         for(i=0; (i < tvec->neig); i++)
1433         {
1434             snew(tvec->vec[i],nr);
1435             scan_edvec(in,nr,tvec->vec[i]);
1436         }
1437     }
1438 }
1439
1440
1441 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1442 static void read_edvecs(FILE *in,int nr,t_edvecs *vecs)
1443 {
1444         gmx_bool bHaveReference = FALSE;
1445
1446
1447     read_edvec(in, nr, &vecs->mon   , FALSE, &bHaveReference);
1448     read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1449     read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1450     read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1451     read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1452     read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1453 }
1454
1455
1456 /* Check if the same atom indices are used for reference and average positions */
1457 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1458 {
1459     int i;
1460
1461
1462     /* If the number of atoms differs between the two structures,
1463      * they cannot be identical */
1464     if (sref.nr != sav.nr)
1465         return FALSE;
1466
1467     /* Now that we know that both stuctures have the same number of atoms,
1468      * check if also the indices are identical */
1469     for (i=0; i < sav.nr; i++)
1470     {
1471         if (sref.anrs[i] != sav.anrs[i])
1472             return FALSE;
1473     }
1474     fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1475
1476     return TRUE;
1477 }
1478
1479
1480 static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int edi_nr, t_commrec *cr)
1481 {
1482     int readmagic;
1483     const int magic=670;
1484     gmx_bool bEOF;
1485
1486     /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1487     gmx_bool bHaveReference = FALSE;
1488
1489
1490     /* the edi file is not free format, so expect problems if the input is corrupt. */
1491
1492     /* check the magic number */
1493     readmagic=read_edint(in,&bEOF);
1494     /* Check whether we have reached the end of the input file */
1495     if (bEOF)
1496         return 0;
1497
1498     if (readmagic != magic)
1499     {
1500         if (readmagic==666 || readmagic==667 || readmagic==668)
1501             gmx_fatal(FARGS,"Wrong magic number: Use newest version of make_edi to produce edi file");
1502         else if (readmagic != 669)
1503             gmx_fatal(FARGS,"Wrong magic number %d in %s",readmagic,ed->edinam);
1504     }
1505
1506     /* check the number of atoms */
1507     edi->nini=read_edint(in,&bEOF);
1508     if (edi->nini != nr_mdatoms)
1509         gmx_fatal(FARGS,"Nr of atoms in %s (%d) does not match nr of md atoms (%d)",
1510                 ed->edinam,edi->nini,nr_mdatoms);
1511
1512     /* Done checking. For the rest we blindly trust the input */
1513     edi->fitmas          = read_checked_edint(in,"FITMAS");
1514     edi->pcamas          = read_checked_edint(in,"ANALYSIS_MAS");
1515     edi->outfrq          = read_checked_edint(in,"OUTFRQ");
1516     edi->maxedsteps      = read_checked_edint(in,"MAXLEN");
1517     edi->slope           = read_checked_edreal(in,"SLOPECRIT");
1518
1519     edi->presteps        = read_checked_edint(in,"PRESTEPS");
1520     edi->flood.deltaF0   = read_checked_edreal(in,"DELTA_F0");
1521     edi->flood.deltaF    = read_checked_edreal(in,"INIT_DELTA_F");
1522     edi->flood.tau       = read_checked_edreal(in,"TAU");
1523     edi->flood.constEfl  = read_checked_edreal(in,"EFL_NULL");
1524     edi->flood.alpha2    = read_checked_edreal(in,"ALPHA2");
1525     edi->flood.kT        = read_checked_edreal(in,"KT");
1526     edi->flood.bHarmonic = read_checked_edint(in,"HARMONIC");
1527     if (readmagic > 669)
1528         edi->flood.bConstForce = read_checked_edint(in,"CONST_FORCE_FLOODING");
1529     else
1530         edi->flood.bConstForce = FALSE;
1531     edi->flood.flood_id  = edi_nr;
1532     edi->sref.nr         = read_checked_edint(in,"NREF");
1533
1534     /* allocate space for reference positions and read them */
1535     snew(edi->sref.anrs,edi->sref.nr);
1536     snew(edi->sref.x   ,edi->sref.nr);
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     snew(edi->sav.x_old,edi->sav.nr);
1546     read_edx(in,edi->sav.nr,edi->sav.anrs,edi->sav.x);
1547
1548     /* Check if the same atom indices are used for reference and average positions */
1549     edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1550
1551     /* eigenvectors */
1552     read_edvecs(in,edi->sav.nr,&edi->vecs);
1553     read_edvec(in,edi->sav.nr,&edi->flood.vecs,edi->flood.bHarmonic, &bHaveReference);
1554
1555     /* target positions */
1556     edi->star.nr=read_edint(in,&bEOF);
1557     if (edi->star.nr > 0)
1558     {
1559         snew(edi->star.anrs,edi->star.nr);
1560         snew(edi->star.x   ,edi->star.nr);
1561         edi->star.sqrtm    =NULL;
1562         read_edx(in,edi->star.nr,edi->star.anrs,edi->star.x);
1563     }
1564
1565     /* positions defining origin of expansion circle */
1566     edi->sori.nr=read_edint(in,&bEOF);
1567     if (edi->sori.nr > 0)
1568     {
1569         if (bHaveReference)
1570         {
1571                 /* Both an -ori structure and a at least one manual reference point have been
1572                  * specified. That's ambiguous and probably not intentional. */
1573                 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1574                                  "    point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1575         }
1576         snew(edi->sori.anrs,edi->sori.nr);
1577         snew(edi->sori.x   ,edi->sori.nr);
1578         edi->sori.sqrtm    =NULL;
1579         read_edx(in,edi->sori.nr,edi->sori.anrs,edi->sori.x);
1580     }
1581
1582     /* all done */
1583     return 1;
1584 }
1585
1586
1587
1588 /* Read in the edi input file. Note that it may contain several ED data sets which were
1589  * achieved by concatenating multiple edi files. The standard case would be a single ED
1590  * data set, though. */
1591 static void read_edi_file(gmx_edsam_t ed, t_edpar *edi, int nr_mdatoms, t_commrec *cr)
1592 {
1593     FILE    *in;
1594     t_edpar *curr_edi,*last_edi;
1595     t_edpar *edi_read;
1596     int     edi_nr = 0;
1597
1598
1599     /* This routine is executed on the master only */
1600
1601     /* Open the .edi parameter input file */
1602     in = gmx_fio_fopen(ed->edinam,"r");
1603     fprintf(stderr, "ED: Reading edi file %s\n", ed->edinam);
1604
1605     /* Now read a sequence of ED input parameter sets from the edi file */
1606     curr_edi=edi;
1607     last_edi=edi;
1608     while( read_edi(in, ed, curr_edi, nr_mdatoms, edi_nr, cr) )
1609     {
1610         edi_nr++;
1611         /* Make shure that the number of atoms in each dataset is the same as in the tpr file */
1612         if (edi->nini != nr_mdatoms)
1613             gmx_fatal(FARGS,"edi file %s (dataset #%d) was made for %d atoms, but the simulation contains %d atoms.",
1614                     ed->edinam, edi_nr, edi->nini, nr_mdatoms);
1615         /* Since we arrived within this while loop we know that there is still another data set to be read in */
1616         /* We need to allocate space for the data: */
1617         snew(edi_read,1);
1618         /* Point the 'next_edi' entry to the next edi: */
1619         curr_edi->next_edi=edi_read;
1620         /* Keep the curr_edi pointer for the case that the next dataset is empty: */
1621         last_edi = curr_edi;
1622         /* Let's prepare to read in the next edi data set: */
1623         curr_edi = edi_read;
1624     }
1625     if (edi_nr == 0)
1626         gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", ed->edinam);
1627
1628     /* Terminate the edi dataset list with a NULL pointer: */
1629     last_edi->next_edi = NULL;
1630
1631     fprintf(stderr, "ED: Found %d ED dataset%s.\n", edi_nr, edi_nr>1? "s" : "");
1632
1633     /* Close the .edi file again */
1634     gmx_fio_fclose(in);
1635 }
1636
1637
1638 struct t_fit_to_ref {
1639     rvec *xcopy;       /* Working copy of the positions in fit_to_reference */
1640 };
1641
1642 /* Fit the current positions to the reference positions
1643  * Do not actually do the fit, just return rotation and translation.
1644  * Note that the COM of the reference structure was already put into
1645  * the origin by init_edi. */
1646 static void fit_to_reference(rvec      *xcoll,    /* The positions to be fitted */
1647                              rvec      transvec,  /* The translation vector */
1648                              matrix    rotmat,    /* The rotation matrix */
1649                              t_edpar   *edi)      /* Just needed for do_edfit */
1650 {
1651     rvec com;          /* center of mass */
1652     int  i;
1653     struct t_fit_to_ref *loc;
1654
1655
1656     /* Allocate memory the first time this routine is called for each edi dataset */
1657     if (NULL == edi->buf->fit_to_ref)
1658     {
1659         snew(edi->buf->fit_to_ref, 1);
1660         snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1661     }
1662     loc = edi->buf->fit_to_ref;
1663
1664     /* We do not touch the original positions but work on a copy. */
1665     for (i=0; i<edi->sref.nr; i++)
1666         copy_rvec(xcoll[i], loc->xcopy[i]);
1667
1668     /* Calculate the center of mass */
1669     get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1670
1671     transvec[XX] = -com[XX];
1672     transvec[YY] = -com[YY];
1673     transvec[ZZ] = -com[ZZ];
1674
1675     /* Subtract the center of mass from the copy */
1676     translate_x(loc->xcopy, edi->sref.nr, transvec);
1677
1678     /* Determine the rotation matrix */
1679     do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1680 }
1681
1682
1683 static void translate_and_rotate(rvec *x,         /* The positions to be translated and rotated */
1684                                  int nat,         /* How many positions are there? */
1685                                  rvec transvec,   /* The translation vector */
1686                                  matrix rotmat)   /* The rotation matrix */
1687 {
1688     /* Translation */
1689     translate_x(x, nat, transvec);
1690
1691     /* Rotation */
1692     rotate_x(x, nat, rotmat);
1693 }
1694
1695
1696 /* Gets the rms deviation of the positions to the structure s */
1697 /* fit_to_structure has to be called before calling this routine! */
1698 static real rmsd_from_structure(rvec           *x,  /* The positions under consideration */
1699                                 struct gmx_edx *s)  /* The structure from which the rmsd shall be computed */
1700 {
1701     real  rmsd=0.0;
1702     int   i;
1703
1704
1705     for (i=0; i < s->nr; i++)
1706         rmsd += distance2(s->x[i], x[i]);
1707
1708     rmsd /= (real) s->nr;
1709     rmsd = sqrt(rmsd);
1710
1711     return rmsd;
1712 }
1713
1714
1715 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1716 {
1717     t_edpar *edi;
1718
1719
1720     if (ed->eEDtype != eEDnone)
1721     {
1722         /* Loop over ED datasets (usually there is just one dataset, though) */
1723         edi=ed->edpar;
1724         while (edi)
1725         {
1726             /* Local atoms of the reference structure (for fitting), need only be assembled
1727              * if their indices differ from the average ones */
1728             if (!edi->bRefEqAv)
1729                 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1730                         &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1731
1732             /* Local atoms of the average structure (on these ED will be performed) */
1733             dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1734                     &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1735
1736             /* Indicate that the ED shift vectors for this structure need to be updated
1737              * at the next call to communicate_group_positions, since obviously we are in a NS step */
1738             edi->buf->do_edsam->bUpdateShifts = TRUE;
1739
1740             /* Set the pointer to the next ED dataset (if any) */
1741             edi=edi->next_edi;
1742         }
1743     }
1744 }
1745
1746
1747 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1748 {
1749     int tx,ty,tz;
1750
1751
1752     tx=is[XX];
1753     ty=is[YY];
1754     tz=is[ZZ];
1755
1756     if(TRICLINIC(box))
1757     {
1758         xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1759         xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1760         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1761     } else
1762     {
1763         xu[XX] = x[XX]-tx*box[XX][XX];
1764         xu[YY] = x[YY]-ty*box[YY][YY];
1765         xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1766     }
1767 }
1768
1769
1770 static void do_linfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
1771 {
1772     int  i, j;
1773     real proj, add;
1774     rvec vec_dum;
1775
1776
1777     /* loop over linfix vectors */
1778     for (i=0; i<edi->vecs.linfix.neig; i++)
1779     {
1780         /* calculate the projection */
1781         proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1782
1783         /* calculate the correction */
1784         add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1785
1786         /* apply the correction */
1787         add /= edi->sav.sqrtm[i];
1788         for (j=0; j<edi->sav.nr; j++)
1789         {
1790             svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1791             rvec_inc(xcoll[j], vec_dum);
1792         }
1793     }
1794 }
1795
1796
1797 static void do_linacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1798 {
1799     int  i, j;
1800     real proj, add;
1801     rvec vec_dum;
1802
1803
1804     /* loop over linacc vectors */
1805     for (i=0; i<edi->vecs.linacc.neig; i++)
1806     {
1807         /* calculate the projection */
1808         proj=projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1809
1810         /* calculate the correction */
1811         add = 0.0;
1812         if (edi->vecs.linacc.stpsz[i] > 0.0)
1813         {
1814             if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1815                 add = edi->vecs.linacc.refproj[i] - proj;
1816         }
1817         if (edi->vecs.linacc.stpsz[i] < 0.0)
1818         {
1819             if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1820                 add = edi->vecs.linacc.refproj[i] - proj;
1821         }
1822
1823         /* apply the correction */
1824         add /= edi->sav.sqrtm[i];
1825         for (j=0; j<edi->sav.nr; j++)
1826         {
1827             svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
1828             rvec_inc(xcoll[j], vec_dum);
1829         }
1830
1831         /* new positions will act as reference */
1832         edi->vecs.linacc.refproj[i] = proj + add;
1833     }
1834 }
1835
1836
1837 static void do_radfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
1838 {
1839     int  i,j;
1840     real *proj, rad=0.0, ratio;
1841     rvec vec_dum;
1842
1843
1844     if (edi->vecs.radfix.neig == 0)
1845         return;
1846
1847     snew(proj, edi->vecs.radfix.neig);
1848
1849     /* loop over radfix vectors */
1850     for (i=0; i<edi->vecs.radfix.neig; i++)
1851     {
1852         /* calculate the projections, radius */
1853         proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
1854         rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
1855     }
1856
1857     rad   = sqrt(rad);
1858     ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
1859     edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
1860
1861     /* loop over radfix vectors */
1862     for (i=0; i<edi->vecs.radfix.neig; i++)
1863     {
1864         proj[i] -= edi->vecs.radfix.refproj[i];
1865
1866         /* apply the correction */
1867         proj[i] /= edi->sav.sqrtm[i];
1868         proj[i] *= ratio;
1869         for (j=0; j<edi->sav.nr; j++) {
1870             svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
1871             rvec_inc(xcoll[j], vec_dum);
1872         }
1873     }
1874
1875     sfree(proj);
1876 }
1877
1878
1879 static void do_radacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1880 {
1881     int  i,j;
1882     real *proj, rad=0.0, ratio=0.0;
1883     rvec vec_dum;
1884
1885
1886     if (edi->vecs.radacc.neig == 0)
1887         return;
1888
1889     snew(proj,edi->vecs.radacc.neig);
1890
1891     /* loop over radacc vectors */
1892     for (i=0; i<edi->vecs.radacc.neig; i++)
1893     {
1894         /* calculate the projections, radius */
1895         proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
1896         rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
1897     }
1898     rad = sqrt(rad);
1899
1900     /* only correct when radius decreased */
1901     if (rad < edi->vecs.radacc.radius)
1902     {
1903         ratio = edi->vecs.radacc.radius/rad - 1.0;
1904         rad   = edi->vecs.radacc.radius;
1905     }
1906     else
1907         edi->vecs.radacc.radius = rad;
1908
1909     /* loop over radacc vectors */
1910     for (i=0; i<edi->vecs.radacc.neig; i++)
1911     {
1912         proj[i] -= edi->vecs.radacc.refproj[i];
1913
1914         /* apply the correction */
1915         proj[i] /= edi->sav.sqrtm[i];
1916         proj[i] *= ratio;
1917         for (j=0; j<edi->sav.nr; j++)
1918         {
1919             svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
1920             rvec_inc(xcoll[j], vec_dum);
1921         }
1922     }
1923     sfree(proj);
1924 }
1925
1926
1927 struct t_do_radcon {
1928     real *proj;
1929 };
1930
1931 static void do_radcon(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1932 {
1933     int  i,j;
1934     real rad=0.0, ratio=0.0;
1935     struct t_do_radcon *loc;
1936     gmx_bool bFirst;
1937     rvec vec_dum;
1938
1939
1940     if(edi->buf->do_radcon != NULL)
1941     {
1942         bFirst = FALSE;
1943         loc    = edi->buf->do_radcon;
1944     }
1945     else
1946     {
1947         bFirst = TRUE;
1948         snew(edi->buf->do_radcon, 1);
1949     }
1950     loc = edi->buf->do_radcon;
1951
1952     if (edi->vecs.radcon.neig == 0)
1953         return;
1954
1955     if (bFirst)
1956         snew(loc->proj, edi->vecs.radcon.neig);
1957
1958     /* loop over radcon vectors */
1959     for (i=0; i<edi->vecs.radcon.neig; i++)
1960     {
1961         /* calculate the projections, radius */
1962         loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
1963         rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
1964     }
1965     rad = sqrt(rad);
1966     /* only correct when radius increased */
1967     if (rad > edi->vecs.radcon.radius)
1968     {
1969         ratio = edi->vecs.radcon.radius/rad - 1.0;
1970
1971         /* loop over radcon vectors */
1972         for (i=0; i<edi->vecs.radcon.neig; i++)
1973         {
1974             /* apply the correction */
1975             loc->proj[i] -= edi->vecs.radcon.refproj[i];
1976             loc->proj[i] /= edi->sav.sqrtm[i];
1977             loc->proj[i] *= ratio;
1978
1979             for (j=0; j<edi->sav.nr; j++)
1980             {
1981                 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
1982                 rvec_inc(xcoll[j], vec_dum);
1983             }
1984         }
1985     }
1986     else
1987         edi->vecs.radcon.radius = rad;
1988
1989     if (rad != edi->vecs.radcon.radius)
1990     {
1991         rad = 0.0;
1992         for (i=0; i<edi->vecs.radcon.neig; i++)
1993         {
1994             /* calculate the projections, radius */
1995             loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
1996             rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
1997         }
1998         rad = sqrt(rad);
1999     }
2000 }
2001
2002
2003 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step, t_commrec *cr)
2004 {
2005     int i;
2006
2007
2008     /* subtract the average positions */
2009     for (i=0; i<edi->sav.nr; i++)
2010         rvec_dec(xcoll[i], edi->sav.x[i]);
2011
2012     /* apply the constraints */
2013     if (step >= 0)
2014         do_linfix(xcoll, edi, step, cr);
2015     do_linacc(xcoll, edi, cr);
2016     if (step >= 0)
2017         do_radfix(xcoll, edi, step, cr);
2018     do_radacc(xcoll, edi, cr);
2019     do_radcon(xcoll, edi, cr);
2020
2021     /* add back the average positions */
2022     for (i=0; i<edi->sav.nr; i++)
2023         rvec_inc(xcoll[i], edi->sav.x[i]);
2024 }
2025
2026
2027 /* Write out the projections onto the eigenvectors */
2028 static void write_edo(int nr_edi, t_edpar *edi, gmx_edsam_t ed, gmx_large_int_t step,real rmsd)
2029 {
2030     int i;
2031     char buf[22];
2032
2033
2034     if (edi->bNeedDoEdsam)
2035     {
2036         if (step == -1)
2037             fprintf(ed->edo, "Initial projections:\n");
2038         else
2039         {
2040             fprintf(ed->edo,"Step %s, ED #%d  ", gmx_step_str(step, buf), nr_edi);
2041             fprintf(ed->edo,"  RMSD %f nm\n",rmsd);
2042         }
2043
2044         if (edi->vecs.mon.neig)
2045         {
2046             fprintf(ed->edo,"  Monitor eigenvectors");
2047             for (i=0; i<edi->vecs.mon.neig; i++)
2048                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.mon.ieig[i],edi->vecs.mon.xproj[i]);
2049             fprintf(ed->edo,"\n");
2050         }
2051         if (edi->vecs.linfix.neig)
2052         {
2053             fprintf(ed->edo,"  Linfix  eigenvectors");
2054             for (i=0; i<edi->vecs.linfix.neig; i++)
2055                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.linfix.ieig[i],edi->vecs.linfix.xproj[i]);
2056             fprintf(ed->edo,"\n");
2057         }
2058         if (edi->vecs.linacc.neig)
2059         {
2060             fprintf(ed->edo,"  Linacc  eigenvectors");
2061             for (i=0; i<edi->vecs.linacc.neig; i++)
2062                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.linacc.ieig[i],edi->vecs.linacc.xproj[i]);
2063             fprintf(ed->edo,"\n");
2064         }
2065         if (edi->vecs.radfix.neig)
2066         {
2067             fprintf(ed->edo,"  Radfix  eigenvectors");
2068             for (i=0; i<edi->vecs.radfix.neig; i++)
2069                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radfix.ieig[i],edi->vecs.radfix.xproj[i]);
2070             fprintf(ed->edo,"\n");
2071             fprintf(ed->edo,"  fixed increment radius = %f\n", calc_radius(&edi->vecs.radfix));
2072         }
2073         if (edi->vecs.radacc.neig)
2074         {
2075             fprintf(ed->edo,"  Radacc  eigenvectors");
2076             for (i=0; i<edi->vecs.radacc.neig; i++)
2077                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radacc.ieig[i],edi->vecs.radacc.xproj[i]);
2078             fprintf(ed->edo,"\n");
2079             fprintf(ed->edo,"  acceptance radius      = %f\n", calc_radius(&edi->vecs.radacc));
2080         }
2081         if (edi->vecs.radcon.neig)
2082         {
2083             fprintf(ed->edo,"  Radcon  eigenvectors");
2084             for (i=0; i<edi->vecs.radcon.neig; i++)
2085                 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radcon.ieig[i],edi->vecs.radcon.xproj[i]);
2086             fprintf(ed->edo,"\n");
2087             fprintf(ed->edo,"  contracting radius     = %f\n", calc_radius(&edi->vecs.radcon));
2088         }
2089     }
2090 }
2091
2092 /* Returns if any constraints are switched on */
2093 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2094 {
2095     if (edtype == eEDedsam || edtype == eEDflood)
2096     {
2097         return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2098                 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2099                 edi->vecs.radcon.neig);
2100     }
2101     return 0;
2102 }
2103
2104
2105 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2106  * umbrella sampling simulations. */
2107 static void copyEvecReference(t_eigvec* floodvecs)
2108 {
2109     int i;
2110
2111
2112     if (NULL==floodvecs->refproj0)
2113         snew(floodvecs->refproj0, floodvecs->neig);
2114
2115     for (i=0; i<floodvecs->neig; i++)
2116     {
2117         floodvecs->refproj0[i] = floodvecs->refproj[i];
2118     }
2119 }
2120
2121
2122 void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
2123                 t_inputrec  *ir,     /* input record                       */
2124                 t_commrec   *cr,     /* communication record               */
2125                 gmx_edsam_t ed,      /* contains all ED data               */
2126                 rvec        x[],     /* positions of the whole MD system   */
2127                 matrix      box)     /* the box                            */
2128 {
2129     t_edpar *edi = NULL;    /* points to a single edi data set */
2130     int     numedis=0;      /* keep track of the number of ED data sets in edi file */
2131     int     i,nr_edi,avindex;
2132     rvec    *x_pbc  = NULL; /* positions of the whole MD system with pbc removed  */
2133     rvec    *xfit   = NULL; /* the positions which will be fitted to the reference structure  */
2134     rvec    *xstart = NULL; /* the positions which are subject to ED sampling */
2135     rvec    fit_transvec;   /* translation ... */
2136     matrix  fit_rotmat;     /* ... and rotation from fit to reference structure */
2137
2138
2139     if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
2140         gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
2141
2142     if (MASTER(cr))
2143         fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2144
2145     /* Needed for initializing radacc radius in do_edsam */
2146     ed->bFirst = 1;
2147
2148     /* The input file is read by the master and the edi structures are
2149      * initialized here. Input is stored in ed->edpar. Then the edi
2150      * structures are transferred to the other nodes */
2151     if (MASTER(cr))
2152     {
2153         snew(ed->edpar,1);
2154         /* Read the whole edi file at once: */
2155         read_edi_file(ed,ed->edpar,mtop->natoms,cr);
2156
2157         /* Initialization for every ED/flooding dataset. Flooding uses one edi dataset per
2158          * flooding vector, Essential dynamics can be applied to more than one structure
2159          * as well, but will be done in the order given in the edi file, so
2160          * expect different results for different order of edi file concatenation! */
2161         edi=ed->edpar;
2162         while(edi != NULL)
2163         {
2164             init_edi(mtop,ir,cr,ed,edi);
2165
2166             /* Init flooding parameters if needed */
2167             init_flood(edi,ed,ir->delta_t,cr);
2168
2169             edi=edi->next_edi;
2170             numedis++;
2171         }
2172     }
2173
2174     /* The master does the work here. The other nodes get the positions
2175      * not before dd_partition_system which is called after init_edsam */
2176     if (MASTER(cr))
2177     {
2178         /* Remove pbc, make molecule whole.
2179          * When ir->bContinuation=TRUE this has already been done, but ok.
2180          */
2181         snew(x_pbc,mtop->natoms);
2182         m_rveccopy(mtop->natoms,x,x_pbc);
2183         do_pbc_first_mtop(NULL,ir->ePBC,box,mtop,x_pbc);
2184
2185         /* Reset pointer to first ED data set which contains the actual ED data */
2186         edi=ed->edpar;
2187
2188         /* Loop over all ED/flooding data sets (usually only one, though) */
2189         for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2190         {
2191             /* We use srenew to allocate memory since the size of the buffers
2192              * is likely to change with every ED dataset */
2193             srenew(xfit  , edi->sref.nr );
2194             srenew(xstart, edi->sav.nr  );
2195
2196             /* Extract the positions of the atoms to which will be fitted */
2197             for (i=0; i < edi->sref.nr; i++)
2198             {
2199                 copy_rvec(x_pbc[edi->sref.anrs[i]], xfit[i]);
2200
2201                 /* Save the sref positions such that in the next time step we can make the ED group whole
2202                  * in case any of the atoms do not have the correct PBC representation */
2203                 copy_rvec(xfit[i], edi->sref.x_old[i]);
2204             }
2205
2206             /* Extract the positions of the atoms subject to ED sampling */
2207             for (i=0; i < edi->sav.nr; i++)
2208             {
2209                 copy_rvec(x_pbc[edi->sav.anrs[i]], xstart[i]);
2210
2211                 /* Save the sav positions such that in the next time step we can make the ED group whole
2212                  * in case any of the atoms do not have the correct PBC representation */
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             /* For the target and origin structure both a reference (fit) and an
2231              * average structure can be provided in make_edi. If both structures
2232              * are the same, make_edi only stores one of them in the .edi file.
2233              * If they differ, first the fit and then the average structure is stored
2234              * in star (or sor), thus the number of entries in star/sor is
2235              * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2236              * the size of the average group. */
2237
2238             /* process target structure, if required */
2239             if (edi->star.nr > 0)
2240             {
2241                 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2242
2243                 /* get translation & rotation for fit of target structure to reference structure */
2244                 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2245                 /* do the fit */
2246                 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2247                 if (edi->star.nr == edi->sav.nr)
2248                 {
2249                     avindex = 0;
2250                 }
2251                 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2252                 {
2253                     /* The last sav.nr indices of the target structure correspond to
2254                      * the average structure, which must be projected */
2255                     avindex = edi->star.nr - edi->sav.nr;
2256                 }
2257                 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon, cr);
2258             } else
2259                 rad_project(edi, xstart, &edi->vecs.radcon, cr);
2260
2261             /* process structure that will serve as origin of expansion circle */
2262             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2263                 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2264
2265             if (edi->sori.nr > 0)
2266             {
2267                 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2268
2269                 /* fit this structure to reference structure */
2270                 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2271                 /* do the fit */
2272                 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2273                 if (edi->sori.nr == edi->sav.nr)
2274                 {
2275                     avindex = 0;
2276                 }
2277                 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2278                 {
2279                     /* For the projection, we need the last sav.nr indices of sori */
2280                     avindex = edi->sori.nr - edi->sav.nr;
2281                 }
2282
2283                 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc, cr);
2284                 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix, cr);
2285                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2286                 {
2287                     fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2288                     /* Set center of flooding potential to the ORIGIN structure */
2289                     rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs, cr);
2290                     /* We already know that no (moving) reference position was provided,
2291                      * therefore we can overwrite refproj[0]*/
2292                     copyEvecReference(&edi->flood.vecs);
2293                 }
2294             }
2295             else /* No origin structure given */
2296             {
2297                 rad_project(edi, xstart, &edi->vecs.radacc, cr);
2298                 rad_project(edi, xstart, &edi->vecs.radfix, cr);
2299                 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2300                 {
2301                     if (edi->flood.bHarmonic)
2302                     {
2303                         fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2304                         for (i=0; i<edi->flood.vecs.neig; i++)
2305                             edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2306                     }
2307                     else
2308                     {
2309                         fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2310                         /* Set center of flooding potential to the center of the covariance matrix,
2311                          * i.e. the average structure, i.e. zero in the projected system */
2312                         for (i=0; i<edi->flood.vecs.neig; i++)
2313                             edi->flood.vecs.refproj[i] = 0.0;
2314                     }
2315                 }
2316             }
2317             /* For convenience, output the center of the flooding potential for the eigenvectors */
2318             if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2319             {
2320                 for (i=0; i<edi->flood.vecs.neig; i++)
2321                 {
2322                     fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", i, edi->flood.vecs.refproj[i]);
2323                     if (edi->flood.bHarmonic)
2324                         fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2325                     fprintf(stdout, "\n");
2326                 }
2327             }
2328
2329             /* set starting projections for linsam */
2330             rad_project(edi, xstart, &edi->vecs.linacc, cr);
2331             rad_project(edi, xstart, &edi->vecs.linfix, cr);
2332
2333             /* Output to file, set the step to -1 so that write_edo knows it was called from init_edsam */
2334             if (ed->edo && !(ed->bStartFromCpt))
2335                 write_edo(nr_edi, edi, ed, -1, 0);
2336
2337             /* Prepare for the next edi data set: */
2338             edi=edi->next_edi;
2339         }
2340         /* Cleaning up on the master node: */
2341         sfree(x_pbc);
2342         sfree(xfit);
2343         sfree(xstart);
2344
2345     } /* end of MASTER only section */
2346
2347     if (PAR(cr))
2348     {
2349         /* First let everybody know how many ED data sets to expect */
2350         gmx_bcast(sizeof(numedis), &numedis, cr);
2351         /* Broadcast the essential dynamics / flooding data to all nodes */
2352         broadcast_ed_data(cr, ed, numedis);
2353     }
2354     else
2355     {
2356         /* In the single-CPU case, point the local atom numbers pointers to the global
2357          * one, so that we can use the same notation in serial and parallel case: */
2358
2359         /* Loop over all ED data sets (usually only one, though) */
2360         edi=ed->edpar;
2361         for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2362         {
2363             edi->sref.anrs_loc = edi->sref.anrs;
2364             edi->sav.anrs_loc  = edi->sav.anrs;
2365             edi->star.anrs_loc = edi->star.anrs;
2366             edi->sori.anrs_loc = edi->sori.anrs;
2367             /* For the same reason as above, make a dummy c_ind array: */
2368             snew(edi->sav.c_ind, edi->sav.nr);
2369             /* Initialize the array */
2370             for (i=0; i<edi->sav.nr; i++)
2371                 edi->sav.c_ind[i] = i;
2372             /* In the general case we will need a different-sized array for the reference indices: */
2373             if (!edi->bRefEqAv)
2374             {
2375                 snew(edi->sref.c_ind, edi->sref.nr);
2376                 for (i=0; i<edi->sref.nr; i++)
2377                     edi->sref.c_ind[i] = i;
2378             }
2379             /* Point to the very same array in case of other structures: */
2380             edi->star.c_ind = edi->sav.c_ind;
2381             edi->sori.c_ind = edi->sav.c_ind;
2382             /* In the serial case, the local number of atoms is the global one: */
2383             edi->sref.nr_loc = edi->sref.nr;
2384             edi->sav.nr_loc  = edi->sav.nr;
2385             edi->star.nr_loc = edi->star.nr;
2386             edi->sori.nr_loc = edi->sori.nr;
2387
2388             /* An on we go to the next edi dataset */
2389             edi=edi->next_edi;
2390         }
2391     }
2392
2393     /* Allocate space for ED buffer variables */
2394     /* Again, loop over ED data sets */
2395     edi=ed->edpar;
2396     for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2397     {
2398         /* Allocate space for ED buffer */
2399         snew(edi->buf, 1);
2400         snew(edi->buf->do_edsam, 1);
2401
2402         /* Space for collective ED buffer variables */
2403
2404         /* Collective positions of atoms with the average indices */
2405         snew(edi->buf->do_edsam->xcoll                  , edi->sav.nr);
2406         snew(edi->buf->do_edsam->shifts_xcoll           , edi->sav.nr); /* buffer for xcoll shifts */
2407         snew(edi->buf->do_edsam->extra_shifts_xcoll     , edi->sav.nr);
2408         /* Collective positions of atoms with the reference indices */
2409         if (!edi->bRefEqAv)
2410         {
2411             snew(edi->buf->do_edsam->xc_ref             , edi->sref.nr);
2412             snew(edi->buf->do_edsam->shifts_xc_ref      , edi->sref.nr); /* To store the shifts in */
2413             snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2414         }
2415
2416         /* Get memory for flooding forces */
2417         snew(edi->flood.forces_cartesian                , edi->sav.nr);
2418
2419 #ifdef DUMPEDI
2420         /* Dump it all into one file per process */
2421         dump_edi(edi, cr, nr_edi);
2422 #endif
2423
2424         /* An on we go to the next edi dataset */
2425         edi=edi->next_edi;
2426     }
2427
2428     /* Flush the edo file so that the user can check some things
2429      * when the simulation has started */
2430     if (ed->edo)
2431         fflush(ed->edo);
2432 }
2433
2434
2435 void do_edsam(t_inputrec  *ir,
2436               gmx_large_int_t step,
2437               t_mdatoms   *md,
2438               t_commrec   *cr,
2439               rvec        xs[],   /* The local current positions on this processor */
2440               rvec        v[],    /* The velocities */
2441               matrix      box,
2442               gmx_edsam_t ed)
2443 {
2444     int     i,edinr,iupdate=500;
2445     matrix  rotmat;         /* rotation matrix */
2446     rvec    transvec;       /* translation vector */
2447     rvec    dv,dx,x_unsh;   /* tmp vectors for velocity, distance, unshifted x coordinate */
2448     real    dt_1;           /* 1/dt */
2449     struct t_do_edsam *buf;
2450     t_edpar *edi;
2451     real    rmsdev=-1;      /* RMSD from reference structure prior to applying the constraints */
2452     gmx_bool bSuppress=FALSE; /* Write .edo file on master? */
2453
2454
2455     /* Check if ED sampling has to be performed */
2456     if ( ed->eEDtype==eEDnone )
2457         return;
2458
2459     /* Suppress output on first call of do_edsam if
2460      * two-step sd2 integrator is used */
2461     if ( (ir->eI==eiSD2) && (v != NULL) )
2462         bSuppress = TRUE;
2463
2464     dt_1 = 1.0/ir->delta_t;
2465
2466     /* Loop over all ED datasets (usually one) */
2467     edi  = ed->edpar;
2468     edinr = 0;
2469     while (edi != NULL)
2470     {
2471         edinr++;
2472         if (edi->bNeedDoEdsam)
2473         {
2474
2475             buf=edi->buf->do_edsam;
2476
2477             if (ed->bFirst)
2478                 /* initialise radacc radius for slope criterion */
2479                 buf->oldrad=calc_radius(&edi->vecs.radacc);
2480
2481             /* Copy the positions into buf->xc* arrays and after ED
2482              * feed back corrections to the official positions */
2483
2484             /* Broadcast the ED positions such that every node has all of them
2485              * Every node contributes its local positions xs and stores it in
2486              * the collective buf->xcoll array. Note that for edinr > 1
2487              * xs could already have been modified by an earlier ED */
2488
2489             communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
2490                     edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old,  box);
2491
2492 #ifdef DEBUG_ED
2493             dump_xcoll(edi, buf, cr, step);
2494 #endif
2495             /* Only assembly reference positions if their indices differ from the average ones */
2496             if (!edi->bRefEqAv)
2497                 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
2498                         edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
2499
2500             /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
2501              * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
2502              * set bUpdateShifts=TRUE in the parallel case. */
2503             buf->bUpdateShifts = FALSE;
2504
2505             /* Now all nodes have all of the ED positions in edi->sav->xcoll,
2506              * as well as the indices in edi->sav.anrs */
2507
2508             /* Fit the reference indices to the reference structure */
2509             if (edi->bRefEqAv)
2510                 fit_to_reference(buf->xcoll , transvec, rotmat, edi);
2511             else
2512                 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
2513
2514             /* Now apply the translation and rotation to the ED structure */
2515             translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
2516
2517             /* Find out how well we fit to the reference (just for output steps) */
2518             if (do_per_step(step,edi->outfrq) && MASTER(cr))
2519             {
2520                 if (edi->bRefEqAv)
2521                 {
2522                     /* Indices of reference and average structures are identical,
2523                      * thus we can calculate the rmsd to SREF using xcoll */
2524                     rmsdev = rmsd_from_structure(buf->xcoll,&edi->sref);
2525                 }
2526                 else
2527                 {
2528                     /* We have to translate & rotate the reference atoms first */
2529                     translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
2530                     rmsdev = rmsd_from_structure(buf->xc_ref,&edi->sref);
2531                 }
2532             }
2533
2534             /* update radsam references, when required */
2535             if (do_per_step(step,edi->maxedsteps) && step >= edi->presteps)
2536             {
2537                 project(buf->xcoll, edi);
2538                 rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
2539                 rad_project(edi, buf->xcoll, &edi->vecs.radfix, cr);
2540                 buf->oldrad=-1.e5;
2541             }
2542
2543             /* update radacc references, when required */
2544             if (do_per_step(step,iupdate) && step >= edi->presteps)
2545             {
2546                 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
2547                 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
2548                 {
2549                     project(buf->xcoll, edi);
2550                     rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
2551                     buf->oldrad = 0.0;
2552                 } else
2553                     buf->oldrad = edi->vecs.radacc.radius;
2554             }
2555
2556             /* apply the constraints */
2557             if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
2558             {
2559                 /* ED constraints should be applied already in the first MD step
2560                  * (which is step 0), therefore we pass step+1 to the routine */
2561                 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step, cr);
2562             }
2563
2564             /* write to edo, when required */
2565             if (do_per_step(step,edi->outfrq))
2566             {
2567                 project(buf->xcoll, edi);
2568                 if (MASTER(cr) && !bSuppress)
2569                     write_edo(edinr, edi, ed, step, rmsdev);
2570             }
2571
2572             /* Copy back the positions unless monitoring only */
2573             if (ed_constraints(ed->eEDtype, edi))
2574             {
2575                 /* remove fitting */
2576                 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
2577
2578                 /* Copy the ED corrected positions into the coordinate array */
2579                 /* Each node copies its local part. In the serial case, nat_loc is the
2580                  * total number of ED atoms */
2581                 for (i=0; i<edi->sav.nr_loc; i++)
2582                 {
2583                     /* Unshift local ED coordinate and store in x_unsh */
2584                     ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
2585                                             buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
2586
2587                     /* dx is the ED correction to the positions: */
2588                     rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
2589
2590                     if (v != NULL)
2591                     {
2592                         /* dv is the ED correction to the velocity: */
2593                         svmul(dt_1, dx, dv);
2594                         /* apply the velocity correction: */
2595                         rvec_inc(v[edi->sav.anrs_loc[i]], dv);
2596                     }
2597                     /* Finally apply the position correction due to ED: */
2598                     copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
2599                 }
2600             }
2601         } /* END of if (edi->bNeedDoEdsam) */
2602
2603         /* Prepare for the next ED dataset */
2604         edi = edi->next_edi;
2605
2606     } /* END of loop over ED datasets */
2607
2608     ed->bFirst = FALSE;
2609 }