Converted the opaque C struct to a pimpl-ed C++ class.
Numerous callers of constraint routines now don't have to pass
parameters that are embedded within the class at setup time,
e.g. for logging, communication, per-atom information,
performance counters.
Some of those parameters have been converted to use const references
per style, which requires various callers and callees to be modified
accordingly. In particular, the mtop utility functions that take const
pointers have been deprecated, and some temporary wrapper functions
used so that we can defer the update of code completely unrelated to
constraints until another time. Similarly, t_commrec is retained as a
pointer, since it also makes sense to change globally.
Made ConstraintVariable an enum class. This generates some compiler
warnings to force us to cover potential bug cases with fatal errors.
Used more complete names for some of the enum members.
Introduced a factory function to continue the design that constr is
nullptr when we're not using constraints.
Added some const correctness where it now became necessary.
Refs #2423
Change-Id: I7a3833489b675f30863ca37c0359cd3e950b5494
comm->bInterCGMultiBody = FALSE;
}
- dd->bInterCGcons = inter_charge_group_constraints(mtop);
- dd->bInterCGsettles = inter_charge_group_settles(mtop);
+ dd->bInterCGcons = inter_charge_group_constraints(*mtop);
+ dd->bInterCGsettles = inter_charge_group_settles(*mtop);
if (ir->rlist == 0)
{
dd_pme_maxshift_x(dd), dd_pme_maxshift_y(dd));
}
+ // TODO constr should probably be a member of the object that
+ // manages work on a PP rank. Or this and similar code should be
+ // structured as callbacks.
if (constr)
{
- set_constraints(constr, top_local, ir, mdatoms, cr);
+ constr->setConstraints(*top_local, *mdatoms);
}
if (ir->bPull)
ilc_local->nr = 0;
if (dd->constraint_comm)
{
- at2con_mt = atom2constraints_moltype(constr);
+ // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
+ GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
+ at2con_mt = constr->atom2constraints_moltype();
ireq = &dd->constraint_comm->ireq[0];
ireq->n = 0;
}
if (dd->bInterCGsettles)
{
- at2settle_mt = atom2settle_moltype(constr);
+ // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
+ GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
+ at2settle_mt = constr->atom2settle_moltype();
ils_local->nr = 0;
}
else
/* Open input and output files, allocate space for ED data structure */
gmx_edsam_t ed = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr);
- saveEdsamPointer(constr, ed);
+ GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
+ constr->saveEdsamPointer(ed);
/* Needed for initializing radacc radius in do_edsam */
ed->bFirst = TRUE;
namespace gmx
{
-class Constraints
+/* \brief Impl class for Constraints
+ *
+ * \todo Members like md, idef are valid only for the lifetime of a
+ * domain, which would be good to make clearer in the structure of the
+ * code. It should not be possible to call apply() if setConstraints()
+ * has not been called. For example, this could be achieved if
+ * setConstraints returned a valid object with such a method. That
+ * still requires that we manage the lifetime of that object
+ * correctly, however. */
+class Constraints::Impl
{
public:
- int ncon_tot; /* The total number of constraints */
- int nflexcon; /* The number of flexible constraints */
- int n_at2con_mt; /* The size of at2con = #moltypes */
- t_blocka *at2con_mt; /* A list of atoms to constraints */
- int n_at2settle_mt; /* The size of at2settle = #moltypes */
- int **at2settle_mt; /* A list of atoms to settles */
- bool bInterCGsettles;
- Lincs *lincsd; /* LINCS data */
- shakedata *shaked; /* SHAKE data */
- settledata *settled; /* SETTLE data */
- int maxwarn; /* The maximum number of warnings */
- int warncount_lincs;
- int warncount_settle;
- gmx_edsam_t ed; /* The essential dynamics data */
-
- /* Thread local working data */
- tensor *vir_r_m_dr_th; /* Thread virial contribution */
- bool *bSettleErrorHasOccurred; /* Did a settle error occur? */
-
- /* Only used for printing warnings */
- const gmx_mtop_t *warn_mtop; /* Pointer to the global topology */
+ Impl(const gmx_mtop_t &mtop_p,
+ const t_inputrec &ir_p,
+ FILE *log_p,
+ const t_mdatoms &md_p,
+ const t_commrec *cr_p,
+ const gmx_multisim_t &ms,
+ t_nrnb *nrnb,
+ gmx_wallcycle *wcycle_p,
+ bool pbcHandlingRequired,
+ int numConstraints,
+ int numSettles);
+ void setConstraints(const gmx_localtop_t &top,
+ const t_mdatoms &md);
+ bool apply(bool bLog,
+ bool bEner,
+ gmx_int64_t step,
+ int delta_step,
+ real step_scaling,
+ rvec *x,
+ rvec *xprime,
+ rvec *min_proj,
+ matrix box,
+ real lambda,
+ real *dvdlambda,
+ rvec *v,
+ tensor *vir,
+ ConstraintVariable econq);
+ //! The total number of constraints.
+ int ncon_tot = 0;
+ //! The number of flexible constraints.
+ int nflexcon = 0;
+ //! The size of at2con = number of moltypes.
+ int n_at2con_mt = 0;
+ //! A list of atoms to constraints.
+ t_blocka *at2con_mt = nullptr;
+ //! The size of at2settle = number of moltypes
+ int n_at2settle_mt = 0;
+ //! A list of atoms to settles.
+ int **at2settle_mt = nullptr;
+ //! Whether any SETTLES cross charge-group boundaries.
+ bool bInterCGsettles = false;
+ //! LINCS data.
+ Lincs *lincsd = nullptr;
+ //! SHAKE data.
+ shakedata *shaked = nullptr;
+ //! SETTLE data.
+ settledata *settled = nullptr;
+ //! The maximum number of warnings.
+ int maxwarn = 0;
+ //! The number of warnings for LINCS.
+ int warncount_lincs = 0;
+ //! The number of warnings for SETTLE.
+ int warncount_settle = 0;
+ //! The essential dynamics data.
+ gmx_edsam_t ed = nullptr;
+
+ //! Thread-local virial contribution.
+ tensor *vir_r_m_dr_th = {0};
+ //! Did a settle error occur?
+ bool *bSettleErrorHasOccurred = nullptr;
+
+ //! Pointer to the global topology - only used for printing warnings.
+ const gmx_mtop_t &mtop;
+ //! Parameters for the interactions in this domain.
+ const t_idef *idef = nullptr;
+ //! Data about atoms in this domain.
+ const t_mdatoms &md;
+ //! Whether we need to do pbc for handling bonds.
+ bool pbcHandlingRequired_ = false;
+
+ //! Logging support.
+ FILE *log = nullptr;
+ //! Communication support.
+ const t_commrec *cr = nullptr;
+ //! Multi-sim support.
+ const gmx_multisim_t &ms;
+ /*!\brief Input options.
+ *
+ * \todo Replace with IMdpOptions */
+ const t_inputrec &ir;
+ //! Flop counting support.
+ t_nrnb *nrnb = nullptr;
+ //! Tracks wallcycle usage.
+ gmx_wallcycle *wcycle;
};
-int n_flexible_constraints(const Constraints *constr)
-{
- int nflexcon;
-
- if (constr)
- {
- nflexcon = constr->nflexcon;
- }
- else
- {
- nflexcon = 0;
- }
+Constraints::~Constraints() = default;
- return nflexcon;
+int Constraints::numFlexibleConstraints() const
+{
+ return impl_->nflexcon;
}
//! Clears constraint quantities for atoms in nonlocal region.
//! Writes out coordinates.
static void write_constr_pdb(const char *fn, const char *title,
- const gmx_mtop_t *mtop,
+ const gmx_mtop_t &mtop,
int start, int homenr, const t_commrec *cr,
- rvec x[], matrix box)
+ const rvec x[], matrix box)
{
char fname[STRLEN];
FILE *out;
}
//! Writes out domain contents to help diagnose crashes.
-static void dump_confs(FILE *fplog, gmx_int64_t step, const gmx_mtop_t *mtop,
+static void dump_confs(FILE *log, gmx_int64_t step, const gmx_mtop_t &mtop,
int start, int homenr, const t_commrec *cr,
- rvec x[], rvec xprime[], matrix box)
+ const rvec x[], rvec xprime[], matrix box)
{
char buf[STRLEN], buf2[22];
sprintf(buf, "step%sc", gmx_step_str(step, buf2));
write_constr_pdb(buf, "coordinates after constraining",
mtop, start, homenr, cr, xprime, box);
- if (fplog)
+ if (log)
{
- fprintf(fplog, "Wrote pdb files with previous and current coordinates\n");
+ fprintf(log, "Wrote pdb files with previous and current coordinates\n");
}
fprintf(stderr, "Wrote pdb files with previous and current coordinates\n");
}
-bool constrain(FILE *fplog, bool bLog, bool bEner,
- Constraints *constr,
- const t_idef *idef, const t_inputrec *ir,
- const t_commrec *cr,
- const gmx_multisim_t *ms,
- gmx_int64_t step, int delta_step,
- real step_scaling,
- const t_mdatoms *md,
- rvec *x, rvec *xprime, rvec *min_proj,
- bool bMolPBC, matrix box,
- real lambda, real *dvdlambda,
- rvec *v, tensor *vir,
- t_nrnb *nrnb, int econq)
+bool
+Constraints::apply(bool bLog,
+ bool bEner,
+ gmx_int64_t step,
+ int delta_step,
+ real step_scaling,
+ rvec *x,
+ rvec *xprime,
+ rvec *min_proj,
+ matrix box,
+ real lambda,
+ real *dvdlambda,
+ rvec *v,
+ tensor *vir,
+ ConstraintVariable econq)
{
- bool bOK, bDump;
- int start, homenr;
- tensor vir_r_m_dr;
- real scaled_delta_t;
- real invdt, vir_fac = 0, t;
- const t_ilist *settle;
- int nsettle;
- t_pbc pbc, *pbc_null;
- char buf[22];
- int nth, th;
-
- if (econq == econqForceDispl && !EI_ENERGY_MINIMIZATION(ir->eI))
+ return impl_->apply(bLog,
+ bEner,
+ step,
+ delta_step,
+ step_scaling,
+ x,
+ xprime,
+ min_proj,
+ box,
+ lambda,
+ dvdlambda,
+ v,
+ vir,
+ econq);
+}
+
+bool
+Constraints::Impl::apply(bool bLog,
+ bool bEner,
+ gmx_int64_t step,
+ int delta_step,
+ real step_scaling,
+ rvec *x,
+ rvec *xprime,
+ rvec *min_proj,
+ matrix box,
+ real lambda,
+ real *dvdlambda,
+ rvec *v,
+ tensor *vir,
+ ConstraintVariable econq)
+{
+ bool bOK, bDump;
+ int start, homenr;
+ tensor vir_r_m_dr;
+ real scaled_delta_t;
+ real invdt, vir_fac = 0, t;
+ int nsettle;
+ t_pbc pbc, *pbc_null;
+ char buf[22];
+ int nth, th;
+
+ wallcycle_start(wcycle, ewcCONSTR);
+
+ if (econq == ConstraintVariable::ForceDispl && !EI_ENERGY_MINIMIZATION(ir.eI))
{
gmx_incons("constrain called for forces displacements while not doing energy minimization, can not do this while the LINCS and SETTLE constraint connection matrices are mass weighted");
}
bDump = FALSE;
start = 0;
- homenr = md->homenr;
+ homenr = md.homenr;
- scaled_delta_t = step_scaling * ir->delta_t;
+ scaled_delta_t = step_scaling * ir.delta_t;
/* Prepare time step for use in constraint implementations, and
- avoid generating inf when ir->delta_t = 0. */
- if (ir->delta_t == 0)
+ avoid generating inf when ir.delta_t = 0. */
+ if (ir.delta_t == 0)
{
invdt = 0.0;
}
invdt = 1.0/scaled_delta_t;
}
- if (ir->efep != efepNO && EI_DYNAMICS(ir->eI))
+ if (ir.efep != efepNO && EI_DYNAMICS(ir.eI))
{
/* Set the constraint lengths for the step at which this configuration
* is meant to be. The invmasses should not be changed.
*/
- lambda += delta_step*ir->fepvals->delta_lambda;
+ lambda += delta_step*ir.fepvals->delta_lambda;
}
if (vir != nullptr)
{
clear_mat(vir_r_m_dr);
}
-
- settle = &idef->il[F_SETTLE];
+ const t_ilist *settle = &idef->il[F_SETTLE];
nsettle = settle->nr/(1+NRAL(F_SETTLE));
if (nsettle > 0)
* Note that PBC for constraints is different from PBC for bondeds.
* For constraints there is both forward and backward communication.
*/
- if (ir->ePBC != epbcNONE &&
- (cr->dd || bMolPBC) && !(cr->dd && cr->dd->constraint_comm == nullptr))
+ if (ir.ePBC != epbcNONE &&
+ (cr->dd || pbcHandlingRequired_) && !(cr->dd && cr->dd->constraint_comm == nullptr))
{
/* With pbc=screw the screw has been changed to a shift
* by the constraint coordinate communication routine,
* so that here we can use normal pbc.
*/
- pbc_null = set_pbc_dd(&pbc, ir->ePBC,
+ pbc_null = set_pbc_dd(&pbc, ir.ePBC,
DOMAINDECOMP(cr) ? cr->dd->nc : nullptr,
FALSE, box);
}
*/
if (cr->dd)
{
- dd_move_x_constraints(cr->dd, box, x, xprime, econq == econqCoord);
+ dd_move_x_constraints(cr->dd, box, x, xprime, econq == ConstraintVariable::Positions);
if (v != nullptr)
{
}
}
- if (constr->lincsd != nullptr)
+ if (lincsd != nullptr)
{
- bOK = constrain_lincs(fplog, bLog, bEner, ir, step, constr->lincsd, md, cr, ms,
+ bOK = constrain_lincs(log, bLog, bEner, ir, step, lincsd, md, cr, ms,
x, xprime, min_proj,
box, pbc_null, lambda, dvdlambda,
invdt, v, vir != nullptr, vir_r_m_dr,
econq, nrnb,
- constr->maxwarn, &constr->warncount_lincs);
- if (!bOK && constr->maxwarn < INT_MAX)
+ maxwarn, &warncount_lincs);
+ if (!bOK && maxwarn < INT_MAX)
{
- if (fplog != nullptr)
+ if (log != nullptr)
{
- fprintf(fplog, "Constraint error in algorithm %s at step %s\n",
+ fprintf(log, "Constraint error in algorithm %s at step %s\n",
econstr_names[econtLINCS], gmx_step_str(step, buf));
}
bDump = TRUE;
}
}
- if (constr->shaked != nullptr)
+ if (shaked != nullptr)
{
- bOK = constrain_shake(fplog, constr->shaked,
- md->invmass,
- idef, ir, x, xprime, min_proj, nrnb,
+ bOK = constrain_shake(log, shaked,
+ md.invmass,
+ *idef, ir, x, xprime, min_proj, nrnb,
lambda, dvdlambda,
invdt, v, vir != nullptr, vir_r_m_dr,
- constr->maxwarn < INT_MAX, econq);
+ maxwarn < INT_MAX, econq);
- if (!bOK && constr->maxwarn < INT_MAX)
+ if (!bOK && maxwarn < INT_MAX)
{
- if (fplog != nullptr)
+ if (log != nullptr)
{
- fprintf(fplog, "Constraint error in algorithm %s at step %s\n",
+ fprintf(log, "Constraint error in algorithm %s at step %s\n",
econstr_names[econtSHAKE], gmx_step_str(step, buf));
}
bDump = TRUE;
if (nsettle > 0)
{
- bool bSettleErrorHasOccurred = false;
+ bool bSettleErrorHasOccurred0 = false;
switch (econq)
{
- case econqCoord:
+ case ConstraintVariable::Positions:
#pragma omp parallel for num_threads(nth) schedule(static)
for (th = 0; th < nth; th++)
{
{
if (th > 0)
{
- clear_mat(constr->vir_r_m_dr_th[th]);
+ clear_mat(vir_r_m_dr_th[th]);
}
- csettle(constr->settled,
+ csettle(settled,
nth, th,
pbc_null,
x[0], xprime[0],
invdt, v ? v[0] : nullptr,
vir != nullptr,
- th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
- th == 0 ? &bSettleErrorHasOccurred : &constr->bSettleErrorHasOccurred[th]);
+ th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th],
+ th == 0 ? &bSettleErrorHasOccurred0 : &bSettleErrorHasOccurred[th]);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
}
inc_nrnb(nrnb, eNR_CONSTR_VIR, nsettle*3);
}
break;
- case econqVeloc:
- case econqDeriv:
- case econqForce:
- case econqForceDispl:
+ case ConstraintVariable::Velocities:
+ case ConstraintVariable::Derivative:
+ case ConstraintVariable::Force:
+ case ConstraintVariable::ForceDispl:
#pragma omp parallel for num_threads(nth) schedule(static)
for (th = 0; th < nth; th++)
{
}
else
{
- calcvir_atom_end = md->homenr;
+ calcvir_atom_end = md.homenr;
}
if (th > 0)
{
- clear_mat(constr->vir_r_m_dr_th[th]);
+ clear_mat(vir_r_m_dr_th[th]);
}
int start_th = (nsettle* th )/nth;
if (start_th >= 0 && end_th - start_th > 0)
{
- settle_proj(constr->settled, econq,
+ settle_proj(settled, econq,
end_th-start_th,
settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
pbc_null,
x,
xprime, min_proj, calcvir_atom_end,
- th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th]);
+ th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]);
}
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
/* This is an overestimate */
inc_nrnb(nrnb, eNR_SETTLE, nsettle);
break;
- case econqDeriv_FlexCon:
+ case ConstraintVariable::Deriv_FlexCon:
/* Nothing to do, since the are no flexible constraints in settles */
break;
default:
/* Reduce the virial contributions over the threads */
for (int th = 1; th < nth; th++)
{
- m_add(vir_r_m_dr, constr->vir_r_m_dr_th[th], vir_r_m_dr);
+ m_add(vir_r_m_dr, vir_r_m_dr_th[th], vir_r_m_dr);
}
}
- if (econq == econqCoord)
+ if (econq == ConstraintVariable::Positions)
{
for (int th = 1; th < nth; th++)
{
- bSettleErrorHasOccurred = bSettleErrorHasOccurred || constr->bSettleErrorHasOccurred[th];
+ bSettleErrorHasOccurred0 = bSettleErrorHasOccurred0 || bSettleErrorHasOccurred[th];
}
- if (bSettleErrorHasOccurred)
+ if (bSettleErrorHasOccurred0)
{
char buf[STRLEN];
sprintf(buf,
"\nstep " "%" GMX_PRId64 ": One or more water molecules can not be settled.\n"
"Check for bad contacts and/or reduce the timestep if appropriate.\n",
step);
- if (fplog)
+ if (log)
{
- fprintf(fplog, "%s", buf);
+ fprintf(log, "%s", buf);
}
fprintf(stderr, "%s", buf);
- constr->warncount_settle++;
- if (constr->warncount_settle > constr->maxwarn)
+ warncount_settle++;
+ if (warncount_settle > maxwarn)
{
- too_many_constraint_warnings(-1, constr->warncount_settle);
+ too_many_constraint_warnings(-1, warncount_settle);
}
bDump = TRUE;
* 0.5 also passes vir = NULL, so cannot reach this
* assertion. This assertion should remain until someone knows
* that this path works for their intended purpose, and then
- * they can use scaled_delta_t instead of ir->delta_t
+ * they can use scaled_delta_t instead of ir.delta_t
* below. */
assert(gmx_within_tol(step_scaling, 1.0, GMX_REAL_EPS));
switch (econq)
{
- case econqCoord:
- vir_fac = 0.5/(ir->delta_t*ir->delta_t);
+ case ConstraintVariable::Positions:
+ vir_fac = 0.5/(ir.delta_t*ir.delta_t);
break;
- case econqVeloc:
- vir_fac = 0.5/ir->delta_t;
+ case ConstraintVariable::Velocities:
+ vir_fac = 0.5/ir.delta_t;
break;
- case econqForce:
- case econqForceDispl:
+ case ConstraintVariable::Force:
+ case ConstraintVariable::ForceDispl:
vir_fac = 0.5;
break;
default:
gmx_incons("Unsupported constraint quantity for virial");
}
- if (EI_VV(ir->eI))
+ if (EI_VV(ir.eI))
{
vir_fac *= 2; /* only constraining over half the distance here */
}
if (bDump)
{
- dump_confs(fplog, step, constr->warn_mtop, start, homenr, cr, x, xprime, box);
+ dump_confs(log, step, mtop, start, homenr, cr, x, xprime, box);
}
- if (econq == econqCoord)
+ if (econq == ConstraintVariable::Positions)
{
- if (ir->bPull && pull_have_constraint(ir->pull_work))
+ if (ir.bPull && pull_have_constraint(ir.pull_work))
{
- if (EI_DYNAMICS(ir->eI))
+ if (EI_DYNAMICS(ir.eI))
{
- t = ir->init_t + (step + delta_step)*ir->delta_t;
+ t = ir.init_t + (step + delta_step)*ir.delta_t;
}
else
{
- t = ir->init_t;
+ t = ir.init_t;
}
- set_pbc(&pbc, ir->ePBC, box);
- pull_constraint(ir->pull_work, md, &pbc, cr, ir->delta_t, t, x, xprime, v, *vir);
+ set_pbc(&pbc, ir.ePBC, box);
+ pull_constraint(ir.pull_work, &md, &pbc, cr, ir.delta_t, t, x, xprime, v, *vir);
}
- if (constr->ed && delta_step > 0)
+ if (ed && delta_step > 0)
{
/* apply the essential dynamics constraints here */
- do_edsam(ir, step, cr, xprime, v, box, constr->ed);
+ do_edsam(&ir, step, cr, xprime, v, box, ed);
}
}
+ wallcycle_stop(wcycle, ewcCONSTR);
return bOK;
}
-real *constr_rmsd_data(Constraints *constr)
+real *Constraints::rmsdData() const
{
- if (constr->lincsd)
+ if (impl_->lincsd)
{
- return lincs_rmsd_data(constr->lincsd);
+ return lincs_rmsd_data(impl_->lincsd);
}
else
{
}
}
-real constr_rmsd(const Constraints *constr)
+real Constraints::rmsd() const
{
- if (constr->lincsd)
+ if (impl_->lincsd)
{
- return lincs_rmsd(constr->lincsd);
+ return lincs_rmsd(impl_->lincsd);
}
else
{
return at2s;
}
-void set_constraints(Constraints *constr,
- gmx_localtop_t *top, const t_inputrec *ir,
- const t_mdatoms *md, const t_commrec *cr)
+void
+Constraints::Impl::setConstraints(const gmx_localtop_t &top,
+ const t_mdatoms &md)
{
- t_idef *idef = &top->idef;
+ idef = &top.idef;
- if (constr->ncon_tot > 0)
+ if (ncon_tot > 0)
{
/* With DD we might also need to call LINCS on a domain no constraints for
* communicating coordinates to other nodes that do have constraints.
*/
- if (ir->eConstrAlg == econtLINCS)
+ if (ir.eConstrAlg == econtLINCS)
{
- set_lincs(idef, md, EI_DYNAMICS(ir->eI), cr, constr->lincsd);
+ set_lincs(top.idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
}
- if (ir->eConstrAlg == econtSHAKE)
+ if (ir.eConstrAlg == econtSHAKE)
{
if (cr->dd)
{
// We are using the local topology, so there are only
// F_CONSTR constraints.
- make_shake_sblock_dd(constr->shaked, &idef->il[F_CONSTR], &top->cgs, cr->dd);
+ make_shake_sblock_dd(shaked, &idef->il[F_CONSTR], &top.cgs, cr->dd);
}
else
{
- make_shake_sblock_serial(constr->shaked, idef, md);
+ make_shake_sblock_serial(shaked, idef, md);
}
}
}
- if (constr->settled)
+ if (settled)
{
- settle_set_constraints(constr->settled,
+ settle_set_constraints(settled,
&idef->il[F_SETTLE], md);
}
/* Make a selection of the local atoms for essential dynamics */
- if (constr->ed && cr->dd)
+ if (ed && cr->dd)
{
- dd_make_local_ed_indices(cr->dd, constr->ed);
+ dd_make_local_ed_indices(cr->dd, ed);
}
}
-Constraints *init_constraints(FILE *fplog,
- const gmx_mtop_t *mtop, const t_inputrec *ir,
- bool doEssentialDynamics,
- const t_commrec *cr)
+void
+Constraints::setConstraints(const gmx_localtop_t &top,
+ const t_mdatoms &md)
{
- int nconstraints =
- gmx_mtop_ftype_count(mtop, F_CONSTR) +
- gmx_mtop_ftype_count(mtop, F_CONSTRNC);
- int nsettles =
- gmx_mtop_ftype_count(mtop, F_SETTLE);
+ impl_->setConstraints(top, md);
+}
- GMX_RELEASE_ASSERT(!ir->bPull || ir->pull_work != nullptr, "init_constraints called with COM pulling before/without initializing the pull code");
+Constraints::Constraints(const gmx_mtop_t &mtop,
+ const t_inputrec &ir,
+ FILE *log,
+ const t_mdatoms &md,
+ const t_commrec *cr,
+ const gmx_multisim_t &ms,
+ t_nrnb *nrnb,
+ gmx_wallcycle *wcycle,
+ bool pbcHandlingRequired,
+ int numConstraints,
+ int numSettles)
+ : impl_(new Impl(mtop,
+ ir,
+ log,
+ md,
+ cr,
+ ms,
+ nrnb,
+ wcycle,
+ pbcHandlingRequired,
+ numConstraints,
+ numSettles))
+{
+}
- if (nconstraints + nsettles == 0 &&
- !(ir->bPull && pull_have_constraint(ir->pull_work)) &&
- !doEssentialDynamics)
+Constraints::Impl::Impl(const gmx_mtop_t &mtop_p,
+ const t_inputrec &ir_p,
+ FILE *log_p,
+ const t_mdatoms &md_p,
+ const t_commrec *cr_p,
+ const gmx_multisim_t &ms_p,
+ t_nrnb *nrnb_p,
+ gmx_wallcycle *wcycle_p,
+ bool pbcHandlingRequired,
+ int numConstraints,
+ int numSettles)
+ : ncon_tot(numConstraints),
+ mtop(mtop_p),
+ md(md_p),
+ pbcHandlingRequired_(pbcHandlingRequired),
+ log(log_p),
+ cr(cr_p),
+ ms(ms_p),
+ ir(ir_p),
+ nrnb(nrnb_p),
+ wcycle(wcycle_p)
+{
+ if (numConstraints + numSettles > 0 && ir.epc == epcMTTK)
{
- return nullptr;
+ gmx_fatal(FARGS, "Constraints are not implemented with MTTK pressure control.");
}
- Constraints *constr;
-
- snew(constr, 1);
-
- constr->ncon_tot = nconstraints;
- constr->nflexcon = 0;
- if (nconstraints > 0)
+ nflexcon = 0;
+ if (numConstraints > 0)
{
- constr->n_at2con_mt = mtop->moltype.size();
- snew(constr->at2con_mt, constr->n_at2con_mt);
- for (int mt = 0; mt < static_cast<int>(mtop->moltype.size()); mt++)
+ n_at2con_mt = mtop.moltype.size();
+ snew(at2con_mt, n_at2con_mt);
+ for (int mt = 0; mt < static_cast<int>(mtop.moltype.size()); mt++)
{
int nflexcon;
- constr->at2con_mt[mt] = make_at2con(0, mtop->moltype[mt].atoms.nr,
- mtop->moltype[mt].ilist,
- mtop->ffparams.iparams,
- EI_DYNAMICS(ir->eI), &nflexcon);
- for (const gmx_molblock_t &molblock : mtop->molblock)
+ at2con_mt[mt] = make_at2con(0, mtop.moltype[mt].atoms.nr,
+ mtop.moltype[mt].ilist,
+ mtop.ffparams.iparams,
+ EI_DYNAMICS(ir.eI), &nflexcon);
+ for (const gmx_molblock_t &molblock : mtop.molblock)
{
if (molblock.type == mt)
{
- constr->nflexcon += molblock.nmol*nflexcon;
+ nflexcon += molblock.nmol*nflexcon;
}
}
}
- if (constr->nflexcon > 0)
+ if (nflexcon > 0)
{
- if (fplog)
+ if (log)
{
- fprintf(fplog, "There are %d flexible constraints\n",
- constr->nflexcon);
- if (ir->fc_stepsize == 0)
+ fprintf(log, "There are %d flexible constraints\n",
+ nflexcon);
+ if (ir.fc_stepsize == 0)
{
- fprintf(fplog, "\n"
+ fprintf(log, "\n"
"WARNING: step size for flexible constraining = 0\n"
" All flexible constraints will be rigid.\n"
" Will try to keep all flexible constraints at their original length,\n"
" but the lengths may exhibit some drift.\n\n");
- constr->nflexcon = 0;
+ nflexcon = 0;
}
}
- if (constr->nflexcon > 0)
+ if (nflexcon > 0)
{
- please_cite(fplog, "Hess2002");
+ please_cite(log, "Hess2002");
}
}
- if (ir->eConstrAlg == econtLINCS)
+ if (ir.eConstrAlg == econtLINCS)
{
- constr->lincsd = init_lincs(fplog, mtop,
- constr->nflexcon, constr->at2con_mt,
- DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
- ir->nLincsIter, ir->nProjOrder);
+ lincsd = init_lincs(log, mtop,
+ nflexcon, at2con_mt,
+ DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
+ ir.nLincsIter, ir.nProjOrder);
}
- if (ir->eConstrAlg == econtSHAKE)
+ if (ir.eConstrAlg == econtSHAKE)
{
if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
{
gmx_fatal(FARGS, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
}
- if (constr->nflexcon)
+ if (nflexcon)
{
gmx_fatal(FARGS, "For this system also velocities and/or forces need to be constrained, this can not be done with SHAKE, you should select LINCS");
}
- please_cite(fplog, "Ryckaert77a");
- if (ir->bShakeSOR)
+ please_cite(log, "Ryckaert77a");
+ if (ir.bShakeSOR)
{
- please_cite(fplog, "Barth95a");
+ please_cite(log, "Barth95a");
}
- constr->shaked = shake_init();
+ shaked = shake_init();
}
}
- if (nsettles > 0)
+ if (numSettles > 0)
{
- please_cite(fplog, "Miyamoto92a");
+ please_cite(log, "Miyamoto92a");
- constr->bInterCGsettles = inter_charge_group_settles(mtop);
+ bInterCGsettles = inter_charge_group_settles(mtop);
- constr->settled = settle_init(mtop);
+ settled = settle_init(mtop);
/* Make an atom to settle index for use in domain decomposition */
- constr->n_at2settle_mt = mtop->moltype.size();
- snew(constr->at2settle_mt, constr->n_at2settle_mt);
- for (size_t mt = 0; mt < mtop->moltype.size(); mt++)
+ n_at2settle_mt = mtop.moltype.size();
+ snew(at2settle_mt, n_at2settle_mt);
+ for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
{
- constr->at2settle_mt[mt] =
- make_at2settle(mtop->moltype[mt].atoms.nr,
- &mtop->moltype[mt].ilist[F_SETTLE]);
+ at2settle_mt[mt] =
+ make_at2settle(mtop.moltype[mt].atoms.nr,
+ &mtop.moltype[mt].ilist[F_SETTLE]);
}
/* Allocate thread-local work arrays */
int nthreads = gmx_omp_nthreads_get(emntSETTLE);
- if (nthreads > 1 && constr->vir_r_m_dr_th == nullptr)
+ if (nthreads > 1 && vir_r_m_dr_th == nullptr)
{
- snew(constr->vir_r_m_dr_th, nthreads);
- snew(constr->bSettleErrorHasOccurred, nthreads);
+ snew(vir_r_m_dr_th, nthreads);
+ snew(bSettleErrorHasOccurred, nthreads);
}
}
- if (nconstraints + nsettles > 0 && ir->epc == epcMTTK)
- {
- gmx_fatal(FARGS, "Constraints are not implemented with MTTK pressure control.");
- }
-
- constr->maxwarn = 999;
+ maxwarn = 999;
char *env = getenv("GMX_MAXCONSTRWARN");
if (env)
{
- constr->maxwarn = 0;
- sscanf(env, "%8d", &constr->maxwarn);
- if (constr->maxwarn < 0)
+ maxwarn = 0;
+ sscanf(env, "%8d", &maxwarn);
+ if (maxwarn < 0)
{
- constr->maxwarn = INT_MAX;
+ maxwarn = INT_MAX;
}
- if (fplog)
+ if (log)
{
- fprintf(fplog,
+ fprintf(log,
"Setting the maximum number of constraint warnings to %d\n",
- constr->maxwarn);
+ maxwarn);
}
if (MASTER(cr))
{
fprintf(stderr,
"Setting the maximum number of constraint warnings to %d\n",
- constr->maxwarn);
+ maxwarn);
}
}
- constr->warncount_lincs = 0;
- constr->warncount_settle = 0;
-
- constr->warn_mtop = mtop;
-
- return constr;
+ warncount_lincs = 0;
+ warncount_settle = 0;
}
-void saveEdsamPointer(Constraints *constr, gmx_edsam *ed)
+void Constraints::saveEdsamPointer(gmx_edsam_t ed)
{
- constr->ed = ed;
+ impl_->ed = ed;
}
-const t_blocka *atom2constraints_moltype(const Constraints *constr)
+const t_blocka *Constraints::atom2constraints_moltype() const
{
- return constr->at2con_mt;
+ return impl_->at2con_mt;
}
-const int **atom2settle_moltype(const Constraints *constr)
+const int **Constraints::atom2settle_moltype() const
{
- return (const int **)constr->at2settle_mt;
+ return (const int **)impl_->at2settle_mt;
}
-bool inter_charge_group_constraints(const gmx_mtop_t *mtop)
+bool inter_charge_group_constraints(const gmx_mtop_t &mtop)
{
const gmx_moltype_t *molt;
const t_block *cgs;
bool bInterCG;
bInterCG = FALSE;
- for (size_t mb = 0; mb < mtop->molblock.size() && !bInterCG; mb++)
+ for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++)
{
- molt = &mtop->moltype[mtop->molblock[mb].type];
+ molt = &mtop.moltype[mtop.molblock[mb].type];
if (molt->ilist[F_CONSTR].nr > 0 ||
molt->ilist[F_CONSTRNC].nr > 0 ||
return bInterCG;
}
-bool inter_charge_group_settles(const gmx_mtop_t *mtop)
+bool inter_charge_group_settles(const gmx_mtop_t &mtop)
{
const gmx_moltype_t *molt;
const t_block *cgs;
bool bInterCG;
bInterCG = FALSE;
- for (size_t mb = 0; mb < mtop->molblock.size() && !bInterCG; mb++)
+ for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++)
{
- molt = &mtop->moltype[mtop->molblock[mb].type];
+ molt = &mtop.moltype[mtop.molblock[mb].type];
if (molt->ilist[F_SETTLE].nr > 0)
{
#include "gromacs/math/vectypes.h"
#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/classhelpers.h"
#include "gromacs/utility/real.h"
struct gmx_edsam;
struct gmx_localtop_t;
struct gmx_mtop_t;
struct gmx_multisim_t;
+struct gmx_wallcycle;
struct t_blocka;
struct t_commrec;
-struct t_idef;
struct t_ilist;
struct t_inputrec;
struct t_mdatoms;
namespace gmx
{
-class Constraints;
-
-enum
+//! Describes supported flavours of constrained updates.
+enum class ConstraintVariable : int
{
- econqCoord, /* Constrain coordinates (mass weighted) */
- econqVeloc, /* Constrain velocities (mass weighted) */
- econqDeriv, /* Constrain a derivative (mass weighted), *
- * for instance velocity or acceleration, *
- * constraint virial can not be calculated. */
- econqDeriv_FlexCon, /* As econqDeriv, but only output flex. con. */
- econqForce, /* Constrain forces (non mass-weighted) */
- econqForceDispl /* Constrain forces (mass-weighted 1/0 for freeze) */
+ Positions, /* Constrain positions (mass weighted) */
+ Velocities, /* Constrain velocities (mass weighted) */
+ Derivative, /* Constrain a derivative (mass weighted), *
+ * for instance velocity or acceleration, *
+ * constraint virial can not be calculated. */
+ Deriv_FlexCon, /* As Derivative, but only output flex. con. */
+ Force, /* Constrain forces (non mass-weighted) */
+ // TODO What does this do? Improve the comment.
+ ForceDispl /* Like Force, but free particles will have mass
+ * 1 and frozen particles mass 0 */
};
-/*! \brief Returns the total number of flexible constraints in the system. */
-int n_flexible_constraints(const Constraints *constr);
+/*! \libinternal
+ * \brief Handles constraints */
+class Constraints
+{
+ private:
+ /*! \brief Constructor
+ *
+ * Private to enforce use of makeConstraints() factory
+ * function. */
+ Constraints(const gmx_mtop_t &mtop,
+ const t_inputrec &ir,
+ FILE *log,
+ const t_mdatoms &md,
+ const t_commrec *cr,
+ const gmx_multisim_t &ms,
+ t_nrnb *nrnb,
+ gmx_wallcycle *wcycle,
+ bool pbcHandlingRequired,
+ int numConstraints,
+ int numSettles);
+ public:
+ /*! \brief This member type helps implement a factory
+ * function, because its objects can access the private
+ * constructor. */
+ struct CreationHelper;
+
+ //! Destructor.
+ ~Constraints();
+
+ /*! \brief Returns the total number of flexible constraints in the system. */
+ int numFlexibleConstraints() const;
+
+ /*! \brief Set up all the local constraints for the domain.
+ *
+ * \todo Make this a callback that is called automatically
+ * once a new domain has been made. */
+ void setConstraints(const gmx_localtop_t &top,
+ const t_mdatoms &md);
+
+ /*! \brief Applies constraints to coordinates.
+ *
+ * When econq=ConstraintVariable::Positions constrains
+ * coordinates xprime using th directions in x, min_proj is
+ * not used.
+ *
+ * When econq=ConstraintVariable::Derivative, calculates the
+ * components xprime in the constraint directions and
+ * subtracts these components from min_proj. So when
+ * min_proj=xprime, the constraint components are projected
+ * out.
+ *
+ * When econq=ConstraintVariable::Deriv_FlexCon, the same is
+ * done as with ConstraintVariable::Derivative, but only the
+ * components of the flexible constraints are stored.
+ *
+ * delta_step is used for determining the constraint reference lengths
+ * when lenA != lenB or will the pull code with a pulling rate.
+ * step + delta_step is the step at which the final configuration
+ * is meant to be; for update delta_step = 1.
+ *
+ * step_scaling can be used to update coordinates based on the time
+ * step multiplied by this factor. Thus, normally 1.0 is passed. The
+ * SD1 integrator uses 0.5 in one of its calls, to correct positions
+ * for half a step of changed velocities.
+ *
+ * If v!=NULL also constrain v by adding the constraint corrections / dt.
+ *
+ * If vir!=NULL calculate the constraint virial.
+ *
+ * Return whether the application of constraints succeeded without error.
+ */
+ bool apply(bool bLog,
+ bool bEner,
+ gmx_int64_t step,
+ int delta_step,
+ real step_scaling,
+ rvec *x,
+ rvec *xprime,
+ rvec *min_proj,
+ matrix box,
+ real lambda,
+ real *dvdlambda,
+ rvec *v,
+ tensor *vir,
+ ConstraintVariable econq);
+ //! Links the essentialdynamics and constraint code.
+ void saveEdsamPointer(gmx_edsam *ed);
+ //! Getter for use by domain decomposition.
+ const t_blocka *atom2constraints_moltype() const;
+ //! Getter for use by domain decomposition.
+ const int **atom2settle_moltype() const;
+
+ /*! \brief Return the data for reduction for determining
+ * constraint RMS relative deviations, or nullptr when not
+ * supported for any active constraints. */
+ real *rmsdData() const;
+ /*! \brief Return the RMSD of the constraints when available. */
+ real rmsd() const;
+
+ private:
+ //! Implementation type.
+ class Impl;
+ //! Implementation object.
+ PrivateImplPointer<Impl> impl_;
+};
/*! \brief Generate a fatal error because of too many LINCS/SETTLE warnings. */
void too_many_constraint_warnings(int eConstrAlg, int warncount);
-/*! \brief Applies constraints to coordinates.
- *
- * When econq=econqCoord constrains coordinates xprime using th
- * directions in x, min_proj is not used.
- *
- * When econq=econqDeriv, calculates the components xprime in
- * the constraint directions and subtracts these components from min_proj.
- * So when min_proj=xprime, the constraint components are projected out.
- *
- * When econq=econqDeriv_FlexCon, the same is done as with econqDeriv,
- * but only the components of the flexible constraints are stored.
- *
- * When bMolPBC=TRUE, assume that molecules might be broken: correct PBC.
- *
- * delta_step is used for determining the constraint reference lengths
- * when lenA != lenB or will the pull code with a pulling rate.
- * step + delta_step is the step at which the final configuration
- * is meant to be; for update delta_step = 1.
- *
- * step_scaling can be used to update coordinates based on the time
- * step multiplied by this factor. Thus, normally 1.0 is passed. The
- * SD1 integrator uses 0.5 in one of its calls, to correct positions
- * for half a step of changed velocities.
- *
- * If v!=NULL also constrain v by adding the constraint corrections / dt.
- *
- * If vir!=NULL calculate the constraint virial.
- *
- * Return TRUE if OK, FALSE in case of shake error
- *
- */
-bool constrain(FILE *log, bool bLog, bool bEner,
- Constraints *constr,
- const t_idef *idef,
- const t_inputrec *ir,
- const t_commrec *cr,
- const gmx_multisim_t *ms,
- gmx_int64_t step, int delta_step,
- real step_scaling,
- const t_mdatoms *md,
- rvec *x, rvec *xprime, rvec *min_proj,
- bool bMolPBC, matrix box,
- real lambda, real *dvdlambda,
- rvec *v, tensor *vir,
- t_nrnb *nrnb, int econq);
-
-/*! \brief Initialize constraints stuff */
-Constraints *init_constraints(FILE *log,
- const gmx_mtop_t *mtop, const t_inputrec *ir,
- bool doEssentialDynamics,
- const t_commrec *cr);
-
-/*! \brief Put a pointer to the essential dynamics constraints into the constr struct. */
-void saveEdsamPointer(Constraints *constr,
- gmx_edsam *ed);
-
-/*! \brief Set up all the local constraints for this rank. */
-void set_constraints(Constraints *constr,
- gmx_localtop_t *top,
- const t_inputrec *ir,
- const t_mdatoms *md,
- const t_commrec *cr);
-
/* The at2con t_blocka struct returned by the routines below
* contains a list of constraints per atom.
* The F_CONSTRNC constraints in this structure number consecutively
const t_ilist *ilist, const t_iparams *iparams,
bool bDynamics, int *nflexiblecons);
-/*! \brief Returns an array of atom to constraints lists for the moltypes */
-const t_blocka *atom2constraints_moltype(const Constraints *constr);
-
-/*! \brief Returns an array of atom to settles lists for the moltypes */
-const int **atom2settle_moltype(const Constraints *constr);
-
/*! \brief Macro for getting the constraint iatoms for a constraint number con
* which comes from a list where F_CONSTR and F_CONSTRNC constraints
* are concatenated. */
#define constr_iatomptr(nconstr, iatom_constr, iatom_constrnc, con) ((con) < (nconstr) ? (iatom_constr)+(con)*3 : (iatom_constrnc)+(con-nconstr)*3)
/*! \brief Returns whether there are inter charge group constraints */
-bool inter_charge_group_constraints(const gmx_mtop_t *mtop);
+bool inter_charge_group_constraints(const gmx_mtop_t &mtop);
/*! \brief Returns whether there are inter charge group settles */
-bool inter_charge_group_settles(const gmx_mtop_t *mtop);
-
-/*! \brief Return the data for determining constraint RMS relative deviations.
- * Returns NULL when LINCS is not used. */
-real *constr_rmsd_data(Constraints *constr);
-
-/*! \brief Return the RMSD of the constraint */
-real constr_rmsd(const Constraints *constr);
+bool inter_charge_group_settles(const gmx_mtop_t &mtop);
} // namespace
#endif // GMX_SIMD_HAVE_REAL
/*! \brief LINCS projection, works on derivatives of the coordinates. */
-static void do_lincsp(rvec *x, rvec *f, rvec *fp, t_pbc *pbc,
+static void do_lincsp(const rvec *x, rvec *f, rvec *fp, t_pbc *pbc,
Lincs *lincsd, int th,
real *invmass,
- int econq, bool bCalcDHDL,
+ ConstraintVariable econq, bool bCalcDHDL,
bool bCalcVir, tensor rmdf)
{
int b0, b1, b;
r = lincsd->tmpv;
blnr = lincsd->blnr;
blbnb = lincsd->blbnb;
- if (econq != econqForce)
+ if (econq != ConstraintVariable::Force)
{
/* Use mass-weighted parameters */
blc = lincsd->blc;
lincs_matrix_expand(lincsd, &lincsd->task[th], blcc, rhs1, rhs2, sol);
/* nrec*(ncons+2*nrtot) flops */
- if (econq == econqDeriv_FlexCon)
+ if (econq == ConstraintVariable::Deriv_FlexCon)
{
/* We only want to constraint the flexible constraints,
* so we mask out the normal ones by setting sol to 0.
* so we pass invmass=NULL, which results in the use of 1 for all atoms.
*/
lincs_update_atoms(lincsd, th, 1.0, sol, r,
- (econq != econqForce) ? invmass : nullptr, fp);
+ (econq != ConstraintVariable::Force) ? invmass : nullptr, fp);
if (bCalcDHDL)
{
#endif // GMX_SIMD_HAVE_REAL
//! Implements LINCS constraining.
-static void do_lincs(rvec *x, rvec *xp, matrix box, t_pbc *pbc,
+static void do_lincs(const rvec *x, rvec *xp, matrix box, t_pbc *pbc,
Lincs *lincsd, int th,
const real *invmass,
const t_commrec *cr,
return (*(int *)a) - (*(int *)b);
}
-Lincs *init_lincs(FILE *fplog, const gmx_mtop_t *mtop,
+Lincs *init_lincs(FILE *fplog, const gmx_mtop_t &mtop,
int nflexcon_global, const t_blocka *at2con,
bool bPLINCS, int nIter, int nProjOrder)
{
li->nOrder = nProjOrder;
li->max_connect = 0;
- for (size_t mt = 0; mt < mtop->moltype.size(); mt++)
+ for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
{
- for (int a = 0; a < mtop->moltype[mt].atoms.nr; a++)
+ for (int a = 0; a < mtop.moltype[mt].atoms.nr; a++)
{
li->max_connect = std::max(li->max_connect,
at2con[mt].index[a + 1] - at2con[mt].index[a]);
li->ncg_triangle = 0;
bMoreThanTwoSeq = FALSE;
- for (const gmx_molblock_t &molb : mtop->molblock)
+ for (const gmx_molblock_t &molb : mtop.molblock)
{
- const gmx_moltype_t &molt = mtop->moltype[molb.type];
+ const gmx_moltype_t &molt = mtop.moltype[molb.type];
li->ncg_triangle +=
molb.nmol*
* to other constraints, and if so add those connected constraints to our task. */
static void check_assign_connected(Lincs *li,
const t_iatom *iatom,
- const t_idef *idef,
+ const t_idef &idef,
int bDynamics,
int a1, int a2,
const t_blocka *at2con)
real lenA, lenB;
type = iatom[cc*3];
- lenA = idef->iparams[type].constr.dA;
- lenB = idef->iparams[type].constr.dB;
+ lenA = idef.iparams[type].constr.dA;
+ lenB = idef.iparams[type].constr.dB;
if (bDynamics || lenA != 0 || lenB != 0)
{
* in the triangle to our task. */
static void check_assign_triangle(Lincs *li,
const t_iatom *iatom,
- const t_idef *idef,
+ const t_idef &idef,
int bDynamics,
int constraint_index,
int a1, int a2,
i = c_triangle[end]*3;
type = iatom[i];
- lenA = idef->iparams[type].constr.dA;
- lenB = idef->iparams[type].constr.dB;
+ lenA = idef.iparams[type].constr.dA;
+ lenB = idef.iparams[type].constr.dB;
if (bDynamics || lenA != 0 || lenB != 0)
{
}
}
-void set_lincs(const t_idef *idef,
- const t_mdatoms *md,
+void set_lincs(const t_idef &idef,
+ const t_mdatoms &md,
bool bDynamics,
const t_commrec *cr,
Lincs *li)
}
/* This is the local topology, so there are only F_CONSTR constraints */
- if (idef->il[F_CONSTR].nr == 0)
+ if (idef.il[F_CONSTR].nr == 0)
{
/* There are no constraints,
* we do not need to fill any data structures.
}
else
{
- natoms = md->homenr;
+ natoms = md.homenr;
}
- at2con = make_at2con(0, natoms, idef->il, idef->iparams, bDynamics,
+ at2con = make_at2con(0, natoms, idef.il, idef.iparams, bDynamics,
&nflexcon);
- ncon_tot = idef->il[F_CONSTR].nr/3;
+ ncon_tot = idef.il[F_CONSTR].nr/3;
/* Ensure we have enough padding for aligned loads for each thread */
if (ncon_tot + li->ntask*simd_width > li->nc_alloc || li->nc_alloc == 0)
resize_real_aligned(&li->mlambda, li->nc_alloc);
}
- iatom = idef->il[F_CONSTR].iatoms;
+ iatom = idef.il[F_CONSTR].iatoms;
ncc_alloc_old = li->ncc_alloc;
li->blnr[0] = li->ncc;
type = iatom[3*con];
a1 = iatom[3*con + 1];
a2 = iatom[3*con + 2];
- lenA = idef->iparams[type].constr.dA;
- lenB = idef->iparams[type].constr.dB;
+ lenA = idef.iparams[type].constr.dA;
+ lenB = idef.iparams[type].constr.dB;
/* Skip the flexible constraints when not doing dynamics */
if (bDynamics || lenA != 0 || lenB != 0)
{
if (li->ntask > 1)
{
- lincs_thread_setup(li, md->nr);
+ lincs_thread_setup(li, md.nr);
}
- set_lincs_matrix(li, md->invmass, md->lambda);
+ set_lincs_matrix(li, md.invmass, md.lambda);
}
//! Issues a warning when LINCS constraints cannot be satisfied.
static void lincs_warning(FILE *fplog,
- gmx_domdec_t *dd, rvec *x, rvec *xprime, t_pbc *pbc,
+ gmx_domdec_t *dd, const rvec *x, rvec *xprime, t_pbc *pbc,
int ncons, int *bla, real *bllen, real wangle,
int maxwarn, int *warncount)
{
}
bool constrain_lincs(FILE *fplog, bool bLog, bool bEner,
- const t_inputrec *ir,
+ const t_inputrec &ir,
gmx_int64_t step,
- Lincs *lincsd, const t_mdatoms *md,
+ Lincs *lincsd, const t_mdatoms &md,
const t_commrec *cr,
- const gmx_multisim_t *ms,
- rvec *x, rvec *xprime, rvec *min_proj,
+ const gmx_multisim_t &ms,
+ const rvec *x, rvec *xprime, rvec *min_proj,
matrix box, t_pbc *pbc,
real lambda, real *dvdlambda,
real invdt, rvec *v,
bool bCalcVir, tensor vir_r_m_dr,
- int econq,
+ ConstraintVariable econq,
t_nrnb *nrnb,
int maxwarn, int *warncount)
{
* We can also easily check if any constraint length is changed,
* if not dH/dlambda=0 and we can also set the boolean to FALSE.
*/
- bCalcDHDL = (ir->efep != efepNO && dvdlambda != nullptr);
+ bCalcDHDL = (ir.efep != efepNO && dvdlambda != nullptr);
if (lincsd->nc == 0 && cr->dd == nullptr)
{
return bOK;
}
- if (econq == econqCoord)
+ if (econq == ConstraintVariable::Positions)
{
/* We can't use bCalcDHDL here, since NULL can be passed for dvdlambda
* also with efep!=fepNO.
*/
- if (ir->efep != efepNO)
+ if (ir.efep != efepNO)
{
- if (md->nMassPerturbed && lincsd->matlam != md->lambda)
+ if (md.nMassPerturbed && lincsd->matlam != md.lambda)
{
- set_lincs_matrix(lincsd, md->invmass, md->lambda);
+ set_lincs_matrix(lincsd, md.invmass, md.lambda);
}
for (i = 0; i < lincsd->nc; i++)
clear_mat(lincsd->task[th].vir_r_m_dr);
do_lincs(x, xprime, box, pbc, lincsd, th,
- md->invmass, cr,
+ md.invmass, cr,
bCalcDHDL,
- ir->LincsWarnAngle, &bWarn,
+ ir.LincsWarnAngle, &bWarn,
invdt, v, bCalcVir,
th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
}
{
cconerr(lincsd, xprime, pbc,
&ncons_loc, &p_ssd, &p_max, &p_imax);
- if (isMultiSim(ms))
+ if (isMultiSim(&ms))
{
- sprintf(buf3, " in simulation %d", ms->sim);
+ sprintf(buf3, " in simulation %d", ms.sim);
}
else
{
sprintf(buf, "\nStep %s, time %g (ps) LINCS WARNING%s\n"
"relative constraint deviation after LINCS:\n"
"rms %.6f, max %.6f (between atoms %d and %d)\n",
- gmx_step_str(step, buf2), ir->init_t+step*ir->delta_t,
+ gmx_step_str(step, buf2), ir.init_t+step*ir.delta_t,
buf3,
std::sqrt(p_ssd/ncons_loc), p_max,
ddglatnr(cr->dd, lincsd->bla[2*p_imax]),
fprintf(stderr, "%s", buf);
lincs_warning(fplog, cr->dd, x, xprime, pbc,
lincsd->nc, lincsd->bla, lincsd->bllen,
- ir->LincsWarnAngle, maxwarn, warncount);
+ ir.LincsWarnAngle, maxwarn, warncount);
}
bOK = (p_max < 0.5);
}
int th = gmx_omp_get_thread_num();
do_lincsp(x, xprime, min_proj, pbc, lincsd, th,
- md->invmass, econq, bCalcDHDL,
+ md.invmass, econq, bCalcDHDL,
bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
{
dhdlambda += lincsd->task[th].dhdlambda;
}
- if (econq == econqCoord)
+ if (econq == ConstraintVariable::Positions)
{
/* dhdlambda contains dH/dlambda*dt^2, correct for this */
/* TODO This should probably use invdt, so that sd integrator scaling works properly */
- dhdlambda /= ir->delta_t*ir->delta_t;
+ dhdlambda /= ir.delta_t*ir.delta_t;
}
*dvdlambda += dhdlambda;
}
namespace gmx
{
+enum class ConstraintVariable : int;
+
/* Abstract type for LINCS that is defined only in the file that uses it */
class Lincs;
real lincs_rmsd(const Lincs *lincsd);
/*! \brief Initializes and returns the lincs data struct. */
-Lincs *init_lincs(FILE *fplog, const gmx_mtop_t *mtop,
+Lincs *init_lincs(FILE *fplog, const gmx_mtop_t &mtop,
int nflexcon_global, const t_blocka *at2con,
bool bPLINCS, int nIter, int nProjOrder);
/*! \brief Initialize lincs stuff */
-void set_lincs(const t_idef *idef, const t_mdatoms *md,
+void set_lincs(const t_idef &idef, const t_mdatoms &md,
bool bDynamics, const t_commrec *cr,
Lincs *li);
*
* \returns true if the constraining succeeded. */
bool
-constrain_lincs(FILE *log, bool bLog, bool bEner,
- const t_inputrec *ir,
+constrain_lincs(FILE *fplog, bool bLog, bool bEner,
+ const t_inputrec &ir,
gmx_int64_t step,
- Lincs *lincsd, const t_mdatoms *md,
+ Lincs *lincsd, const t_mdatoms &md,
const t_commrec *cr,
- const gmx_multisim_t *ms,
- rvec *x, rvec *xprime, rvec *min_proj,
+ const gmx_multisim_t &ms,
+ const rvec *x, rvec *xprime, rvec *min_proj,
matrix box, t_pbc *pbc,
real lambda, real *dvdlambda,
real invdt, rvec *v,
bool bCalcVir, tensor vir_r_m_dr,
- int econ,
+ ConstraintVariable econq,
t_nrnb *nrnb,
int maxwarn, int *warncount);
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \libinternal \file
+ * \brief Declares and implements factory function for Constraints.
+ *
+ * \author Berk Hess <hess@kth.se>
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \ingroup module_mdlib
+ * \inlibraryapi
+ */
+
+#ifndef GMX_MDLIB_MAKECONSTRAINTS_H
+#define GMX_MDLIB_MAKECONSTRAINTS_H
+
+#include <memory>
+
+#include "gromacs/compat/make_unique.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/pulling/pull.h"
+#include "gromacs/topology/ifunc.h"
+#include "gromacs/topology/mtop_util.h"
+
+struct gmx_mtop_t;
+
+namespace gmx
+{
+
+/*! \internal
+ * \brief Support type to help implement makeConstraints().
+ *
+ * This member type of Constraints also inherits from it, so that it
+ * can access the private constructor of Constraints to support the
+ * implementation of the factory function. This approach avoids having
+ * to declare makeConstraints() as a template friend function. */
+struct Constraints::CreationHelper : public Constraints
+{
+ public:
+ /*! \brief Constructor that can call the private constructor
+ * of Constraints.
+ *
+ * The parameter pack insulates this helper type from changes
+ * to the arguments to the constructor. */
+ template<typename ... Args>
+ CreationHelper(Args && ... args)
+ : Constraints(std::forward<Args>(args) ...) {}
+};
+
+/*! \brief Factory function for Constraints.
+ *
+ * We only want an object to manage computing constraints when the
+ * simulation requires one. Checking for whether the object was made
+ * adds overhead to simulations that use constraints, while avoiding
+ * overhead on those that do not, so is a design trade-off we might
+ * reconsider some time.
+ *
+ * Using a private constructor and a factory function ensures that we
+ * can only make a Constraints object when the prerequisites are
+ * satisfied, ie. that something needs them and if necessary has
+ * already been initialized.
+ *
+ * Using the parameter pack insulates the factory function from
+ * changes to the type signature of the constructor that don't
+ * affect the logic here. */
+template<typename ... Args>
+std::unique_ptr<Constraints> makeConstraints(const gmx_mtop_t &mtop,
+ const t_inputrec &ir,
+ bool doEssentialDynamics,
+ Args && ... args)
+{
+ int numConstraints = (gmx_mtop_ftype_count(mtop, F_CONSTR) +
+ gmx_mtop_ftype_count(mtop, F_CONSTRNC));
+ int numSettles = gmx_mtop_ftype_count(mtop, F_SETTLE);
+ GMX_RELEASE_ASSERT(!ir.bPull || ir.pull_work != nullptr,
+ "When COM pulling is active, it must be initialized before constraints are initialized");
+ bool doPullingWithConstraints = ir.bPull && pull_have_constraint(ir.pull_work);
+ if (numConstraints + numSettles == 0 &&
+ !doPullingWithConstraints && !doEssentialDynamics)
+ {
+ // No work, so don't make a Constraints object.
+ return nullptr;
+ }
+ return compat::make_unique<Constraints::CreationHelper>
+ (mtop, ir, std::forward<Args>(args) ..., numConstraints, numSettles);
+}
+
+} // namespace
+
+#endif
add_ebin(md->ebin, md->ie, md->f_nre, ecopy, bSum);
if (md->nCrmsd)
{
- crmsd[0] = gmx::constr_rmsd(constr);
+ crmsd[0] = constr->rmsd();
add_ebin(md->ebin, md->iconrmsd, md->nCrmsd, crmsd, FALSE);
}
if (md->bDynBox)
}
}
-settledata *settle_init(const gmx_mtop_t *mtop)
+settledata *settle_init(const gmx_mtop_t &mtop)
{
/* Check that we have only one settle type */
int settle_type = -1;
*/
settled->massw.mO = -1;
- real dOH = mtop->ffparams.iparams[settle_type].settle.doh;
- real dHH = mtop->ffparams.iparams[settle_type].settle.dhh;
+ real dOH = mtop.ffparams.iparams[settle_type].settle.doh;
+ real dHH = mtop.ffparams.iparams[settle_type].settle.dhh;
settleparam_init(&settled->mass1, 1.0, 1.0, 1.0, 1.0, dOH, dHH);
settled->ow1 = nullptr;
void settle_set_constraints(settledata *settled,
const t_ilist *il_settle,
- const t_mdatoms *mdatoms)
+ const t_mdatoms &mdatoms)
{
#if GMX_SIMD_HAVE_REAL
const int pack_size = GMX_SIMD_REAL_WIDTH;
int firstO = iatoms[1];
int firstH = iatoms[2];
settleparam_init(&settled->massw,
- mdatoms->massT[firstO],
- mdatoms->massT[firstH],
- mdatoms->invmass[firstO],
- mdatoms->invmass[firstH],
+ mdatoms.massT[firstO],
+ mdatoms.massT[firstH],
+ mdatoms.invmass[firstO],
+ mdatoms.invmass[firstH],
settled->mass1.dOH,
settled->mass1.dHH);
}
* SETTLEs that appear in multiple DD domains, so we only count
* the contribution on the home range of the oxygen atom.
*/
- settled->virfac[i] = (iatoms[i*nral1 + 1] < mdatoms->homenr ? 1 : 0);
+ settled->virfac[i] = (iatoms[i*nral1 + 1] < mdatoms.homenr ? 1 : 0);
}
/* Pack the index array to the full SIMD width with copies from
}
}
-void settle_proj(settledata *settled, int econq,
+void settle_proj(settledata *settled, ConstraintVariable econq,
int nsettle, t_iatom iatoms[],
const t_pbc *pbc,
- rvec x[],
+ const rvec x[],
rvec *der, rvec *derp,
int calcvir_atom_end, tensor vir_r_m_dder)
{
calcvir_atom_end *= DIM;
- if (econq == econqForce)
+ if (econq == ConstraintVariable::Force)
{
p = &settled->mass1;
}
namespace gmx
{
+enum class ConstraintVariable : int;
+
/* Abstract type for SETTLE that is defined only in the file that uses it */
struct settledata;
/*! \brief Initializes and returns a structure with SETTLE parameters */
-settledata *settle_init(const gmx_mtop_t *mtop);
+settledata *settle_init(const gmx_mtop_t &mtop);
//! Cleans up.
void settle_free(settledata *settled);
/*! \brief Set up the indices for the settle constraints */
void settle_set_constraints(settledata *settled,
const t_ilist *il_settle,
- const t_mdatoms *mdatoms);
+ const t_mdatoms &mdatoms);
/*! \brief Constrain coordinates using SETTLE.
* Can be called on any number of threads.
/*! \brief Analytical algorithm to subtract the components of derivatives
* of coordinates working on settle type constraint.
*/
-void settle_proj(settledata *settled, int econq,
+void settle_proj(settledata *settled, ConstraintVariable econq,
int nsettle, t_iatom iatoms[],
const t_pbc *pbc, /* PBC data pointer, can be NULL */
- rvec x[],
+ const rvec x[],
rvec *der, rvec *derp,
int CalcVirAtomEnd, tensor vir_r_m_dder);
void
make_shake_sblock_serial(shakedata *shaked,
- const t_idef *idef, const t_mdatoms *md)
+ const t_idef *idef, const t_mdatoms &md)
{
int i, j, m, ncons;
int bstart, bnr;
ncons = idef->il[F_CONSTR].nr/3;
init_blocka(&sblocks);
- gen_sblocks(nullptr, 0, md->homenr, idef, &sblocks, FALSE);
+ gen_sblocks(nullptr, 0, md.homenr, idef, &sblocks, FALSE);
/*
bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
}
/* Calculate block number for each atom */
- inv_sblock = make_invblocka(&sblocks, md->nr);
+ inv_sblock = make_invblocka(&sblocks, md.nr);
done_blocka(&sblocks);
static int vec_shakef(FILE *fplog, shakedata *shaked,
const real invmass[], int ncon,
t_iparams ip[], t_iatom *iatom,
- real tol, rvec x[], rvec prime[], real omega,
+ real tol, const rvec x[], rvec prime[], real omega,
bool bFEP, real lambda, real scaled_lagrange_multiplier[],
real invdt, rvec *v,
- bool bCalcVir, tensor vir_r_m_dr, int econq)
+ bool bCalcVir, tensor vir_r_m_dr, ConstraintVariable econq)
{
rvec *rij;
real *half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared;
switch (econq)
{
- case econqCoord:
+ case ConstraintVariable::Positions:
cshake(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0], half_of_reduced_mass, omega, invmass, distance_squared_tolerance, scaled_lagrange_multiplier, &error);
break;
- case econqVeloc:
+ case ConstraintVariable::Velocities:
crattle(iatom, ncon, &nit, maxnit, constraint_distance_squared, prime[0], rij[0], half_of_reduced_mass, omega, invmass, distance_squared_tolerance, scaled_lagrange_multiplier, &error, invdt);
break;
+ default:
+ gmx_incons("Unknown constraint quantity for SHAKE");
}
if (nit >= maxnit)
i = ia[1];
j = ia[2];
- if ((econq == econqCoord) && v != nullptr)
+ if ((econq == ConstraintVariable::Positions) && v != nullptr)
{
/* Correct the velocities */
mm = scaled_lagrange_multiplier[ll]*invmass[i]*invdt;
}
//! Check that constraints are satisfied.
-static void check_cons(FILE *log, int nc, rvec x[], rvec prime[], rvec v[],
+static void check_cons(FILE *log, int nc, const rvec x[], rvec prime[], rvec v[],
t_iparams ip[], t_iatom *iatom,
- const real invmass[], int econq)
+ const real invmass[], ConstraintVariable econq)
{
t_iatom *ia;
int ai, aj;
switch (econq)
{
- case econqCoord:
+ case ConstraintVariable::Positions:
rvec_sub(prime[ai], prime[aj], dx);
dp = norm(dx);
fprintf(log, "%5d %5.2f %5d %5.2f %10.5f %10.5f %10.5f\n",
ai+1, 1.0/invmass[ai],
aj+1, 1.0/invmass[aj], d, dp, ip[ia[0]].constr.dA);
break;
- case econqVeloc:
+ case ConstraintVariable::Velocities:
rvec_sub(v[ai], v[aj], dv);
d = iprod(dx, dv);
rvec_sub(prime[ai], prime[aj], dv);
ai+1, 1.0/invmass[ai],
aj+1, 1.0/invmass[aj], d, dp, 0.);
break;
+ default:
+ gmx_incons("Unknown constraint quantity for SHAKE");
}
}
}
static bool
bshakef(FILE *log, shakedata *shaked,
const real invmass[],
- const t_idef *idef, const t_inputrec *ir, rvec x_s[], rvec prime[],
+ const t_idef &idef, const t_inputrec &ir, const rvec x_s[], rvec prime[],
t_nrnb *nrnb, real lambda, real *dvdlambda,
real invdt, rvec *v, bool bCalcVir, tensor vir_r_m_dr,
- bool bDumpOnError, int econq)
+ bool bDumpOnError, ConstraintVariable econq)
{
t_iatom *iatoms;
real *lam, dt_2, dvdl;
int i, n0, ncon, blen, type, ll;
int tnit = 0, trij = 0;
- ncon = idef->il[F_CONSTR].nr/3;
+ ncon = idef.il[F_CONSTR].nr/3;
for (ll = 0; ll < ncon; ll++)
{
// TODO Rewrite this block so that it is obvious that i, iatoms
// and lam are all iteration variables. Is this easier if the
// sblock data structure is organized differently?
- iatoms = &(idef->il[F_CONSTR].iatoms[shaked->sblock[0]]);
+ iatoms = &(idef.il[F_CONSTR].iatoms[shaked->sblock[0]]);
lam = shaked->scaled_lagrange_multiplier;
for (i = 0; (i < shaked->nblocks); )
{
blen = (shaked->sblock[i+1]-shaked->sblock[i]);
blen /= 3;
- n0 = vec_shakef(log, shaked, invmass, blen, idef->iparams,
- iatoms, ir->shake_tol, x_s, prime, shaked->omega,
- ir->efep != efepNO, lambda, lam, invdt, v, bCalcVir, vir_r_m_dr,
+ n0 = vec_shakef(log, shaked, invmass, blen, idef.iparams,
+ iatoms, ir.shake_tol, x_s, prime, shaked->omega,
+ ir.efep != efepNO, lambda, lam, invdt, v, bCalcVir, vir_r_m_dr,
econq);
if (n0 == 0)
if (bDumpOnError && log)
{
{
- check_cons(log, blen, x_s, prime, v, idef->iparams, iatoms, invmass, econq);
+ check_cons(log, blen, x_s, prime, v, idef.iparams, iatoms, invmass, econq);
}
}
return FALSE;
i++;
}
/* only for position part? */
- if (econq == econqCoord)
+ if (econq == ConstraintVariable::Positions)
{
- if (ir->efep != efepNO)
+ if (ir.efep != efepNO)
{
real bondA, bondB;
/* TODO This should probably use invdt, so that sd integrator scaling works properly */
- dt_2 = 1/gmx::square(ir->delta_t);
+ dt_2 = 1/gmx::square(ir.delta_t);
dvdl = 0;
for (ll = 0; ll < ncon; ll++)
{
- type = idef->il[F_CONSTR].iatoms[3*ll];
+ type = idef.il[F_CONSTR].iatoms[3*ll];
/* Per equations in the manual, dv/dl = -2 \sum_ll lagrangian_ll * r_ll * (d_B - d_A) */
/* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll (eta_ll is the
estimate of the Langrangian, definition on page 336 of Ryckaert et al 1977),
so the pre-factors are already present. */
- bondA = idef->iparams[type].constr.dA;
- bondB = idef->iparams[type].constr.dB;
+ bondA = idef.iparams[type].constr.dA;
+ bondB = idef.iparams[type].constr.dB;
dvdl += shaked->scaled_lagrange_multiplier[ll] * dt_2 * (bondB - bondA);
}
*dvdlambda += dvdl;
}
}
- if (ir->bShakeSOR)
+ if (ir.bShakeSOR)
{
if (tnit > shaked->gamma)
{
}
bool
-constrain_shake(FILE *log,
- shakedata *shaked,
- const real invmass[],
- const t_idef *idef,
- const t_inputrec *ir,
- rvec x_s[],
- rvec xprime[],
- rvec vprime[],
- t_nrnb *nrnb,
- real lambda,
- real *dvdlambda,
- real invdt,
- rvec *v,
- bool bCalcVir,
- tensor vir_r_m_dr,
- bool bDumpOnError,
- int econq)
+constrain_shake(FILE *log,
+ shakedata *shaked,
+ const real invmass[],
+ const t_idef &idef,
+ const t_inputrec &ir,
+ const rvec x_s[],
+ rvec xprime[],
+ rvec vprime[],
+ t_nrnb *nrnb,
+ real lambda,
+ real *dvdlambda,
+ real invdt,
+ rvec *v,
+ bool bCalcVir,
+ tensor vir_r_m_dr,
+ bool bDumpOnError,
+ ConstraintVariable econq)
{
if (shaked->nblocks == 0)
{
bool bOK;
switch (econq)
{
- case (econqCoord):
+ case (ConstraintVariable::Positions):
bOK = bshakef(log, shaked,
invmass,
idef, ir, x_s, xprime, nrnb,
invdt, v, bCalcVir, vir_r_m_dr,
bDumpOnError, econq);
break;
- case (econqVeloc):
+ case (ConstraintVariable::Velocities):
bOK = bshakef(log, shaked,
invmass,
idef, ir, x_s, vprime, nrnb,
namespace gmx
{
+enum class ConstraintVariable : int;
+
/* Abstract type for SHAKE that is defined only in the file that uses it */
struct shakedata;
//! Make SHAKE blocks when not using DD.
void
make_shake_sblock_serial(shakedata *shaked,
- const t_idef *idef, const t_mdatoms *md);
+ const t_idef *idef, const t_mdatoms &md);
//! Make SHAKE blocks when using DD.
void
* Return TRUE when OK, FALSE when shake-error
*/
bool
-constrain_shake(FILE *log, /* Log file */
- shakedata *shaked, /* Total number of atoms */
- const real invmass[], /* Atomic masses */
- const t_idef *idef, /* The interaction def */
- const t_inputrec *ir, /* Input record */
- rvec x_s[], /* Coords before update */
- rvec xprime[], /* Output coords when constraining x */
- rvec vprime[], /* Output coords when constraining v */
- t_nrnb *nrnb, /* Performance measure */
- real lambda, /* FEP lambda */
- real *dvdlambda, /* FEP force */
- real invdt, /* 1/delta_t */
- rvec *v, /* Also constrain v if v!=NULL */
- bool bCalcVir, /* Calculate r x m delta_r */
- tensor vir_r_m_dr, /* sum r x m delta_r */
- bool bDumpOnError, /* Dump debugging stuff on error*/
- int econq); /* which type of constraint is occurring */
+constrain_shake(FILE *log, /* Log file */
+ shakedata *shaked, /* Total number of atoms */
+ const real invmass[], /* Atomic masses */
+ const t_idef &idef, /* The interaction def */
+ const t_inputrec &ir, /* Input record */
+ const rvec x_s[], /* Coords before update */
+ rvec xprime[], /* Output coords when constraining x */
+ rvec vprime[], /* Output coords when constraining v */
+ t_nrnb *nrnb, /* Performance measure */
+ real lambda, /* FEP lambda */
+ real *dvdlambda, /* FEP force */
+ real invdt, /* 1/delta_t */
+ rvec *v, /* Also constrain v if v!=NULL */
+ bool bCalcVir, /* Calculate r x m delta_r */
+ tensor vir_r_m_dr, /* sum r x m delta_r */
+ bool bDumpOnError, /* Dump debugging stuff on error*/
+ ConstraintVariable econq); /* which type of constraint is occurring */
/*! \brief Regular iterative shake */
void cshake(const int iatom[], int ncon, int *nnit, int maxnit,
}
}
-static void init_adir(FILE *log,
- gmx_shellfc_t *shfc,
+static void init_adir(gmx_shellfc_t *shfc,
gmx::Constraints *constr,
- const t_idef *idef,
const t_inputrec *ir,
const t_commrec *cr,
- const gmx_multisim_t *ms,
int dd_ac1,
gmx_int64_t step,
const t_mdatoms *md,
rvec *x,
rvec *f,
rvec *acc_dir,
- gmx_bool bMolPBC,
matrix box,
gmx::ArrayRef<const real> lambda,
- real *dvdlambda,
- t_nrnb *nrnb)
+ real *dvdlambda)
{
rvec *xnold, *xnew;
double dt, w_dt;
}
}
}
- constrain(log, FALSE, FALSE, constr, idef, ir, cr, ms, step, 0, 1.0, md,
- x, xnold, nullptr, bMolPBC, box,
- lambda[efptBONDED], &(dvdlambda[efptBONDED]),
- nullptr, nullptr, nrnb, gmx::econqCoord);
- constrain(log, FALSE, FALSE, constr, idef, ir, cr, ms, step, 0, 1.0, md,
- x, xnew, nullptr, bMolPBC, box,
- lambda[efptBONDED], &(dvdlambda[efptBONDED]),
- nullptr, nullptr, nrnb, gmx::econqCoord);
+ constr->apply(FALSE, FALSE, step, 0, 1.0,
+ x, xnold, nullptr, box,
+ lambda[efptBONDED], &(dvdlambda[efptBONDED]),
+ nullptr, nullptr, gmx::ConstraintVariable::Positions);
+ constr->apply(FALSE, FALSE, step, 0, 1.0,
+ x, xnew, nullptr, box,
+ lambda[efptBONDED], &(dvdlambda[efptBONDED]),
+ nullptr, nullptr, gmx::ConstraintVariable::Positions);
for (n = 0; n < end; n++)
{
}
/* Project the acceleration on the old bond directions */
- constrain(log, FALSE, FALSE, constr, idef, ir, cr, ms, step, 0, 1.0, md,
- x_old, xnew, acc_dir, bMolPBC, box,
- lambda[efptBONDED], &(dvdlambda[efptBONDED]),
- nullptr, nullptr, nrnb, gmx::econqDeriv_FlexCon);
+ constr->apply(FALSE, FALSE, step, 0, 1.0,
+ x_old, xnew, acc_dir, box,
+ lambda[efptBONDED], &(dvdlambda[efptBONDED]),
+ nullptr, nullptr, gmx::ConstraintVariable::Deriv_FlexCon);
}
void relax_shell_flexcon(FILE *fplog,
sf_dir = 0;
if (nflexcon)
{
- init_adir(fplog, shfc,
- constr, idef, inputrec, cr, ms, dd_ac1, mdstep, md, end,
+ init_adir(shfc,
+ constr, inputrec, cr, dd_ac1, mdstep, md, end,
shfc->x_old, as_rvec_array(state->x.data()), as_rvec_array(state->x.data()), as_rvec_array(force[Min].data()),
shfc->acc_dir,
- fr->bMolPBC, state->box, state->lambda, &dum, nrnb);
+ state->box, state->lambda, &dum);
for (i = 0; i < end; i++)
{
if (nflexcon)
{
- init_adir(fplog, shfc,
- constr, idef, inputrec, cr, ms, dd_ac1, mdstep, md, end,
+ init_adir(shfc,
+ constr, inputrec, cr, dd_ac1, mdstep, md, end,
x_old, as_rvec_array(state->x.data()), as_rvec_array(pos[Min].data()), as_rvec_array(force[Min].data()), acc_dir,
- fr->bMolPBC, state->box, state->lambda, &dum, nrnb);
+ state->box, state->lambda, &dum);
directional_sd(pos[Min], pos[Try], acc_dir, end, fr->fc_stepsize);
}
sf_dir = 0;
if (nflexcon)
{
- init_adir(fplog, shfc,
- constr, idef, inputrec, cr, ms, dd_ac1, mdstep, md, end,
+ init_adir(shfc,
+ constr, inputrec, cr, dd_ac1, mdstep, md, end,
x_old, as_rvec_array(state->x.data()), as_rvec_array(pos[Try].data()), as_rvec_array(force[Try].data()), acc_dir,
- fr->bMolPBC, state->box, state->lambda, &dum, nrnb);
+ state->box, state->lambda, &dum);
for (i = 0; i < end; i++)
{
void do_constrain_first(FILE *fplog, Constraints *constr,
t_inputrec *ir, t_mdatoms *md,
- t_state *state, const t_commrec *cr,
- const gmx_multisim_t *ms,
- t_nrnb *nrnb,
- t_forcerec *fr, gmx_localtop_t *top)
+ t_state *state)
{
int i, m, start, end;
gmx_int64_t step;
dvdl_dum = 0;
/* constrain the current position */
- constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
- ir, cr, ms, step, 0, 1.0, md,
- as_rvec_array(state->x.data()), as_rvec_array(state->x.data()), nullptr,
- fr->bMolPBC, state->box,
- state->lambda[efptBONDED], &dvdl_dum,
- nullptr, nullptr, nrnb, econqCoord);
+ constr->apply(TRUE, FALSE,
+ step, 0, 1.0,
+ as_rvec_array(state->x.data()), as_rvec_array(state->x.data()), nullptr,
+ state->box,
+ state->lambda[efptBONDED], &dvdl_dum,
+ nullptr, nullptr, ConstraintVariable::Positions);
if (EI_VV(ir->eI))
{
/* constrain the inital velocity, and save it */
/* also may be useful if we need the ekin from the halfstep for velocity verlet */
- constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
- ir, cr, ms, step, 0, 1.0, md,
- as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
- fr->bMolPBC, state->box,
- state->lambda[efptBONDED], &dvdl_dum,
- nullptr, nullptr, nrnb, econqVeloc);
+ constr->apply(TRUE, FALSE,
+ step, 0, 1.0,
+ as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
+ state->box,
+ state->lambda[efptBONDED], &dvdl_dum,
+ nullptr, nullptr, ConstraintVariable::Velocities);
}
/* constrain the inital velocities at t-dt/2 */
if (EI_STATE_VELOCITY(ir->eI) && ir->eI != eiVV)
gmx_step_str(step, buf));
}
dvdl_dum = 0;
- constrain(nullptr, TRUE, FALSE, constr, &(top->idef),
- ir, cr, ms, step, -1, 1.0, md,
- as_rvec_array(state->x.data()), savex, nullptr,
- fr->bMolPBC, state->box,
- state->lambda[efptBONDED], &dvdl_dum,
- as_rvec_array(state->v.data()), nullptr, nrnb, econqCoord);
+ constr->apply(TRUE, FALSE,
+ step, -1, 1.0,
+ as_rvec_array(state->x.data()), savex, nullptr,
+ state->box,
+ state->lambda[efptBONDED], &dvdl_dum,
+ as_rvec_array(state->v.data()), nullptr, ConstraintVariable::Positions);
for (i = start; i < end; i++)
{
#include "gromacs/timing/walltime_accounting.h"
#include "gromacs/utility/arrayref.h"
-struct gmx_localtop_t;
-struct gmx_multisim_t;
struct gmx_output_env_t;
struct gmx_update_t;
struct MdrunOptions;
void do_constrain_first(FILE *log, gmx::Constraints *constr,
t_inputrec *inputrec, t_mdatoms *md,
- t_state *state, const t_commrec *cr,
- const gmx_multisim_t *ms,
- t_nrnb *nrnb,
- t_forcerec *fr, gmx_localtop_t *top);
+ t_state *state);
void init_md(FILE *fplog,
const t_commrec *cr, gmx::IMDOutputProvider *outputProvider,
ie = add_binr(rb, nener, copyenerd);
if (constr)
{
- rmsd_data = gmx::constr_rmsd_data(constr);
+ rmsd_data = constr->rmsdData();
if (rmsd_data)
{
irmsd = add_binr(rb, 2, rmsd_data);
mdatoms.homenr = numSettles * atomsPerSettle;
// Finally make the settle data structures
- settledata *settled = settle_init(&mtop);
- settle_set_constraints(settled, &mtop.moltype[0].ilist[F_SETTLE], &mdatoms);
+ settledata *settled = settle_init(mtop);
+ settle_set_constraints(settled, &mtop.moltype[0].ilist[F_SETTLE], mdatoms);
// Copy the original positions from the array of doubles to a vector of reals
std::vector<real> startingPositions(std::begin(g_positions), std::end(g_positions));
void constrain_velocities(gmx_int64_t step,
real *dvdlambda, /* the contribution to be added to the bonded interactions */
- const t_inputrec *inputrec, /* input record and box stuff */
- t_mdatoms *md,
t_state *state,
- gmx_bool bMolPBC,
- t_idef *idef,
tensor vir_part,
- const t_commrec *cr,
- const gmx_multisim_t *ms,
- t_nrnb *nrnb,
gmx_wallcycle_t wcycle,
gmx::Constraints *constr,
gmx_bool bCalcVir,
/* Constrain the coordinates upd->xp */
wallcycle_start(wcycle, ewcCONSTR);
{
- constrain(nullptr, do_log, do_ene, constr, idef,
- inputrec, cr, ms, step, 1, 1.0, md,
- as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
- bMolPBC, state->box,
- state->lambda[efptBONDED], dvdlambda,
- nullptr, bCalcVir ? &vir_con : nullptr, nrnb, econqVeloc);
+ constr->apply(do_log, do_ene,
+ step, 1, 1.0,
+ as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
+ state->box,
+ state->lambda[efptBONDED], dvdlambda,
+ nullptr, bCalcVir ? &vir_con : nullptr, ConstraintVariable::Velocities);
}
wallcycle_stop(wcycle, ewcCONSTR);
void constrain_coordinates(gmx_int64_t step,
real *dvdlambda, /* the contribution to be added to the bonded interactions */
- const t_inputrec *inputrec, /* input record and box stuff */
- t_mdatoms *md,
t_state *state,
- gmx_bool bMolPBC,
- t_idef *idef,
tensor vir_part,
- const t_commrec *cr,
- const gmx_multisim_t *ms,
- t_nrnb *nrnb,
gmx_wallcycle_t wcycle,
gmx_update_t *upd,
gmx::Constraints *constr,
/* Constrain the coordinates upd->xp */
wallcycle_start(wcycle, ewcCONSTR);
{
- constrain(nullptr, do_log, do_ene, constr, idef,
- inputrec, cr, ms, step, 1, 1.0, md,
- as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), nullptr,
- bMolPBC, state->box,
- state->lambda[efptBONDED], dvdlambda,
- as_rvec_array(state->v.data()), bCalcVir ? &vir_con : nullptr, nrnb, econqCoord);
+ constr->apply(do_log, do_ene,
+ step, 1, 1.0,
+ as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), nullptr,
+ state->box,
+ state->lambda[efptBONDED], dvdlambda,
+ as_rvec_array(state->v.data()), bCalcVir ? &vir_con : nullptr, ConstraintVariable::Positions);
}
wallcycle_stop(wcycle, ewcCONSTR);
const t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
- gmx_bool bMolPBC,
- t_idef *idef,
const t_commrec *cr,
- const gmx_multisim_t *ms,
t_nrnb *nrnb,
gmx_wallcycle_t wcycle,
gmx_update_t *upd,
{
/* Constrain the coordinates upd->xp for half a time step */
wallcycle_start(wcycle, ewcCONSTR);
- constrain(nullptr, do_log, do_ene, constr, idef,
- inputrec, cr, ms, step, 1, 0.5, md,
- as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), nullptr,
- bMolPBC, state->box,
- state->lambda[efptBONDED], dvdlambda,
- as_rvec_array(state->v.data()), nullptr, nrnb, econqCoord);
-
+ constr->apply(do_log, do_ene,
+ step, 1, 0.5,
+ as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), nullptr,
+ state->box,
+ state->lambda[efptBONDED], dvdlambda,
+ as_rvec_array(state->v.data()), nullptr, ConstraintVariable::Positions);
wallcycle_stop(wcycle, ewcCONSTR);
}
}
class ekinstate_t;
struct gmx_ekindata_t;
struct gmx_enerdata_t;
-struct gmx_multisim_t;
struct t_extmass;
struct t_fcdata;
struct t_graph;
struct t_grpopts;
-struct t_idef;
struct t_inputrec;
struct t_mdatoms;
struct t_nrnb;
void constrain_velocities(gmx_int64_t step,
real *dvdlambda, /* the contribution to be added to the bonded interactions */
- const t_inputrec *inputrec, /* input record and box stuff */
- t_mdatoms *md,
t_state *state,
- gmx_bool bMolPBC,
- t_idef *idef,
tensor vir_part,
- const t_commrec *cr,
- const gmx_multisim_t *ms,
- t_nrnb *nrnb,
gmx_wallcycle_t wcycle,
gmx::Constraints *constr,
gmx_bool bCalcVir,
void constrain_coordinates(gmx_int64_t step,
real *dvdlambda, /* the contribution to be added to the bonded interactions */
- const t_inputrec *inputrec, /* input record and box stuff */
- t_mdatoms *md,
t_state *state,
- gmx_bool bMolPBC,
- t_idef *idef,
tensor vir_part,
- const t_commrec *cr,
- const gmx_multisim_t *ms,
- t_nrnb *nrnb,
gmx_wallcycle_t wcycle,
gmx_update_t *upd,
gmx::Constraints *constr,
const t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
- gmx_bool bMolPBC,
- t_idef *idef,
const t_commrec *cr,
- const gmx_multisim_t *ms,
t_nrnb *nrnb,
gmx_wallcycle_t wcycle,
gmx_update_t *upd,
/* Check for polarizable models and flexible constraints */
shellfc = init_shell_flexcon(fplog,
- top_global, n_flexible_constraints(constr),
+ top_global, constr ? constr->numFlexibleConstraints() : 0,
ir->nstcalcenergy, DOMAINDECOMP(cr));
if (inputrecDeform(ir))
update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
}
- /* Initialize constraints */
+ /* Initialize constraints (when DD is inactive, this is only done once). */
if (constr && !DOMAINDECOMP(cr))
{
- set_constraints(constr, top, ir, mdatoms, cr);
+ constr->setConstraints(*top, *mdatoms);
}
// TODO: Remove this by converting AWH into a ForceProvider
if (constr)
{
/* Constrain the initial coordinates and velocities */
- do_constrain_first(fplog, constr, ir, mdatoms, state,
- cr, ms, nrnb, fr, top);
+ do_constrain_first(fplog, constr, ir, mdatoms, state);
}
if (vsite)
{
{
fprintf(fplog,
"RMS relative constraint deviation after constraining: %.2e\n",
- constr_rmsd(constr));
+ constr->rmsd());
}
if (EI_STATE_VELOCITY(ir->eI))
{
if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
{
wallcycle_stop(wcycle, ewcUPDATE);
- constrain_velocities(step, nullptr, ir, mdatoms,
- state, fr->bMolPBC,
- &top->idef, shake_vir,
- cr, ms, nrnb, wcycle, constr,
+ constrain_velocities(step, nullptr,
+ state,
+ shake_vir,
+ wcycle, constr,
bCalcVir, do_log, do_ene);
wallcycle_start(wcycle, ewcUPDATE);
}
/* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
if (constr && bIfRandomize)
{
- constrain_velocities(step, nullptr, ir, mdatoms,
- state, fr->bMolPBC,
- &top->idef, tmp_vir,
- cr, ms, nrnb, wcycle, constr,
+ constrain_velocities(step, nullptr,
+ state,
+ tmp_vir,
+ wcycle, constr,
bCalcVir, do_log, do_ene);
}
}
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
- constrain_coordinates(step, &dvdl_constr, ir, mdatoms, state,
- fr->bMolPBC,
- &top->idef, shake_vir,
- cr, ms, nrnb, wcycle, upd, constr,
+ constrain_coordinates(step, &dvdl_constr, state,
+ shake_vir,
+ wcycle, upd, constr,
bCalcVir, do_log, do_ene);
update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
- fr->bMolPBC, &top->idef, cr, ms,
- nrnb, wcycle, upd, constr, do_log, do_ene);
+ cr, nrnb, wcycle, upd, constr, do_log, do_ene);
finish_update(ir, mdatoms,
state, graph,
nrnb, wcycle, upd, constr);
*shellfc = init_shell_flexcon(stdout,
top_global,
- gmx::n_flexible_constraints(constr),
+ constr ? constr->numFlexibleConstraints() : 0,
ir->nstcalcenergy,
DOMAINDECOMP(cr));
}
if (constr)
{
+ // TODO how should this cross-module support dependency be managed?
if (ir->eConstrAlg == econtSHAKE &&
gmx_mtop_ftype_count(top_global, F_CONSTR) > 0)
{
if (!DOMAINDECOMP(cr))
{
- set_constraints(constr, *top, ir, mdatoms, cr);
+ constr->setConstraints(**top, *mdatoms);
}
if (!ir->bContinuation)
{
/* Constrain the starting coordinates */
dvdl_constr = 0;
- constrain(PAR(cr) ? nullptr : fplog, TRUE, TRUE, constr, &(*top)->idef,
- ir, cr, ms, -1, 0, 1.0, mdatoms,
- as_rvec_array(ems->s.x.data()),
- as_rvec_array(ems->s.x.data()),
- nullptr,
- fr->bMolPBC, ems->s.box,
- ems->s.lambda[efptFEP], &dvdl_constr,
- nullptr, nullptr, nrnb, gmx::econqCoord);
+ constr->apply(TRUE, TRUE,
+ -1, 0, 1.0,
+ as_rvec_array(ems->s.x.data()),
+ as_rvec_array(ems->s.x.data()),
+ nullptr,
+ ems->s.box,
+ ems->s.lambda[efptFEP], &dvdl_constr,
+ nullptr, nullptr, gmx::ConstraintVariable::Positions);
}
}
//
// \returns true when the step succeeded, false when a constraint error occurred
static bool do_em_step(const t_commrec *cr,
- const gmx_multisim_t *ms,
t_inputrec *ir, t_mdatoms *md,
- gmx_bool bMolPBC,
em_state_t *ems1, real a, const PaddedRVecVector *force,
em_state_t *ems2,
- gmx::Constraints *constr, gmx_localtop_t *top,
- t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ gmx::Constraints *constr,
gmx_int64_t count)
{
if (constr)
{
- wallcycle_start(wcycle, ewcCONSTR);
dvdl_constr = 0;
validStep =
- constrain(nullptr, TRUE, TRUE, constr, &top->idef,
- ir, cr, ms, count, 0, 1.0, md,
- as_rvec_array(s1->x.data()), as_rvec_array(s2->x.data()),
- nullptr, bMolPBC, s2->box,
- s2->lambda[efptBONDED], &dvdl_constr,
- nullptr, nullptr, nrnb, gmx::econqCoord);
- wallcycle_stop(wcycle, ewcCONSTR);
+ constr->apply(TRUE, TRUE,
+ count, 0, 1.0,
+ as_rvec_array(s1->x.data()), as_rvec_array(s2->x.data()),
+ nullptr, s2->box,
+ s2->lambda[efptBONDED], &dvdl_constr,
+ nullptr, nullptr, gmx::ConstraintVariable::Positions);
// We should move this check to the different minimizers
if (!validStep && ir->eI != eiSteep)
if (constr)
{
/* Project out the constraint components of the force */
- wallcycle_start(wcycle, ewcCONSTR);
dvdl_constr = 0;
rvec *f_rvec = as_rvec_array(ems->f.data());
- constrain(nullptr, FALSE, FALSE, constr, &top->idef,
- inputrec, cr, ms, count, 0, 1.0, mdAtoms->mdatoms(),
- as_rvec_array(ems->s.x.data()), f_rvec, f_rvec,
- fr->bMolPBC, ems->s.box,
- ems->s.lambda[efptBONDED], &dvdl_constr,
- nullptr, &shake_vir, nrnb, gmx::econqForceDispl);
+ constr->apply(FALSE, FALSE,
+ count, 0, 1.0,
+ as_rvec_array(ems->s.x.data()), f_rvec, f_rvec,
+ ems->s.box,
+ ems->s.lambda[efptBONDED], &dvdl_constr,
+ nullptr, &shake_vir, gmx::ConstraintVariable::ForceDispl);
enerd->term[F_DVDL_CONSTR] += dvdl_constr;
m_add(force_vir, shake_vir, vir);
- wallcycle_stop(wcycle, ewcCONSTR);
}
else
{
}
/* Take a trial step (new coords in s_c) */
- do_em_step(cr, ms, inputrec, mdatoms, fr->bMolPBC, s_min, c, &s_min->s.cg_p, s_c,
- constr, top, nrnb, wcycle, -1);
+ do_em_step(cr, inputrec, mdatoms, s_min, c, &s_min->s.cg_p, s_c,
+ constr, -1);
neval++;
/* Calculate energy for the trial step */
}
/* Take a trial step to this new point - new coords in s_b */
- do_em_step(cr, ms, inputrec, mdatoms, fr->bMolPBC, s_min, b, &s_min->s.cg_p, s_b,
- constr, top, nrnb, wcycle, -1);
+ do_em_step(cr, inputrec, mdatoms, s_min, b, &s_min->s.cg_p, s_b,
+ constr, -1);
neval++;
/* Calculate energy for the trial step */
if (count > 0)
{
validStep =
- do_em_step(cr, ms, inputrec, mdatoms, fr->bMolPBC,
+ do_em_step(cr, inputrec, mdatoms,
s_min, stepsize, &s_min->f, s_try,
- constr, top, nrnb, wcycle, count);
+ constr, count);
}
if (validStep)
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/calc_verletbuf.h"
-#include "gromacs/mdlib/constr.h"
#include "gromacs/mdlib/deform.h"
#include "gromacs/mdlib/forcerec.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/main.h"
+#include "gromacs/mdlib/makeconstraints.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdrun.h"
init_rot(fplog, inputrec, nfile, fnm, cr, globalState.get(), &mtop, oenv, mdrunOptions);
}
- /* Let init_constraints know whether we have essential dynamics constraints.
+ /* Let makeConstraints know whether we have essential dynamics constraints.
* TODO: inputrec should tell us whether we use an algorithm, not a file option or the checkpoint
*/
- bool doEdsam = (opt2fn_null("-ei", nfile, fnm) != nullptr || observablesHistory.edsamHistory);
-
- Constraints *constr = init_constraints(fplog, &mtop, inputrec, doEdsam, cr);
+ bool doEssentialDynamics = (opt2fn_null("-ei", nfile, fnm) != nullptr || observablesHistory.edsamHistory);
+ auto constr = makeConstraints(mtop, *inputrec, doEssentialDynamics,
+ fplog, *mdAtoms->mdatoms(),
+ cr, *ms, nrnb, wcycle, fr->bMolPBC);
if (DOMAINDECOMP(cr))
{
fplog, cr, ms, mdlog, nfile, fnm,
oenv,
mdrunOptions,
- vsite, constr,
+ vsite, constr.get(),
mdModules->outputProvider(),
inputrec, &mtop,
fcd,
struct t_atom;
+// TODO All of the functions taking a const gmx_mtop * are deprecated
+// and should be replaced by versions taking const gmx_mtop & when
+// their callers are refactored similarly.
+
/*! \brief Look up the molecule block and other indices of a global atom index
*
* The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
}
}
+//! \copydoc mtopGetAtomAndResidueName()
+static inline void
+mtopGetAtomAndResidueName(const gmx_mtop_t &mtop,
+ int globalAtomIndex,
+ int *moleculeBlock,
+ const char **atomName,
+ int *residueNumber,
+ const char **residueName,
+ int *globalResidueIndex)
+{
+ mtopGetAtomAndResidueName(&mtop, globalAtomIndex, moleculeBlock,
+ atomName, residueNumber, residueName, globalResidueIndex);
+}
+
/*! \brief Returns residue information for an atom based on global atom index
*
* The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
return iloop;
}
+gmx_mtop_ilistloop_t
+gmx_mtop_ilistloop_init(const gmx_mtop_t &mtop)
+{
+ return gmx_mtop_ilistloop_init(&mtop);
+}
+
static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop)
{
sfree(iloop);
return n;
}
+int gmx_mtop_ftype_count(const gmx_mtop_t &mtop, int ftype)
+{
+ return gmx_mtop_ftype_count(&mtop, ftype);
+}
+
t_block gmx_mtop_global_cgs(const gmx_mtop_t *mtop)
{
t_block cgs_gl;
struct t_ilist;
struct t_symtab;
+// TODO All of the functions taking a const gmx_mtop * are deprecated
+// and should be replaced by versions taking const gmx_mtop & when
+// their callers are refactored similarly.
+
/* Should be called after generating or reading mtop,
* to set some compute intesive variables to avoid
* N^2 operations later on.
gmx_mtop_ilistloop_t
gmx_mtop_ilistloop_init(const gmx_mtop_t *mtop);
+/* Initialize an ilist loop over all molecule types in the system. */
+gmx_mtop_ilistloop_t
+gmx_mtop_ilistloop_init(const gmx_mtop_t &mtop);
/* Loop to the next molecule,
* When not at the end:
int
gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype);
+/* Returns the total number of interactions in the system of type ftype */
+int
+gmx_mtop_ftype_count(const gmx_mtop_t &mtop, int ftype);
/* Returns a charge group index for the whole system */
t_block
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff