src/mdlib/edsam.c

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