#include <algorithm>
-#include "gromacs/fileio/filenm.h"
+#include "gromacs/commandline/filenm.h"
+#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/domdec/ga2la.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/xvgr.h"
-#include "gromacs/legacyheaders/copyrite.h"
-#include "gromacs/legacyheaders/gmx_ga2la.h"
-#include "gromacs/legacyheaders/macros.h"
-#include "gromacs/legacyheaders/mdrun.h"
-#include "gromacs/legacyheaders/names.h"
-#include "gromacs/legacyheaders/network.h"
-#include "gromacs/legacyheaders/typedefs.h"
-#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/math/functions.h"
#include "gromacs/math/vec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/mdlib/mdrun.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull_internal.h"
#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/pleasecite.h"
+#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
static std::string append_before_extension(std::string pathname,
}
}
-static FILE *open_pull_out(const char *fn, struct pull_t *pull, const output_env_t oenv,
+static FILE *open_pull_out(const char *fn, struct pull_t *pull,
+ const gmx_output_env_t *oenv,
gmx_bool bCoord, unsigned long Flags)
{
FILE *fp;
}
/* Apply forces in a mass weighted fashion */
-static void apply_forces(struct pull_t * pull, t_mdatoms * md, rvec *f)
+static void apply_forces_coord(struct pull_t * pull, int coord,
+ const t_mdatoms * md,
+ rvec *f)
{
- int c;
+ const pull_coord_work_t *pcrd;
- for (c = 0; c < pull->ncoord; c++)
+ pcrd = &pull->coord[coord];
+
+ if (pcrd->params.eGeom == epullgCYL)
{
- const pull_coord_work_t *pcrd;
+ dvec f_tot;
+ int m;
- pcrd = &pull->coord[c];
+ apply_forces_cyl_grp(&pull->dyna[coord], pcrd->cyl_dev, md,
+ pcrd->f, pcrd->f_scal, -1, f);
- if (pcrd->params.eGeom == epullgCYL)
+ /* Sum the force along the vector and the radial force */
+ for (m = 0; m < DIM; m++)
{
- dvec f_tot;
- int m;
-
- apply_forces_cyl_grp(&pull->dyna[c], pcrd->cyl_dev, md,
- pcrd->f, pcrd->f_scal, -1, f);
-
- /* Sum the force along the vector and the radial force */
- for (m = 0; m < DIM; m++)
- {
- f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
- }
- apply_forces_grp(&pull->group[pcrd->params.group[1]], md, f_tot, 1, f);
+ f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
}
- else
+ apply_forces_grp(&pull->group[pcrd->params.group[1]], md, f_tot, 1, f);
+ }
+ else
+ {
+ if (pcrd->params.eGeom == epullgDIRRELATIVE)
{
- if (pcrd->params.eGeom == epullgDIRRELATIVE)
- {
- /* We need to apply the torque forces to the pull groups
- * that define the pull vector.
- */
- apply_forces_vec_torque(pull, pcrd, md, f);
- }
+ /* We need to apply the torque forces to the pull groups
+ * that define the pull vector.
+ */
+ apply_forces_vec_torque(pull, pcrd, md, f);
+ }
- if (pull->group[pcrd->params.group[0]].params.nat > 0)
- {
- apply_forces_grp(&pull->group[pcrd->params.group[0]], md, pcrd->f, -1, f);
- }
- apply_forces_grp(&pull->group[pcrd->params.group[1]], md, pcrd->f, 1, f);
+ if (pull->group[pcrd->params.group[0]].params.nat > 0)
+ {
+ apply_forces_grp(&pull->group[pcrd->params.group[0]], md, pcrd->f, -1, f);
}
+ apply_forces_grp(&pull->group[pcrd->params.group[1]], md, pcrd->f, 1, f);
}
}
{
gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f).\nYou might want to consider using \"pull-geometry = direction-periodic\" instead.\n",
pcrd->params.group[0], pcrd->params.group[1],
- sqrt(dr2), sqrt(max_dist2));
+ sqrt(dr2), sqrt(0.98*0.98*max_dist2));
}
if (pcrd->params.eGeom == epullgDIRPBC)
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(pcrd->value_ref);
+ c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
if (c_b < 0)
{
- q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
+ q = -0.5*(c_b - std::sqrt(c_b*c_b - 4*c_a*c_c));
lambda = -q/c_a;
}
else
{
- q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
+ q = -0.5*(c_b + std::sqrt(c_b*c_b - 4*c_a*c_c));
lambda = -c_c/q;
}
sfree(rnew);
}
-/* Pulling with a harmonic umbrella potential or constant force */
-static void do_pull_pot(struct pull_t *pull, t_pbc *pbc, double t, real lambda,
- real *V, tensor vir, real *dVdl)
+static void calc_pull_coord_force(pull_coord_work_t *pcrd,
+ double dev, real lambda,
+ real *V, tensor vir, real *dVdl)
{
- int c, j, m;
- double dev, ndr, invdr = 0;
+ int j, m;
+ double ndr, invdr = 0;
real k, dkdl;
- /* loop over the pull coordinates */
- *V = 0;
- *dVdl = 0;
- for (c = 0; c < pull->ncoord; c++)
- {
- pull_coord_work_t *pcrd;
-
- pcrd = &pull->coord[c];
+ k = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
+ dkdl = pcrd->params.kB - pcrd->params.k;
- if (pcrd->params.eType == epullCONSTRAINT)
+ if (pcrd->params.eGeom == epullgDIST)
+ {
+ ndr = dnorm(pcrd->dr);
+ if (ndr > 0)
{
- continue;
+ invdr = 1/ndr;
+ }
+ else
+ {
+ /* With an harmonic umbrella, the force is 0 at r=0,
+ * so we can set invdr to any value.
+ * With a constant force, the force at r=0 is not defined,
+ * so we zero it (this is anyhow a very rare event).
+ */
+ invdr = 0;
+ }
+ }
+ else
+ {
+ ndr = 0;
+ for (m = 0; m < DIM; m++)
+ {
+ ndr += pcrd->vec[m]*pcrd->dr[m];
}
+ }
- dev = get_pull_coord_deviation(pull, c, pbc, t);
+ switch (pcrd->params.eType)
+ {
+ case epullUMBRELLA:
+ case epullFLATBOTTOM:
+ /* The only difference between an umbrella and a flat-bottom
+ * potential is that a flat-bottom is zero below zero.
+ */
+ if (pcrd->params.eType == epullFLATBOTTOM && dev < 0)
+ {
+ dev = 0;
+ }
- k = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
- dkdl = pcrd->params.kB - pcrd->params.k;
+ pcrd->f_scal = -k*dev;
+ *V += 0.5* k*gmx::square(dev);
+ *dVdl += 0.5*dkdl*gmx::square(dev);
+ break;
+ case epullCONST_F:
+ pcrd->f_scal = -k;
+ *V += k*ndr;
+ *dVdl += dkdl*ndr;
+ break;
+ default:
+ gmx_incons("Unsupported pull type in do_pull_pot");
+ }
- if (pcrd->params.eGeom == epullgDIST)
+ if (pcrd->params.eGeom == epullgDIST)
+ {
+ for (m = 0; m < DIM; m++)
{
- ndr = dnorm(pcrd->dr);
- if (ndr > 0)
- {
- invdr = 1/ndr;
- }
- else
- {
- /* With an harmonic umbrella, the force is 0 at r=0,
- * so we can set invdr to any value.
- * With a constant force, the force at r=0 is not defined,
- * so we zero it (this is anyhow a very rare event).
- */
- invdr = 0;
- }
+ pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
}
- else
+ }
+ else
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
+ }
+ }
+
+ if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC)
+ {
+ /* Add the pull contribution to the virial */
+ for (j = 0; j < DIM; j++)
{
- ndr = 0;
for (m = 0; m < DIM; m++)
{
- ndr += pcrd->vec[m]*pcrd->dr[m];
+ vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
}
}
+ }
+}
- switch (pcrd->params.eType)
- {
- case epullUMBRELLA:
- case epullFLATBOTTOM:
- /* The only difference between an umbrella and a flat-bottom
- * potential is that a flat-bottom is zero below zero.
- */
- if (pcrd->params.eType == epullFLATBOTTOM && dev < 0)
- {
- dev = 0;
- }
+void set_pull_coord_reference_value(struct pull_t *pull,
+ int coord_ind, real value_ref,
+ const struct t_pbc *pbc,
+ const t_mdatoms *md,
+ real lambda,
+ gmx_bool bUpdateForce, rvec *f, tensor vir)
+{
+ pull_coord_work_t *pcrd;
- pcrd->f_scal = -k*dev;
- *V += 0.5* k*dsqr(dev);
- *dVdl += 0.5*dkdl*dsqr(dev);
- break;
- case epullCONST_F:
- pcrd->f_scal = -k;
- *V += k*ndr;
- *dVdl += dkdl*ndr;
- break;
- default:
- gmx_incons("Unsupported pull type in do_pull_pot");
- }
+ pcrd = &pull->coord[coord_ind];
+
+ if (pcrd->params.rate != 0)
+ {
+ gmx_incons("Can not update the reference value for pull coordinates with rate!=0");
+ }
+
+ /* Update the reference value */
+ pcrd->value_ref = value_ref;
+
+ if (bUpdateForce)
+ {
+ real V = 0, dVdl = 0;
+ double f_scal_old;
+ dvec f_old;
+ double dev;
+ int j, m;
- if (pcrd->params.eGeom == epullgDIST)
+ if (pcrd->params.eType == epullCONSTRAINT)
{
- for (m = 0; m < DIM; m++)
- {
- pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
- }
+ gmx_incons("Can not update the pull force for constraint coordinates");
}
- else
+
+ /* The time parameter is not used here, so we pass 0.0 */
+ dev = get_pull_coord_deviation(pull, coord_ind, pbc, 0.0);
+
+ f_scal_old = pcrd->f_scal;
+ copy_dvec(pcrd->f, f_old);
+
+ /* Calculate the new forces, ingnore V, vir and dVdl */
+ calc_pull_coord_force(pcrd, dev, lambda, &V, NULL, &dVdl);
+
+ /* Here we determine the force correction:
+ * substract the half step contribution we added already,
+ * add the half step contribution for the new force.
+ * Note that the printing of f_scal is done in pull_potential, which
+ * might be called before this function is called, so the output might
+ * contain values without the correction below.
+ */
+ pcrd->f_scal = 0.5*(-f_scal_old + pcrd->f_scal);
+ for (m = 0; m < DIM; m++)
{
- for (m = 0; m < DIM; m++)
- {
- pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
- }
+ pcrd->f[m] = 0.5*(-f_old[m] + pcrd->f[m]);
}
if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC)
}
}
}
+
+ apply_forces_coord(pull, coord_ind, md, f);
}
}
+/* Calculate the pull potential and scalar force for a pull coordinate */
+static void do_pull_pot_coord(struct pull_t *pull, int coord_ind, t_pbc *pbc,
+ double t, real lambda,
+ real *V, tensor vir, real *dVdl)
+{
+ pull_coord_work_t *pcrd;
+ double dev;
+
+ pcrd = &pull->coord[coord_ind];
+
+ assert(pcrd->params.eType != epullCONSTRAINT);
+
+ dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
+
+ calc_pull_coord_force(pcrd, dev, lambda, V, vir, dVdl);
+}
+
real pull_potential(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
t_commrec *cr, double t, real lambda,
rvec *x, rvec *f, tensor vir, real *dvdlambda)
{
- real V = 0, dVdl;
+ real V = 0;
assert(pull != NULL);
if (pull->comm.bParticipate)
{
+ real dVdl = 0;
+
pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
- do_pull_pot(pull, pbc, t, lambda,
- &V, MASTER(cr) ? vir : NULL, &dVdl);
+ for (int c = 0; c < pull->ncoord; c++)
+ {
+ if (pull->coord[c].params.eType == epullCONSTRAINT)
+ {
+ continue;
+ }
+
+ do_pull_pot_coord(pull, c, pbc, t, lambda,
+ &V, MASTER(cr) ? vir : NULL, &dVdl);
- /* Distribute forces over the pull groups */
- apply_forces(pull, md, f);
+ /* Distribute the force over the atoms in the pulled groups */
+ apply_forces_coord(pull, c, md, f);
+ }
if (MASTER(cr))
{
}
}
-static void make_local_pull_group(gmx_ga2la_t ga2la,
+static void make_local_pull_group(gmx_ga2la_t *ga2la,
pull_group_work_t *pg, int start, int end)
{
int i, ii;
void dd_make_local_pull_groups(t_commrec *cr, struct pull_t *pull, t_mdatoms *md)
{
- gmx_domdec_t *dd;
- pull_comm_t *comm;
- gmx_ga2la_t ga2la;
- gmx_bool bMustParticipate;
- int g;
+ gmx_domdec_t *dd;
+ pull_comm_t *comm;
+ gmx_ga2la_t *ga2la;
+ gmx_bool bMustParticipate;
+ int g;
dd = cr->dd;
struct pull_t *
init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
int nfile, const t_filenm fnm[],
- gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
+ gmx_mtop_t *mtop, t_commrec *cr,
+ const gmx_output_env_t *oenv, real lambda,
gmx_bool bOutFile, unsigned long Flags)
{
struct pull_t *pull;
copy_rvec(pull_params->coord[c].vec, pcrd->vec);
break;
default:
- gmx_incons("Pull geometry not handled");
+ /* We allow reading of newer tpx files with new pull geometries,
+ * but with the same tpx format, with old code. A new geometry
+ * only adds a new enum value, which will be out of range for
+ * old code. The first place we need to generate an error is
+ * here, since the pull code can't handle this.
+ * The enum can be used before arriving here only for printing
+ * the string corresponding to the geometry, which will result
+ * in printing "UNDEFINED".
+ */
+ gmx_fatal(FARGS, "Pull geometry not supported for pull coordinate %d. The geometry enum %d in the input is larger than that supported by the code (up to %d). You are probably reading a tpr file generated with a newer version of Gromacs with an binary from an older version of Gromacs.",
+ c + 1, pcrd->params.eGeom, epullgNR - 1);
}
if (pcrd->params.eType == epullCONSTRAINT)
* when it is not only used as a cylinder group.
*/
calc_com_start = (pcrd->params.eGeom == epullgCYL ? 1 : 0);
- calc_com_end = (pcrd->params.eGeom == epullgDIRRELATIVE ? 4 : 2);
+ calc_com_end = pcrd->params.ngroup;
for (g = calc_com_start; g <= calc_com_end; g++)
{
clear_ivec(pulldim_con);
for (c = 0; c < pull->ncoord; c++)
{
- if (pull->coord[c].params.group[0] == g ||
- pull->coord[c].params.group[1] == g ||
- (pull->coord[c].params.eGeom == epullgDIRRELATIVE &&
- (pull->coord[c].params.group[2] == g ||
- pull->coord[c].params.group[3] == g)))
+ pull_coord_work_t *pcrd;
+ int gi;
+ gmx_bool bGroupUsed;
+
+ pcrd = &pull->coord[c];
+
+ bGroupUsed = FALSE;
+ for (gi = 0; gi < pcrd->params.ngroup; gi++)
+ {
+ if (pcrd->params.group[gi] == g)
+ {
+ bGroupUsed = TRUE;
+ }
+ }
+
+ if (bGroupUsed)
{
for (m = 0; m < DIM; m++)
{
- if (pull->coord[c].params.dim[m] == 1)
+ if (pcrd->params.dim[m] == 1)
{
pulldim[m] = 1;
- if (pull->coord[c].params.eType == epullCONSTRAINT)
+ if (pcrd->params.eType == epullCONSTRAINT)
{
bConstraint = TRUE;
pulldim_con[m] = 1;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff