src/mdlib/edsam.c

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