src/mdlib/edsam.c

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