src/gromacs/mdlib/edsam.c

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