2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
59 #include "mtop_util.h"
61 #include "mpelogging.h"
63 #include "groupcoord.h"
66 /* We use the same defines as in mvdata.c here */
67 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d),(cr))
68 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
69 #define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
72 /* enum to identify the type of ED: none, normal ED, flooding */
73 enum {eEDnone, eEDedsam, eEDflood, eEDnr};
75 /* enum to identify operations on reference, average, origin, target structures */
76 enum {eedREF, eedAV, eedORI, eedTAR, eedNR};
81 int neig; /* nr of eigenvectors */
82 int *ieig; /* index nrs of eigenvectors */
83 real *stpsz; /* stepsizes (per eigenvector) */
84 rvec **vec; /* eigenvector components */
85 real *xproj; /* instantaneous x projections */
86 real *fproj; /* instantaneous f projections */
87 real radius; /* instantaneous radius */
88 real *refproj; /* starting or target projecions */
89 /* When using flooding as harmonic restraint: The current reference projection
90 * is at each step calculated from the initial refproj0 and the slope. */
91 real *refproj0,*refprojslope;
97 t_eigvec mon; /* only monitored, no constraints */
98 t_eigvec linfix; /* fixed linear constraints */
99 t_eigvec linacc; /* acceptance linear constraints */
100 t_eigvec radfix; /* fixed radial constraints (exp) */
101 t_eigvec radacc; /* acceptance radial constraints (exp) */
102 t_eigvec radcon; /* acceptance rad. contraction constr. */
109 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
111 gmx_bool bConstForce; /* Do not calculate a flooding potential,
112 instead flood with a constant force */
122 rvec *forces_cartesian;
123 t_eigvec vecs; /* use flooding for these */
127 /* This type is for the average, reference, target, and origin structure */
128 typedef struct gmx_edx
130 int nr; /* number of atoms this structure contains */
131 int nr_loc; /* number of atoms on local node */
132 int *anrs; /* atom index numbers */
133 int *anrs_loc; /* local atom index numbers */
134 int nalloc_loc; /* allocation size of anrs_loc */
135 int *c_ind; /* at which position of the whole anrs
136 * array is a local atom?, i.e.
137 * c_ind[0...nr_loc-1] gives the atom index
138 * with respect to the collective
139 * anrs[0...nr-1] array */
140 rvec *x; /* positions for this structure */
141 rvec *x_old; /* Last positions which have the correct PBC
142 representation of the ED group. In
143 combination with keeping track of the
144 shift vectors, the ED group can always
146 real *m; /* masses */
147 real mtot; /* total mass (only used in sref) */
148 real *sqrtm; /* sqrt of the masses used for mass-
149 * weighting of analysis (only used in sav) */
155 int nini; /* total Nr of atoms */
156 gmx_bool fitmas; /* true if trans fit with cm */
157 gmx_bool pcamas; /* true if mass-weighted PCA */
158 int presteps; /* number of steps to run without any
159 * perturbations ... just monitoring */
160 int outfrq; /* freq (in steps) of writing to edo */
161 int maxedsteps; /* max nr of steps per cycle */
163 /* all gmx_edx datasets are copied to all nodes in the parallel case */
164 struct gmx_edx sref; /* reference positions, to these fitting
166 gmx_bool bRefEqAv; /* If true, reference & average indices
167 * are the same. Used for optimization */
168 struct gmx_edx sav; /* average positions */
169 struct gmx_edx star; /* target positions */
170 struct gmx_edx sori; /* origin positions */
172 t_edvecs vecs; /* eigenvectors */
173 real slope; /* minimal slope in acceptance radexp */
175 gmx_bool bNeedDoEdsam; /* if any of the options mon, linfix, ...
176 * is used (i.e. apart from flooding) */
177 t_edflood flood; /* parameters especially for flooding */
178 struct t_ed_buffer *buf; /* handle to local buffers */
179 struct edpar *next_edi; /* Pointer to another ed dataset */
183 typedef struct gmx_edsam
185 int eEDtype; /* Type of ED: see enums above */
186 const char *edinam; /* name of ED sampling input file */
187 const char *edonam; /* output */
188 FILE *edo; /* output file pointer */
198 rvec old_transvec,older_transvec,transvec_compact;
199 rvec *xcoll; /* Positions from all nodes, this is the
200 collective set we work on.
201 These are the positions of atoms with
202 average structure indices */
203 rvec *xc_ref; /* same but with reference structure indices */
204 ivec *shifts_xcoll; /* Shifts for xcoll */
205 ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */
206 ivec *shifts_xc_ref; /* Shifts for xc_ref */
207 ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
208 gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
209 ED shifts for this ED dataset need to
214 /* definition of ED buffer structure */
217 struct t_fit_to_ref * fit_to_ref;
218 struct t_do_edfit * do_edfit;
219 struct t_do_edsam * do_edsam;
220 struct t_do_radcon * do_radcon;
224 /* Function declarations */
225 static void fit_to_reference(rvec *xcoll,rvec transvec,matrix rotmat,t_edpar *edi);
227 static void translate_and_rotate(rvec *x,int nat,rvec transvec,matrix rotmat);
228 /* End function declarations */
231 /* Does not subtract average positions, projection on single eigenvector is returned
232 * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
233 * Average position is subtracted in ed_apply_constraints prior to calling projectx
235 static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec)
241 for (i=0; i<edi->sav.nr; i++)
243 proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
250 /* Specialized: projection is stored in vec->refproj
251 * -> used for radacc, radfix, radcon and center of flooding potential
252 * subtracts average positions, projects vector x */
253 static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec, t_commrec *cr)
258 /* Subtract average positions */
259 for (i = 0; i < edi->sav.nr; i++)
261 rvec_dec(x[i], edi->sav.x[i]);
264 for (i = 0; i < vec->neig; i++)
266 vec->refproj[i] = projectx(edi,x,vec->vec[i]);
267 rad += pow((vec->refproj[i]-vec->xproj[i]),2);
269 vec->radius=sqrt(rad);
271 /* Add average positions */
272 for (i = 0; i < edi->sav.nr; i++)
274 rvec_inc(x[i], edi->sav.x[i]);
279 /* Project vector x, subtract average positions prior to projection and add
280 * them afterwards to retain the unchanged vector. Store in xproj. Mass-weighting
282 static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
283 t_eigvec *vec, /* The eigenvectors */
289 if (!vec->neig) return;
291 /* Subtract average positions */
292 for (i=0; i<edi->sav.nr; i++)
294 rvec_dec(x[i], edi->sav.x[i]);
297 for (i=0; i<vec->neig; i++)
299 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
302 /* Add average positions */
303 for (i=0; i<edi->sav.nr; i++)
305 rvec_inc(x[i], edi->sav.x[i]);
310 /* Project vector x onto all edi->vecs (mon, linfix,...) */
311 static void project(rvec *x, /* positions to project */
312 t_edpar *edi) /* edi data set */
314 /* It is not more work to subtract the average position in every
315 * subroutine again, because these routines are rarely used simultanely */
316 project_to_eigvectors(x, &edi->vecs.mon , edi);
317 project_to_eigvectors(x, &edi->vecs.linfix, edi);
318 project_to_eigvectors(x, &edi->vecs.linacc, edi);
319 project_to_eigvectors(x, &edi->vecs.radfix, edi);
320 project_to_eigvectors(x, &edi->vecs.radacc, edi);
321 project_to_eigvectors(x, &edi->vecs.radcon, edi);
325 static real calc_radius(t_eigvec *vec)
331 for (i=0; i<vec->neig; i++)
333 rad += pow((vec->refproj[i]-vec->xproj[i]),2);
336 return rad=sqrt(rad);
342 static void dump_xcoll(t_edpar *edi, struct t_do_edsam *buf, t_commrec *cr,
349 ivec *shifts, *eshifts;
356 shifts = buf->shifts_xcoll;
357 eshifts = buf->extra_shifts_xcoll;
359 sprintf(fn, "xcolldump_step%d.txt", step);
362 for (i=0; i<edi->sav.nr; i++)
364 fprintf(fp, "%d %9.5f %9.5f %9.5f %d %d %d %d %d %d\n",
366 xcoll[i][XX] , xcoll[i][YY] , xcoll[i][ZZ],
367 shifts[i][XX] , shifts[i][YY] , shifts[i][ZZ],
368 eshifts[i][XX], eshifts[i][YY], eshifts[i][ZZ]);
376 static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
381 fprintf(out, "#%s positions:\n%d\n", name, s->nr);
387 fprintf(out, "#index, x, y, z");
390 fprintf(out, ", sqrt(m)");
392 for (i=0; i<s->nr; i++)
394 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]);
397 fprintf(out,"%9.3f",s->sqrtm[i]);
405 static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
406 const char name[], int length)
411 fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig);
412 /* Dump the data for every eigenvector: */
413 for (i=0; i<ev->neig; i++)
415 fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
416 ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
417 for (j=0; j<length; j++)
419 fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
426 static void dump_edi(t_edpar *edpars, t_commrec *cr, int nr_edi)
432 sprintf(fn, "EDdump_node%d_edi%d", cr->nodeid, nr_edi);
433 out = ffopen(fn, "w");
435 fprintf(out,"#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n",
436 edpars->nini,edpars->fitmas,edpars->pcamas);
437 fprintf(out,"#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n",
438 edpars->outfrq,edpars->maxedsteps,edpars->slope);
439 fprintf(out,"#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n",
440 edpars->presteps,edpars->flood.deltaF0,edpars->flood.tau,
441 edpars->flood.constEfl,edpars->flood.alpha2);
443 /* Dump reference, average, target, origin positions */
444 dump_edi_positions(out, &edpars->sref, "REFERENCE");
445 dump_edi_positions(out, &edpars->sav , "AVERAGE" );
446 dump_edi_positions(out, &edpars->star, "TARGET" );
447 dump_edi_positions(out, &edpars->sori, "ORIGIN" );
449 /* Dump eigenvectors */
450 dump_edi_eigenvecs(out, &edpars->vecs.mon , "MONITORED", edpars->sav.nr);
451 dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX" , edpars->sav.nr);
452 dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC" , edpars->sav.nr);
453 dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX" , edpars->sav.nr);
454 dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC" , edpars->sav.nr);
455 dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON" , edpars->sav.nr);
457 /* Dump flooding eigenvectors */
458 dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING" , edpars->sav.nr);
460 /* Dump ed local buffer */
461 fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit );
462 fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam );
463 fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon );
470 static void dump_rotmat(FILE* out,matrix rotmat)
472 fprintf(out,"ROTMAT: %12.8f %12.8f %12.8f\n",rotmat[XX][XX],rotmat[XX][YY],rotmat[XX][ZZ]);
473 fprintf(out,"ROTMAT: %12.8f %12.8f %12.8f\n",rotmat[YY][XX],rotmat[YY][YY],rotmat[YY][ZZ]);
474 fprintf(out,"ROTMAT: %12.8f %12.8f %12.8f\n",rotmat[ZZ][XX],rotmat[ZZ][YY],rotmat[ZZ][ZZ]);
479 static void dump_rvec(FILE *out, int dim, rvec *x)
484 for (i=0; i<dim; i++)
486 fprintf(out,"%4d %f %f %f\n",i,x[i][XX],x[i][YY],x[i][ZZ]);
492 static void dump_mat(FILE* out, int dim, double** mat)
497 fprintf(out,"MATRIX:\n");
502 fprintf(out,"%f ",mat[i][j]);
515 static void do_edfit(int natoms,rvec *xp,rvec *x,matrix R,t_edpar *edi)
517 /* this is a copy of do_fit with some modifications */
524 struct t_do_edfit *loc;
527 if(edi->buf->do_edfit != NULL)
534 snew(edi->buf->do_edfit,1);
536 loc = edi->buf->do_edfit;
540 snew(loc->omega,2*DIM);
542 for(i=0; i<2*DIM; i++)
544 snew(loc->omega[i],2*DIM);
545 snew(loc->om[i],2*DIM);
559 /* calculate the matrix U */
561 for(n=0;(n<natoms);n++)
563 for(c=0; (c<DIM); c++)
566 for(r=0; (r<DIM); r++)
574 /* construct loc->omega */
575 /* loc->omega is symmetric -> loc->omega==loc->omega' */
582 loc->omega[r][c]=u[r-3][c];
583 loc->omega[c][r]=u[r-3][c];
593 /* determine h and k */
597 dump_mat(stderr,2*DIM,loc->omega);
600 fprintf(stderr,"d[%d] = %f\n",i,d[i]);
604 jacobi(loc->omega,6,d,loc->om,&irot);
608 fprintf(stderr,"IROT=0\n");
611 index=0; /* For the compiler only */
627 vh[j][i]=M_SQRT2*loc->om[i][index];
628 vk[j][i]=M_SQRT2*loc->om[i+DIM][index];
637 R[c][r]=vk[0][r]*vh[0][c]+
648 R[c][r]=vk[0][r]*vh[0][c]+
657 static void rmfit(int nat, rvec *xcoll, rvec transvec, matrix rotmat)
664 * The inverse rotation is described by the transposed rotation matrix */
665 transpose(rotmat,tmat);
666 rotate_x(xcoll, nat, tmat);
668 /* Remove translation */
669 vec[XX]=-transvec[XX];
670 vec[YY]=-transvec[YY];
671 vec[ZZ]=-transvec[ZZ];
672 translate_x(xcoll, nat, vec);
676 /**********************************************************************************
677 ******************** FLOODING ****************************************************
678 **********************************************************************************
680 The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose.
681 The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are
682 read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc).
684 do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in
685 the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure.
687 since the flooding acts on forces do_flood is called from the function force() (force.c), while the other
688 edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions.
690 do_flood makes a copy of the positions,
691 fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the
692 space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the
693 fit. Then do_flood adds these forces to the forcefield-forces
694 (given as parameter) and updates the adaptive flooding parameters Efl and deltaF.
696 To center the flooding potential at a different location one can use the -ori option in make_edi. The ori
697 structure is projected to the system of eigenvectors and then this position in the subspace is used as
698 center of the flooding potential. If the option is not used, the center will be zero in the subspace,
699 i.e. the average structure as given in the make_edi file.
701 To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted
702 signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of
703 Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly
704 so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no
705 further adaption is applied, Efl will stay constant at zero.
707 To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are
708 used as spring constants for the harmonic potential.
709 Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \
710 as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant.
712 To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi)
713 the routine read_edi_file reads all of theses flooding files.
714 The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list
715 calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one
716 edi there is no interdependence whatsoever. The forces are added together.
718 To write energies into the .edr file, call the function
719 get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln)
721 get_flood_energies(real Vfl[],int nnames);
724 - 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.
726 Maybe one should give a range of atoms for which to remove motion, so that motion is removed with
727 two edsam files from two peptide chains
730 static void write_edo_flood(t_edpar *edi, FILE *fp, gmx_large_int_t step)
734 gmx_bool bOutputRef=FALSE;
737 fprintf(fp,"%d.th FL: %s %12.5e %12.5e %12.5e\n",
738 edi->flood.flood_id, gmx_step_str(step,buf),
739 edi->flood.Efl, edi->flood.Vfl, edi->flood.deltaF);
742 /* Check whether any of the references changes with time (this can happen
743 * in case flooding is used as harmonic restraint). If so, output all the
744 * current reference projections. */
745 if (edi->flood.bHarmonic)
747 for (i = 0; i < edi->flood.vecs.neig; i++)
749 if (edi->flood.vecs.refprojslope[i] != 0.0)
756 fprintf(fp, "Ref. projs.: ");
757 for (i = 0; i < edi->flood.vecs.neig; i++)
759 fprintf(fp, "%12.5e ", edi->flood.vecs.refproj[i]);
764 fprintf(fp,"FL_FORCES: ");
766 for (i=0; i<edi->flood.vecs.neig; i++)
768 fprintf(fp," %12.5e",edi->flood.vecs.fproj[i]);
775 /* From flood.xproj compute the Vfl(x) at this point */
776 static real flood_energy(t_edpar *edi, gmx_large_int_t step)
778 /* compute flooding energy Vfl
779 Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } )
780 \lambda_i is the reciprocal eigenvalue 1/\sigma_i
781 it is already computed by make_edi and stored in stpsz[i]
783 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
790 /* Each time this routine is called (i.e. each time step), we add a small
791 * value to the reference projection. This way a harmonic restraint towards
792 * a moving reference is realized. If no value for the additive constant
793 * is provided in the edi file, the reference will not change. */
794 if (edi->flood.bHarmonic)
796 for (i=0; i<edi->flood.vecs.neig; i++)
798 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
803 /* Compute sum which will be the exponent of the exponential */
804 for (i=0; i<edi->flood.vecs.neig; i++)
806 /* stpsz stores the reciprocal eigenvalue 1/sigma_i */
807 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]);
810 /* Compute the Gauss function*/
811 if (edi->flood.bHarmonic)
813 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
817 Vfl = edi->flood.Efl!=0 ? edi->flood.Efl*exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) :0;
824 /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */
825 static void flood_forces(t_edpar *edi)
827 /* compute the forces in the subspace of the flooding eigenvectors
828 * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */
831 real energy=edi->flood.Vfl;
834 if (edi->flood.bHarmonic)
836 for (i=0; i<edi->flood.vecs.neig; i++)
838 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
843 for (i=0; i<edi->flood.vecs.neig; i++)
845 /* if Efl is zero the forces are zero if not use the formula */
846 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;
852 /* Raise forces from subspace into cartesian space */
853 static void flood_blowup(t_edpar *edi, rvec *forces_cart)
855 /* this function lifts the forces from the subspace to the cartesian space
856 all the values not contained in the subspace are assumed to be zero and then
857 a coordinate transformation from eigenvector to cartesian vectors is performed
858 The nonexistent values don't have to be set to zero explicitly, they would occur
859 as zero valued summands, hence we just stop to compute this part of the sum.
861 for every atom we add all the contributions to this atom from all the different eigenvectors.
863 NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the
864 field forces_cart prior the computation, but we compute the forces separately
865 to have them accessible for diagnostics
872 forces_sub = edi->flood.vecs.fproj;
875 /* Calculate the cartesian forces for the local atoms */
877 /* Clear forces first */
878 for (j=0; j<edi->sav.nr_loc; j++)
880 clear_rvec(forces_cart[j]);
883 /* Now compute atomwise */
884 for (j=0; j<edi->sav.nr_loc; j++)
886 /* Compute forces_cart[edi->sav.anrs[j]] */
887 for (eig=0; eig<edi->flood.vecs.neig; eig++)
889 /* Force vector is force * eigenvector (compute only atom j) */
890 svmul(forces_sub[eig],edi->flood.vecs.vec[eig][edi->sav.c_ind[j]],dum);
891 /* Add this vector to the cartesian forces */
892 rvec_inc(forces_cart[j],dum);
898 /* Update the values of Efl, deltaF depending on tau and Vfl */
899 static void update_adaption(t_edpar *edi)
901 /* this function updates the parameter Efl and deltaF according to the rules given in
902 * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al,
905 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
907 edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
908 /* check if restrain (inverted flooding) -> don't let EFL become positive */
909 if (edi->flood.alpha2<0 && edi->flood.Efl>-0.00000001)
914 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
919 static void do_single_flood(
924 gmx_large_int_t step,
927 gmx_bool bNS) /* Are we in a neighbor searching step? */
930 matrix rotmat; /* rotation matrix */
931 matrix tmat; /* inverse rotation */
932 rvec transvec; /* translation vector */
933 struct t_do_edsam *buf;
936 buf=edi->buf->do_edsam;
939 /* Broadcast the positions of the AVERAGE structure such that they are known on
940 * every processor. Each node contributes its local positions x and stores them in
941 * the collective ED array buf->xcoll */
942 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
943 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
945 /* Only assembly REFERENCE positions if their indices differ from the average ones */
948 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
949 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
952 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
953 * We do not need to update the shifts until the next NS step */
954 buf->bUpdateShifts = FALSE;
956 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
957 * as well as the indices in edi->sav.anrs */
959 /* Fit the reference indices to the reference structure */
962 fit_to_reference(buf->xcoll , transvec, rotmat, edi);
966 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
969 /* Now apply the translation and rotation to the ED structure */
970 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
972 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
973 project_to_eigvectors(buf->xcoll,&edi->flood.vecs,edi);
975 if (FALSE == edi->flood.bConstForce)
977 /* Compute Vfl(x) from flood.xproj */
978 edi->flood.Vfl = flood_energy(edi, step);
980 update_adaption(edi);
982 /* Compute the flooding forces */
986 /* Translate them into cartesian positions */
987 flood_blowup(edi, edi->flood.forces_cartesian);
989 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
990 /* Each node rotates back its local forces */
991 transpose(rotmat,tmat);
992 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
994 /* Finally add forces to the main force variable */
995 for (i=0; i<edi->sav.nr_loc; i++)
997 rvec_inc(force[edi->sav.anrs_loc[i]],edi->flood.forces_cartesian[i]);
1000 /* Output is written by the master process */
1001 if (do_per_step(step,edi->outfrq) && MASTER(cr))
1003 write_edo_flood(edi,edo,step);
1008 /* Main flooding routine, called from do_force */
1009 extern void do_flood(
1010 FILE *log, /* md.log file */
1011 t_commrec *cr, /* Communication record */
1012 rvec x[], /* Positions on the local processor */
1013 rvec force[], /* forcefield forces, to these the flooding forces are added */
1014 gmx_edsam_t ed, /* ed data structure contains all ED and flooding datasets */
1015 matrix box, /* the box */
1016 gmx_large_int_t step, /* The relative time step since ir->init_step is already subtracted */
1017 gmx_bool bNS) /* Are we in a neighbor searching step? */
1022 if (ed->eEDtype != eEDflood)
1028 /* Call flooding for one matrix */
1029 if (edi->flood.vecs.neig)
1031 do_single_flood(ed->edo,x,force,edi,step,box,cr,bNS);
1033 edi = edi->next_edi;
1038 /* Called by init_edi, configure some flooding related variables and structures,
1039 * print headers to output files */
1040 static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt, t_commrec *cr, gmx_bool bPrintheader)
1045 edi->flood.Efl = edi->flood.constEfl;
1049 if (edi->flood.vecs.neig)
1051 /* If in any of the datasets we find a flooding vector, flooding is turned on */
1052 ed->eEDtype = eEDflood;
1054 fprintf(stderr,"ED: Flooding of matrix %d is switched on.\n", edi->flood.flood_id);
1056 if (edi->flood.bConstForce)
1058 /* We have used stpsz as a vehicle to carry the fproj values for constant
1059 * force flooding. Now we copy that to flood.vecs.fproj. Note that
1060 * in const force flooding, fproj is never changed. */
1061 for (i=0; i<edi->flood.vecs.neig; i++)
1063 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
1065 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
1066 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
1072 fprintf(ed->edo,"FL_HEADER: Flooding of matrix %d is switched on! The flooding output will have the following format:\n",
1073 edi->flood.flood_id);
1074 fprintf(ed->edo,"FL_HEADER: Step Efl Vfl deltaF\n");
1081 /*********** Energy book keeping ******/
1082 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1091 srenew(names,count);
1092 sprintf(buf,"Vfl_%d",count);
1093 names[count-1]=strdup(buf);
1094 actual=actual->next_edi;
1101 static void get_flood_energies(t_edpar *edi, real Vfl[],int nnames)
1103 /*fl has to be big enough to capture nnames-many entries*/
1112 Vfl[count-1]=actual->flood.Vfl;
1113 actual=actual->next_edi;
1116 if (nnames!=count-1)
1118 gmx_fatal(FARGS,"Number of energies is not consistent with t_edi structure");
1121 /************* END of FLOODING IMPLEMENTATION ****************************/
1125 gmx_edsam_t ed_open(int nfile,const t_filenm fnm[],unsigned long Flags,t_commrec *cr)
1130 /* Allocate space for the ED data structure */
1133 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1134 ed->eEDtype = eEDedsam;
1138 /* Open .edi input file: */
1139 ed->edinam=ftp2fn(efEDI,nfile,fnm);
1140 /* The master opens the .edo output file */
1141 fprintf(stderr,"ED sampling will be performed!\n");
1142 ed->edonam = ftp2fn(efEDO,nfile,fnm);
1143 ed->edo = gmx_fio_fopen(ed->edonam,(Flags & MD_APPENDFILES)? "a+" : "w+");
1149 /* Broadcasts the structure data */
1150 static void bc_ed_positions(t_commrec *cr, struct gmx_edx *s, int stype)
1152 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1153 snew_bc(cr, s->x , s->nr ); /* Positions */
1154 nblock_bc(cr, s->nr, s->anrs );
1155 nblock_bc(cr, s->nr, s->x );
1157 /* For the average & reference structures we need an array for the collective indices,
1158 * and we need to broadcast the masses as well */
1159 if (stype == eedAV || stype == eedREF)
1161 /* We need these additional variables in the parallel case: */
1162 snew(s->c_ind , s->nr ); /* Collective indices */
1163 /* Local atom indices get assigned in dd_make_local_group_indices.
1164 * There, also memory is allocated */
1165 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1166 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1167 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1170 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1171 if (stype == eedREF)
1173 snew_bc(cr, s->m, s->nr);
1174 nblock_bc(cr, s->nr, s->m);
1177 /* For the average structure we might need the masses for mass-weighting */
1180 snew_bc(cr, s->sqrtm, s->nr);
1181 nblock_bc(cr, s->nr, s->sqrtm);
1182 snew_bc(cr, s->m, s->nr);
1183 nblock_bc(cr, s->nr, s->m);
1188 /* Broadcasts the eigenvector data */
1189 static void bc_ed_vecs(t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
1193 snew_bc(cr, ev->ieig , ev->neig); /* index numbers of eigenvector */
1194 snew_bc(cr, ev->stpsz , ev->neig); /* stepsizes per eigenvector */
1195 snew_bc(cr, ev->xproj , ev->neig); /* instantaneous x projection */
1196 snew_bc(cr, ev->fproj , ev->neig); /* instantaneous f projection */
1197 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1199 nblock_bc(cr, ev->neig, ev->ieig );
1200 nblock_bc(cr, ev->neig, ev->stpsz );
1201 nblock_bc(cr, ev->neig, ev->xproj );
1202 nblock_bc(cr, ev->neig, ev->fproj );
1203 nblock_bc(cr, ev->neig, ev->refproj);
1205 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1206 for (i=0; i<ev->neig; i++)
1208 snew_bc(cr, ev->vec[i], length);
1209 nblock_bc(cr, length, ev->vec[i]);
1212 /* For harmonic restraints the reference projections can change with time */
1215 snew_bc(cr, ev->refproj0 , ev->neig);
1216 snew_bc(cr, ev->refprojslope, ev->neig);
1217 nblock_bc(cr, ev->neig, ev->refproj0 );
1218 nblock_bc(cr, ev->neig, ev->refprojslope);
1223 /* Broadcasts the ED / flooding data to other nodes
1224 * and allocates memory where needed */
1225 static void broadcast_ed_data(t_commrec *cr, gmx_edsam_t ed, int numedis)
1231 /* Master lets the other nodes know if its ED only or also flooding */
1232 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1234 snew_bc(cr, ed->edpar,1);
1235 /* Now transfer the ED data set(s) */
1237 for (nr=0; nr<numedis; nr++)
1239 /* Broadcast a single ED data set */
1242 /* Broadcast positions */
1243 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1244 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1245 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1246 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1248 /* Broadcast eigenvectors */
1249 bc_ed_vecs(cr, &edi->vecs.mon , edi->sav.nr, FALSE);
1250 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
1251 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
1252 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
1253 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
1254 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
1255 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1256 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
1258 /* Set the pointer to the next ED dataset */
1261 snew_bc(cr, edi->next_edi, 1);
1262 edi = edi->next_edi;
1268 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1269 static void init_edi(gmx_mtop_t *mtop,t_inputrec *ir,
1270 t_commrec *cr,gmx_edsam_t ed,t_edpar *edi)
1273 real totalmass = 0.0;
1275 gmx_mtop_atomlookup_t alook=NULL;
1278 /* NOTE Init_edi is executed on the master process only
1279 * The initialized data sets are then transmitted to the
1280 * other nodes in broadcast_ed_data */
1282 edi->bNeedDoEdsam = edi->vecs.mon.neig
1283 || edi->vecs.linfix.neig
1284 || edi->vecs.linacc.neig
1285 || edi->vecs.radfix.neig
1286 || edi->vecs.radacc.neig
1287 || edi->vecs.radcon.neig;
1289 alook = gmx_mtop_atomlookup_init(mtop);
1291 /* evaluate masses (reference structure) */
1292 snew(edi->sref.m, edi->sref.nr);
1293 for (i = 0; i < edi->sref.nr; i++)
1297 gmx_mtop_atomnr_to_atom(alook,edi->sref.anrs[i],&atom);
1298 edi->sref.m[i] = atom->m;
1302 edi->sref.m[i] = 1.0;
1305 /* Check that every m > 0. Bad things will happen otherwise. */
1306 if (edi->sref.m[i] <= 0.0)
1308 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1309 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1310 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1311 "atoms from the reference structure by creating a proper index group.\n",
1312 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1315 totalmass += edi->sref.m[i];
1317 edi->sref.mtot = totalmass;
1319 /* Masses m and sqrt(m) for the average structure. Note that m
1320 * is needed if forces have to be evaluated in do_edsam */
1321 snew(edi->sav.sqrtm, edi->sav.nr );
1322 snew(edi->sav.m , edi->sav.nr );
1323 for (i = 0; i < edi->sav.nr; i++)
1325 gmx_mtop_atomnr_to_atom(alook,edi->sav.anrs[i],&atom);
1326 edi->sav.m[i] = atom->m;
1329 edi->sav.sqrtm[i] = sqrt(atom->m);
1333 edi->sav.sqrtm[i] = 1.0;
1336 /* Check that every m > 0. Bad things will happen otherwise. */
1337 if (edi->sav.sqrtm[i] <= 0.0)
1339 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1340 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1341 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1342 "atoms from the average structure by creating a proper index group.\n",
1343 i, edi->sav.anrs[i]+1, atom->m);
1347 gmx_mtop_atomlookup_destroy(alook);
1349 /* put reference structure in origin */
1350 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1354 translate_x(edi->sref.x, edi->sref.nr, com);
1356 /* Init ED buffer */
1361 static void check(const char *line, const char *label)
1363 if (!strstr(line,label))
1365 gmx_fatal(FARGS,"Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s",label,line);
1370 static int read_checked_edint(FILE *file,const char *label)
1372 char line[STRLEN+1];
1376 fgets2 (line,STRLEN,file);
1378 fgets2 (line,STRLEN,file);
1379 sscanf (line,"%d",&idum);
1384 static int read_edint(FILE *file,gmx_bool *bEOF)
1386 char line[STRLEN+1];
1391 eof=fgets2 (line,STRLEN,file);
1397 eof=fgets2 (line,STRLEN,file);
1403 sscanf (line,"%d",&idum);
1409 static real read_checked_edreal(FILE *file,const char *label)
1411 char line[STRLEN+1];
1415 fgets2 (line,STRLEN,file);
1417 fgets2 (line,STRLEN,file);
1418 sscanf (line,"%lf",&rdum);
1419 return (real) rdum; /* always read as double and convert to single */
1423 static void read_edx(FILE *file,int number,int *anrs,rvec *x)
1426 char line[STRLEN+1];
1430 for(i=0; i<number; i++)
1432 fgets2 (line,STRLEN,file);
1433 sscanf (line,"%d%lf%lf%lf",&anrs[i],&d[0],&d[1],&d[2]);
1434 anrs[i]--; /* we are reading FORTRAN indices */
1437 x[i][j]=d[j]; /* always read as double and convert to single */
1443 static void scan_edvec(FILE *in,int nr,rvec *vec)
1445 char line[STRLEN+1];
1450 for(i=0; (i < nr); i++)
1452 fgets2 (line,STRLEN,in);
1453 sscanf (line,"%le%le%le",&x,&y,&z);
1461 static void read_edvec(FILE *in,int nr,t_eigvec *tvec,gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1464 double rdum,refproj_dum=0.0,refprojslope_dum=0.0;
1465 char line[STRLEN+1];
1468 tvec->neig=read_checked_edint(in,"NUMBER OF EIGENVECTORS");
1471 snew(tvec->ieig ,tvec->neig);
1472 snew(tvec->stpsz ,tvec->neig);
1473 snew(tvec->vec ,tvec->neig);
1474 snew(tvec->xproj ,tvec->neig);
1475 snew(tvec->fproj ,tvec->neig);
1476 snew(tvec->refproj,tvec->neig);
1479 snew(tvec->refproj0 ,tvec->neig);
1480 snew(tvec->refprojslope,tvec->neig);
1483 for(i=0; (i < tvec->neig); i++)
1485 fgets2 (line,STRLEN,in);
1486 if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1488 nscan = sscanf(line,"%d%lf%lf%lf",&idum,&rdum,&refproj_dum,&refprojslope_dum);
1489 /* Zero out values which were not scanned */
1493 /* Every 4 values read, including reference position */
1494 *bHaveReference = TRUE;
1497 /* A reference position is provided */
1498 *bHaveReference = TRUE;
1499 /* No value for slope, set to 0 */
1500 refprojslope_dum = 0.0;
1503 /* No values for reference projection and slope, set to 0 */
1505 refprojslope_dum = 0.0;
1508 gmx_fatal(FARGS,"Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1511 tvec->refproj[i]=refproj_dum;
1512 tvec->refproj0[i]=refproj_dum;
1513 tvec->refprojslope[i]=refprojslope_dum;
1515 else /* Normal flooding */
1517 nscan = sscanf(line,"%d%lf",&idum,&rdum);
1520 gmx_fatal(FARGS,"Expected 2 values for flooding vec: <nr> <stpsz>\n");
1524 tvec->stpsz[i]=rdum;
1525 } /* end of loop over eigenvectors */
1527 for(i=0; (i < tvec->neig); i++)
1529 snew(tvec->vec[i],nr);
1530 scan_edvec(in,nr,tvec->vec[i]);
1536 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1537 static void read_edvecs(FILE *in,int nr,t_edvecs *vecs)
1539 gmx_bool bHaveReference = FALSE;
1542 read_edvec(in, nr, &vecs->mon , FALSE, &bHaveReference);
1543 read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1544 read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1545 read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1546 read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1547 read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1551 /* Check if the same atom indices are used for reference and average positions */
1552 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1557 /* If the number of atoms differs between the two structures,
1558 * they cannot be identical */
1559 if (sref.nr != sav.nr)
1564 /* Now that we know that both stuctures have the same number of atoms,
1565 * check if also the indices are identical */
1566 for (i=0; i < sav.nr; i++)
1568 if (sref.anrs[i] != sav.anrs[i])
1573 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1579 static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int edi_nr, t_commrec *cr)
1582 const int magic=670;
1585 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1586 gmx_bool bHaveReference = FALSE;
1589 /* the edi file is not free format, so expect problems if the input is corrupt. */
1591 /* check the magic number */
1592 readmagic=read_edint(in,&bEOF);
1593 /* Check whether we have reached the end of the input file */
1599 if (readmagic != magic)
1601 if (readmagic==666 || readmagic==667 || readmagic==668)
1603 gmx_fatal(FARGS,"Wrong magic number: Use newest version of make_edi to produce edi file");
1605 else if (readmagic != 669)
1607 gmx_fatal(FARGS,"Wrong magic number %d in %s",readmagic,ed->edinam);
1611 /* check the number of atoms */
1612 edi->nini=read_edint(in,&bEOF);
1613 if (edi->nini != nr_mdatoms)
1615 gmx_fatal(FARGS,"Nr of atoms in %s (%d) does not match nr of md atoms (%d)",
1616 ed->edinam,edi->nini,nr_mdatoms);
1619 /* Done checking. For the rest we blindly trust the input */
1620 edi->fitmas = read_checked_edint(in,"FITMAS");
1621 edi->pcamas = read_checked_edint(in,"ANALYSIS_MAS");
1622 edi->outfrq = read_checked_edint(in,"OUTFRQ");
1623 edi->maxedsteps = read_checked_edint(in,"MAXLEN");
1624 edi->slope = read_checked_edreal(in,"SLOPECRIT");
1626 edi->presteps = read_checked_edint(in,"PRESTEPS");
1627 edi->flood.deltaF0 = read_checked_edreal(in,"DELTA_F0");
1628 edi->flood.deltaF = read_checked_edreal(in,"INIT_DELTA_F");
1629 edi->flood.tau = read_checked_edreal(in,"TAU");
1630 edi->flood.constEfl = read_checked_edreal(in,"EFL_NULL");
1631 edi->flood.alpha2 = read_checked_edreal(in,"ALPHA2");
1632 edi->flood.kT = read_checked_edreal(in,"KT");
1633 edi->flood.bHarmonic = read_checked_edint(in,"HARMONIC");
1634 if (readmagic > 669)
1636 edi->flood.bConstForce = read_checked_edint(in,"CONST_FORCE_FLOODING");
1640 edi->flood.bConstForce = FALSE;
1642 edi->flood.flood_id = edi_nr;
1643 edi->sref.nr = read_checked_edint(in,"NREF");
1645 /* allocate space for reference positions and read them */
1646 snew(edi->sref.anrs,edi->sref.nr);
1647 snew(edi->sref.x ,edi->sref.nr);
1648 snew(edi->sref.x_old,edi->sref.nr);
1649 edi->sref.sqrtm =NULL;
1650 read_edx(in,edi->sref.nr,edi->sref.anrs,edi->sref.x);
1652 /* average positions. they define which atoms will be used for ED sampling */
1653 edi->sav.nr=read_checked_edint(in,"NAV");
1654 snew(edi->sav.anrs,edi->sav.nr);
1655 snew(edi->sav.x ,edi->sav.nr);
1656 snew(edi->sav.x_old,edi->sav.nr);
1657 read_edx(in,edi->sav.nr,edi->sav.anrs,edi->sav.x);
1659 /* Check if the same atom indices are used for reference and average positions */
1660 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1663 read_edvecs(in,edi->sav.nr,&edi->vecs);
1664 read_edvec(in,edi->sav.nr,&edi->flood.vecs,edi->flood.bHarmonic, &bHaveReference);
1666 /* target positions */
1667 edi->star.nr=read_edint(in,&bEOF);
1668 if (edi->star.nr > 0)
1670 snew(edi->star.anrs,edi->star.nr);
1671 snew(edi->star.x ,edi->star.nr);
1672 edi->star.sqrtm =NULL;
1673 read_edx(in,edi->star.nr,edi->star.anrs,edi->star.x);
1676 /* positions defining origin of expansion circle */
1677 edi->sori.nr=read_edint(in,&bEOF);
1678 if (edi->sori.nr > 0)
1682 /* Both an -ori structure and a at least one manual reference point have been
1683 * specified. That's ambiguous and probably not intentional. */
1684 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1685 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1687 snew(edi->sori.anrs,edi->sori.nr);
1688 snew(edi->sori.x ,edi->sori.nr);
1689 edi->sori.sqrtm =NULL;
1690 read_edx(in,edi->sori.nr,edi->sori.anrs,edi->sori.x);
1699 /* Read in the edi input file. Note that it may contain several ED data sets which were
1700 * achieved by concatenating multiple edi files. The standard case would be a single ED
1701 * data set, though. */
1702 static int read_edi_file(gmx_edsam_t ed, t_edpar *edi, int nr_mdatoms, t_commrec *cr)
1705 t_edpar *curr_edi,*last_edi;
1710 /* This routine is executed on the master only */
1712 /* Open the .edi parameter input file */
1713 in = gmx_fio_fopen(ed->edinam,"r");
1714 fprintf(stderr, "ED: Reading edi file %s\n", ed->edinam);
1716 /* Now read a sequence of ED input parameter sets from the edi file */
1719 while( read_edi(in, ed, curr_edi, nr_mdatoms, edi_nr, cr) )
1722 /* Make shure that the number of atoms in each dataset is the same as in the tpr file */
1723 if (edi->nini != nr_mdatoms)
1725 gmx_fatal(FARGS,"edi file %s (dataset #%d) was made for %d atoms, but the simulation contains %d atoms.",
1726 ed->edinam, edi_nr, edi->nini, nr_mdatoms);
1728 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1729 /* We need to allocate space for the data: */
1731 /* Point the 'next_edi' entry to the next edi: */
1732 curr_edi->next_edi=edi_read;
1733 /* Keep the curr_edi pointer for the case that the next dataset is empty: */
1734 last_edi = curr_edi;
1735 /* Let's prepare to read in the next edi data set: */
1736 curr_edi = edi_read;
1740 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", ed->edinam);
1743 /* Terminate the edi dataset list with a NULL pointer: */
1744 last_edi->next_edi = NULL;
1746 fprintf(stderr, "ED: Found %d ED dataset%s.\n", edi_nr, edi_nr>1? "s" : "");
1748 /* Close the .edi file again */
1755 struct t_fit_to_ref {
1756 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1759 /* Fit the current positions to the reference positions
1760 * Do not actually do the fit, just return rotation and translation.
1761 * Note that the COM of the reference structure was already put into
1762 * the origin by init_edi. */
1763 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1764 rvec transvec, /* The translation vector */
1765 matrix rotmat, /* The rotation matrix */
1766 t_edpar *edi) /* Just needed for do_edfit */
1768 rvec com; /* center of mass */
1770 struct t_fit_to_ref *loc;
1773 GMX_MPE_LOG(ev_fit_to_reference_start);
1775 /* Allocate memory the first time this routine is called for each edi dataset */
1776 if (NULL == edi->buf->fit_to_ref)
1778 snew(edi->buf->fit_to_ref, 1);
1779 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1781 loc = edi->buf->fit_to_ref;
1783 /* We do not touch the original positions but work on a copy. */
1784 for (i=0; i<edi->sref.nr; i++)
1786 copy_rvec(xcoll[i], loc->xcopy[i]);
1789 /* Calculate the center of mass */
1790 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1792 transvec[XX] = -com[XX];
1793 transvec[YY] = -com[YY];
1794 transvec[ZZ] = -com[ZZ];
1796 /* Subtract the center of mass from the copy */
1797 translate_x(loc->xcopy, edi->sref.nr, transvec);
1799 /* Determine the rotation matrix */
1800 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1802 GMX_MPE_LOG(ev_fit_to_reference_finish);
1806 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1807 int nat, /* How many positions are there? */
1808 rvec transvec, /* The translation vector */
1809 matrix rotmat) /* The rotation matrix */
1812 translate_x(x, nat, transvec);
1815 rotate_x(x, nat, rotmat);
1819 /* Gets the rms deviation of the positions to the structure s */
1820 /* fit_to_structure has to be called before calling this routine! */
1821 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1822 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1828 for (i=0; i < s->nr; i++)
1830 rmsd += distance2(s->x[i], x[i]);
1833 rmsd /= (real) s->nr;
1840 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1845 if (ed->eEDtype != eEDnone)
1847 /* Loop over ED datasets (usually there is just one dataset, though) */
1851 /* Local atoms of the reference structure (for fitting), need only be assembled
1852 * if their indices differ from the average ones */
1855 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1856 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1859 /* Local atoms of the average structure (on these ED will be performed) */
1860 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1861 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1863 /* Indicate that the ED shift vectors for this structure need to be updated
1864 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1865 edi->buf->do_edsam->bUpdateShifts = TRUE;
1867 /* Set the pointer to the next ED dataset (if any) */
1874 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1879 GMX_MPE_LOG(ev_unshift_start);
1887 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1888 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1889 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1893 xu[XX] = x[XX]-tx*box[XX][XX];
1894 xu[YY] = x[YY]-ty*box[YY][YY];
1895 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1898 GMX_MPE_LOG(ev_unshift_finish);
1902 static void do_linfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
1909 /* loop over linfix vectors */
1910 for (i=0; i<edi->vecs.linfix.neig; i++)
1912 /* calculate the projection */
1913 proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1915 /* calculate the correction */
1916 add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1918 /* apply the correction */
1919 add /= edi->sav.sqrtm[i];
1920 for (j=0; j<edi->sav.nr; j++)
1922 svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1923 rvec_inc(xcoll[j], vec_dum);
1929 static void do_linacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1936 /* loop over linacc vectors */
1937 for (i=0; i<edi->vecs.linacc.neig; i++)
1939 /* calculate the projection */
1940 proj=projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1942 /* calculate the correction */
1944 if (edi->vecs.linacc.stpsz[i] > 0.0)
1946 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1948 add = edi->vecs.linacc.refproj[i] - proj;
1951 if (edi->vecs.linacc.stpsz[i] < 0.0)
1953 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1955 add = edi->vecs.linacc.refproj[i] - proj;
1959 /* apply the correction */
1960 add /= edi->sav.sqrtm[i];
1961 for (j=0; j<edi->sav.nr; j++)
1963 svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
1964 rvec_inc(xcoll[j], vec_dum);
1967 /* new positions will act as reference */
1968 edi->vecs.linacc.refproj[i] = proj + add;
1973 static void do_radfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
1976 real *proj, rad=0.0, ratio;
1980 if (edi->vecs.radfix.neig == 0)
1983 snew(proj, edi->vecs.radfix.neig);
1985 /* loop over radfix vectors */
1986 for (i=0; i<edi->vecs.radfix.neig; i++)
1988 /* calculate the projections, radius */
1989 proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
1990 rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
1994 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
1995 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
1997 /* loop over radfix vectors */
1998 for (i=0; i<edi->vecs.radfix.neig; i++)
2000 proj[i] -= edi->vecs.radfix.refproj[i];
2002 /* apply the correction */
2003 proj[i] /= edi->sav.sqrtm[i];
2005 for (j=0; j<edi->sav.nr; j++)
2007 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
2008 rvec_inc(xcoll[j], vec_dum);
2016 static void do_radacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
2019 real *proj, rad=0.0, ratio=0.0;
2023 if (edi->vecs.radacc.neig == 0)
2026 snew(proj,edi->vecs.radacc.neig);
2028 /* loop over radacc vectors */
2029 for (i=0; i<edi->vecs.radacc.neig; i++)
2031 /* calculate the projections, radius */
2032 proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
2033 rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
2037 /* only correct when radius decreased */
2038 if (rad < edi->vecs.radacc.radius)
2040 ratio = edi->vecs.radacc.radius/rad - 1.0;
2041 rad = edi->vecs.radacc.radius;
2044 edi->vecs.radacc.radius = rad;
2046 /* loop over radacc vectors */
2047 for (i=0; i<edi->vecs.radacc.neig; i++)
2049 proj[i] -= edi->vecs.radacc.refproj[i];
2051 /* apply the correction */
2052 proj[i] /= edi->sav.sqrtm[i];
2054 for (j=0; j<edi->sav.nr; j++)
2056 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
2057 rvec_inc(xcoll[j], vec_dum);
2064 struct t_do_radcon {
2068 static void do_radcon(rvec *xcoll, t_edpar *edi, t_commrec *cr)
2071 real rad=0.0, ratio=0.0;
2072 struct t_do_radcon *loc;
2077 if(edi->buf->do_radcon != NULL)
2080 loc = edi->buf->do_radcon;
2085 snew(edi->buf->do_radcon, 1);
2087 loc = edi->buf->do_radcon;
2089 if (edi->vecs.radcon.neig == 0)
2093 snew(loc->proj, edi->vecs.radcon.neig);
2095 /* loop over radcon vectors */
2096 for (i=0; i<edi->vecs.radcon.neig; i++)
2098 /* calculate the projections, radius */
2099 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2100 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2103 /* only correct when radius increased */
2104 if (rad > edi->vecs.radcon.radius)
2106 ratio = edi->vecs.radcon.radius/rad - 1.0;
2108 /* loop over radcon vectors */
2109 for (i=0; i<edi->vecs.radcon.neig; i++)
2111 /* apply the correction */
2112 loc->proj[i] -= edi->vecs.radcon.refproj[i];
2113 loc->proj[i] /= edi->sav.sqrtm[i];
2114 loc->proj[i] *= ratio;
2116 for (j=0; j<edi->sav.nr; j++)
2118 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
2119 rvec_inc(xcoll[j], vec_dum);
2124 edi->vecs.radcon.radius = rad;
2126 if (rad != edi->vecs.radcon.radius)
2129 for (i=0; i<edi->vecs.radcon.neig; i++)
2131 /* calculate the projections, radius */
2132 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2133 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2140 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step, t_commrec *cr)
2145 GMX_MPE_LOG(ev_ed_apply_cons_start);
2147 /* subtract the average positions */
2148 for (i=0; i<edi->sav.nr; i++)
2150 rvec_dec(xcoll[i], edi->sav.x[i]);
2153 /* apply the constraints */
2156 do_linfix(xcoll, edi, step, cr);
2158 do_linacc(xcoll, edi, cr);
2161 do_radfix(xcoll, edi, step, cr);
2163 do_radacc(xcoll, edi, cr);
2164 do_radcon(xcoll, edi, cr);
2166 /* add back the average positions */
2167 for (i=0; i<edi->sav.nr; i++)
2169 rvec_inc(xcoll[i], edi->sav.x[i]);
2172 GMX_MPE_LOG(ev_ed_apply_cons_finish);
2176 /* Write out the projections onto the eigenvectors */
2177 static void write_edo(int nr_edi, t_edpar *edi, gmx_edsam_t ed, gmx_large_int_t step,real rmsd)
2183 if (edi->bNeedDoEdsam)
2187 fprintf(ed->edo, "Initial projections:\n");
2191 fprintf(ed->edo,"Step %s, ED #%d ", gmx_step_str(step, buf), nr_edi);
2192 fprintf(ed->edo," RMSD %f nm\n",rmsd);
2195 if (edi->vecs.mon.neig)
2197 fprintf(ed->edo," Monitor eigenvectors");
2198 for (i=0; i<edi->vecs.mon.neig; i++)
2200 fprintf(ed->edo," %d: %12.5e ",edi->vecs.mon.ieig[i],edi->vecs.mon.xproj[i]);
2202 fprintf(ed->edo,"\n");
2204 if (edi->vecs.linfix.neig)
2206 fprintf(ed->edo," Linfix eigenvectors");
2207 for (i=0; i<edi->vecs.linfix.neig; i++)
2209 fprintf(ed->edo," %d: %12.5e ",edi->vecs.linfix.ieig[i],edi->vecs.linfix.xproj[i]);
2211 fprintf(ed->edo,"\n");
2213 if (edi->vecs.linacc.neig)
2215 fprintf(ed->edo," Linacc eigenvectors");
2216 for (i=0; i<edi->vecs.linacc.neig; i++)
2218 fprintf(ed->edo," %d: %12.5e ",edi->vecs.linacc.ieig[i],edi->vecs.linacc.xproj[i]);
2220 fprintf(ed->edo,"\n");
2222 if (edi->vecs.radfix.neig)
2224 fprintf(ed->edo," Radfix eigenvectors");
2225 for (i=0; i<edi->vecs.radfix.neig; i++)
2227 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radfix.ieig[i],edi->vecs.radfix.xproj[i]);
2229 fprintf(ed->edo,"\n");
2230 fprintf(ed->edo," fixed increment radius = %f\n", calc_radius(&edi->vecs.radfix));
2232 if (edi->vecs.radacc.neig)
2234 fprintf(ed->edo," Radacc eigenvectors");
2235 for (i=0; i<edi->vecs.radacc.neig; i++)
2237 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radacc.ieig[i],edi->vecs.radacc.xproj[i]);
2239 fprintf(ed->edo,"\n");
2240 fprintf(ed->edo," acceptance radius = %f\n", calc_radius(&edi->vecs.radacc));
2242 if (edi->vecs.radcon.neig)
2244 fprintf(ed->edo," Radcon eigenvectors");
2245 for (i=0; i<edi->vecs.radcon.neig; i++)
2247 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radcon.ieig[i],edi->vecs.radcon.xproj[i]);
2249 fprintf(ed->edo,"\n");
2250 fprintf(ed->edo," contracting radius = %f\n", calc_radius(&edi->vecs.radcon));
2255 /* Returns if any constraints are switched on */
2256 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2258 if (edtype == eEDedsam || edtype == eEDflood)
2260 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2261 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2262 edi->vecs.radcon.neig);
2268 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2269 * umbrella sampling simulations. */
2270 static void copyEvecReference(t_eigvec* floodvecs)
2275 if (NULL==floodvecs->refproj0)
2277 snew(floodvecs->refproj0, floodvecs->neig);
2280 for (i=0; i<floodvecs->neig; i++)
2282 floodvecs->refproj0[i] = floodvecs->refproj[i];
2287 /* Call on MASTER only. Check whether the essential dynamics / flooding
2288 * datasets of the checkpoint file are consistent with the provided .edi file. */
2289 static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate)
2291 t_edpar *edi = NULL; /* points to a single edi data set */
2295 if (NULL == EDstate->nref || NULL == EDstate->nav)
2297 gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
2298 "start of a new simulation. If a simulation runs with/without ED constraints,\n"
2299 "it must also continue with/without ED constraints when checkpointing.\n"
2300 "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
2301 "from without a checkpoint.\n");
2308 /* Check number of atoms in the reference and average structures */
2309 if (EDstate->nref[edinum] != edi->sref.nr)
2311 gmx_fatal(FARGS, "The number of reference structure atoms in ED dataset #%d is\n"
2312 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2313 edinum+1, EDstate->nref[edinum], edi->sref.nr);
2315 if (EDstate->nav[edinum] != edi->sav.nr)
2317 gmx_fatal(FARGS, "The number of average structure atoms in ED dataset #%d is\n"
2318 "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
2319 edinum+1, EDstate->nav[edinum], edi->sav.nr);
2325 if (edinum != EDstate->nED)
2327 gmx_fatal(FARGS, "The number of essential dynamics / flooding datasets is not consistent.\n"
2328 "There are %d ED datasets in .cpt file, but %d in .edi file!\n"
2329 "Are you shure this is the correct .edi file?\n", EDstate->nED, edinum);
2334 /* The edsamstate struct stores the information we need to make the ED group
2335 * whole again after restarts from a checkpoint file. Here we do the following:
2336 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
2337 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
2338 * c) in any case, for subsequent checkpoint writing, we set the pointers in
2339 * edsamstate to the x_old arrays, which contain the correct PBC representation of
2340 * all ED structures at the last time step. */
2341 static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate)
2347 snew(EDstate->old_sref_p, EDstate->nED);
2348 snew(EDstate->old_sav_p , EDstate->nED);
2350 /* If we did not read in a .cpt file, these arrays are not yet allocated */
2351 if (!EDstate->bFromCpt)
2353 snew(EDstate->nref, EDstate->nED);
2354 snew(EDstate->nav , EDstate->nED);
2357 /* Loop over all ED/flooding data sets (usually only one, though) */
2359 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2361 /* We always need the last reference and average positions such that
2362 * in the next time step we can make the ED group whole again
2363 * if the atoms do not have the correct PBC representation */
2364 if (EDstate->bFromCpt)
2366 /* Copy the last whole positions of reference and average group from .cpt */
2367 for (i=0; i<edi->sref.nr; i++)
2369 copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
2371 for (i=0; i<edi->sav.nr ; i++)
2373 copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
2378 EDstate->nref[nr_edi-1] = edi->sref.nr;
2379 EDstate->nav [nr_edi-1] = edi->sav.nr;
2382 /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
2383 EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
2384 EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old ;
2386 edi = edi->next_edi;
2391 void init_edsam(gmx_mtop_t *mtop, /* global topology */
2392 t_inputrec *ir, /* input record */
2393 t_commrec *cr, /* communication record */
2394 gmx_edsam_t ed, /* contains all ED data */
2395 rvec x[], /* positions of the whole MD system */
2396 matrix box, /* the box */
2397 edsamstate_t *EDstate)
2399 t_edpar *edi = NULL; /* points to a single edi data set */
2400 int i,nr_edi,avindex;
2401 rvec *x_pbc = NULL; /* positions of the whole MD system with pbc removed */
2402 rvec *xfit=NULL, *xstart=NULL; /* dummy arrays to determine initial RMSDs */
2403 rvec fit_transvec; /* translation ... */
2404 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2407 if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
2409 gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
2412 GMX_MPE_LOG(ev_edsam_start);
2416 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2420 gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n"
2421 "flooding simulation. Please also provide the correct .edi file with -ei.\n");
2425 /* Needed for initializing radacc radius in do_edsam */
2428 /* The input file is read by the master and the edi structures are
2429 * initialized here. Input is stored in ed->edpar. Then the edi
2430 * structures are transferred to the other nodes */
2434 /* Read the whole edi file at once: */
2435 EDstate->nED = read_edi_file(ed,ed->edpar,mtop->natoms,cr);
2437 /* Make shure the checkpoint was produced in a run using this .edi file */
2438 if (EDstate->bFromCpt)
2440 crosscheck_edi_file_vs_checkpoint(ed, EDstate);
2442 init_edsamstate(ed, EDstate);
2444 /* Initialization for every ED/flooding dataset. Flooding uses one edi dataset per
2445 * flooding vector, Essential dynamics can be applied to more than one structure
2446 * as well, but will be done in the order given in the edi file, so
2447 * expect different results for different order of edi file concatenation! */
2451 init_edi(mtop,ir,cr,ed,edi);
2453 /* Init flooding parameters if needed */
2454 init_flood(edi,ed,ir->delta_t,cr,!EDstate->bFromCpt);
2460 /* The master does the work here. The other nodes get the positions
2461 * not before dd_partition_system which is called after init_edsam */
2464 /* Remove pbc, make molecule whole.
2465 * When ir->bContinuation=TRUE this has already been done, but ok.
2467 snew(x_pbc,mtop->natoms);
2468 m_rveccopy(mtop->natoms,x,x_pbc);
2469 do_pbc_first_mtop(NULL,ir->ePBC,box,mtop,x_pbc);
2471 /* Reset pointer to first ED data set which contains the actual ED data */
2474 /* Loop over all ED/flooding data sets (usually only one, though) */
2475 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2477 /* Extract the initial reference and average positions. When starting
2478 * from .cpt, these have already been read into sref.x_old
2479 * in init_edsamstate() */
2480 if (!EDstate->bFromCpt)
2482 /* If this is the first run (i.e. no checkpoint present) we assume
2483 * that the starting positions give us the correct PBC representation */
2484 for (i=0; i < edi->sref.nr; i++)
2486 copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
2489 for (i=0; i < edi->sav.nr; i++)
2491 copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
2495 /* Now we have the PBC-correct start positions of the reference and
2496 average structure. We copy that over to dummy arrays on which we
2497 can apply fitting to print out the RMSD. We srenew the memory since
2498 the size of the buffers is likely different for every ED dataset */
2499 srenew(xfit , edi->sref.nr );
2500 srenew(xstart, edi->sav.nr );
2501 copy_rvecn(edi->sref.x_old, xfit, 0, edi->sref.nr);
2502 copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
2504 /* Make the fit to the REFERENCE structure, get translation and rotation */
2505 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2507 /* Output how well we fit to the reference at the start */
2508 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2509 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm (dataset #%d)\n",
2510 rmsd_from_structure(xfit, &edi->sref), nr_edi);
2512 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2513 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2515 /* calculate initial projections */
2516 project(xstart, edi);
2518 /* For the target and origin structure both a reference (fit) and an
2519 * average structure can be provided in make_edi. If both structures
2520 * are the same, make_edi only stores one of them in the .edi file.
2521 * If they differ, first the fit and then the average structure is stored
2522 * in star (or sor), thus the number of entries in star/sor is
2523 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2524 * the size of the average group. */
2526 /* process target structure, if required */
2527 if (edi->star.nr > 0)
2529 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2531 /* get translation & rotation for fit of target structure to reference structure */
2532 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2534 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2535 if (edi->star.nr == edi->sav.nr)
2539 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2541 /* The last sav.nr indices of the target structure correspond to
2542 * the average structure, which must be projected */
2543 avindex = edi->star.nr - edi->sav.nr;
2545 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon, cr);
2549 rad_project(edi, xstart, &edi->vecs.radcon, cr);
2552 /* process structure that will serve as origin of expansion circle */
2553 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2555 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2558 if (edi->sori.nr > 0)
2560 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2562 /* fit this structure to reference structure */
2563 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2565 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2566 if (edi->sori.nr == edi->sav.nr)
2570 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2572 /* For the projection, we need the last sav.nr indices of sori */
2573 avindex = edi->sori.nr - edi->sav.nr;
2576 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc, cr);
2577 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix, cr);
2578 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2580 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2581 /* Set center of flooding potential to the ORIGIN structure */
2582 rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs, cr);
2583 /* We already know that no (moving) reference position was provided,
2584 * therefore we can overwrite refproj[0]*/
2585 copyEvecReference(&edi->flood.vecs);
2588 else /* No origin structure given */
2590 rad_project(edi, xstart, &edi->vecs.radacc, cr);
2591 rad_project(edi, xstart, &edi->vecs.radfix, cr);
2592 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2594 if (edi->flood.bHarmonic)
2596 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2597 for (i=0; i<edi->flood.vecs.neig; i++)
2599 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2604 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2605 /* Set center of flooding potential to the center of the covariance matrix,
2606 * i.e. the average structure, i.e. zero in the projected system */
2607 for (i=0; i<edi->flood.vecs.neig; i++)
2609 edi->flood.vecs.refproj[i] = 0.0;
2614 /* For convenience, output the center of the flooding potential for the eigenvectors */
2615 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2617 for (i=0; i<edi->flood.vecs.neig; i++)
2619 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
2620 if (edi->flood.bHarmonic)
2622 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2624 fprintf(stdout, "\n");
2628 /* set starting projections for linsam */
2629 rad_project(edi, xstart, &edi->vecs.linacc, cr);
2630 rad_project(edi, xstart, &edi->vecs.linfix, cr);
2632 /* Output to file, set the step to -1 so that write_edo knows it was called from init_edsam */
2633 if (ed->edo && !(EDstate->bFromCpt))
2635 write_edo(nr_edi, edi, ed, -1, 0);
2638 /* Prepare for the next edi data set: */
2641 /* Cleaning up on the master node: */
2646 } /* end of MASTER only section */
2650 /* First let everybody know how many ED data sets to expect */
2651 gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr);
2652 /* Broadcast the essential dynamics / flooding data to all nodes */
2653 broadcast_ed_data(cr, ed, EDstate->nED);
2657 /* In the single-CPU case, point the local atom numbers pointers to the global
2658 * one, so that we can use the same notation in serial and parallel case: */
2660 /* Loop over all ED data sets (usually only one, though) */
2662 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2664 edi->sref.anrs_loc = edi->sref.anrs;
2665 edi->sav.anrs_loc = edi->sav.anrs;
2666 edi->star.anrs_loc = edi->star.anrs;
2667 edi->sori.anrs_loc = edi->sori.anrs;
2668 /* For the same reason as above, make a dummy c_ind array: */
2669 snew(edi->sav.c_ind, edi->sav.nr);
2670 /* Initialize the array */
2671 for (i=0; i<edi->sav.nr; i++)
2673 edi->sav.c_ind[i] = i;
2675 /* In the general case we will need a different-sized array for the reference indices: */
2678 snew(edi->sref.c_ind, edi->sref.nr);
2679 for (i=0; i<edi->sref.nr; i++)
2681 edi->sref.c_ind[i] = i;
2684 /* Point to the very same array in case of other structures: */
2685 edi->star.c_ind = edi->sav.c_ind;
2686 edi->sori.c_ind = edi->sav.c_ind;
2687 /* In the serial case, the local number of atoms is the global one: */
2688 edi->sref.nr_loc = edi->sref.nr;
2689 edi->sav.nr_loc = edi->sav.nr;
2690 edi->star.nr_loc = edi->star.nr;
2691 edi->sori.nr_loc = edi->sori.nr;
2693 /* An on we go to the next edi dataset */
2698 /* Allocate space for ED buffer variables */
2699 /* Again, loop over ED data sets */
2701 for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
2703 /* Allocate space for ED buffer */
2705 snew(edi->buf->do_edsam, 1);
2707 /* Space for collective ED buffer variables */
2709 /* Collective positions of atoms with the average indices */
2710 snew(edi->buf->do_edsam->xcoll , edi->sav.nr);
2711 snew(edi->buf->do_edsam->shifts_xcoll , edi->sav.nr); /* buffer for xcoll shifts */
2712 snew(edi->buf->do_edsam->extra_shifts_xcoll , edi->sav.nr);
2713 /* Collective positions of atoms with the reference indices */
2716 snew(edi->buf->do_edsam->xc_ref , edi->sref.nr);
2717 snew(edi->buf->do_edsam->shifts_xc_ref , edi->sref.nr); /* To store the shifts in */
2718 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2721 /* Get memory for flooding forces */
2722 snew(edi->flood.forces_cartesian , edi->sav.nr);
2725 /* Dump it all into one file per process */
2726 dump_edi(edi, cr, nr_edi);
2729 /* An on we go to the next edi dataset */
2733 /* Flush the edo file so that the user can check some things
2734 * when the simulation has started */
2740 GMX_MPE_LOG(ev_edsam_finish);
2744 void do_edsam(t_inputrec *ir,
2745 gmx_large_int_t step,
2748 rvec xs[], /* The local current positions on this processor */
2749 rvec v[], /* The velocities */
2753 int i,edinr,iupdate=500;
2754 matrix rotmat; /* rotation matrix */
2755 rvec transvec; /* translation vector */
2756 rvec dv,dx,x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
2757 real dt_1; /* 1/dt */
2758 struct t_do_edsam *buf;
2760 real rmsdev=-1; /* RMSD from reference structure prior to applying the constraints */
2761 gmx_bool bSuppress=FALSE; /* Write .edo file on master? */
2764 /* Check if ED sampling has to be performed */
2765 if ( ed->eEDtype==eEDnone )
2770 /* Suppress output on first call of do_edsam if
2771 * two-step sd2 integrator is used */
2772 if ( (ir->eI==eiSD2) && (v != NULL) )
2777 dt_1 = 1.0/ir->delta_t;
2779 /* Loop over all ED datasets (usually one) */
2785 if (edi->bNeedDoEdsam)
2788 buf=edi->buf->do_edsam;
2792 /* initialise radacc radius for slope criterion */
2793 buf->oldrad=calc_radius(&edi->vecs.radacc);
2796 /* Copy the positions into buf->xc* arrays and after ED
2797 * feed back corrections to the official positions */
2799 /* Broadcast the ED positions such that every node has all of them
2800 * Every node contributes its local positions xs and stores it in
2801 * the collective buf->xcoll array. Note that for edinr > 1
2802 * xs could already have been modified by an earlier ED */
2804 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
2805 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
2808 dump_xcoll(edi, buf, cr, step);
2810 /* Only assembly reference positions if their indices differ from the average ones */
2813 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
2814 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
2817 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
2818 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
2819 * set bUpdateShifts=TRUE in the parallel case. */
2820 buf->bUpdateShifts = FALSE;
2822 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
2823 * as well as the indices in edi->sav.anrs */
2825 /* Fit the reference indices to the reference structure */
2828 fit_to_reference(buf->xcoll , transvec, rotmat, edi);
2832 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
2835 /* Now apply the translation and rotation to the ED structure */
2836 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
2838 /* Find out how well we fit to the reference (just for output steps) */
2839 if (do_per_step(step,edi->outfrq) && MASTER(cr))
2843 /* Indices of reference and average structures are identical,
2844 * thus we can calculate the rmsd to SREF using xcoll */
2845 rmsdev = rmsd_from_structure(buf->xcoll,&edi->sref);
2849 /* We have to translate & rotate the reference atoms first */
2850 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
2851 rmsdev = rmsd_from_structure(buf->xc_ref,&edi->sref);
2855 /* update radsam references, when required */
2856 if (do_per_step(step,edi->maxedsteps) && step >= edi->presteps)
2858 project(buf->xcoll, edi);
2859 rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
2860 rad_project(edi, buf->xcoll, &edi->vecs.radfix, cr);
2864 /* update radacc references, when required */
2865 if (do_per_step(step,iupdate) && step >= edi->presteps)
2867 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
2868 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
2870 project(buf->xcoll, edi);
2871 rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
2876 buf->oldrad = edi->vecs.radacc.radius;
2880 /* apply the constraints */
2881 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
2883 /* ED constraints should be applied already in the first MD step
2884 * (which is step 0), therefore we pass step+1 to the routine */
2885 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step, cr);
2888 /* write to edo, when required */
2889 if (do_per_step(step,edi->outfrq))
2891 project(buf->xcoll, edi);
2892 if (MASTER(cr) && !bSuppress)
2894 write_edo(edinr, edi, ed, step, rmsdev);
2898 /* Copy back the positions unless monitoring only */
2899 if (ed_constraints(ed->eEDtype, edi))
2901 /* remove fitting */
2902 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
2904 /* Copy the ED corrected positions into the coordinate array */
2905 /* Each node copies its local part. In the serial case, nat_loc is the
2906 * total number of ED atoms */
2907 for (i=0; i<edi->sav.nr_loc; i++)
2909 /* Unshift local ED coordinate and store in x_unsh */
2910 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
2911 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
2913 /* dx is the ED correction to the positions: */
2914 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
2918 /* dv is the ED correction to the velocity: */
2919 svmul(dt_1, dx, dv);
2920 /* apply the velocity correction: */
2921 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
2923 /* Finally apply the position correction due to ED: */
2924 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
2927 } /* END of if (edi->bNeedDoEdsam) */
2929 /* Prepare for the next ED dataset */
2930 edi = edi->next_edi;
2932 } /* END of loop over ED datasets */
2936 GMX_MPE_LOG(ev_edsam_finish);