#include "mtop_util.h"
#include "edsam.h"
#include "gmxfio.h"
+#include "xvgr.h"
#include "groupcoord.h"
#define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
#define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
+/* These macros determine the column width in the output file */
+#define EDcol_sfmt "%17s"
+#define EDcol_efmt "%17.5e"
+#define EDcol_ffmt "%17f"
/* enum to identify the type of ED: none, normal ED, flooding */
enum {eEDnone, eEDedsam, eEDflood, eEDnr};
real dt;
real constEfl;
real alpha2;
- int flood_id;
rvec *forces_cartesian;
t_eigvec vecs; /* use flooding for these */
} t_edflood;
* with respect to the collective
* anrs[0...nr-1] array */
rvec *x; /* positions for this structure */
- rvec *x_old; /* used to keep track of the shift vectors
- such that the ED molecule can always be
- made whole in the parallel case */
+ rvec *x_old; /* Last positions which have the correct PBC
+ representation of the ED group. In
+ combination with keeping track of the
+ shift vectors, the ED group can always
+ be made whole */
real *m; /* masses */
real mtot; /* total mass (only used in sref) */
real *sqrtm; /* sqrt of the masses used for mass-
* is used (i.e. apart from flooding) */
t_edflood flood; /* parameters especially for flooding */
struct t_ed_buffer *buf; /* handle to local buffers */
- struct edpar *next_edi; /* Pointer to another ed dataset */
+ struct edpar *next_edi; /* Pointer to another ED group */
} t_edpar;
typedef struct gmx_edsam
{
int eEDtype; /* Type of ED: see enums above */
- const char *edinam; /* name of ED sampling input file */
- const char *edonam; /* output */
FILE *edo; /* output file pointer */
t_edpar *edpar;
gmx_bool bFirst;
- gmx_bool bStartFromCpt;
} t_gmx_edsam;
ivec *shifts_xc_ref; /* Shifts for xc_ref */
ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */
gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the
- ED shifts for this ED dataset need to
+ ED shifts for this ED group need to
be updated */
};
/* Function declarations */
static void fit_to_reference(rvec *xcoll,rvec transvec,matrix rotmat,t_edpar *edi);
-
static void translate_and_rotate(rvec *x,int nat,rvec transvec,matrix rotmat);
+static real rmsd_from_structure(rvec *x, struct gmx_edx *s);
+static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms);
+static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate);
+static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate);
+static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv);
/* End function declarations */
+/* Multiple ED groups will be labeled with letters instead of numbers
+ * to avoid confusion with eigenvector indices */
+static char get_EDgroupChar(int nr_edi, int nED)
+{
+ if (nED == 1)
+ {
+ return ' ';
+ }
+
+ /* nr_edi = 1 -> A
+ * nr_edi = 2 -> B ...
+ */
+ return 'A' + nr_edi - 1;
+}
+
+
/* Does not subtract average positions, projection on single eigenvector is returned
* used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon
* Average position is subtracted in ed_apply_constraints prior to calling projectx
for (i=0; i<edi->sav.nr; i++)
+ {
proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]);
+ }
return proj;
}
/* Specialized: projection is stored in vec->refproj
* -> used for radacc, radfix, radcon and center of flooding potential
* subtracts average positions, projects vector x */
-static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec, t_commrec *cr)
+static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec)
{
int i;
real rad=0.0;
/* Subtract average positions */
for (i = 0; i < edi->sav.nr; i++)
+ {
rvec_dec(x[i], edi->sav.x[i]);
+ }
for (i = 0; i < vec->neig; i++)
{
/* Add average positions */
for (i = 0; i < edi->sav.nr; i++)
+ {
rvec_inc(x[i], edi->sav.x[i]);
+ }
}
/* Subtract average positions */
for (i=0; i<edi->sav.nr; i++)
+ {
rvec_dec(x[i], edi->sav.x[i]);
+ }
for (i=0; i<vec->neig; i++)
+ {
vec->xproj[i] = projectx(edi, x, vec->vec[i]);
+ }
/* Add average positions */
for (i=0; i<edi->sav.nr; i++)
+ {
rvec_inc(x[i], edi->sav.x[i]);
+ }
}
for (i=0; i<vec->neig; i++)
+ {
rad += pow((vec->refproj[i]-vec->xproj[i]),2);
+ }
return rad=sqrt(rad);
}
fp = fopen(fn, "w");
for (i=0; i<edi->sav.nr; i++)
+ {
fprintf(fp, "%d %9.5f %9.5f %9.5f %d %d %d %d %d %d\n",
edi->sav.anrs[i]+1,
xcoll[i][XX] , xcoll[i][YY] , xcoll[i][ZZ],
shifts[i][XX] , shifts[i][YY] , shifts[i][ZZ],
eshifts[i][XX], eshifts[i][YY], eshifts[i][ZZ]);
+ }
fclose(fp);
}
fprintf(out, "#%s positions:\n%d\n", name, s->nr);
if (s->nr == 0)
+ {
return;
+ }
fprintf(out, "#index, x, y, z");
if (s->sqrtm)
+ {
fprintf(out, ", sqrt(m)");
+ }
for (i=0; i<s->nr; i++)
{
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]);
if (s->sqrtm)
+ {
fprintf(out,"%9.3f",s->sqrtm[i]);
+ }
}
fprintf(out, "\n");
}
fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n",
ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius);
for (j=0; j<length; j++)
+ {
fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX], ev->vec[i][j][YY], ev->vec[i][j][ZZ]);
+ }
}
}
for (i=0; i<dim; i++)
+ {
fprintf(out,"%4d %f %f %f\n",i,x[i][XX],x[i][YY],x[i][ZZ]);
+ }
}
for (i=0;i<dim;i++)
{
for (j=0;j<dim;j++)
+ {
fprintf(out,"%f ",mat[i][j]);
+ }
fprintf(out,"\n");
}
}
gmx_bool bFirst;
if(edi->buf->do_edfit != NULL)
+ {
bFirst = FALSE;
+ }
else
{
bFirst = TRUE;
/* construct loc->omega */
/* loc->omega is symmetric -> loc->omega==loc->omega' */
for(r=0;(r<6);r++)
+ {
for(c=0;(c<=r);c++)
+ {
if ((r>=3) && (c<3))
{
loc->omega[r][c]=u[r-3][c];
loc->omega[r][c]=0;
loc->omega[c][r]=0;
}
+ }
+ }
/* determine h and k */
#ifdef DEBUG
int i;
dump_mat(stderr,2*DIM,loc->omega);
for (i=0; i<6; i++)
+ {
fprintf(stderr,"d[%d] = %f\n",i,d[i]);
+ }
}
#endif
jacobi(loc->omega,6,d,loc->om,&irot);
if (irot==0)
+ {
fprintf(stderr,"IROT=0\n");
+ }
index=0; /* For the compiler only */
{
max_d=-1000;
for(i=0;(i<6);i++)
+ {
if (d[i]>max_d)
{
max_d=d[i];
index=i;
}
+ }
d[index]=-10000;
for(i=0;(i<3);i++)
{
/* determine R */
for(c=0;(c<3);c++)
+ {
for(r=0;(r<3);r++)
+ {
R[c][r]=vk[0][r]*vh[0][c]+
- vk[1][r]*vh[1][c]+
- vk[2][r]*vh[2][c];
+ vk[1][r]*vh[1][c]+
+ vk[2][r]*vh[2][c];
+ }
+ }
if (det(R) < 0)
+ {
for(c=0;(c<3);c++)
+ {
for(r=0;(r<3);r++)
+ {
R[c][r]=vk[0][r]*vh[0][c]+
- vk[1][r]*vh[1][c]-
- vk[2][r]*vh[2][c];
+ vk[1][r]*vh[1][c]-
+ vk[2][r]*vh[2][c];
+ }
+ }
+ }
}
two edsam files from two peptide chains
*/
-static void write_edo_flood(t_edpar *edi, FILE *fp, gmx_large_int_t step)
+static void write_edo_flood(t_edpar *edi, FILE *fp, real rmsd)
{
int i;
- char buf[22];
- gmx_bool bOutputRef=FALSE;
- fprintf(fp,"%d.th FL: %s %12.5e %12.5e %12.5e\n",
- edi->flood.flood_id, gmx_step_str(step,buf),
- edi->flood.Efl, edi->flood.Vfl, edi->flood.deltaF);
+ /* Output how well we fit to the reference structure */
+ fprintf(fp, EDcol_ffmt, rmsd);
-
- /* Check whether any of the references changes with time (this can happen
- * in case flooding is used as harmonic restraint). If so, output all the
- * current reference projections. */
- if (edi->flood.bHarmonic)
+ for (i=0; i<edi->flood.vecs.neig; i++)
{
- for (i = 0; i < edi->flood.vecs.neig; i++)
+ fprintf(fp, EDcol_efmt, edi->flood.vecs.xproj[i]);
+
+ /* Check whether the reference projection changes with time (this can happen
+ * in case flooding is used as harmonic restraint). If so, output the
+ * current reference projection */
+ if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
{
- if (edi->flood.vecs.refprojslope[i] != 0.0)
- bOutputRef=TRUE;
+ fprintf(fp, EDcol_efmt, edi->flood.vecs.refproj[i]);
}
- if (bOutputRef)
+
+ /* Output Efl if we are doing adaptive flooding */
+ if (0 != edi->flood.tau)
{
- fprintf(fp, "Ref. projs.: ");
- for (i = 0; i < edi->flood.vecs.neig; i++)
- {
- fprintf(fp, "%12.5e ", edi->flood.vecs.refproj[i]);
- }
- fprintf(fp, "\n");
+ fprintf(fp, EDcol_efmt, edi->flood.Efl);
}
- }
- fprintf(fp,"FL_FORCES: ");
-
- for (i=0; i<edi->flood.vecs.neig; i++)
- fprintf(fp," %12.5e",edi->flood.vecs.fproj[i]);
+ fprintf(fp, EDcol_efmt, edi->flood.Vfl);
- fprintf(fp,"\n");
+ /* Output deltaF if we are doing adaptive flooding */
+ if (0 != edi->flood.tau)
+ {
+ fprintf(fp, EDcol_efmt, edi->flood.deltaF);
+ }
+ fprintf(fp, EDcol_efmt, edi->flood.vecs.fproj[i]);
+ }
}
if (edi->flood.bHarmonic)
+ {
for (i=0; i<edi->flood.vecs.neig; i++)
{
edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]);
}
+ }
else
+ {
for (i=0; i<edi->flood.vecs.neig; i++)
{
/* if Efl is zero the forces are zero if not use the formula */
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;
}
+ }
}
/* Clear forces first */
for (j=0; j<edi->sav.nr_loc; j++)
+ {
clear_rvec(forces_cart[j]);
+ }
/* Now compute atomwise */
for (j=0; j<edi->sav.nr_loc; j++)
edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF);
/* check if restrain (inverted flooding) -> don't let EFL become positive */
if (edi->flood.alpha2<0 && edi->flood.Efl>-0.00000001)
+ {
edi->flood.Efl = 0;
+ }
edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl;
}
matrix rotmat; /* rotation matrix */
matrix tmat; /* inverse rotation */
rvec transvec; /* translation vector */
+ real rmsdev;
struct t_do_edsam *buf;
/* Only assembly REFERENCE positions if their indices differ from the average ones */
if (!edi->bRefEqAv)
+ {
communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
+ }
/* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
* We do not need to update the shifts until the next NS step */
/* Fit the reference indices to the reference structure */
if (edi->bRefEqAv)
+ {
fit_to_reference(buf->xcoll , transvec, rotmat, edi);
+ }
else
+ {
fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
+ }
/* Now apply the translation and rotation to the ED structure */
translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
/* Finally add forces to the main force variable */
for (i=0; i<edi->sav.nr_loc; i++)
+ {
rvec_inc(force[edi->sav.anrs_loc[i]],edi->flood.forces_cartesian[i]);
+ }
/* Output is written by the master process */
if (do_per_step(step,edi->outfrq) && MASTER(cr))
- write_edo_flood(edi,edo,step);
+ {
+ /* Output how well we fit to the reference */
+ if (edi->bRefEqAv)
+ {
+ /* Indices of reference and average structures are identical,
+ * thus we can calculate the rmsd to SREF using xcoll */
+ rmsdev = rmsd_from_structure(buf->xcoll,&edi->sref);
+ }
+ else
+ {
+ /* We have to translate & rotate the reference atoms first */
+ translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat);
+ rmsdev = rmsd_from_structure(buf->xc_ref,&edi->sref);
+ }
+
+ write_edo_flood(edi,edo,rmsdev);
+ }
}
/* Main flooding routine, called from do_force */
extern void do_flood(
- FILE *log, /* md.log file */
t_commrec *cr, /* Communication record */
+ t_inputrec *ir, /* Input record */
rvec x[], /* Positions on the local processor */
rvec force[], /* forcefield forces, to these the flooding forces are added */
- gmx_edsam_t ed, /* ed data structure contains all ED and flooding datasets */
+ gmx_edsam_t ed, /* ed data structure contains all ED and flooding groups */
matrix box, /* the box */
gmx_large_int_t step, /* The relative time step since ir->init_step is already subtracted */
gmx_bool bNS) /* Are we in a neighbor searching step? */
t_edpar *edi;
+ edi = ed->edpar;
+
+ /* Write time to edo, when required. Output the time anyhow since we need
+ * it in the output file for ED constraints. */
+ if (MASTER(cr) && do_per_step(step,edi->outfrq))
+ {
+ fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t);
+ }
+
if (ed->eEDtype != eEDflood)
+ {
return;
+ }
- edi = ed->edpar;
while (edi)
{
/* Call flooding for one matrix */
if (edi->flood.vecs.neig)
+ {
do_single_flood(ed->edo,x,force,edi,step,box,cr,bNS);
+ }
edi = edi->next_edi;
}
}
/* Called by init_edi, configure some flooding related variables and structures,
* print headers to output files */
-static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt, t_commrec *cr)
+static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt)
{
int i;
if (edi->flood.vecs.neig)
{
- /* If in any of the datasets we find a flooding vector, flooding is turned on */
+ /* If in any of the ED groups we find a flooding vector, flooding is turned on */
ed->eEDtype = eEDflood;
- fprintf(stderr,"ED: Flooding of matrix %d is switched on.\n", edi->flood.flood_id);
+ fprintf(stderr,"ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s":"");
if (edi->flood.bConstForce)
{
edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]);
}
}
- fprintf(ed->edo,"FL_HEADER: Flooding of matrix %d is switched on! The flooding output will have the following format:\n",
- edi->flood.flood_id);
- fprintf(ed->edo,"FL_HEADER: Step Efl Vfl deltaF\n");
}
}
count++;
}
if (nnames!=count-1)
+ {
gmx_fatal(FARGS,"Number of energies is not consistent with t_edi structure");
+ }
}
/************* END of FLOODING IMPLEMENTATION ****************************/
#endif
-gmx_edsam_t ed_open(int nfile,const t_filenm fnm[],unsigned long Flags,t_commrec *cr)
+gmx_edsam_t ed_open(int natoms, edsamstate_t *EDstate, int nfile,const t_filenm fnm[],unsigned long Flags, const output_env_t oenv, t_commrec *cr)
{
gmx_edsam_t ed;
+ int nED;
/* Allocate space for the ED data structure */
if (MASTER(cr))
{
- /* Open .edi input file: */
- ed->edinam=ftp2fn(efEDI,nfile,fnm);
- /* The master opens the .edo output file */
fprintf(stderr,"ED sampling will be performed!\n");
- ed->edonam = ftp2fn(efEDO,nfile,fnm);
- ed->edo = gmx_fio_fopen(ed->edonam,(Flags & MD_APPENDFILES)? "a+" : "w+");
- ed->bStartFromCpt = Flags & MD_STARTFROMCPT;
+ snew(ed->edpar,1);
+
+ /* Read the edi input file: */
+ nED = read_edi_file(ftp2fn(efEDI,nfile,fnm),ed->edpar,natoms);
+
+ /* Make sure the checkpoint was produced in a run using this .edi file */
+ if (EDstate->bFromCpt)
+ {
+ crosscheck_edi_file_vs_checkpoint(ed, EDstate);
+ }
+ else
+ {
+ EDstate->nED = nED;
+ }
+ init_edsamstate(ed, EDstate);
+
+ /* The master opens the ED output file */
+ if (Flags & MD_APPENDFILES)
+ {
+ ed->edo = gmx_fio_fopen(opt2fn("-eo",nfile,fnm),"a+");
+ }
+ else
+ {
+ ed->edo = xvgropen(opt2fn("-eo",nfile,fnm),
+ "Essential dynamics / flooding output",
+ "Time (ps)",
+ "RMSDs (nm), projections on EVs (nm), ...", oenv);
+
+ /* Make a descriptive legend */
+ write_edo_legend(ed, EDstate->nED, oenv);
+ }
}
return ed;
}
/* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
- /* Set the pointer to the next ED dataset */
+ /* Set the pointer to the next ED group */
if (edi->next_edi)
{
snew_bc(cr, edi->next_edi, 1);
/* init-routine called for every *.edi-cycle, initialises t_edpar structure */
-static void init_edi(gmx_mtop_t *mtop,t_inputrec *ir,
- t_commrec *cr,gmx_edsam_t ed,t_edpar *edi)
+static void init_edi(gmx_mtop_t *mtop,t_edpar *edi)
{
int i;
real totalmass = 0.0;
static void check(const char *line, const char *label)
{
if (!strstr(line,label))
+ {
gmx_fatal(FARGS,"Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s",label,line);
+ }
}
sscanf (line,"%d%lf%lf%lf",&anrs[i],&d[0],&d[1],&d[2]);
anrs[i]--; /* we are reading FORTRAN indices */
for(j=0; j<3; j++)
+ {
x[i][j]=d[j]; /* always read as double and convert to single */
+ }
}
}
{
nscan = sscanf(line,"%d%lf",&idum,&rdum);
if (nscan != 2)
+ {
gmx_fatal(FARGS,"Expected 2 values for flooding vec: <nr> <stpsz>\n");
+ }
}
tvec->ieig[i]=idum;
tvec->stpsz[i]=rdum;
/* If the number of atoms differs between the two structures,
* they cannot be identical */
if (sref.nr != sav.nr)
+ {
return FALSE;
+ }
/* Now that we know that both stuctures have the same number of atoms,
* check if also the indices are identical */
for (i=0; i < sav.nr; i++)
{
if (sref.anrs[i] != sav.anrs[i])
+ {
return FALSE;
+ }
}
fprintf(stderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n");
}
-static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int edi_nr, t_commrec *cr)
+static int read_edi(FILE* in,t_edpar *edi,int nr_mdatoms, const char *fn)
{
int readmagic;
const int magic=670;
readmagic=read_edint(in,&bEOF);
/* Check whether we have reached the end of the input file */
if (bEOF)
+ {
return 0;
+ }
if (readmagic != magic)
{
if (readmagic==666 || readmagic==667 || readmagic==668)
+ {
gmx_fatal(FARGS,"Wrong magic number: Use newest version of make_edi to produce edi file");
+ }
else if (readmagic != 669)
- gmx_fatal(FARGS,"Wrong magic number %d in %s",readmagic,ed->edinam);
+ {
+ gmx_fatal(FARGS,"Wrong magic number %d in %s",readmagic,fn);
+ }
}
/* check the number of atoms */
edi->nini=read_edint(in,&bEOF);
if (edi->nini != nr_mdatoms)
- gmx_fatal(FARGS,"Nr of atoms in %s (%d) does not match nr of md atoms (%d)",
- ed->edinam,edi->nini,nr_mdatoms);
+ {
+ gmx_fatal(FARGS,"Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn,edi->nini,nr_mdatoms);
+ }
/* Done checking. For the rest we blindly trust the input */
edi->fitmas = read_checked_edint(in,"FITMAS");
edi->flood.kT = read_checked_edreal(in,"KT");
edi->flood.bHarmonic = read_checked_edint(in,"HARMONIC");
if (readmagic > 669)
+ {
edi->flood.bConstForce = read_checked_edint(in,"CONST_FORCE_FLOODING");
+ }
else
+ {
edi->flood.bConstForce = FALSE;
- edi->flood.flood_id = edi_nr;
+ }
edi->sref.nr = read_checked_edint(in,"NREF");
/* allocate space for reference positions and read them */
edi->sori.nr=read_edint(in,&bEOF);
if (edi->sori.nr > 0)
{
- if (bHaveReference)
- {
- /* Both an -ori structure and a at least one manual reference point have been
- * specified. That's ambiguous and probably not intentional. */
- gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
- " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
- }
+ if (bHaveReference)
+ {
+ /* Both an -ori structure and a at least one manual reference point have been
+ * specified. That's ambiguous and probably not intentional. */
+ gmx_fatal(FARGS, "ED: An origin structure has been provided and a at least one (moving) reference\n"
+ " point was manually specified in the edi file. That is ambiguous. Aborting.\n");
+ }
snew(edi->sori.anrs,edi->sori.nr);
snew(edi->sori.x ,edi->sori.nr);
edi->sori.sqrtm =NULL;
/* Read in the edi input file. Note that it may contain several ED data sets which were
* achieved by concatenating multiple edi files. The standard case would be a single ED
* data set, though. */
-static void read_edi_file(gmx_edsam_t ed, t_edpar *edi, int nr_mdatoms, t_commrec *cr)
+static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms)
{
FILE *in;
t_edpar *curr_edi,*last_edi;
/* This routine is executed on the master only */
/* Open the .edi parameter input file */
- in = gmx_fio_fopen(ed->edinam,"r");
- fprintf(stderr, "ED: Reading edi file %s\n", ed->edinam);
+ in = gmx_fio_fopen(fn,"r");
+ fprintf(stderr, "ED: Reading edi file %s\n", fn);
/* Now read a sequence of ED input parameter sets from the edi file */
curr_edi=edi;
last_edi=edi;
- while( read_edi(in, ed, curr_edi, nr_mdatoms, edi_nr, cr) )
+ while( read_edi(in, curr_edi, nr_mdatoms, fn) )
{
edi_nr++;
- /* Make shure that the number of atoms in each dataset is the same as in the tpr file */
- if (edi->nini != nr_mdatoms)
- gmx_fatal(FARGS,"edi file %s (dataset #%d) was made for %d atoms, but the simulation contains %d atoms.",
- ed->edinam, edi_nr, edi->nini, nr_mdatoms);
+
/* Since we arrived within this while loop we know that there is still another data set to be read in */
/* We need to allocate space for the data: */
snew(edi_read,1);
/* Point the 'next_edi' entry to the next edi: */
curr_edi->next_edi=edi_read;
- /* Keep the curr_edi pointer for the case that the next dataset is empty: */
+ /* Keep the curr_edi pointer for the case that the next group is empty: */
last_edi = curr_edi;
/* Let's prepare to read in the next edi data set: */
curr_edi = edi_read;
}
if (edi_nr == 0)
- gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", ed->edinam);
+ {
+ gmx_fatal(FARGS, "No complete ED data set found in edi file %s.", fn);
+ }
- /* Terminate the edi dataset list with a NULL pointer: */
+ /* Terminate the edi group list with a NULL pointer: */
last_edi->next_edi = NULL;
- fprintf(stderr, "ED: Found %d ED dataset%s.\n", edi_nr, edi_nr>1? "s" : "");
+ fprintf(stderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr>1? "s" : "");
/* Close the .edi file again */
gmx_fio_fclose(in);
+
+ return edi_nr;
}
struct t_fit_to_ref *loc;
- /* Allocate memory the first time this routine is called for each edi dataset */
+ /* Allocate memory the first time this routine is called for each edi group */
if (NULL == edi->buf->fit_to_ref)
{
snew(edi->buf->fit_to_ref, 1);
/* We do not touch the original positions but work on a copy. */
for (i=0; i<edi->sref.nr; i++)
+ {
copy_rvec(xcoll[i], loc->xcopy[i]);
+ }
/* Calculate the center of mass */
get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com);
for (i=0; i < s->nr; i++)
+ {
rmsd += distance2(s->x[i], x[i]);
+ }
rmsd /= (real) s->nr;
rmsd = sqrt(rmsd);
if (ed->eEDtype != eEDnone)
{
- /* Loop over ED datasets (usually there is just one dataset, though) */
+ /* Loop over ED groups */
edi=ed->edpar;
while (edi)
{
/* Local atoms of the reference structure (for fitting), need only be assembled
* if their indices differ from the average ones */
if (!edi->bRefEqAv)
+ {
dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs,
&edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind);
+ }
/* Local atoms of the average structure (on these ED will be performed) */
dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs,
* at the next call to communicate_group_positions, since obviously we are in a NS step */
edi->buf->do_edsam->bUpdateShifts = TRUE;
- /* Set the pointer to the next ED dataset (if any) */
+ /* Set the pointer to the next ED group (if any) */
edi=edi->next_edi;
}
}
xu[XX] = x[XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
xu[YY] = x[YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
xu[ZZ] = x[ZZ]-tz*box[ZZ][ZZ];
- } else
+ }
+ else
{
xu[XX] = x[XX]-tx*box[XX][XX];
xu[YY] = x[YY]-ty*box[YY][YY];
}
-static void do_linfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
+static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_large_int_t step)
{
int i, j;
real proj, add;
}
-static void do_linacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
+static void do_linacc(rvec *xcoll, t_edpar *edi)
{
int i, j;
real proj, add;
if (edi->vecs.linacc.stpsz[i] > 0.0)
{
if ((proj-edi->vecs.linacc.refproj[i]) < 0.0)
+ {
add = edi->vecs.linacc.refproj[i] - proj;
+ }
}
if (edi->vecs.linacc.stpsz[i] < 0.0)
{
if ((proj-edi->vecs.linacc.refproj[i]) > 0.0)
+ {
add = edi->vecs.linacc.refproj[i] - proj;
+ }
}
/* apply the correction */
}
-static void do_radfix(rvec *xcoll, t_edpar *edi, int step, t_commrec *cr)
+static void do_radfix(rvec *xcoll, t_edpar *edi)
{
int i,j;
real *proj, rad=0.0, ratio;
/* apply the correction */
proj[i] /= edi->sav.sqrtm[i];
proj[i] *= ratio;
- for (j=0; j<edi->sav.nr; j++) {
+ for (j=0; j<edi->sav.nr; j++)
+ {
svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum);
rvec_inc(xcoll[j], vec_dum);
}
}
-static void do_radacc(rvec *xcoll, t_edpar *edi, t_commrec *cr)
+static void do_radacc(rvec *xcoll, t_edpar *edi)
{
int i,j;
real *proj, rad=0.0, ratio=0.0;
real *proj;
};
-static void do_radcon(rvec *xcoll, t_edpar *edi, t_commrec *cr)
+static void do_radcon(rvec *xcoll, t_edpar *edi)
{
int i,j;
real rad=0.0, ratio=0.0;
loc = edi->buf->do_radcon;
if (edi->vecs.radcon.neig == 0)
+ {
return;
-
+ }
+
if (bFirst)
+ {
snew(loc->proj, edi->vecs.radcon.neig);
+ }
/* loop over radcon vectors */
for (i=0; i<edi->vecs.radcon.neig; i++)
}
-static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step, t_commrec *cr)
+static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_large_int_t step)
{
int i;
/* subtract the average positions */
for (i=0; i<edi->sav.nr; i++)
+ {
rvec_dec(xcoll[i], edi->sav.x[i]);
+ }
/* apply the constraints */
if (step >= 0)
- do_linfix(xcoll, edi, step, cr);
- do_linacc(xcoll, edi, cr);
+ {
+ do_linfix(xcoll, edi, step);
+ }
+ do_linacc(xcoll, edi);
if (step >= 0)
- do_radfix(xcoll, edi, step, cr);
- do_radacc(xcoll, edi, cr);
- do_radcon(xcoll, edi, cr);
+ {
+ do_radfix(xcoll, edi);
+ }
+ do_radacc(xcoll, edi);
+ do_radcon(xcoll, edi);
/* add back the average positions */
for (i=0; i<edi->sav.nr; i++)
+ {
rvec_inc(xcoll[i], edi->sav.x[i]);
+ }
}
-/* Write out the projections onto the eigenvectors */
-static void write_edo(int nr_edi, t_edpar *edi, gmx_edsam_t ed, gmx_large_int_t step,real rmsd)
+/* Write out the projections onto the eigenvectors. The order of output
+ * corresponds to ed_output_legend() */
+static void write_edo(t_edpar *edi, FILE *fp,real rmsd)
{
int i;
- char buf[22];
- if (edi->bNeedDoEdsam)
+ /* Output how well we fit to the reference structure */
+ fprintf(fp, EDcol_ffmt, rmsd);
+
+ for (i=0; i<edi->vecs.mon.neig; i++)
{
- if (step == -1)
- fprintf(ed->edo, "Initial projections:\n");
- else
- {
- fprintf(ed->edo,"Step %s, ED #%d ", gmx_step_str(step, buf), nr_edi);
- fprintf(ed->edo," RMSD %f nm\n",rmsd);
- }
+ fprintf(fp, EDcol_efmt, edi->vecs.mon.xproj[i]);
+ }
- if (edi->vecs.mon.neig)
- {
- fprintf(ed->edo," Monitor eigenvectors");
- for (i=0; i<edi->vecs.mon.neig; i++)
- fprintf(ed->edo," %d: %12.5e ",edi->vecs.mon.ieig[i],edi->vecs.mon.xproj[i]);
- fprintf(ed->edo,"\n");
- }
- if (edi->vecs.linfix.neig)
- {
- fprintf(ed->edo," Linfix eigenvectors");
- for (i=0; i<edi->vecs.linfix.neig; i++)
- fprintf(ed->edo," %d: %12.5e ",edi->vecs.linfix.ieig[i],edi->vecs.linfix.xproj[i]);
- fprintf(ed->edo,"\n");
- }
- if (edi->vecs.linacc.neig)
- {
- fprintf(ed->edo," Linacc eigenvectors");
- for (i=0; i<edi->vecs.linacc.neig; i++)
- fprintf(ed->edo," %d: %12.5e ",edi->vecs.linacc.ieig[i],edi->vecs.linacc.xproj[i]);
- fprintf(ed->edo,"\n");
- }
- if (edi->vecs.radfix.neig)
- {
- fprintf(ed->edo," Radfix eigenvectors");
- for (i=0; i<edi->vecs.radfix.neig; i++)
- fprintf(ed->edo," %d: %12.5e ",edi->vecs.radfix.ieig[i],edi->vecs.radfix.xproj[i]);
- fprintf(ed->edo,"\n");
- fprintf(ed->edo," fixed increment radius = %f\n", calc_radius(&edi->vecs.radfix));
- }
- if (edi->vecs.radacc.neig)
- {
- fprintf(ed->edo," Radacc eigenvectors");
- for (i=0; i<edi->vecs.radacc.neig; i++)
- fprintf(ed->edo," %d: %12.5e ",edi->vecs.radacc.ieig[i],edi->vecs.radacc.xproj[i]);
- fprintf(ed->edo,"\n");
- fprintf(ed->edo," acceptance radius = %f\n", calc_radius(&edi->vecs.radacc));
- }
- if (edi->vecs.radcon.neig)
- {
- fprintf(ed->edo," Radcon eigenvectors");
- for (i=0; i<edi->vecs.radcon.neig; i++)
- fprintf(ed->edo," %d: %12.5e ",edi->vecs.radcon.ieig[i],edi->vecs.radcon.xproj[i]);
- fprintf(ed->edo,"\n");
- fprintf(ed->edo," contracting radius = %f\n", calc_radius(&edi->vecs.radcon));
- }
+ for (i=0; i<edi->vecs.linfix.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.linfix.xproj[i]);
+ }
+
+ for (i=0; i<edi->vecs.linacc.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.linacc.xproj[i]);
+ }
+
+ for (i=0; i<edi->vecs.radfix.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.radfix.xproj[i]);
+ }
+ if (edi->vecs.radfix.neig)
+ {
+ fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radfix)); /* fixed increment radius */
+ }
+
+ for (i=0; i<edi->vecs.radacc.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.radacc.xproj[i]);
+ }
+ if (edi->vecs.radacc.neig)
+ {
+ fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radacc)); /* acceptance radius */
+ }
+
+ for (i=0; i<edi->vecs.radcon.neig; i++)
+ {
+ fprintf(fp, EDcol_efmt, edi->vecs.radcon.xproj[i]);
+ }
+ if (edi->vecs.radcon.neig)
+ {
+ fprintf(fp, EDcol_ffmt, calc_radius(&edi->vecs.radcon)); /* contracting radius */
}
}
if (NULL==floodvecs->refproj0)
+ {
snew(floodvecs->refproj0, floodvecs->neig);
+ }
for (i=0; i<floodvecs->neig; i++)
{
}
+/* Call on MASTER only. Check whether the essential dynamics / flooding
+ * groups of the checkpoint file are consistent with the provided .edi file. */
+static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate)
+{
+ t_edpar *edi = NULL; /* points to a single edi data set */
+ int edinum;
+
+
+ if (NULL == EDstate->nref || NULL == EDstate->nav)
+ {
+ gmx_fatal(FARGS, "Essential dynamics and flooding can only be switched on (or off) at the\n"
+ "start of a new simulation. If a simulation runs with/without ED constraints,\n"
+ "it must also continue with/without ED constraints when checkpointing.\n"
+ "To switch on (or off) ED constraints, please prepare a new .tpr to start\n"
+ "from without a checkpoint.\n");
+ }
+
+ edi=ed->edpar;
+ edinum = 0;
+ while(edi != NULL)
+ {
+ /* Check number of atoms in the reference and average structures */
+ if (EDstate->nref[edinum] != edi->sref.nr)
+ {
+ gmx_fatal(FARGS, "The number of reference structure atoms in ED group %c is\n"
+ "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
+ get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr);
+ }
+ if (EDstate->nav[edinum] != edi->sav.nr)
+ {
+ gmx_fatal(FARGS, "The number of average structure atoms in ED group %c is\n"
+ "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n",
+ get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr);
+ }
+ edi=edi->next_edi;
+ edinum++;
+ }
+
+ if (edinum != EDstate->nED)
+ {
+ gmx_fatal(FARGS, "The number of essential dynamics / flooding groups is not consistent.\n"
+ "There are %d ED groups in the .cpt file, but %d in the .edi file!\n"
+ "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum);
+ }
+}
+
+
+/* The edsamstate struct stores the information we need to make the ED group
+ * whole again after restarts from a checkpoint file. Here we do the following:
+ * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
+ * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
+ * c) in any case, for subsequent checkpoint writing, we set the pointers in
+ * edsamstate to the x_old arrays, which contain the correct PBC representation of
+ * all ED structures at the last time step. */
+static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate)
+{
+ int i, nr_edi;
+ t_edpar *edi;
+
+
+ snew(EDstate->old_sref_p, EDstate->nED);
+ snew(EDstate->old_sav_p , EDstate->nED);
+
+ /* If we did not read in a .cpt file, these arrays are not yet allocated */
+ if (!EDstate->bFromCpt)
+ {
+ snew(EDstate->nref, EDstate->nED);
+ snew(EDstate->nav , EDstate->nED);
+ }
+
+ /* Loop over all ED/flooding data sets (usually only one, though) */
+ edi = ed->edpar;
+ for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
+ {
+ /* We always need the last reference and average positions such that
+ * in the next time step we can make the ED group whole again
+ * if the atoms do not have the correct PBC representation */
+ if (EDstate->bFromCpt)
+ {
+ /* Copy the last whole positions of reference and average group from .cpt */
+ for (i=0; i<edi->sref.nr; i++)
+ {
+ copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]);
+ }
+ for (i=0; i<edi->sav.nr ; i++)
+ {
+ copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]);
+ }
+ }
+ else
+ {
+ EDstate->nref[nr_edi-1] = edi->sref.nr;
+ EDstate->nav [nr_edi-1] = edi->sav.nr;
+ }
+
+ /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */
+ EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old;
+ EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old ;
+
+ edi = edi->next_edi;
+ }
+}
+
+
+/* Adds 'buf' to 'str' */
+static void add_to_string(char **str, char *buf)
+{
+ int len;
+
+
+ len = strlen(*str) + strlen(buf) + 1;
+ srenew(*str, len);
+ strcat(*str, buf);
+}
+
+
+static void add_to_string_aligned(char **str, char *buf)
+{
+ char buf_aligned[STRLEN];
+
+ sprintf(buf_aligned, EDcol_sfmt, buf);
+ add_to_string(str, buf_aligned);
+}
+
+
+static void nice_legend(const char ***setname, int *nsets, char **LegendStr, char *value, char *unit, char EDgroupchar)
+{
+ char tmp[STRLEN], tmp2[STRLEN];
+
+
+ sprintf(tmp, "%c %s", EDgroupchar, value);
+ add_to_string_aligned(LegendStr, tmp);
+ sprintf(tmp2, "%s (%s)", tmp, unit);
+ (*setname)[*nsets] = strdup(tmp2);
+ (*nsets)++;
+}
+
+
+static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype)
+{
+ int i;
+ char tmp[STRLEN];
+
+
+ for (i=0; i<evec->neig; i++)
+ {
+ sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype);
+ nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar);
+ }
+}
+
+
+/* Makes a legend for the xvg output file. Call on MASTER only! */
+static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv)
+{
+ t_edpar *edi = NULL;
+ int i;
+ int nr_edi, nsets, n_flood, n_edsam;
+ const char **setname;
+ char buf[STRLEN];
+ char *LegendStr=NULL;
+
+
+ edi = ed->edpar;
+
+ fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s":"");
+
+ for (nr_edi = 1; nr_edi <= nED; nr_edi++)
+ {
+ fprintf(ed->edo, "#\n");
+ fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED));
+ fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr);
+ fprintf(ed->edo, "# monitor : %d vec%s\n" , edi->vecs.mon.neig , edi->vecs.mon.neig != 1 ? "s":"");
+ fprintf(ed->edo, "# LINFIX : %d vec%s\n" , edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s":"");
+ fprintf(ed->edo, "# LINACC : %d vec%s\n" , edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s":"");
+ fprintf(ed->edo, "# RADFIX : %d vec%s\n" , edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s":"");
+ fprintf(ed->edo, "# RADACC : %d vec%s\n" , edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s":"");
+ fprintf(ed->edo, "# RADCON : %d vec%s\n" , edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s":"");
+ fprintf(ed->edo, "# FLOODING : %d vec%s " , edi->flood.vecs.neig , edi->flood.vecs.neig != 1 ? "s":"");
+
+ if (edi->flood.vecs.neig)
+ {
+ /* If in any of the groups we find a flooding vector, flooding is turned on */
+ ed->eEDtype = eEDflood;
+
+ /* Print what flavor of flooding we will do */
+ if (0 == edi->flood.tau) /* constant flooding strength */
+ {
+ fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl);
+ if (edi->flood.bHarmonic)
+ {
+ fprintf(ed->edo, ", harmonic");
+ }
+ }
+ else /* adaptive flooding */
+ {
+ fprintf(ed->edo, ", adaptive");
+ }
+ }
+ fprintf(ed->edo, "\n");
+
+ edi = edi->next_edi;
+ }
+
+ /* Print a nice legend */
+ snew(LegendStr, 1);
+ LegendStr[0] = '\0';
+ sprintf(buf, "# %6s", "time");
+ add_to_string(&LegendStr, buf);
+
+ /* Calculate the maximum number of columns we could end up with */
+ edi = ed->edpar;
+ nsets = 0;
+ for (nr_edi = 1; nr_edi <= nED; nr_edi++)
+ {
+ nsets += 5 +edi->vecs.mon.neig
+ +edi->vecs.linfix.neig
+ +edi->vecs.linacc.neig
+ +edi->vecs.radfix.neig
+ +edi->vecs.radacc.neig
+ +edi->vecs.radcon.neig
+ + 6*edi->flood.vecs.neig;
+ edi = edi->next_edi;
+ }
+ snew(setname, nsets);
+
+ /* In the mdrun time step in a first function call (do_flood()) the flooding
+ * forces are calculated and in a second function call (do_edsam()) the
+ * ED constraints. To get a corresponding legend, we need to loop twice
+ * over the edi groups and output first the flooding, then the ED part */
+
+ /* The flooding-related legend entries, if flooding is done */
+ nsets = 0;
+ if (eEDflood == ed->eEDtype)
+ {
+ edi = ed->edpar;
+ for (nr_edi = 1; nr_edi <= nED; nr_edi++)
+ {
+ /* Always write out the projection on the flooding EVs. Of course, this can also
+ * be achieved with the monitoring option in do_edsam() (if switched on by the
+ * user), but in that case the positions need to be communicated in do_edsam(),
+ * which is not necessary when doing flooding only. */
+ nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
+
+ for (i=0; i<edi->flood.vecs.neig; i++)
+ {
+ sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
+
+ /* Output the current reference projection if it changes with time;
+ * this can happen when flooding is used as harmonic restraint */
+ if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0)
+ {
+ sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED));
+ }
+
+ /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */
+ if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */
+ {
+ sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
+ }
+
+ sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
+
+ if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */
+ {
+ sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED));
+ }
+
+ sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]);
+ nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED));
+ }
+
+ edi = edi->next_edi;
+ } /* End of flooding-related legend entries */
+ }
+ n_flood = nsets;
+
+ /* Now the ED-related entries, if essential dynamics is done */
+ edi = ed->edpar;
+ for (nr_edi = 1; nr_edi <= nED; nr_edi++)
+ {
+ nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) );
+
+ /* Essential dynamics, projections on eigenvectors */
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon , get_EDgroupChar(nr_edi, nED), "MON" );
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX");
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC");
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX");
+ if (edi->vecs.radfix.neig)
+ {
+ nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED));
+ }
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC");
+ if (edi->vecs.radacc.neig)
+ {
+ nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED));
+ }
+ nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON");
+ if (edi->vecs.radcon.neig)
+ {
+ nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED));
+ }
+
+ edi = edi->next_edi;
+ } /* end of 'pure' essential dynamics legend entries */
+ n_edsam = nsets - n_flood;
+
+ xvgr_legend(ed->edo, nsets, setname, oenv);
+ sfree(setname);
+
+ fprintf(ed->edo, "#\n"
+ "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n",
+ n_flood, 1 == n_flood ? "":"s",
+ n_edsam, 1 == n_edsam ? "":"s");
+ fprintf(ed->edo, "%s", LegendStr);
+ sfree(LegendStr);
+
+ fflush(ed->edo);
+}
+
+
void init_edsam(gmx_mtop_t *mtop, /* global topology */
t_inputrec *ir, /* input record */
t_commrec *cr, /* communication record */
gmx_edsam_t ed, /* contains all ED data */
rvec x[], /* positions of the whole MD system */
- matrix box) /* the box */
+ matrix box, /* the box */
+ edsamstate_t *EDstate)
{
t_edpar *edi = NULL; /* points to a single edi data set */
- int numedis=0; /* keep track of the number of ED data sets in edi file */
int i,nr_edi,avindex;
rvec *x_pbc = NULL; /* positions of the whole MD system with pbc removed */
- rvec *xfit = NULL; /* the positions which will be fitted to the reference structure */
- rvec *xstart = NULL; /* the positions which are subject to ED sampling */
+ rvec *xfit=NULL, *xstart=NULL; /* dummy arrays to determine initial RMSDs */
rvec fit_transvec; /* translation ... */
matrix fit_rotmat; /* ... and rotation from fit to reference structure */
if (!DOMAINDECOMP(cr) && PAR(cr) && MASTER(cr))
+ {
gmx_fatal(FARGS, "Please switch on domain decomposition to use essential dynamics in parallel.");
+ }
if (MASTER(cr))
+ {
fprintf(stderr, "ED: Initializing essential dynamics constraints.\n");
+ if (NULL == ed)
+ {
+ gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n"
+ "flooding simulation. Please also provide the correct .edi file with -ei.\n");
+ }
+ }
+
/* Needed for initializing radacc radius in do_edsam */
- ed->bFirst = 1;
+ ed->bFirst = TRUE;
/* The input file is read by the master and the edi structures are
* initialized here. Input is stored in ed->edpar. Then the edi
* structures are transferred to the other nodes */
if (MASTER(cr))
{
- snew(ed->edpar,1);
- /* Read the whole edi file at once: */
- read_edi_file(ed,ed->edpar,mtop->natoms,cr);
-
- /* Initialization for every ED/flooding dataset. Flooding uses one edi dataset per
+ /* Initialization for every ED/flooding group. Flooding uses one edi group per
* flooding vector, Essential dynamics can be applied to more than one structure
* as well, but will be done in the order given in the edi file, so
* expect different results for different order of edi file concatenation! */
edi=ed->edpar;
while(edi != NULL)
{
- init_edi(mtop,ir,cr,ed,edi);
-
- /* Init flooding parameters if needed */
- init_flood(edi,ed,ir->delta_t,cr);
-
+ init_edi(mtop,edi);
+ init_flood(edi,ed,ir->delta_t);
edi=edi->next_edi;
- numedis++;
}
}
/* Reset pointer to first ED data set which contains the actual ED data */
edi=ed->edpar;
-
/* Loop over all ED/flooding data sets (usually only one, though) */
- for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
+ for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
{
- /* We use srenew to allocate memory since the size of the buffers
- * is likely to change with every ED dataset */
- srenew(xfit , edi->sref.nr );
- srenew(xstart, edi->sav.nr );
-
- /* Extract the positions of the atoms to which will be fitted */
- for (i=0; i < edi->sref.nr; i++)
+ /* For multiple ED groups we use the output frequency that was specified
+ * in the first set */
+ if (nr_edi > 1)
{
- copy_rvec(x_pbc[edi->sref.anrs[i]], xfit[i]);
-
- /* Save the sref positions such that in the next time step we can make the ED group whole
- * in case any of the atoms do not have the correct PBC representation */
- copy_rvec(xfit[i], edi->sref.x_old[i]);
+ edi->outfrq = ed->edpar->outfrq;
}
- /* Extract the positions of the atoms subject to ED sampling */
- for (i=0; i < edi->sav.nr; i++)
+ /* Extract the initial reference and average positions. When starting
+ * from .cpt, these have already been read into sref.x_old
+ * in init_edsamstate() */
+ if (!EDstate->bFromCpt)
{
- copy_rvec(x_pbc[edi->sav.anrs[i]], xstart[i]);
+ /* If this is the first run (i.e. no checkpoint present) we assume
+ * that the starting positions give us the correct PBC representation */
+ for (i=0; i < edi->sref.nr; i++)
+ {
+ copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]);
+ }
- /* Save the sav positions such that in the next time step we can make the ED group whole
- * in case any of the atoms do not have the correct PBC representation */
- copy_rvec(xstart[i], edi->sav.x_old[i]);
+ for (i=0; i < edi->sav.nr; i++)
+ {
+ copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]);
+ }
}
+ /* Now we have the PBC-correct start positions of the reference and
+ average structure. We copy that over to dummy arrays on which we
+ can apply fitting to print out the RMSD. We srenew the memory since
+ the size of the buffers is likely different for every ED group */
+ srenew(xfit , edi->sref.nr );
+ srenew(xstart, edi->sav.nr );
+ copy_rvecn(edi->sref.x_old, xfit, 0, edi->sref.nr);
+ copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr);
+
/* Make the fit to the REFERENCE structure, get translation and rotation */
fit_to_reference(xfit, fit_transvec, fit_rotmat, edi);
/* Output how well we fit to the reference at the start */
translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat);
- fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm (dataset #%d)\n",
- rmsd_from_structure(xfit, &edi->sref), nr_edi);
+ fprintf(stderr, "ED: Initial RMSD from reference after fit = %f nm",
+ rmsd_from_structure(xfit, &edi->sref));
+ if (EDstate->nED > 1)
+ {
+ fprintf(stderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED));
+ }
+ fprintf(stderr, "\n");
/* Now apply the translation and rotation to the atoms on which ED sampling will be performed */
translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat);
* the average structure, which must be projected */
avindex = edi->star.nr - edi->sav.nr;
}
- rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon, cr);
- } else
- rad_project(edi, xstart, &edi->vecs.radcon, cr);
+ rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon);
+ }
+ else
+ {
+ rad_project(edi, xstart, &edi->vecs.radcon);
+ }
/* process structure that will serve as origin of expansion circle */
if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
+ {
fprintf(stderr, "ED: Setting center of flooding potential (0 = average structure)\n");
+ }
if (edi->sori.nr > 0)
{
avindex = edi->sori.nr - edi->sav.nr;
}
- rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc, cr);
- rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix, cr);
+ rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc);
+ rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix);
if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
{
fprintf(stderr, "ED: The ORIGIN structure will define the flooding potential center.\n");
/* Set center of flooding potential to the ORIGIN structure */
- rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs, cr);
+ rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs);
/* We already know that no (moving) reference position was provided,
* therefore we can overwrite refproj[0]*/
copyEvecReference(&edi->flood.vecs);
}
else /* No origin structure given */
{
- rad_project(edi, xstart, &edi->vecs.radacc, cr);
- rad_project(edi, xstart, &edi->vecs.radfix, cr);
+ rad_project(edi, xstart, &edi->vecs.radacc);
+ rad_project(edi, xstart, &edi->vecs.radfix);
if ( (eEDflood == ed->eEDtype) && (FALSE == edi->flood.bConstForce) )
{
if (edi->flood.bHarmonic)
{
fprintf(stderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n");
for (i=0; i<edi->flood.vecs.neig; i++)
+ {
edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i];
+ }
}
else
{
/* Set center of flooding potential to the center of the covariance matrix,
* i.e. the average structure, i.e. zero in the projected system */
for (i=0; i<edi->flood.vecs.neig; i++)
+ {
edi->flood.vecs.refproj[i] = 0.0;
+ }
}
}
}
{
for (i=0; i<edi->flood.vecs.neig; i++)
{
- fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", i, edi->flood.vecs.refproj[i]);
+ fprintf(stdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]);
if (edi->flood.bHarmonic)
+ {
fprintf(stdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]);
+ }
fprintf(stdout, "\n");
}
}
/* set starting projections for linsam */
- rad_project(edi, xstart, &edi->vecs.linacc, cr);
- rad_project(edi, xstart, &edi->vecs.linfix, cr);
-
- /* Output to file, set the step to -1 so that write_edo knows it was called from init_edsam */
- if (ed->edo && !(ed->bStartFromCpt))
- write_edo(nr_edi, edi, ed, -1, 0);
+ rad_project(edi, xstart, &edi->vecs.linacc);
+ rad_project(edi, xstart, &edi->vecs.linfix);
/* Prepare for the next edi data set: */
edi=edi->next_edi;
if (PAR(cr))
{
/* First let everybody know how many ED data sets to expect */
- gmx_bcast(sizeof(numedis), &numedis, cr);
+ gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr);
/* Broadcast the essential dynamics / flooding data to all nodes */
- broadcast_ed_data(cr, ed, numedis);
+ broadcast_ed_data(cr, ed, EDstate->nED);
}
else
{
/* Loop over all ED data sets (usually only one, though) */
edi=ed->edpar;
- for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
+ for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
{
edi->sref.anrs_loc = edi->sref.anrs;
edi->sav.anrs_loc = edi->sav.anrs;
snew(edi->sav.c_ind, edi->sav.nr);
/* Initialize the array */
for (i=0; i<edi->sav.nr; i++)
+ {
edi->sav.c_ind[i] = i;
+ }
/* In the general case we will need a different-sized array for the reference indices: */
if (!edi->bRefEqAv)
{
snew(edi->sref.c_ind, edi->sref.nr);
for (i=0; i<edi->sref.nr; i++)
+ {
edi->sref.c_ind[i] = i;
+ }
}
/* Point to the very same array in case of other structures: */
edi->star.c_ind = edi->sav.c_ind;
edi->star.nr_loc = edi->star.nr;
edi->sori.nr_loc = edi->sori.nr;
- /* An on we go to the next edi dataset */
+ /* An on we go to the next ED group */
edi=edi->next_edi;
}
}
/* Allocate space for ED buffer variables */
/* Again, loop over ED data sets */
edi=ed->edpar;
- for (nr_edi = 1; nr_edi <= numedis; nr_edi++)
+ for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++)
{
/* Allocate space for ED buffer */
snew(edi->buf, 1);
dump_edi(edi, cr, nr_edi);
#endif
- /* An on we go to the next edi dataset */
+ /* Next ED group */
edi=edi->next_edi;
}
/* Flush the edo file so that the user can check some things
* when the simulation has started */
if (ed->edo)
+ {
fflush(ed->edo);
+ }
}
void do_edsam(t_inputrec *ir,
gmx_large_int_t step,
- t_mdatoms *md,
t_commrec *cr,
rvec xs[], /* The local current positions on this processor */
rvec v[], /* The velocities */
struct t_do_edsam *buf;
t_edpar *edi;
real rmsdev=-1; /* RMSD from reference structure prior to applying the constraints */
- gmx_bool bSuppress=FALSE; /* Write .edo file on master? */
+ gmx_bool bSuppress=FALSE; /* Write .xvg output file on master? */
/* Check if ED sampling has to be performed */
if ( ed->eEDtype==eEDnone )
+ {
return;
+ }
/* Suppress output on first call of do_edsam if
* two-step sd2 integrator is used */
if ( (ir->eI==eiSD2) && (v != NULL) )
+ {
bSuppress = TRUE;
+ }
dt_1 = 1.0/ir->delta_t;
- /* Loop over all ED datasets (usually one) */
+ /* Loop over all ED groups (usually one) */
edi = ed->edpar;
edinr = 0;
while (edi != NULL)
buf=edi->buf->do_edsam;
if (ed->bFirst)
+ {
/* initialise radacc radius for slope criterion */
buf->oldrad=calc_radius(&edi->vecs.radacc);
+ }
/* Copy the positions into buf->xc* arrays and after ED
* feed back corrections to the official positions */
communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
-#ifdef DEBUG_ED
- dump_xcoll(edi, buf, cr, step);
-#endif
/* Only assembly reference positions if their indices differ from the average ones */
if (!edi->bRefEqAv)
+ {
communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
+ }
/* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
* We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
/* Fit the reference indices to the reference structure */
if (edi->bRefEqAv)
+ {
fit_to_reference(buf->xcoll , transvec, rotmat, edi);
+ }
else
+ {
fit_to_reference(buf->xc_ref, transvec, rotmat, edi);
+ }
/* Now apply the translation and rotation to the ED structure */
translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat);
if (do_per_step(step,edi->maxedsteps) && step >= edi->presteps)
{
project(buf->xcoll, edi);
- rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
- rad_project(edi, buf->xcoll, &edi->vecs.radfix, cr);
+ rad_project(edi, buf->xcoll, &edi->vecs.radacc);
+ rad_project(edi, buf->xcoll, &edi->vecs.radfix);
buf->oldrad=-1.e5;
}
if (edi->vecs.radacc.radius - buf->oldrad < edi->slope)
{
project(buf->xcoll, edi);
- rad_project(edi, buf->xcoll, &edi->vecs.radacc, cr);
+ rad_project(edi, buf->xcoll, &edi->vecs.radacc);
buf->oldrad = 0.0;
- } else
+ }
+ else
+ {
buf->oldrad = edi->vecs.radacc.radius;
+ }
}
/* apply the constraints */
{
/* ED constraints should be applied already in the first MD step
* (which is step 0), therefore we pass step+1 to the routine */
- ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step, cr);
+ ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step);
}
/* write to edo, when required */
{
project(buf->xcoll, edi);
if (MASTER(cr) && !bSuppress)
- write_edo(edinr, edi, ed, step, rmsdev);
+ {
+ write_edo(edi, ed->edo, rmsdev);
+ }
}
/* Copy back the positions unless monitoring only */
}
} /* END of if (edi->bNeedDoEdsam) */
- /* Prepare for the next ED dataset */
+ /* Prepare for the next ED group */
edi = edi->next_edi;
- } /* END of loop over ED datasets */
+ } /* END of loop over ED groups */
ed->bFirst = FALSE;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff