}
static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd,
- gmx_bool bPrintRefValue, double t,
+ gmx_bool bPrintRefValue,
gmx_bool bPrintComponents)
{
- double distance, distance2;
- int m;
-
- if (pcrd->eGeom == epullgDIST)
- {
- /* Geometry=distance: print the length of the distance vector */
- distance2 = 0;
- for (m = 0; m < DIM; m++)
- {
- if (pcrd->dim[m])
- {
- distance2 += pcrd->dr[m]*pcrd->dr[m];
- }
- }
- distance = sqrt(distance2);
- }
- else
- {
- /* Geometry=direction-like: print distance along the pull vector */
- distance = 0;
- for (m = 0; m < DIM; m++)
- {
- if (pcrd->dim[m])
- {
- distance += pcrd->dr[m]*pcrd->vec[m];
- }
- }
- }
- fprintf(out, "\t%g", distance);
+ fprintf(out, "\t%g", pcrd->value);
if (bPrintRefValue)
{
- fprintf(out, "\t%g", pcrd->init + pcrd->rate*t);
+ fprintf(out, "\t%g", pcrd->value_ref);
}
if (bPrintComponents)
{
+ int m;
+
for (m = 0; m < DIM; m++)
{
if (pcrd->dim[m])
pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[1]]);
}
pull_print_coord_dr(out, pcrd,
- pull->bPrintRefValue, t, pull->bPrintComp);
+ pull->bPrintRefValue, pull->bPrintComp);
}
fprintf(out, "\n");
}
*/
static void low_get_pull_coord_dr(const t_pull *pull,
const t_pull_coord *pcrd,
- const t_pbc *pbc, double t,
+ const t_pbc *pbc,
dvec xg, dvec xref, double max_dist2,
dvec dr)
{
{
for (m = 0; m < DIM; m++)
{
- dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m];
+ dref[m] = pcrd->value_ref*pcrd->vec[m];
}
/* Add the reference position, so we use the correct periodic image */
dvec_inc(xrefr, dref);
/* This function returns the distance based on the contents of the pull struct.
* pull->coord[coord_ind].dr, and potentially vector, are updated.
*/
-static void get_pull_coord_dr(const t_pull *pull,
- int coord_ind,
- const t_pbc *pbc, double t)
+static void get_pull_coord_dr(t_pull *pull,
+ int coord_ind,
+ const t_pbc *pbc)
{
double md2;
t_pull_coord *pcrd;
}
}
- low_get_pull_coord_dr(pull, pcrd, pbc, t,
+ low_get_pull_coord_dr(pull, pcrd, pbc,
pull->group[pcrd->group[1]].x,
pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
md2,
pcrd->dr);
}
-void get_pull_coord_distance(const t_pull *pull,
- int coord_ind,
- const t_pbc *pbc, double t,
- double *dev)
+static void update_pull_coord_reference_value(t_pull_coord *pcrd, double t)
{
- static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
- but is fairly benign */
- const t_pull_coord *pcrd;
- int m;
- double ref, drs, inpr;
+ /* With zero rate the reference value is set initially and doesn't change */
+ if (pcrd->rate != 0)
+ {
+ pcrd->value_ref = pcrd->init + pcrd->rate*t;
+ }
+}
+
+/* Calculates pull->coord[coord_ind].value.
+ * This function also updates pull->coord[coord_ind].dr.
+ */
+static void get_pull_coord_distance(t_pull *pull,
+ int coord_ind,
+ const t_pbc *pbc)
+{
+ t_pull_coord *pcrd;
+ int m;
pcrd = &pull->coord[coord_ind];
- get_pull_coord_dr(pull, coord_ind, pbc, t);
+ get_pull_coord_dr(pull, coord_ind, pbc);
- ref = pcrd->init + pcrd->rate*t;
+ switch (pcrd->eGeom)
+ {
+ case epullgDIST:
+ /* Pull along the vector between the com's */
+ pcrd->value = dnorm(pcrd->dr);
+ break;
+ case epullgDIR:
+ case epullgDIRPBC:
+ case epullgDIRRELATIVE:
+ case epullgCYL:
+ /* Pull along vec */
+ pcrd->value = 0;
+ for (m = 0; m < DIM; m++)
+ {
+ pcrd->value += pcrd->vec[m]*pcrd->dr[m];
+ }
+ break;
+ }
+}
+
+/* Returns the deviation from the reference value.
+ * Updates pull->coord[coord_ind].dr, .value and .value_ref.
+ */
+static double get_pull_coord_deviation(t_pull *pull,
+ int coord_ind,
+ const t_pbc *pbc,
+ double t)
+{
+ static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
+ but is fairly benign */
+ t_pull_coord *pcrd;
+ double dev = 0;
+
+ pcrd = &pull->coord[coord_ind];
+
+ get_pull_coord_distance(pull, coord_ind, pbc);
+
+ update_pull_coord_reference_value(pcrd, t);
switch (pcrd->eGeom)
{
case epullgDIST:
/* Pull along the vector between the com's */
- if (ref < 0 && !bWarned)
+ if (pcrd->value_ref < 0 && !bWarned)
{
- fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref);
+ fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, pcrd->value_ref);
bWarned = TRUE;
}
- drs = dnorm(pcrd->dr);
- if (drs == 0)
+ if (pcrd->value == 0)
{
/* With no vector we can not determine the direction for the force,
* so we set the force to zero.
*/
- *dev = 0;
+ dev = 0;
}
else
{
/* Determine the deviation */
- *dev = drs - ref;
+ dev = pcrd->value - pcrd->value_ref;
}
break;
case epullgDIR:
case epullgDIRRELATIVE:
case epullgCYL:
/* Pull along vec */
- inpr = 0;
- for (m = 0; m < DIM; m++)
- {
- inpr += pcrd->vec[m]*pcrd->dr[m];
- }
- *dev = inpr - ref;
+ dev = pcrd->value - pcrd->value_ref;
break;
}
+
+ return dev;
+}
+
+void get_pull_coord_value(t_pull *pull,
+ int coord_ind,
+ const t_pbc *pbc,
+ double *value)
+{
+ get_pull_coord_distance(pull, coord_ind, pbc);
+
+ *value = pull->coord[coord_ind].value;
}
void clear_pull_forces(t_pull *pull)
dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
dvec *rnew; /* current 'new' positions of the groups */
double *dr_tot; /* the total update of the coords */
- double ref;
dvec vec;
double inpr;
double lambda, rm, invdt = 0;
continue;
}
- get_pull_coord_dr(pull, c, pbc, t);
- /* Store the difference vector at time t for printing */
+ /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
+ * We don't modify dr and value anymore, so these values are also used
+ * for printing.
+ */
+ get_pull_coord_distance(pull, c, pbc);
+
if (debug)
{
fprintf(debug, "Pull coord %d dr %f %f %f\n",
continue;
}
+ update_pull_coord_reference_value(pcrd, t);
+
pgrp0 = &pull->group[pcrd->group[0]];
pgrp1 = &pull->group[pcrd->group[1]];
/* Get the current difference vector */
- low_get_pull_coord_dr(pull, pcrd, pbc, t,
+ low_get_pull_coord_dr(pull, pcrd, pbc,
rnew[pcrd->group[1]],
rnew[pcrd->group[0]],
-1, unc_ij);
- ref = pcrd->init + pcrd->rate*t;
-
if (debug)
{
fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
switch (pcrd->eGeom)
{
case epullgDIST:
- if (ref <= 0)
+ if (pcrd->value_ref <= 0)
{
- gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref);
+ gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, pcrd->value_ref);
}
{
c_a = diprod(r_ij[c], r_ij[c]);
c_b = diprod(unc_ij, r_ij[c])*2;
- c_c = diprod(unc_ij, unc_ij) - dsqr(ref);
+ c_c = diprod(unc_ij, unc_ij) - dsqr(pcrd->value_ref);
if (c_b < 0)
{
}
/* Select only the component along the vector */
dsvmul(a, vec, unc_ij);
- lambda = a - ref;
+ lambda = a - pcrd->value_ref;
if (debug)
{
fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
g0 = pcrd->group[0];
g1 = pcrd->group[1];
- low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp);
- low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3);
+ low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
+ low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
fprintf(debug,
"Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
rnew[g0][0], rnew[g0][1], rnew[g0][2],
rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
fprintf(debug,
"Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
- "", "", "", "", "", "", ref);
+ "", "", "", "", "", "", pcrd->value_ref);
fprintf(debug,
"Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
dr0[0], dr0[1], dr0[2],
continue;
}
- ref = pcrd->init + pcrd->rate*t;
-
- low_get_pull_coord_dr(pull, pcrd, pbc, t,
+ low_get_pull_coord_dr(pull, pcrd, pbc,
rnew[pcrd->group[1]],
rnew[pcrd->group[0]],
-1, unc_ij);
switch (pcrd->eGeom)
{
case epullgDIST:
- bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol;
+ bConverged =
+ fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->constr_tol;
break;
case epullgDIR:
case epullgDIRPBC:
inpr = diprod(unc_ij, vec);
dsvmul(inpr, vec, unc_ij);
bConverged =
- fabs(diprod(unc_ij, vec) - ref) < pull->constr_tol;
+ fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->constr_tol;
break;
}
{
fprintf(debug, "NOT CONVERGED YET: Group %d:"
"d_ref = %f, current d = %f\n",
- g, ref, dnorm(unc_ij));
+ g, pcrd->value_ref, dnorm(unc_ij));
}
bConverged_all = FALSE;
continue;
}
- get_pull_coord_distance(pull, c, pbc, t, &dev);
+ dev = get_pull_coord_deviation(pull, c, pbc, t);
k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
dkdl = pcrd->kB - pcrd->k;
{
pull->group[pcrd->group[g]].bCalcCOM = TRUE;
}
+
+ /* With non-zero rate the reference value is set at every step */
+ if (pcrd->rate == 0)
+ {
+ /* Initialize the constant reference value */
+ pcrd->value_ref = pcrd->init;
+ }
}
pull->ePBC = ir->ePBC;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff