real *fproj; /* instantaneous f projections */
real radius; /* instantaneous radius */
real *refproj; /* starting or target projections */
- /* When using flooding as harmonic restraint: The current reference projection
- * is at each step calculated from the initial refproj0 and the slope. */
- real *refproj0, *refprojslope;
} t_eigvec;
real alpha2;
rvec *forces_cartesian;
t_eigvec vecs; /* use flooding for these */
+ /* When using flooding as harmonic restraint: The current reference projection
+ * is at each step calculated from the initial initialReferenceProjection and the slope. */
+ real *initialReferenceProjection;
+ real *referenceProjectionSlope;
} t_edflood;
/* 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.bHarmonic && edi->flood.referenceProjectionSlope[i] != 0.0)
{
fprintf(fp, EDcol_efmt, edi->flood.vecs.refproj[i]);
}
{
for (i = 0; i < edi->flood.vecs.neig; i++)
{
- edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i];
+ edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i] + step * edi->flood.referenceProjectionSlope[i];
}
}
/* Broadcasts the eigenvector data */
-static void bc_ed_vecs(const t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic)
+static void bc_ed_vecs(const t_commrec *cr, t_eigvec *ev, int length)
{
int i;
nblock_bc(cr, length, ev->vec[i]);
}
- /* For harmonic restraints the reference projections can change with time */
- if (bHarmonic)
- {
- snew_bc(cr, ev->refproj0, ev->neig);
- snew_bc(cr, ev->refprojslope, ev->neig);
- nblock_bc(cr, ev->neig, ev->refproj0 );
- nblock_bc(cr, ev->neig, ev->refprojslope);
- }
}
bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */
/* Broadcast eigenvectors */
- bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr, FALSE);
- bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE);
- bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE);
- bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE);
- bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE);
- bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE);
+ bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr);
+ bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr);
+ bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr);
+ bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr);
+ bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr);
+ bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr);
/* Broadcast flooding eigenvectors and, if needed, values for the moving reference */
- bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic);
+ bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr);
+
+ /* For harmonic restraints the reference projections can change with time */
+ if (edi->flood.bHarmonic)
+ {
+ snew_bc(cr, edi->flood.initialReferenceProjection, edi->flood.vecs.neig);
+ snew_bc(cr, edi->flood.referenceProjectionSlope, edi->flood.vecs.neig);
+ nblock_bc(cr, edi->flood.vecs.neig, edi->flood.initialReferenceProjection);
+ nblock_bc(cr, edi->flood.vecs.neig, edi->flood.referenceProjectionSlope);
+ }
+
/* Set the pointer to the next ED group */
if (edi->next_edi)
}
}
-
-static void scan_edvec(FILE *in, int nr, rvec *vec)
+namespace
{
- char line[STRLEN+1];
- int i;
- double x, y, z;
-
-
- for (i = 0; (i < nr); i++)
+/*!\brief Read eigenvector coordinates for multiple eigenvectors from a file.
+ * \param[in] in the file to read from
+ * \param[in] numAtoms the number of atoms/coordinates for the eigenvector
+ * \param[out] vec the eigen-vectors
+ * \param[in] nEig the number of eigenvectors
+ */
+void scan_edvec(FILE *in, int numAtoms, rvec ***vec, int nEig)
+{
+ snew(*vec, nEig);
+ for (int iEigenvector = 0; iEigenvector < nEig; iEigenvector++)
{
- fgets2 (line, STRLEN, in);
- sscanf (line, max_ev_fmt_lelele, &x, &y, &z);
- vec[i][XX] = x;
- vec[i][YY] = y;
- vec[i][ZZ] = z;
+ snew((*vec)[iEigenvector], numAtoms);
+ for (int iAtom = 0; iAtom < numAtoms; iAtom++)
+ {
+ char line[STRLEN+1];
+ double x, y, z;
+ fgets2(line, STRLEN, in);
+ sscanf(line, max_ev_fmt_lelele, &x, &y, &z);
+ (*vec)[iEigenvector][iAtom][XX] = x;
+ (*vec)[iEigenvector][iAtom][YY] = y;
+ (*vec)[iEigenvector][iAtom][ZZ] = z;
+ }
}
-}
-
-static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference)
+}
+/*!\brief Read an essential dynamics eigenvector and corresponding step size.
+ * \param[in] in input file
+ * \param[in] nrAtoms number of atoms/coordinates
+ * \param[out] tvec the eigenvector
+ */
+void read_edvec(FILE *in, int nrAtoms, t_eigvec *tvec)
{
- int i, idum, nscan;
- double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0;
- char line[STRLEN+1];
-
-
tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
- if (tvec->neig > 0)
+ if (tvec->neig <= 0)
{
- snew(tvec->ieig, tvec->neig);
- snew(tvec->stpsz, tvec->neig);
- snew(tvec->vec, tvec->neig);
- snew(tvec->xproj, tvec->neig);
- snew(tvec->fproj, tvec->neig);
- snew(tvec->refproj, tvec->neig);
- if (bReadRefproj)
- {
- snew(tvec->refproj0, tvec->neig);
- snew(tvec->refprojslope, tvec->neig);
- }
-
- for (i = 0; (i < tvec->neig); i++)
- {
- fgets2 (line, STRLEN, in);
- if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */
- {
- nscan = sscanf(line, max_ev_fmt_dlflflf, &idum, &rdum, &refproj_dum, &refprojslope_dum);
- /* Zero out values which were not scanned */
- switch (nscan)
- {
- case 4:
- /* Every 4 values read, including reference position */
- *bHaveReference = TRUE;
- break;
- case 3:
- /* A reference position is provided */
- *bHaveReference = TRUE;
- /* No value for slope, set to 0 */
- refprojslope_dum = 0.0;
- break;
- case 2:
- /* No values for reference projection and slope, set to 0 */
- refproj_dum = 0.0;
- refprojslope_dum = 0.0;
- break;
- default:
- gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
- }
- tvec->refproj[i] = refproj_dum;
- tvec->refproj0[i] = refproj_dum;
- tvec->refprojslope[i] = refprojslope_dum;
- }
- else /* Normal flooding */
- {
- nscan = sscanf(line, max_ev_fmt_dlf, &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;
- } /* end of loop over eigenvectors */
+ return;
+ }
- for (i = 0; (i < tvec->neig); i++)
+ snew(tvec->ieig, tvec->neig);
+ snew(tvec->stpsz, tvec->neig);
+ for (int i = 0; i < tvec->neig; i++)
+ {
+ char line[STRLEN+1];
+ fgets2(line, STRLEN, in);
+ int numEigenVectors;
+ double stepSize;
+ const int nscan = sscanf(line, max_ev_fmt_dlf, &numEigenVectors, &stepSize);
+ if (nscan != 2)
{
- snew(tvec->vec[i], nr);
- scan_edvec(in, nr, tvec->vec[i]);
+ gmx_fatal(FARGS, "Expected 2 values for flooding vec: <nr> <stpsz>\n");
}
- }
+ tvec->ieig[i] = numEigenVectors;
+ tvec->stpsz[i] = stepSize;
+ } /* end of loop over eigenvectors */
+
+ scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
}
+/*!\brief Read an essential dynamics eigenvector and intial reference projections and slopes if available.
+ *
+ * Eigenvector and its intial reference and slope are stored on the
+ * same line in the .edi format. To avoid re-winding the .edi file,
+ * the reference eigen-vector and reference data are read in one go.
+ *
+ * \param[in] in input file
+ * \param[in] nrAtoms number of atoms/coordinates
+ * \param[out] tvec the eigenvector
+ * \param[out] initialReferenceProjection The projection onto eigenvectors as read from file.
+ * \param[out] referenceProjectionSlope The slope of the reference projections.
+ */
+bool readEdVecWithReferenceProjection(FILE *in, int nrAtoms, t_eigvec *tvec, real ** initialReferenceProjection, real ** referenceProjectionSlope)
+{
+ bool providesReference = false;
+ tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS");
+ if (tvec->neig <= 0)
+ {
+ return false;
+ }
+ snew(tvec->ieig, tvec->neig);
+ snew(tvec->stpsz, tvec->neig);
+ snew(*initialReferenceProjection, tvec->neig);
+ snew(*referenceProjectionSlope, tvec->neig);
+ for (int i = 0; i < tvec->neig; i++)
+ {
+ char line[STRLEN+1];
+ fgets2 (line, STRLEN, in);
+ int numEigenVectors;
+ double stepSize = 0.;
+ double referenceProjection = 0.;
+ double referenceSlope = 0.;
+ const int nscan = sscanf(line, max_ev_fmt_dlflflf, &numEigenVectors, &stepSize, &referenceProjection, &referenceSlope);
+ /* Zero out values which were not scanned */
+ switch (nscan)
+ {
+ case 4:
+ /* Every 4 values read, including reference position */
+ providesReference = true;
+ break;
+ case 3:
+ /* A reference position is provided */
+ providesReference = true;
+ /* No value for slope, set to 0 */
+ referenceSlope = 0.0;
+ break;
+ case 2:
+ /* No values for reference projection and slope, set to 0 */
+ referenceProjection = 0.0;
+ referenceSlope = 0.0;
+ break;
+ default:
+ gmx_fatal(FARGS, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan);
+ }
+ (*initialReferenceProjection)[i] = referenceProjection;
+ (*referenceProjectionSlope)[i] = referenceSlope;
-/* calls read_edvec for the vector groups, only for flooding there is an extra call */
-static void read_edvecs(FILE *in, int nr, t_edvecs *vecs)
-{
- gmx_bool bHaveReference = FALSE;
+ tvec->ieig[i] = numEigenVectors;
+ tvec->stpsz[i] = stepSize;
+ } /* end of loop over eigenvectors */
+ scan_edvec(in, nrAtoms, &tvec->vec, tvec->neig);
+ return providesReference;
+}
- read_edvec(in, nr, &vecs->mon, FALSE, &bHaveReference);
- read_edvec(in, nr, &vecs->linfix, FALSE, &bHaveReference);
- read_edvec(in, nr, &vecs->linacc, FALSE, &bHaveReference);
- read_edvec(in, nr, &vecs->radfix, FALSE, &bHaveReference);
- read_edvec(in, nr, &vecs->radacc, FALSE, &bHaveReference);
- read_edvec(in, nr, &vecs->radcon, FALSE, &bHaveReference);
+/*!\brief Allocate working memory for the eigenvectors.
+ * \param[in,out] tvec eigenvector for which memory will be allocated
+ */
+void setup_edvec(t_eigvec *tvec)
+{
+ snew(tvec->xproj, tvec->neig);
+ snew(tvec->fproj, tvec->neig);
+ snew(tvec->refproj, tvec->neig);
+}
}
{
edi->sori.sqrtm = nullptr;
}
+ setup_edvec(&edi->vecs.linacc);
+ setup_edvec(&edi->vecs.mon);
+ setup_edvec(&edi->vecs.linfix);
+ setup_edvec(&edi->vecs.linacc);
+ setup_edvec(&edi->vecs.radfix);
+ setup_edvec(&edi->vecs.radacc);
+ setup_edvec(&edi->vecs.radcon);
+ setup_edvec(&edi->flood.vecs);
}
/*!\brief Check if edi file version supports CONST_FORCE_FLOODING.
*/
void read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, bool hasConstForceFlooding, const char *fn)
{
- /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
- gmx_bool bHaveReference = FALSE;
-
bool bEOF;
/* check the number of atoms */
edi->nini = read_edint(in, &bEOF);
edi->bRefEqAv = check_if_same(edi->sref, edi->sav);
/* eigenvectors */
- read_edvecs(in, edi->sav.nr, &edi->vecs);
- read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference);
+
+ read_edvec(in, edi->sav.nr, &edi->vecs.mon);
+ read_edvec(in, edi->sav.nr, &edi->vecs.linfix);
+ read_edvec(in, edi->sav.nr, &edi->vecs.linacc);
+ read_edvec(in, edi->sav.nr, &edi->vecs.radfix);
+ read_edvec(in, edi->sav.nr, &edi->vecs.radacc);
+ read_edvec(in, edi->sav.nr, &edi->vecs.radcon);
+
+ /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */
+ bool bHaveReference = false;
+ if (edi->flood.bHarmonic)
+ {
+ bHaveReference = readEdVecWithReferenceProjection(in, edi->sav.nr, &edi->flood.vecs, &edi->flood.initialReferenceProjection, &edi->flood.referenceProjectionSlope);
+ }
+ else
+ {
+ read_edvec(in, edi->sav.nr, &edi->flood.vecs);
+ }
/* target positions */
edi->star.nr = read_edint(in, &bEOF);
}
-/* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/
+/* Copies reference projection 'refproj' to fixed 'initialReferenceProjection' variable for flooding/
* umbrella sampling simulations. */
-static void copyEvecReference(t_eigvec* floodvecs)
+static void copyEvecReference(t_eigvec* floodvecs, real * initialReferenceProjection)
{
int i;
- if (nullptr == floodvecs->refproj0)
+ if (nullptr == initialReferenceProjection)
{
- snew(floodvecs->refproj0, floodvecs->neig);
+ snew(initialReferenceProjection, floodvecs->neig);
}
for (i = 0; i < floodvecs->neig; i++)
{
- floodvecs->refproj0[i] = floodvecs->refproj[i];
+ initialReferenceProjection[i] = floodvecs->refproj[i];
}
}
/* 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)
+ if (edi->flood.bHarmonic && edi->flood.referenceProjectionSlope[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));
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);
+ copyEvecReference(&edi->flood.vecs, edi->flood.initialReferenceProjection);
}
}
else /* No origin structure given */
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];
+ edi->flood.vecs.refproj[i] = edi->flood.initialReferenceProjection[i];
}
}
else
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, " (adding %11.4e/timestep)", edi->flood.referenceProjectionSlope[i]);
}
fprintf(stdout, "\n");
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff