1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
56 #include "mtop_util.h"
58 #include "mpelogging.h"
60 #include "groupcoord.h"
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)); }
69 /* enum to identify the type of ED: none, normal ED, flooding */
70 enum {eEDnone, eEDedsam, eEDflood, eEDnr};
72 /* enum to identify operations on reference, average, origin, target structures */
73 enum {eedREF, eedAV, eedORI, eedTAR, eedNR};
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;
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. */
106 gmx_bool bHarmonic; /* Use flooding for harmonic restraint on
108 gmx_bool bConstForce; /* Do not calculate a flooding potential,
109 instead flood with a constant force */
119 rvec *forces_cartesian;
120 t_eigvec vecs; /* use flooding for these */
124 /* This type is for the average, reference, target, and origin structure */
125 typedef struct gmx_edx
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) */
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 */
158 /* all gmx_edx datasets are copied to all nodes in the parallel case */
159 struct gmx_edx sref; /* reference positions, to these fitting
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 */
167 t_edvecs vecs; /* eigenvectors */
168 real slope; /* minimal slope in acceptance radexp */
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 */
178 typedef struct gmx_edsam
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 */
186 gmx_bool bStartFromCpt;
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
210 /* definition of ED buffer structure */
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;
220 /* Function declarations */
221 static void fit_to_reference(rvec *xcoll,rvec transvec,matrix rotmat,t_edpar *edi);
223 static void translate_and_rotate(rvec *x,int nat,rvec transvec,matrix rotmat);
224 /* End function declarations */
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
231 static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec)
237 for (i=0; i<edi->sav.nr; i++)
238 proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
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)
252 /* Subtract average positions */
253 for (i = 0; i < edi->sav.nr; i++)
254 rvec_dec(x[i], edi->sav.x[i]);
256 for (i = 0; i < vec->neig; i++)
258 vec->refproj[i] = projectx(edi,x,vec->vec[i]);
259 rad += pow((vec->refproj[i]-vec->xproj[i]),2);
261 vec->radius=sqrt(rad);
263 /* Add average positions */
264 for (i = 0; i < edi->sav.nr; i++)
265 rvec_inc(x[i], edi->sav.x[i]);
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
272 static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */
273 t_eigvec *vec, /* The eigenvectors */
279 if (!vec->neig) return;
281 /* Subtract average positions */
282 for (i=0; i<edi->sav.nr; i++)
283 rvec_dec(x[i], edi->sav.x[i]);
285 for (i=0; i<vec->neig; i++)
286 vec->xproj[i] = projectx(edi, x, vec->vec[i]);
288 /* Add average positions */
289 for (i=0; i<edi->sav.nr; i++)
290 rvec_inc(x[i], edi->sav.x[i]);
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 */
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);
309 static real calc_radius(t_eigvec *vec)
315 for (i=0; i<vec->neig; i++)
316 rad += pow((vec->refproj[i]-vec->xproj[i]),2);
318 return rad=sqrt(rad);
324 static void dump_xcoll(t_edpar *edi, struct t_do_edsam *buf, t_commrec *cr,
331 ivec *shifts, *eshifts;
338 shifts = buf->shifts_xcoll;
339 eshifts = buf->extra_shifts_xcoll;
341 sprintf(fn, "xcolldump_step%d.txt", step);
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",
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]);
356 static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[])
361 fprintf(out, "#%s positions:\n%d\n", name, s->nr);
365 fprintf(out, "#index, x, y, z");
367 fprintf(out, ", sqrt(m)");
368 for (i=0; i<s->nr; i++)
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]);
372 fprintf(out,"%9.3f",s->sqrtm[i]);
379 static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev,
380 const char name[], int length)
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++)
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]);
398 static void dump_edi(t_edpar *edpars, t_commrec *cr, int nr_edi)
404 sprintf(fn, "EDdump_node%d_edi%d", cr->nodeid, nr_edi);
405 out = ffopen(fn, "w");
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);
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" );
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);
429 /* Dump flooding eigenvectors */
430 dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING" , edpars->sav.nr);
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 );
442 static void dump_rotmat(FILE* out,matrix rotmat)
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]);
451 static void dump_rvec(FILE *out, int dim, rvec *x)
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]);
462 static void dump_mat(FILE* out, int dim, double** mat)
467 fprintf(out,"MATRIX:\n");
471 fprintf(out,"%f ",mat[i][j]);
483 static void do_edfit(int natoms,rvec *xp,rvec *x,matrix R,t_edpar *edi)
485 /* this is a copy of do_fit with some modifications */
492 struct t_do_edfit *loc;
495 if(edi->buf->do_edfit != NULL)
500 snew(edi->buf->do_edfit,1);
502 loc = edi->buf->do_edfit;
506 snew(loc->omega,2*DIM);
508 for(i=0; i<2*DIM; i++)
510 snew(loc->omega[i],2*DIM);
511 snew(loc->om[i],2*DIM);
525 /* calculate the matrix U */
527 for(n=0;(n<natoms);n++)
529 for(c=0; (c<DIM); c++)
532 for(r=0; (r<DIM); r++)
540 /* construct loc->omega */
541 /* loc->omega is symmetric -> loc->omega==loc->omega' */
546 loc->omega[r][c]=u[r-3][c];
547 loc->omega[c][r]=u[r-3][c];
555 /* determine h and k */
559 dump_mat(stderr,2*DIM,loc->omega);
561 fprintf(stderr,"d[%d] = %f\n",i,d[i]);
564 jacobi(loc->omega,6,d,loc->om,&irot);
567 fprintf(stderr,"IROT=0\n");
569 index=0; /* For the compiler only */
583 vh[j][i]=M_SQRT2*loc->om[i][index];
584 vk[j][i]=M_SQRT2*loc->om[i+DIM][index];
591 R[c][r]=vk[0][r]*vh[0][c]+
597 R[c][r]=vk[0][r]*vh[0][c]+
603 static void rmfit(int nat, rvec *xcoll, rvec transvec, matrix rotmat)
610 * The inverse rotation is described by the transposed rotation matrix */
611 transpose(rotmat,tmat);
612 rotate_x(xcoll, nat, tmat);
614 /* Remove translation */
615 vec[XX]=-transvec[XX];
616 vec[YY]=-transvec[YY];
617 vec[ZZ]=-transvec[ZZ];
618 translate_x(xcoll, nat, vec);
622 /**********************************************************************************
623 ******************** FLOODING ****************************************************
624 **********************************************************************************
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).
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.
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.
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.
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.
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.
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.
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.
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)
667 get_flood_energies(real Vfl[],int nnames);
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.
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
676 static void write_edo_flood(t_edpar *edi, FILE *fp, gmx_large_int_t step)
680 gmx_bool bOutputRef=FALSE;
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);
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)
693 for (i = 0; i < edi->flood.vecs.neig; i++)
695 if (edi->flood.vecs.refprojslope[i] != 0.0)
700 fprintf(fp, "Ref. projs.: ");
701 for (i = 0; i < edi->flood.vecs.neig; i++)
703 fprintf(fp, "%12.5e ", edi->flood.vecs.refproj[i]);
708 fprintf(fp,"FL_FORCES: ");
710 for (i=0; i<edi->flood.vecs.neig; i++)
711 fprintf(fp," %12.5e",edi->flood.vecs.fproj[i]);
717 /* From flood.xproj compute the Vfl(x) at this point */
718 static real flood_energy(t_edpar *edi, gmx_large_int_t step)
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]
725 Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2})
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)
738 for (i=0; i<edi->flood.vecs.neig; i++)
740 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
745 /* Compute sum which will be the exponent of the exponential */
746 for (i=0; i<edi->flood.vecs.neig; i++)
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]);
752 /* Compute the Gauss function*/
753 if (edi->flood.bHarmonic)
755 Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */
759 Vfl = edi->flood.Efl!=0 ? edi->flood.Efl*exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) :0;
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)
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 */
773 real energy=edi->flood.Vfl;
776 if (edi->flood.bHarmonic)
777 for (i=0; i<edi->flood.vecs.neig; i++)
779 edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
782 for (i=0; i<edi->flood.vecs.neig; i++)
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;
790 /* Raise forces from subspace into cartesian space */
791 static void flood_blowup(t_edpar *edi, rvec *forces_cart)
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.
799 for every atom we add all the contributions to this atom from all the different eigenvectors.
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
810 forces_sub = edi->flood.vecs.fproj;
813 /* Calculate the cartesian forces for the local atoms */
815 /* Clear forces first */
816 for (j=0; j<edi->sav.nr_loc; j++)
817 clear_rvec(forces_cart[j]);
819 /* Now compute atomwise */
820 for (j=0; j<edi->sav.nr_loc; j++)
822 /* Compute forces_cart[edi->sav.anrs[j]] */
823 for (eig=0; eig<edi->flood.vecs.neig; eig++)
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);
834 /* Update the values of Efl, deltaF depending on tau and Vfl */
835 static void update_adaption(t_edpar *edi)
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,
841 if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001)
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)
848 edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
853 static void do_single_flood(
858 gmx_large_int_t step,
861 gmx_bool bNS) /* Are we in a neighbor searching step? */
864 matrix rotmat; /* rotation matrix */
865 matrix tmat; /* inverse rotation */
866 rvec transvec; /* translation vector */
867 struct t_do_edsam *buf;
870 buf=edi->buf->do_edsam;
873 /* Broadcast the positions of the AVERAGE structure such that they are known on
874 * every processor. Each node contributes its local positions x and stores them in
875 * the collective ED array buf->xcoll */
876 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
877 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
879 /* Only assembly REFERENCE positions if their indices differ from the average ones */
881 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
882 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
884 /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
885 * We do not need to update the shifts until the next NS step */
886 buf->bUpdateShifts = FALSE;
888 /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll,
889 * as well as the indices in edi->sav.anrs */
891 /* Fit the reference indices to the reference structure */
893 fit_to_reference(buf->xcoll , transvec, rotmat, edi);
895 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
897 /* Now apply the translation and rotation to the ED structure */
898 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
900 /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */
901 project_to_eigvectors(buf->xcoll,&edi->flood.vecs,edi);
903 if (FALSE == edi->flood.bConstForce)
905 /* Compute Vfl(x) from flood.xproj */
906 edi->flood.Vfl = flood_energy(edi, step);
908 update_adaption(edi);
910 /* Compute the flooding forces */
914 /* Translate them into cartesian positions */
915 flood_blowup(edi, edi->flood.forces_cartesian);
917 /* Rotate forces back so that they correspond to the given structure and not to the fitted one */
918 /* Each node rotates back its local forces */
919 transpose(rotmat,tmat);
920 rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat);
922 /* Finally add forces to the main force variable */
923 for (i=0; i<edi->sav.nr_loc; i++)
924 rvec_inc(force[edi->sav.anrs_loc[i]],edi->flood.forces_cartesian[i]);
926 /* Output is written by the master process */
927 if (do_per_step(step,edi->outfrq) && MASTER(cr))
928 write_edo_flood(edi,edo,step);
932 /* Main flooding routine, called from do_force */
933 extern void do_flood(
934 FILE *log, /* md.log file */
935 t_commrec *cr, /* Communication record */
936 rvec x[], /* Positions on the local processor */
937 rvec force[], /* forcefield forces, to these the flooding forces are added */
938 gmx_edsam_t ed, /* ed data structure contains all ED and flooding datasets */
939 matrix box, /* the box */
940 gmx_large_int_t step, /* The relative time step since ir->init_step is already subtracted */
941 gmx_bool bNS) /* Are we in a neighbor searching step? */
946 if (ed->eEDtype != eEDflood)
952 /* Call flooding for one matrix */
953 if (edi->flood.vecs.neig)
954 do_single_flood(ed->edo,x,force,edi,step,box,cr,bNS);
960 /* Called by init_edi, configure some flooding related variables and structures,
961 * print headers to output files */
962 static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt, t_commrec *cr)
967 edi->flood.Efl = edi->flood.constEfl;
971 if (edi->flood.vecs.neig)
973 /* If in any of the datasets we find a flooding vector, flooding is turned on */
974 ed->eEDtype = eEDflood;
976 fprintf(stderr,"ED: Flooding of matrix %d is switched on.\n", edi->flood.flood_id);
978 if (edi->flood.bConstForce)
980 /* We have used stpsz as a vehicle to carry the fproj values for constant
981 * force flooding. Now we copy that to flood.vecs.fproj. Note that
982 * in const force flooding, fproj is never changed. */
983 for (i=0; i<edi->flood.vecs.neig; i++)
985 edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i];
987 fprintf(stderr, "ED: applying on eigenvector %d a constant force of %g\n",
988 edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
991 fprintf(ed->edo,"FL_HEADER: Flooding of matrix %d is switched on! The flooding output will have the following format:\n",
992 edi->flood.flood_id);
993 fprintf(ed->edo,"FL_HEADER: Step Efl Vfl deltaF\n");
999 /*********** Energy book keeping ******/
1000 static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */
1009 srenew(names,count);
1010 sprintf(buf,"Vfl_%d",count);
1011 names[count-1]=strdup(buf);
1012 actual=actual->next_edi;
1019 static void get_flood_energies(t_edpar *edi, real Vfl[],int nnames)
1021 /*fl has to be big enough to capture nnames-many entries*/
1030 Vfl[count-1]=actual->flood.Vfl;
1031 actual=actual->next_edi;
1034 if (nnames!=count-1)
1035 gmx_fatal(FARGS,"Number of energies is not consistent with t_edi structure");
1037 /************* END of FLOODING IMPLEMENTATION ****************************/
1041 gmx_edsam_t ed_open(int nfile,const t_filenm fnm[],unsigned long Flags,t_commrec *cr)
1046 /* Allocate space for the ED data structure */
1049 /* We want to perform ED (this switch might later be upgraded to eEDflood) */
1050 ed->eEDtype = eEDedsam;
1054 /* Open .edi input file: */
1055 ed->edinam=ftp2fn(efEDI,nfile,fnm);
1056 /* The master opens the .edo output file */
1057 fprintf(stderr,"ED sampling will be performed!\n");
1058 ed->edonam = ftp2fn(efEDO,nfile,fnm);
1059 ed->edo = gmx_fio_fopen(ed->edonam,(Flags & MD_APPENDFILES)? "a+" : "w+");
1060 ed->bStartFromCpt = Flags & MD_STARTFROMCPT;
1066 /* Broadcasts the structure data */
1067 static void bc_ed_positions(t_commrec *cr, struct gmx_edx *s, int stype)
1069 snew_bc(cr, s->anrs, s->nr ); /* Index numbers */
1070 snew_bc(cr, s->x , s->nr ); /* Positions */
1071 nblock_bc(cr, s->nr, s->anrs );
1072 nblock_bc(cr, s->nr, s->x );
1074 /* For the average & reference structures we need an array for the collective indices,
1075 * and we need to broadcast the masses as well */
1076 if (stype == eedAV || stype == eedREF)
1078 /* We need these additional variables in the parallel case: */
1079 snew(s->c_ind , s->nr ); /* Collective indices */
1080 /* Local atom indices get assigned in dd_make_local_group_indices.
1081 * There, also memory is allocated */
1082 s->nalloc_loc = 0; /* allocation size of s->anrs_loc */
1083 snew_bc(cr, s->x_old, s->nr); /* To be able to always make the ED molecule whole, ... */
1084 nblock_bc(cr, s->nr, s->x_old); /* ... keep track of shift changes with the help of old coords */
1087 /* broadcast masses for the reference structure (for mass-weighted fitting) */
1088 if (stype == eedREF)
1090 snew_bc(cr, s->m, s->nr);
1091 nblock_bc(cr, s->nr, s->m);
1094 /* For the average structure we might need the masses for mass-weighting */
1097 snew_bc(cr, s->sqrtm, s->nr);
1098 nblock_bc(cr, s->nr, s->sqrtm);
1099 snew_bc(cr, s->m, s->nr);
1100 nblock_bc(cr, s->nr, s->m);
1105 /* Broadcasts the eigenvector data */
1106 static void bc_ed_vecs(t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
1110 snew_bc(cr, ev->ieig , ev->neig); /* index numbers of eigenvector */
1111 snew_bc(cr, ev->stpsz , ev->neig); /* stepsizes per eigenvector */
1112 snew_bc(cr, ev->xproj , ev->neig); /* instantaneous x projection */
1113 snew_bc(cr, ev->fproj , ev->neig); /* instantaneous f projection */
1114 snew_bc(cr, ev->refproj, ev->neig); /* starting or target projection */
1116 nblock_bc(cr, ev->neig, ev->ieig );
1117 nblock_bc(cr, ev->neig, ev->stpsz );
1118 nblock_bc(cr, ev->neig, ev->xproj );
1119 nblock_bc(cr, ev->neig, ev->fproj );
1120 nblock_bc(cr, ev->neig, ev->refproj);
1122 snew_bc(cr, ev->vec, ev->neig); /* Eigenvector components */
1123 for (i=0; i<ev->neig; i++)
1125 snew_bc(cr, ev->vec[i], length);
1126 nblock_bc(cr, length, ev->vec[i]);
1129 /* For harmonic restraints the reference projections can change with time */
1132 snew_bc(cr, ev->refproj0 , ev->neig);
1133 snew_bc(cr, ev->refprojslope, ev->neig);
1134 nblock_bc(cr, ev->neig, ev->refproj0 );
1135 nblock_bc(cr, ev->neig, ev->refprojslope);
1140 /* Broadcasts the ED / flooding data to other nodes
1141 * and allocates memory where needed */
1142 static void broadcast_ed_data(t_commrec *cr, gmx_edsam_t ed, int numedis)
1148 /* Master lets the other nodes know if its ED only or also flooding */
1149 gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr);
1151 snew_bc(cr, ed->edpar,1);
1152 /* Now transfer the ED data set(s) */
1154 for (nr=0; nr<numedis; nr++)
1156 /* Broadcast a single ED data set */
1159 /* Broadcast positions */
1160 bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */
1161 bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */
1162 bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */
1163 bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
1165 /* Broadcast eigenvectors */
1166 bc_ed_vecs(cr, &edi->vecs.mon , edi->sav.nr, FALSE);
1167 bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
1168 bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
1169 bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
1170 bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
1171 bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
1172 /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
1173 bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
1175 /* Set the pointer to the next ED dataset */
1178 snew_bc(cr, edi->next_edi, 1);
1179 edi = edi->next_edi;
1185 /* init-routine called for every *.edi-cycle, initialises t_edpar structure */
1186 static void init_edi(gmx_mtop_t *mtop,t_inputrec *ir,
1187 t_commrec *cr,gmx_edsam_t ed,t_edpar *edi)
1190 real totalmass = 0.0;
1194 /* NOTE Init_edi is executed on the master process only
1195 * The initialized data sets are then transmitted to the
1196 * other nodes in broadcast_ed_data */
1198 edi->bNeedDoEdsam = edi->vecs.mon.neig
1199 || edi->vecs.linfix.neig
1200 || edi->vecs.linacc.neig
1201 || edi->vecs.radfix.neig
1202 || edi->vecs.radacc.neig
1203 || edi->vecs.radcon.neig;
1205 /* evaluate masses (reference structure) */
1206 snew(edi->sref.m, edi->sref.nr);
1207 for (i = 0; i < edi->sref.nr; i++)
1211 gmx_mtop_atomnr_to_atom(mtop,edi->sref.anrs[i],&atom);
1212 edi->sref.m[i] = atom->m;
1216 edi->sref.m[i] = 1.0;
1219 /* Check that every m > 0. Bad things will happen otherwise. */
1220 if (edi->sref.m[i] <= 0.0)
1222 gmx_fatal(FARGS, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n"
1223 "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n"
1224 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1225 "atoms from the reference structure by creating a proper index group.\n",
1226 i, edi->sref.anrs[i]+1, edi->sref.m[i]);
1229 totalmass += edi->sref.m[i];
1231 edi->sref.mtot = totalmass;
1233 /* Masses m and sqrt(m) for the average structure. Note that m
1234 * is needed if forces have to be evaluated in do_edsam */
1235 snew(edi->sav.sqrtm, edi->sav.nr );
1236 snew(edi->sav.m , edi->sav.nr );
1237 for (i = 0; i < edi->sav.nr; i++)
1239 gmx_mtop_atomnr_to_atom(mtop,edi->sav.anrs[i],&atom);
1240 edi->sav.m[i] = atom->m;
1243 edi->sav.sqrtm[i] = sqrt(atom->m);
1247 edi->sav.sqrtm[i] = 1.0;
1250 /* Check that every m > 0. Bad things will happen otherwise. */
1251 if (edi->sav.sqrtm[i] <= 0.0)
1253 gmx_fatal(FARGS, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n"
1254 "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n"
1255 "Either make the covariance analysis non-mass-weighted, or exclude massless\n"
1256 "atoms from the average structure by creating a proper index group.\n",
1257 i, edi->sav.anrs[i]+1, atom->m);
1261 /* put reference structure in origin */
1262 get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com);
1266 translate_x(edi->sref.x, edi->sref.nr, com);
1268 /* Init ED buffer */
1273 static void check(const char *line, const char *label)
1275 if (!strstr(line,label))
1276 gmx_fatal(FARGS,"Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s",label,line);
1280 static int read_checked_edint(FILE *file,const char *label)
1282 char line[STRLEN+1];
1286 fgets2 (line,STRLEN,file);
1288 fgets2 (line,STRLEN,file);
1289 sscanf (line,"%d",&idum);
1294 static int read_edint(FILE *file,gmx_bool *bEOF)
1296 char line[STRLEN+1];
1301 eof=fgets2 (line,STRLEN,file);
1307 eof=fgets2 (line,STRLEN,file);
1313 sscanf (line,"%d",&idum);
1319 static real read_checked_edreal(FILE *file,const char *label)
1321 char line[STRLEN+1];
1325 fgets2 (line,STRLEN,file);
1327 fgets2 (line,STRLEN,file);
1328 sscanf (line,"%lf",&rdum);
1329 return (real) rdum; /* always read as double and convert to single */
1333 static void read_edx(FILE *file,int number,int *anrs,rvec *x)
1336 char line[STRLEN+1];
1340 for(i=0; i<number; i++)
1342 fgets2 (line,STRLEN,file);
1343 sscanf (line,"%d%lf%lf%lf",&anrs[i],&d[0],&d[1],&d[2]);
1344 anrs[i]--; /* we are reading FORTRAN indices */
1346 x[i][j]=d[j]; /* always read as double and convert to single */
1351 static void scan_edvec(FILE *in,int nr,rvec *vec)
1353 char line[STRLEN+1];
1358 for(i=0; (i < nr); i++)
1360 fgets2 (line,STRLEN,in);
1361 sscanf (line,"%le%le%le",&x,&y,&z);
1369 static void read_edvec(FILE *in,int nr,t_eigvec *tvec,gmx_bool bReadRefproj, gmx_bool *bHaveReference)
1372 double rdum,refproj_dum=0.0,refprojslope_dum=0.0;
1373 char line[STRLEN+1];
1376 tvec->neig=read_checked_edint(in,"NUMBER OF EIGENVECTORS");
1379 snew(tvec->ieig ,tvec->neig);
1380 snew(tvec->stpsz ,tvec->neig);
1381 snew(tvec->vec ,tvec->neig);
1382 snew(tvec->xproj ,tvec->neig);
1383 snew(tvec->fproj ,tvec->neig);
1384 snew(tvec->refproj,tvec->neig);
1387 snew(tvec->refproj0 ,tvec->neig);
1388 snew(tvec->refprojslope,tvec->neig);
1391 for(i=0; (i < tvec->neig); i++)
1393 fgets2 (line,STRLEN,in);
1394 if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
1396 nscan = sscanf(line,"%d%lf%lf%lf",&idum,&rdum,&refproj_dum,&refprojslope_dum);
1397 /* Zero out values which were not scanned */
1401 /* Every 4 values read, including reference position */
1402 *bHaveReference = TRUE;
1405 /* A reference position is provided */
1406 *bHaveReference = TRUE;
1407 /* No value for slope, set to 0 */
1408 refprojslope_dum = 0.0;
1411 /* No values for reference projection and slope, set to 0 */
1413 refprojslope_dum = 0.0;
1416 gmx_fatal(FARGS,"Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
1419 tvec->refproj[i]=refproj_dum;
1420 tvec->refproj0[i]=refproj_dum;
1421 tvec->refprojslope[i]=refprojslope_dum;
1423 else /* Normal flooding */
1425 nscan = sscanf(line,"%d%lf",&idum,&rdum);
1427 gmx_fatal(FARGS,"Expected 2 values for flooding vec: <nr> <stpsz>\n");
1430 tvec->stpsz[i]=rdum;
1431 } /* end of loop over eigenvectors */
1433 for(i=0; (i < tvec->neig); i++)
1435 snew(tvec->vec[i],nr);
1436 scan_edvec(in,nr,tvec->vec[i]);
1442 /* calls read_edvec for the vector groups, only for flooding there is an extra call */
1443 static void read_edvecs(FILE *in,int nr,t_edvecs *vecs)
1445 gmx_bool bHaveReference = FALSE;
1448 read_edvec(in, nr, &vecs->mon , FALSE, &bHaveReference);
1449 read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
1450 read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
1451 read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
1452 read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
1453 read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
1457 /* Check if the same atom indices are used for reference and average positions */
1458 static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav)
1463 /* If the number of atoms differs between the two structures,
1464 * they cannot be identical */
1465 if (sref.nr != sav.nr)
1468 /* Now that we know that both stuctures have the same number of atoms,
1469 * check if also the indices are identical */
1470 for (i=0; i < sav.nr; i++)
1472 if (sref.anrs[i] != sav.anrs[i])
1475 fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
1481 static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int edi_nr, t_commrec *cr)
1484 const int magic=670;
1487 /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
1488 gmx_bool bHaveReference = FALSE;
1491 /* the edi file is not free format, so expect problems if the input is corrupt. */
1493 /* check the magic number */
1494 readmagic=read_edint(in,&bEOF);
1495 /* Check whether we have reached the end of the input file */
1499 if (readmagic != magic)
1501 if (readmagic==666 || readmagic==667 || readmagic==668)
1502 gmx_fatal(FARGS,"Wrong magic number: Use newest version of make_edi to produce edi file");
1503 else if (readmagic == 669)
1506 gmx_fatal(FARGS,"Wrong magic number %d in %s",readmagic,ed->edinam);
1509 /* check the number of atoms */
1510 edi->nini=read_edint(in,&bEOF);
1511 if (edi->nini != nr_mdatoms)
1512 gmx_fatal(FARGS,"Nr of atoms in %s (%d) does not match nr of md atoms (%d)",
1513 ed->edinam,edi->nini,nr_mdatoms);
1515 /* Done checking. For the rest we blindly trust the input */
1516 edi->fitmas = read_checked_edint(in,"FITMAS");
1517 edi->pcamas = read_checked_edint(in,"ANALYSIS_MAS");
1518 edi->outfrq = read_checked_edint(in,"OUTFRQ");
1519 edi->maxedsteps = read_checked_edint(in,"MAXLEN");
1520 edi->slope = read_checked_edreal(in,"SLOPECRIT");
1522 edi->presteps = read_checked_edint(in,"PRESTEPS");
1523 edi->flood.deltaF0 = read_checked_edreal(in,"DELTA_F0");
1524 edi->flood.deltaF = read_checked_edreal(in,"INIT_DELTA_F");
1525 edi->flood.tau = read_checked_edreal(in,"TAU");
1526 edi->flood.constEfl = read_checked_edreal(in,"EFL_NULL");
1527 edi->flood.alpha2 = read_checked_edreal(in,"ALPHA2");
1528 edi->flood.kT = read_checked_edreal(in,"KT");
1529 edi->flood.bHarmonic = read_checked_edint(in,"HARMONIC");
1530 if (readmagic > 669)
1531 edi->flood.bConstForce = read_checked_edint(in,"CONST_FORCE_FLOODING");
1533 edi->flood.bConstForce = FALSE;
1534 edi->flood.flood_id = edi_nr;
1535 edi->sref.nr = read_checked_edint(in,"NREF");
1537 /* allocate space for reference positions and read them */
1538 snew(edi->sref.anrs,edi->sref.nr);
1539 snew(edi->sref.x ,edi->sref.nr);
1540 snew(edi->sref.x_old,edi->sref.nr);
1541 edi->sref.sqrtm =NULL;
1542 read_edx(in,edi->sref.nr,edi->sref.anrs,edi->sref.x);
1544 /* average positions. they define which atoms will be used for ED sampling */
1545 edi->sav.nr=read_checked_edint(in,"NAV");
1546 snew(edi->sav.anrs,edi->sav.nr);
1547 snew(edi->sav.x ,edi->sav.nr);
1548 snew(edi->sav.x_old,edi->sav.nr);
1549 read_edx(in,edi->sav.nr,edi->sav.anrs,edi->sav.x);
1551 /* Check if the same atom indices are used for reference and average positions */
1552 edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
1555 read_edvecs(in,edi->sav.nr,&edi->vecs);
1556 read_edvec(in,edi->sav.nr,&edi->flood.vecs,edi->flood.bHarmonic, &bHaveReference);
1558 /* target positions */
1559 edi->star.nr=read_edint(in,&bEOF);
1560 if (edi->star.nr > 0)
1562 snew(edi->star.anrs,edi->star.nr);
1563 snew(edi->star.x ,edi->star.nr);
1564 edi->star.sqrtm =NULL;
1565 read_edx(in,edi->star.nr,edi->star.anrs,edi->star.x);
1568 /* positions defining origin of expansion circle */
1569 edi->sori.nr=read_edint(in,&bEOF);
1570 if (edi->sori.nr > 0)
1574 /* Both an -ori structure and a at least one manual reference point have been
1575 * specified. That's ambiguous and probably not intentional. */
1576 gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
1577 " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
1579 snew(edi->sori.anrs,edi->sori.nr);
1580 snew(edi->sori.x ,edi->sori.nr);
1581 edi->sori.sqrtm =NULL;
1582 read_edx(in,edi->sori.nr,edi->sori.anrs,edi->sori.x);
1591 /* Read in the edi input file. Note that it may contain several ED data sets which were
1592 * achieved by concatenating multiple edi files. The standard case would be a single ED
1593 * data set, though. */
1594 static void read_edi_file(gmx_edsam_t ed, t_edpar *edi, int nr_mdatoms, t_commrec *cr)
1597 t_edpar *curr_edi,*last_edi;
1602 /* This routine is executed on the master only */
1604 /* Open the .edi parameter input file */
1605 in = gmx_fio_fopen(ed->edinam,"r");
1606 fprintf(stderr, "ED: Reading edi file %s\n", ed->edinam);
1608 /* Now read a sequence of ED input parameter sets from the edi file */
1611 while( read_edi(in, ed, curr_edi, nr_mdatoms, edi_nr, cr) )
1614 /* Make shure that the number of atoms in each dataset is the same as in the tpr file */
1615 if (edi->nini != nr_mdatoms)
1616 gmx_fatal(FARGS,"edi file %s (dataset #%d) was made for %d atoms, but the simulation contains %d atoms.",
1617 ed->edinam, edi_nr, edi->nini, nr_mdatoms);
1618 /* Since we arrived within this while loop we know that there is still another data set to be read in */
1619 /* We need to allocate space for the data: */
1621 /* Point the 'next_edi' entry to the next edi: */
1622 curr_edi->next_edi=edi_read;
1623 /* Keep the curr_edi pointer for the case that the next dataset is empty: */
1624 last_edi = curr_edi;
1625 /* Let's prepare to read in the next edi data set: */
1626 curr_edi = edi_read;
1629 gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", ed->edinam);
1631 /* Terminate the edi dataset list with a NULL pointer: */
1632 last_edi->next_edi = NULL;
1634 fprintf(stderr, "ED: Found %d ED dataset%s.\n", edi_nr, edi_nr>1? "s" : "");
1636 /* Close the .edi file again */
1641 struct t_fit_to_ref {
1642 rvec *xcopy; /* Working copy of the positions in fit_to_reference */
1645 /* Fit the current positions to the reference positions
1646 * Do not actually do the fit, just return rotation and translation.
1647 * Note that the COM of the reference structure was already put into
1648 * the origin by init_edi. */
1649 static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */
1650 rvec transvec, /* The translation vector */
1651 matrix rotmat, /* The rotation matrix */
1652 t_edpar *edi) /* Just needed for do_edfit */
1654 rvec com; /* center of mass */
1656 struct t_fit_to_ref *loc;
1659 GMX_MPE_LOG(ev_fit_to_reference_start);
1661 /* Allocate memory the first time this routine is called for each edi dataset */
1662 if (NULL == edi->buf->fit_to_ref)
1664 snew(edi->buf->fit_to_ref, 1);
1665 snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr);
1667 loc = edi->buf->fit_to_ref;
1669 /* We do not touch the original positions but work on a copy. */
1670 for (i=0; i<edi->sref.nr; i++)
1671 copy_rvec(xcoll[i], loc->xcopy[i]);
1673 /* Calculate the center of mass */
1674 get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
1676 transvec[XX] = -com[XX];
1677 transvec[YY] = -com[YY];
1678 transvec[ZZ] = -com[ZZ];
1680 /* Subtract the center of mass from the copy */
1681 translate_x(loc->xcopy, edi->sref.nr, transvec);
1683 /* Determine the rotation matrix */
1684 do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi);
1686 GMX_MPE_LOG(ev_fit_to_reference_finish);
1690 static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */
1691 int nat, /* How many positions are there? */
1692 rvec transvec, /* The translation vector */
1693 matrix rotmat) /* The rotation matrix */
1696 translate_x(x, nat, transvec);
1699 rotate_x(x, nat, rotmat);
1703 /* Gets the rms deviation of the positions to the structure s */
1704 /* fit_to_structure has to be called before calling this routine! */
1705 static real rmsd_from_structure(rvec *x, /* The positions under consideration */
1706 struct gmx_edx *s) /* The structure from which the rmsd shall be computed */
1712 for (i=0; i < s->nr; i++)
1713 rmsd += distance2(s->x[i], x[i]);
1715 rmsd /= (real) s->nr;
1722 void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed)
1727 if (ed->eEDtype != eEDnone)
1729 /* Loop over ED datasets (usually there is just one dataset, though) */
1733 /* Local atoms of the reference structure (for fitting), need only be assembled
1734 * if their indices differ from the average ones */
1736 dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
1737 &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
1739 /* Local atoms of the average structure (on these ED will be performed) */
1740 dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
1741 &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind);
1743 /* Indicate that the ED shift vectors for this structure need to be updated
1744 * at the next call to communicate_group_positions, since obviously we are in a NS step */
1745 edi->buf->do_edsam->bUpdateShifts = TRUE;
1747 /* Set the pointer to the next ED dataset (if any) */
1754 static inline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu)
1759 GMX_MPE_LOG(ev_unshift_start);
1767 xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1768 xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1769 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1772 xu[XX] = x[XX]-tx*box[XX][XX];
1773 xu[YY] = x[YY]-ty*box[YY][YY];
1774 xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
1777 GMX_MPE_LOG(ev_unshift_finish);
1781 static void do_linfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
1788 /* loop over linfix vectors */
1789 for (i=0; i<edi->vecs.linfix.neig; i++)
1791 /* calculate the projection */
1792 proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]);
1794 /* calculate the correction */
1795 add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj;
1797 /* apply the correction */
1798 add /= edi->sav.sqrtm[i];
1799 for (j=0; j<edi->sav.nr; j++)
1801 svmul(add, edi->vecs.linfix.vec[i][j], vec_dum);
1802 rvec_inc(xcoll[j], vec_dum);
1808 static void do_linacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1815 /* loop over linacc vectors */
1816 for (i=0; i<edi->vecs.linacc.neig; i++)
1818 /* calculate the projection */
1819 proj=projectx(edi, xcoll, edi->vecs.linacc.vec[i]);
1821 /* calculate the correction */
1823 if (edi->vecs.linacc.stpsz[i] > 0.0)
1825 if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
1826 add = edi->vecs.linacc.refproj[i] - proj;
1828 if (edi->vecs.linacc.stpsz[i] < 0.0)
1830 if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
1831 add = edi->vecs.linacc.refproj[i] - proj;
1834 /* apply the correction */
1835 add /= edi->sav.sqrtm[i];
1836 for (j=0; j<edi->sav.nr; j++)
1838 svmul(add, edi->vecs.linacc.vec[i][j], vec_dum);
1839 rvec_inc(xcoll[j], vec_dum);
1842 /* new positions will act as reference */
1843 edi->vecs.linacc.refproj[i] = proj + add;
1848 static void do_radfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
1851 real *proj, rad=0.0, ratio;
1855 if (edi->vecs.radfix.neig == 0)
1858 snew(proj, edi->vecs.radfix.neig);
1860 /* loop over radfix vectors */
1861 for (i=0; i<edi->vecs.radfix.neig; i++)
1863 /* calculate the projections, radius */
1864 proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]);
1865 rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2);
1869 ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0;
1870 edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0];
1872 /* loop over radfix vectors */
1873 for (i=0; i<edi->vecs.radfix.neig; i++)
1875 proj[i] -= edi->vecs.radfix.refproj[i];
1877 /* apply the correction */
1878 proj[i] /= edi->sav.sqrtm[i];
1880 for (j=0; j<edi->sav.nr; j++) {
1881 svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
1882 rvec_inc(xcoll[j], vec_dum);
1890 static void do_radacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1893 real *proj, rad=0.0, ratio=0.0;
1897 if (edi->vecs.radacc.neig == 0)
1900 snew(proj,edi->vecs.radacc.neig);
1902 /* loop over radacc vectors */
1903 for (i=0; i<edi->vecs.radacc.neig; i++)
1905 /* calculate the projections, radius */
1906 proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]);
1907 rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2);
1911 /* only correct when radius decreased */
1912 if (rad < edi->vecs.radacc.radius)
1914 ratio = edi->vecs.radacc.radius/rad - 1.0;
1915 rad = edi->vecs.radacc.radius;
1918 edi->vecs.radacc.radius = rad;
1920 /* loop over radacc vectors */
1921 for (i=0; i<edi->vecs.radacc.neig; i++)
1923 proj[i] -= edi->vecs.radacc.refproj[i];
1925 /* apply the correction */
1926 proj[i] /= edi->sav.sqrtm[i];
1928 for (j=0; j<edi->sav.nr; j++)
1930 svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum);
1931 rvec_inc(xcoll[j], vec_dum);
1938 struct t_do_radcon {
1942 static void do_radcon(rvec *xcoll, t_edpar *edi, t_commrec *cr)
1945 real rad=0.0, ratio=0.0;
1946 struct t_do_radcon *loc;
1951 if(edi->buf->do_radcon != NULL)
1954 loc = edi->buf->do_radcon;
1959 snew(edi->buf->do_radcon, 1);
1961 loc = edi->buf->do_radcon;
1963 if (edi->vecs.radcon.neig == 0)
1967 snew(loc->proj, edi->vecs.radcon.neig);
1969 /* loop over radcon vectors */
1970 for (i=0; i<edi->vecs.radcon.neig; i++)
1972 /* calculate the projections, radius */
1973 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
1974 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
1977 /* only correct when radius increased */
1978 if (rad > edi->vecs.radcon.radius)
1980 ratio = edi->vecs.radcon.radius/rad - 1.0;
1982 /* loop over radcon vectors */
1983 for (i=0; i<edi->vecs.radcon.neig; i++)
1985 /* apply the correction */
1986 loc->proj[i] -= edi->vecs.radcon.refproj[i];
1987 loc->proj[i] /= edi->sav.sqrtm[i];
1988 loc->proj[i] *= ratio;
1990 for (j=0; j<edi->sav.nr; j++)
1992 svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum);
1993 rvec_inc(xcoll[j], vec_dum);
1998 edi->vecs.radcon.radius = rad;
2000 if (rad != edi->vecs.radcon.radius)
2003 for (i=0; i<edi->vecs.radcon.neig; i++)
2005 /* calculate the projections, radius */
2006 loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]);
2007 rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2);
2014 static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step, t_commrec *cr)
2019 GMX_MPE_LOG(ev_ed_apply_cons_start);
2021 /* subtract the average positions */
2022 for (i=0; i<edi->sav.nr; i++)
2023 rvec_dec(xcoll[i], edi->sav.x[i]);
2025 /* apply the constraints */
2027 do_linfix(xcoll, edi, step, cr);
2028 do_linacc(xcoll, edi, cr);
2030 do_radfix(xcoll, edi, step, cr);
2031 do_radacc(xcoll, edi, cr);
2032 do_radcon(xcoll, edi, cr);
2034 /* add back the average positions */
2035 for (i=0; i<edi->sav.nr; i++)
2036 rvec_inc(xcoll[i], edi->sav.x[i]);
2038 GMX_MPE_LOG(ev_ed_apply_cons_finish);
2042 /* Write out the projections onto the eigenvectors */
2043 static void write_edo(int nr_edi, t_edpar *edi, gmx_edsam_t ed, gmx_large_int_t step,real rmsd)
2049 if (edi->bNeedDoEdsam)
2052 fprintf(ed->edo, "Initial projections:\n");
2055 fprintf(ed->edo,"Step %s, ED #%d ", gmx_step_str(step, buf), nr_edi);
2056 fprintf(ed->edo," RMSD %f nm\n",rmsd);
2059 if (edi->vecs.mon.neig)
2061 fprintf(ed->edo," Monitor eigenvectors");
2062 for (i=0; i<edi->vecs.mon.neig; i++)
2063 fprintf(ed->edo," %d: %12.5e ",edi->vecs.mon.ieig[i],edi->vecs.mon.xproj[i]);
2064 fprintf(ed->edo,"\n");
2066 if (edi->vecs.linfix.neig)
2068 fprintf(ed->edo," Linfix eigenvectors");
2069 for (i=0; i<edi->vecs.linfix.neig; i++)
2070 fprintf(ed->edo," %d: %12.5e ",edi->vecs.linfix.ieig[i],edi->vecs.linfix.xproj[i]);
2071 fprintf(ed->edo,"\n");
2073 if (edi->vecs.linacc.neig)
2075 fprintf(ed->edo," Linacc eigenvectors");
2076 for (i=0; i<edi->vecs.linacc.neig; i++)
2077 fprintf(ed->edo," %d: %12.5e ",edi->vecs.linacc.ieig[i],edi->vecs.linacc.xproj[i]);
2078 fprintf(ed->edo,"\n");
2080 if (edi->vecs.radfix.neig)
2082 fprintf(ed->edo," Radfix eigenvectors");
2083 for (i=0; i<edi->vecs.radfix.neig; i++)
2084 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radfix.ieig[i],edi->vecs.radfix.xproj[i]);
2085 fprintf(ed->edo,"\n");
2086 fprintf(ed->edo," fixed increment radius = %f\n", calc_radius(&edi->vecs.radfix));
2088 if (edi->vecs.radacc.neig)
2090 fprintf(ed->edo," Radacc eigenvectors");
2091 for (i=0; i<edi->vecs.radacc.neig; i++)
2092 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radacc.ieig[i],edi->vecs.radacc.xproj[i]);
2093 fprintf(ed->edo,"\n");
2094 fprintf(ed->edo," acceptance radius = %f\n", calc_radius(&edi->vecs.radacc));
2096 if (edi->vecs.radcon.neig)
2098 fprintf(ed->edo," Radcon eigenvectors");
2099 for (i=0; i<edi->vecs.radcon.neig; i++)
2100 fprintf(ed->edo," %d: %12.5e ",edi->vecs.radcon.ieig[i],edi->vecs.radcon.xproj[i]);
2101 fprintf(ed->edo,"\n");
2102 fprintf(ed->edo," contracting radius = %f\n", calc_radius(&edi->vecs.radcon));
2107 /* Returns if any constraints are switched on */
2108 static int ed_constraints(gmx_bool edtype, t_edpar *edi)
2110 if (edtype == eEDedsam || edtype == eEDflood)
2112 return (edi->vecs.linfix.neig || edi->vecs.linacc.neig ||
2113 edi->vecs.radfix.neig || edi->vecs.radacc.neig ||
2114 edi->vecs.radcon.neig);
2120 /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
2121 * umbrella sampling simulations. */
2122 static void copyEvecReference(t_eigvec* floodvecs)
2127 if (NULL==floodvecs->refproj0)
2128 snew(floodvecs->refproj0, floodvecs->neig);
2130 for (i=0; i<floodvecs->neig; i++)
2132 floodvecs->refproj0[i] = floodvecs->refproj[i];
2137 void init_edsam(gmx_mtop_t *mtop, /* global topology */
2138 t_inputrec *ir, /* input record */
2139 t_commrec *cr, /* communication record */
2140 gmx_edsam_t ed, /* contains all ED data */
2141 rvec x[], /* positions of the whole MD system */
2142 matrix box) /* the box */
2144 t_edpar *edi = NULL; /* points to a single edi data set */
2145 int numedis=0; /* keep track of the number of ED data sets in edi file */
2146 int i,nr_edi,avindex;
2147 rvec *x_pbc = NULL; /* positions of the whole MD system with pbc removed */
2148 rvec *xfit = NULL; /* the positions which will be fitted to the reference structure */
2149 rvec *xstart = NULL; /* the positions which are subject to ED sampling */
2150 rvec fit_transvec; /* translation ... */
2151 matrix fit_rotmat; /* ... and rotation from fit to reference structure */
2154 if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
2155 gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
2157 GMX_MPE_LOG(ev_edsam_start);
2160 fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
2162 /* Needed for initializing radacc radius in do_edsam */
2165 /* The input file is read by the master and the edi structures are
2166 * initialized here. Input is stored in ed->edpar. Then the edi
2167 * structures are transferred to the other nodes */
2171 /* Read the whole edi file at once: */
2172 read_edi_file(ed,ed->edpar,mtop->natoms,cr);
2174 /* Initialization for every ED/flooding dataset. Flooding uses one edi dataset per
2175 * flooding vector, Essential dynamics can be applied to more than one structure
2176 * as well, but will be done in the order given in the edi file, so
2177 * expect different results for different order of edi file concatenation! */
2181 init_edi(mtop,ir,cr,ed,edi);
2183 /* Init flooding parameters if needed */
2184 init_flood(edi,ed,ir->delta_t,cr);
2191 /* The master does the work here. The other nodes get the positions
2192 * not before dd_partition_system which is called after init_edsam */
2195 /* Remove pbc, make molecule whole.
2196 * When ir->bContinuation=TRUE this has already been done, but ok.
2198 snew(x_pbc,mtop->natoms);
2199 m_rveccopy(mtop->natoms,x,x_pbc);
2200 do_pbc_first_mtop(NULL,ir->ePBC,box,mtop,x_pbc);
2202 /* Reset pointer to first ED data set which contains the actual ED data */
2205 /* Loop over all ED/flooding data sets (usually only one, though) */
2206 for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2208 /* We use srenew to allocate memory since the size of the buffers
2209 * is likely to change with every ED dataset */
2210 srenew(xfit , edi->sref.nr );
2211 srenew(xstart, edi->sav.nr );
2213 /* Extract the positions of the atoms to which will be fitted */
2214 for (i=0; i < edi->sref.nr; i++)
2216 copy_rvec(x_pbc[edi->sref.anrs[i]], xfit[i]);
2218 /* Save the sref positions such that in the next time step we can make the ED group whole
2219 * in case any of the atoms do not have the correct PBC representation */
2220 copy_rvec(xfit[i], edi->sref.x_old[i]);
2223 /* Extract the positions of the atoms subject to ED sampling */
2224 for (i=0; i < edi->sav.nr; i++)
2226 copy_rvec(x_pbc[edi->sav.anrs[i]], xstart[i]);
2228 /* Save the sav positions such that in the next time step we can make the ED group whole
2229 * in case any of the atoms do not have the correct PBC representation */
2230 copy_rvec(xstart[i], edi->sav.x_old[i]);
2233 /* Make the fit to the REFERENCE structure, get translation and rotation */
2234 fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
2236 /* Output how well we fit to the reference at the start */
2237 translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
2238 fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm (dataset #%d)\n",
2239 rmsd_from_structure(xfit, &edi->sref), nr_edi);
2241 /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
2242 translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
2244 /* calculate initial projections */
2245 project(xstart, edi);
2247 /* For the target and origin structure both a reference (fit) and an
2248 * average structure can be provided in make_edi. If both structures
2249 * are the same, make_edi only stores one of them in the .edi file.
2250 * If they differ, first the fit and then the average structure is stored
2251 * in star (or sor), thus the number of entries in star/sor is
2252 * (n_fit + n_av) with n_fit the size of the fitting group and n_av
2253 * the size of the average group. */
2255 /* process target structure, if required */
2256 if (edi->star.nr > 0)
2258 fprintf(stderr, "ED: Fitting target structure to reference structure\n");
2260 /* get translation & rotation for fit of target structure to reference structure */
2261 fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi);
2263 translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat);
2264 if (edi->star.nr == edi->sav.nr)
2268 else /* edi->star.nr = edi->sref.nr + edi->sav.nr */
2270 /* The last sav.nr indices of the target structure correspond to
2271 * the average structure, which must be projected */
2272 avindex = edi->star.nr - edi->sav.nr;
2274 rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon, cr);
2276 rad_project(edi, xstart, &edi->vecs.radcon, cr);
2278 /* process structure that will serve as origin of expansion circle */
2279 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2280 fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
2282 if (edi->sori.nr > 0)
2284 fprintf(stderr, "ED: Fitting origin structure to reference structure\n");
2286 /* fit this structure to reference structure */
2287 fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi);
2289 translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat);
2290 if (edi->sori.nr == edi->sav.nr)
2294 else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */
2296 /* For the projection, we need the last sav.nr indices of sori */
2297 avindex = edi->sori.nr - edi->sav.nr;
2300 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc, cr);
2301 rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix, cr);
2302 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2304 fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
2305 /* Set center of flooding potential to the ORIGIN structure */
2306 rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs, cr);
2307 /* We already know that no (moving) reference position was provided,
2308 * therefore we can overwrite refproj[0]*/
2309 copyEvecReference(&edi->flood.vecs);
2312 else /* No origin structure given */
2314 rad_project(edi, xstart, &edi->vecs.radacc, cr);
2315 rad_project(edi, xstart, &edi->vecs.radfix, cr);
2316 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2318 if (edi->flood.bHarmonic)
2320 fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
2321 for (i=0; i<edi->flood.vecs.neig; i++)
2322 edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
2326 fprintf(stderr, "ED: The AVERAGE structure will define the flooding potential center.\n");
2327 /* Set center of flooding potential to the center of the covariance matrix,
2328 * i.e. the average structure, i.e. zero in the projected system */
2329 for (i=0; i<edi->flood.vecs.neig; i++)
2330 edi->flood.vecs.refproj[i] = 0.0;
2334 /* For convenience, output the center of the flooding potential for the eigenvectors */
2335 if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
2337 for (i=0; i<edi->flood.vecs.neig; i++)
2339 fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", i, edi->flood.vecs.refproj[i]);
2340 if (edi->flood.bHarmonic)
2341 fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
2342 fprintf(stdout, "\n");
2346 /* set starting projections for linsam */
2347 rad_project(edi, xstart, &edi->vecs.linacc, cr);
2348 rad_project(edi, xstart, &edi->vecs.linfix, cr);
2350 /* Output to file, set the step to -1 so that write_edo knows it was called from init_edsam */
2351 if (ed->edo && !(ed->bStartFromCpt))
2352 write_edo(nr_edi, edi, ed, -1, 0);
2354 /* Prepare for the next edi data set: */
2357 /* Cleaning up on the master node: */
2362 } /* end of MASTER only section */
2366 /* First let everybody know how many ED data sets to expect */
2367 gmx_bcast(sizeof(numedis), &numedis, cr);
2368 /* Broadcast the essential dynamics / flooding data to all nodes */
2369 broadcast_ed_data(cr, ed, numedis);
2373 /* In the single-CPU case, point the local atom numbers pointers to the global
2374 * one, so that we can use the same notation in serial and parallel case: */
2376 /* Loop over all ED data sets (usually only one, though) */
2378 for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2380 edi->sref.anrs_loc = edi->sref.anrs;
2381 edi->sav.anrs_loc = edi->sav.anrs;
2382 edi->star.anrs_loc = edi->star.anrs;
2383 edi->sori.anrs_loc = edi->sori.anrs;
2384 /* For the same reason as above, make a dummy c_ind array: */
2385 snew(edi->sav.c_ind, edi->sav.nr);
2386 /* Initialize the array */
2387 for (i=0; i<edi->sav.nr; i++)
2388 edi->sav.c_ind[i] = i;
2389 /* In the general case we will need a different-sized array for the reference indices: */
2392 snew(edi->sref.c_ind, edi->sref.nr);
2393 for (i=0; i<edi->sref.nr; i++)
2394 edi->sref.c_ind[i] = i;
2396 /* Point to the very same array in case of other structures: */
2397 edi->star.c_ind = edi->sav.c_ind;
2398 edi->sori.c_ind = edi->sav.c_ind;
2399 /* In the serial case, the local number of atoms is the global one: */
2400 edi->sref.nr_loc = edi->sref.nr;
2401 edi->sav.nr_loc = edi->sav.nr;
2402 edi->star.nr_loc = edi->star.nr;
2403 edi->sori.nr_loc = edi->sori.nr;
2405 /* An on we go to the next edi dataset */
2410 /* Allocate space for ED buffer variables */
2411 /* Again, loop over ED data sets */
2413 for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
2415 /* Allocate space for ED buffer */
2417 snew(edi->buf->do_edsam, 1);
2419 /* Space for collective ED buffer variables */
2421 /* Collective positions of atoms with the average indices */
2422 snew(edi->buf->do_edsam->xcoll , edi->sav.nr);
2423 snew(edi->buf->do_edsam->shifts_xcoll , edi->sav.nr); /* buffer for xcoll shifts */
2424 snew(edi->buf->do_edsam->extra_shifts_xcoll , edi->sav.nr);
2425 /* Collective positions of atoms with the reference indices */
2428 snew(edi->buf->do_edsam->xc_ref , edi->sref.nr);
2429 snew(edi->buf->do_edsam->shifts_xc_ref , edi->sref.nr); /* To store the shifts in */
2430 snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr);
2433 /* Get memory for flooding forces */
2434 snew(edi->flood.forces_cartesian , edi->sav.nr);
2437 /* Dump it all into one file per process */
2438 dump_edi(edi, cr, nr_edi);
2441 /* An on we go to the next edi dataset */
2445 /* Flush the edo file so that the user can check some things
2446 * when the simulation has started */
2450 GMX_MPE_LOG(ev_edsam_finish);
2454 void do_edsam(t_inputrec *ir,
2455 gmx_large_int_t step,
2458 rvec xs[], /* The local current positions on this processor */
2459 rvec v[], /* The velocities */
2463 int i,edinr,iupdate=500;
2464 matrix rotmat; /* rotation matrix */
2465 rvec transvec; /* translation vector */
2466 rvec dv,dx,x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */
2467 real dt_1; /* 1/dt */
2468 struct t_do_edsam *buf;
2470 real rmsdev=-1; /* RMSD from reference structure prior to applying the constraints */
2471 gmx_bool bSuppress=FALSE; /* Write .edo file on master? */
2474 /* Check if ED sampling has to be performed */
2475 if ( ed->eEDtype==eEDnone )
2478 /* Suppress output on first call of do_edsam if
2479 * two-step sd2 integrator is used */
2480 if ( (ir->eI==eiSD2) && (v != NULL) )
2483 dt_1 = 1.0/ir->delta_t;
2485 /* Loop over all ED datasets (usually one) */
2491 if (edi->bNeedDoEdsam)
2494 buf=edi->buf->do_edsam;
2497 /* initialise radacc radius for slope criterion */
2498 buf->oldrad=calc_radius(&edi->vecs.radacc);
2500 /* Copy the positions into buf->xc* arrays and after ED
2501 * feed back corrections to the official positions */
2503 /* Broadcast the ED positions such that every node has all of them
2504 * Every node contributes its local positions xs and stores it in
2505 * the collective buf->xcoll array. Note that for edinr > 1
2506 * xs could already have been modified by an earlier ED */
2508 communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
2509 edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
2512 dump_xcoll(edi, buf, cr, step);
2514 /* Only assembly reference positions if their indices differ from the average ones */
2516 communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
2517 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
2519 /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
2520 * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
2521 * set bUpdateShifts=TRUE in the parallel case. */
2522 buf->bUpdateShifts = FALSE;
2524 /* Now all nodes have all of the ED positions in edi->sav->xcoll,
2525 * as well as the indices in edi->sav.anrs */
2527 /* Fit the reference indices to the reference structure */
2529 fit_to_reference(buf->xcoll , transvec, rotmat, edi);
2531 fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
2533 /* Now apply the translation and rotation to the ED structure */
2534 translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
2536 /* Find out how well we fit to the reference (just for output steps) */
2537 if (do_per_step(step,edi->outfrq) && MASTER(cr))
2541 /* Indices of reference and average structures are identical,
2542 * thus we can calculate the rmsd to SREF using xcoll */
2543 rmsdev = rmsd_from_structure(buf->xcoll,&edi->sref);
2547 /* We have to translate & rotate the reference atoms first */
2548 translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
2549 rmsdev = rmsd_from_structure(buf->xc_ref,&edi->sref);
2553 /* update radsam references, when required */
2554 if (do_per_step(step,edi->maxedsteps) && step >= edi->presteps)
2556 project(buf->xcoll, edi);
2557 rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
2558 rad_project(edi, buf->xcoll, &edi->vecs.radfix, cr);
2562 /* update radacc references, when required */
2563 if (do_per_step(step,iupdate) && step >= edi->presteps)
2565 edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc);
2566 if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
2568 project(buf->xcoll, edi);
2569 rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
2572 buf->oldrad = edi->vecs.radacc.radius;
2575 /* apply the constraints */
2576 if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi))
2578 /* ED constraints should be applied already in the first MD step
2579 * (which is step 0), therefore we pass step+1 to the routine */
2580 ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step, cr);
2583 /* write to edo, when required */
2584 if (do_per_step(step,edi->outfrq))
2586 project(buf->xcoll, edi);
2587 if (MASTER(cr) && !bSuppress)
2588 write_edo(edinr, edi, ed, step, rmsdev);
2591 /* Copy back the positions unless monitoring only */
2592 if (ed_constraints(ed->eEDtype, edi))
2594 /* remove fitting */
2595 rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat);
2597 /* Copy the ED corrected positions into the coordinate array */
2598 /* Each node copies its local part. In the serial case, nat_loc is the
2599 * total number of ED atoms */
2600 for (i=0; i<edi->sav.nr_loc; i++)
2602 /* Unshift local ED coordinate and store in x_unsh */
2603 ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]],
2604 buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh);
2606 /* dx is the ED correction to the positions: */
2607 rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx);
2611 /* dv is the ED correction to the velocity: */
2612 svmul(dt_1, dx, dv);
2613 /* apply the velocity correction: */
2614 rvec_inc(v[edi->sav.anrs_loc[i]], dv);
2616 /* Finally apply the position correction due to ED: */
2617 copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]);
2620 } /* END of if (edi->bNeedDoEdsam) */
2622 /* Prepare for the next ED dataset */
2623 edi = edi->next_edi;
2625 } /* END of loop over ED datasets */
2629 GMX_MPE_LOG(ev_edsam_finish);