#include <cstdlib>
#include <algorithm>
+#include <memory>
#include "gromacs/commandline/filenm.h"
-#include "gromacs/compat/make_unique.h"
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/domdec/localatomset.h"
#include "gromacs/domdec/localatomsetmanager.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
-#include "gromacs/mdlib/mdrun.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/forceoutput.h"
#include "gromacs/mdtypes/inputrec.h"
#include "pull_internal.h"
-static int groupPbcFromParams(const t_pull_group ¶ms, bool setPbcRefToPrevStepCOM)
+namespace gmx
+{
+extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
+} // namespace gmx
+
+static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
{
if (params.nat <= 1)
{
* should not end before that of this object.
* This will be fixed when the pointers are replacted by std::vector.
*/
-pull_group_work_t::pull_group_work_t(const t_pull_group ¶ms,
+pull_group_work_t::pull_group_work_t(const t_pull_group& params,
gmx::LocalAtomSet atomSet,
bool bSetPbcRefToPrevStepCOM) :
params(params),
clear_dvec(xp);
};
-bool pull_coordinate_is_angletype(const t_pull_coord *pcrd)
+bool pull_coordinate_is_angletype(const t_pull_coord* pcrd)
{
- return (pcrd->eGeom == epullgANGLE ||
- pcrd->eGeom == epullgDIHEDRAL ||
- pcrd->eGeom == epullgANGLEAXIS);
+ return (pcrd->eGeom == epullgANGLE || pcrd->eGeom == epullgDIHEDRAL || pcrd->eGeom == epullgANGLEAXIS);
}
-static bool pull_coordinate_is_directional(const t_pull_coord *pcrd)
+static bool pull_coordinate_is_directional(const t_pull_coord* pcrd)
{
- return (pcrd->eGeom == epullgDIR ||
- pcrd->eGeom == epullgDIRPBC ||
- pcrd->eGeom == epullgDIRRELATIVE ||
- pcrd->eGeom == epullgCYL);
+ return (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC
+ || pcrd->eGeom == epullgDIRRELATIVE || pcrd->eGeom == epullgCYL);
}
-const char *pull_coordinate_units(const t_pull_coord *pcrd)
+const char* pull_coordinate_units(const t_pull_coord* pcrd)
{
- return pull_coordinate_is_angletype(pcrd) ? "deg" : "nm";
+ return pull_coordinate_is_angletype(pcrd) ? "deg" : "nm";
}
-double pull_conversion_factor_userinput2internal(const t_pull_coord *pcrd)
+double pull_conversion_factor_userinput2internal(const t_pull_coord* pcrd)
{
if (pull_coordinate_is_angletype(pcrd))
{
}
}
-double pull_conversion_factor_internal2userinput(const t_pull_coord *pcrd)
+double pull_conversion_factor_internal2userinput(const t_pull_coord* pcrd)
{
if (pull_coordinate_is_angletype(pcrd))
{
}
/* Apply forces in a mass weighted fashion for part of the pull group */
-static void apply_forces_grp_part(const pull_group_work_t *pgrp,
- int ind_start, int ind_end,
- const t_mdatoms *md,
- const dvec f_pull, int sign, rvec *f)
+static void apply_forces_grp_part(const pull_group_work_t* pgrp,
+ int ind_start,
+ int ind_end,
+ const t_mdatoms* md,
+ const dvec f_pull,
+ int sign,
+ rvec* f)
{
double inv_wm = pgrp->mwscale;
- auto localAtomIndices = pgrp->atomSet.localIndex();
+ auto localAtomIndices = pgrp->atomSet.localIndex();
for (int i = ind_start; i < ind_end; i++)
{
int ii = localAtomIndices[i];
}
/* Apply forces in a mass weighted fashion */
-static void apply_forces_grp(const pull_group_work_t *pgrp,
- const t_mdatoms *md,
- const dvec f_pull, int sign, rvec *f,
- int nthreads)
+static void apply_forces_grp(const pull_group_work_t* pgrp,
+ const t_mdatoms* md,
+ const dvec f_pull,
+ int sign,
+ rvec* f,
+ int nthreads)
{
auto localAtomIndices = pgrp->atomSet.localIndex();
*/
for (int d = 0; d < DIM; d++)
{
- f[localAtomIndices[0]][d] += sign*f_pull[d];
+ f[localAtomIndices[0]][d] += sign * f_pull[d];
}
}
else
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (int th = 0; th < nthreads; th++)
{
- int ind_start = (localAtomIndices.size()*(th + 0))/nthreads;
- int ind_end = (localAtomIndices.size()*(th + 1))/nthreads;
- apply_forces_grp_part(pgrp, ind_start, ind_end,
- md, f_pull, sign, f);
+ int ind_start = (localAtomIndices.size() * (th + 0)) / nthreads;
+ int ind_end = (localAtomIndices.size() * (th + 1)) / nthreads;
+ apply_forces_grp_part(pgrp, ind_start, ind_end, md, f_pull, sign, f);
}
}
}
}
/* Apply forces in a mass weighted fashion to a cylinder group */
-static void apply_forces_cyl_grp(const pull_group_work_t *pgrp,
- const double dv_corr,
- const t_mdatoms *md,
- const dvec f_pull, double f_scal,
- int sign, rvec *f,
+static void apply_forces_cyl_grp(const pull_group_work_t* pgrp,
+ const double dv_corr,
+ const t_mdatoms* md,
+ const dvec f_pull,
+ double f_scal,
+ int sign,
+ rvec* f,
int gmx_unused nthreads)
{
double inv_wm = pgrp->mwscale;
- auto localAtomIndices = pgrp->atomSet.localIndex();
+ auto localAtomIndices = pgrp->atomSet.localIndex();
/* The cylinder group is always a slab in the system, thus large.
* Therefore we always thread-parallelize this group.
{
continue;
}
- int ii = localAtomIndices[i];
- double mass = md->massT[ii];
+ int ii = localAtomIndices[i];
+ double mass = md->massT[ii];
/* The stored axial distance from the cylinder center (dv) needs
* to be corrected for an offset (dv_corr), which was unknown when
* we calculated dv.
*/
for (int m = 0; m < DIM; m++)
{
- f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
- pgrp->mdw[i][m]*dv_com*f_scal);
+ f[ii][m] += sign * inv_wm * (mass * weight * f_pull[m] + pgrp->mdw[i][m] * dv_com * f_scal);
}
}
}
/* Apply torque forces in a mass weighted fashion to the groups that define
* the pull vector direction for pull coordinate pcrd.
*/
-static void apply_forces_vec_torque(const struct pull_t *pull,
- const pull_coord_work_t *pcrd,
- const t_mdatoms *md,
- rvec *f)
+static void apply_forces_vec_torque(const struct pull_t* pull,
+ const pull_coord_work_t* pcrd,
+ const t_mdatoms* md,
+ rvec* f)
{
- const PullCoordSpatialData &spatialData = pcrd->spatialData;
+ const PullCoordSpatialData& spatialData = pcrd->spatialData;
/* The component inpr along the pull vector is accounted for in the usual
* way. Here we account for the component perpendicular to vec.
double inpr = 0;
for (int m = 0; m < DIM; m++)
{
- inpr += spatialData.dr01[m]*spatialData.vec[m];
+ inpr += spatialData.dr01[m] * spatialData.vec[m];
}
/* The torque force works along the component of the distance vector
* of between the two "usual" pull groups that is perpendicular to
dvec f_perp;
for (int m = 0; m < DIM; m++)
{
- f_perp[m] = (spatialData.dr01[m] - inpr*spatialData.vec[m])/spatialData.vec_len*pcrd->scalarForce;
+ f_perp[m] = (spatialData.dr01[m] - inpr * spatialData.vec[m]) / spatialData.vec_len * pcrd->scalarForce;
}
/* Apply the force to the groups defining the vector using opposite signs */
- apply_forces_grp(&pull->group[pcrd->params.group[2]], md,
- f_perp, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd->params.group[3]], md,
- f_perp, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd->params.group[2]], md, f_perp, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd->params.group[3]], md, f_perp, 1, f, pull->nthreads);
}
/* Apply forces in a mass weighted fashion */
-static void apply_forces_coord(struct pull_t * pull, int coord,
- const PullCoordVectorForces &forces,
- const t_mdatoms * md,
- rvec *f)
+static void apply_forces_coord(struct pull_t* pull,
+ int coord,
+ const PullCoordVectorForces& forces,
+ const t_mdatoms* md,
+ rvec* f)
{
/* Here it would be more efficient to use one large thread-parallel
* region instead of potential parallel regions within apply_forces_grp.
* to data races.
*/
- const pull_coord_work_t &pcrd = pull->coord[coord];
+ const pull_coord_work_t& pcrd = pull->coord[coord];
if (pcrd.params.eGeom == epullgCYL)
{
- apply_forces_cyl_grp(&pull->dyna[coord], pcrd.spatialData.cyl_dev, md,
- forces.force01, pcrd.scalarForce, -1, f,
- pull->nthreads);
+ apply_forces_cyl_grp(&pull->dyna[coord], pcrd.spatialData.cyl_dev, md, forces.force01,
+ pcrd.scalarForce, -1, f, pull->nthreads);
/* Sum the force along the vector and the radial force */
dvec f_tot;
for (int m = 0; m < DIM; m++)
{
- f_tot[m] = forces.force01[m] + pcrd.scalarForce*pcrd.spatialData.ffrad[m];
+ f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
}
- apply_forces_grp(&pull->group[pcrd.params.group[1]], md,
- f_tot, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[1]], md, f_tot, 1, f, pull->nthreads);
}
else
{
if (pull->group[pcrd.params.group[0]].params.nat > 0)
{
- apply_forces_grp(&pull->group[pcrd.params.group[0]], md,
- forces.force01, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[0]], md, forces.force01, -1, f, pull->nthreads);
}
- apply_forces_grp(&pull->group[pcrd.params.group[1]], md,
- forces.force01, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[1]], md, forces.force01, 1, f, pull->nthreads);
if (pcrd.params.ngroup >= 4)
{
- apply_forces_grp(&pull->group[pcrd.params.group[2]], md,
- forces.force23, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd.params.group[3]], md,
- forces.force23, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[2]], md, forces.force23, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[3]], md, forces.force23, 1, f, pull->nthreads);
}
if (pcrd.params.ngroup >= 6)
{
- apply_forces_grp(&pull->group[pcrd.params.group[4]], md,
- forces.force45, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd.params.group[5]], md,
- forces.force45, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[4]], md, forces.force45, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[5]], md, forces.force45, 1, f, pull->nthreads);
}
}
}
-real max_pull_distance2(const pull_coord_work_t *pcrd,
- const t_pbc *pbc)
+real max_pull_distance2(const pull_coord_work_t* pcrd, const t_pbc* pbc)
{
/* Note that this maximum distance calculation is more complex than
* most other cases in GROMACS, since here we have to take care of
}
}
- return 0.25*max_d2;
+ return 0.25 * max_d2;
}
/* This function returns the distance based on coordinates xg and xref.
* Note that the pull coordinate struct pcrd is not modified.
*/
-static void low_get_pull_coord_dr(const struct pull_t *pull,
- const pull_coord_work_t *pcrd,
- const t_pbc *pbc,
- dvec xg, dvec xref, double max_dist2,
- dvec dr)
+static void low_get_pull_coord_dr(const struct pull_t* pull,
+ const pull_coord_work_t* pcrd,
+ const t_pbc* pbc,
+ dvec xg,
+ dvec xref,
+ double max_dist2,
+ dvec dr)
{
- const pull_group_work_t *pgrp0 = &pull->group[pcrd->params.group[0]];
+ const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
/* Only the first group can be an absolute reference, in that case nat=0 */
if (pgrp0->params.nat == 0)
dvec xrefr;
copy_dvec(xref, xrefr);
- dvec dref = {0, 0, 0};
+ dvec dref = { 0, 0, 0 };
if (pcrd->params.eGeom == epullgDIRPBC)
{
for (int m = 0; m < DIM; m++)
{
- dref[m] = pcrd->value_ref*pcrd->spatialData.vec[m];
+ dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
}
/* Add the reference position, so we use the correct periodic image */
dvec_inc(xrefr, dref);
dr[m] *= pcrd->params.dim[m];
if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
{
- dr2 += dr[m]*dr[m];
+ dr2 += dr[m] * dr[m];
}
}
/* Check if we are close to switching to another periodic image */
- if (max_dist2 > 0 && dr2 > 0.98*0.98*max_dist2)
+ if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
{
/* Note that technically there is no issue with switching periodic
* image, as pbc_dx_d returns the distance to the closest periodic
* image. However in all cases where periodic image switches occur,
* the pull results will be useless in practice.
*/
- gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f).\n%s",
- pcrd->params.group[0], pcrd->params.group[1],
- sqrt(dr2), sqrt(0.98*0.98*max_dist2),
- pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
+ gmx_fatal(FARGS,
+ "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
+ "box size (%f).\n%s",
+ pcrd->params.group[0], pcrd->params.group[1], sqrt(dr2),
+ sqrt(0.98 * 0.98 * max_dist2), pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
}
if (pcrd->params.eGeom == epullgDIRPBC)
/* This function returns the distance based on the contents of the pull struct.
* pull->coord[coord_ind].dr, and potentially vector, are updated.
*/
-static void get_pull_coord_dr(struct pull_t *pull,
- int coord_ind,
- const t_pbc *pbc)
+static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
{
- pull_coord_work_t *pcrd = &pull->coord[coord_ind];
- PullCoordSpatialData &spatialData = pcrd->spatialData;
+ pull_coord_work_t* pcrd = &pull->coord[coord_ind];
+ PullCoordSpatialData& spatialData = pcrd->spatialData;
- double md2;
+ double md2;
/* With AWH pulling we allow for periodic pulling with geometry=direction.
* TODO: Store a periodicity flag instead of checking for external pull provider.
*/
- if (pcrd->params.eGeom == epullgDIRPBC || (pcrd->params.eGeom == epullgDIR &&
- pcrd->params.eType == epullEXTERNAL))
+ if (pcrd->params.eGeom == epullgDIRPBC
+ || (pcrd->params.eGeom == epullgDIR && pcrd->params.eType == epullEXTERNAL))
{
md2 = -1;
}
if (pcrd->params.eGeom == epullgDIRRELATIVE)
{
/* We need to determine the pull vector */
- dvec vec;
- int m;
+ dvec vec;
+ int m;
- const pull_group_work_t &pgrp2 = pull->group[pcrd->params.group[2]];
- const pull_group_work_t &pgrp3 = pull->group[pcrd->params.group[3]];
+ const pull_group_work_t& pgrp2 = pull->group[pcrd->params.group[2]];
+ const pull_group_work_t& pgrp3 = pull->group[pcrd->params.group[3]];
pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
spatialData.vec_len = dnorm(vec);
for (m = 0; m < DIM; m++)
{
- spatialData.vec[m] = vec[m]/spatialData.vec_len;
+ spatialData.vec[m] = vec[m] / spatialData.vec_len;
}
if (debug)
{
fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
- coord_ind,
- vec[XX], vec[YY], vec[ZZ],
- spatialData.vec[XX], spatialData.vec[YY], spatialData.vec[ZZ]);
+ coord_ind, vec[XX], vec[YY], vec[ZZ], spatialData.vec[XX], spatialData.vec[YY],
+ spatialData.vec[ZZ]);
}
}
- pull_group_work_t *pgrp0 = &pull->group[pcrd->params.group[0]];
- pull_group_work_t *pgrp1 = &pull->group[pcrd->params.group[1]];
+ pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
+ pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
- low_get_pull_coord_dr(pull, pcrd, pbc,
- pgrp1->x,
- pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x,
- md2,
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp1->x,
+ pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, md2,
spatialData.dr01);
if (pcrd->params.ngroup >= 4)
pgrp2 = &pull->group[pcrd->params.group[2]];
pgrp3 = &pull->group[pcrd->params.group[3]];
- low_get_pull_coord_dr(pull, pcrd, pbc,
- pgrp3->x,
- pgrp2->x,
- md2,
- spatialData.dr23);
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, md2, spatialData.dr23);
}
if (pcrd->params.ngroup >= 6)
{
pgrp4 = &pull->group[pcrd->params.group[4]];
pgrp5 = &pull->group[pcrd->params.group[5]];
- low_get_pull_coord_dr(pull, pcrd, pbc,
- pgrp5->x,
- pgrp4->x,
- md2,
- spatialData.dr45);
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, md2, spatialData.dr45);
}
}
* It is assumed that x is in [-3pi, 3pi) so that x
* needs to be shifted by at most one period.
*/
-static void make_periodic_2pi(double *x)
+static void make_periodic_2pi(double* x)
{
if (*x >= M_PI)
{
}
/* This function should always be used to modify pcrd->value_ref */
-static void low_set_pull_coord_reference_value(pull_coord_work_t *pcrd,
- int coord_ind,
- double value_ref)
+static void low_set_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double value_ref)
{
- GMX_ASSERT(pcrd->params.eType != epullEXTERNAL, "The pull coord reference value should not be used with type external-potential");
+ GMX_ASSERT(pcrd->params.eType != epullEXTERNAL,
+ "The pull coord reference value should not be used with type external-potential");
if (pcrd->params.eGeom == epullgDIST)
{
if (value_ref < 0)
{
- gmx_fatal(FARGS, "Pull reference distance for coordinate %d (%f) needs to be non-negative", coord_ind + 1, value_ref);
+ gmx_fatal(FARGS,
+ "Pull reference distance for coordinate %d (%f) needs to be non-negative",
+ coord_ind + 1, value_ref);
}
}
else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgANGLEAXIS)
{
if (value_ref < 0 || value_ref > M_PI)
{
- gmx_fatal(FARGS, "Pull reference angle for coordinate %d (%f) needs to be in the allowed interval [0,180] deg", coord_ind + 1, value_ref*pull_conversion_factor_internal2userinput(&pcrd->params));
+ gmx_fatal(FARGS,
+ "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
+ "interval [0,180] deg",
+ coord_ind + 1,
+ value_ref * pull_conversion_factor_internal2userinput(&pcrd->params));
}
}
else if (pcrd->params.eGeom == epullgDIHEDRAL)
pcrd->value_ref = value_ref;
}
-static void update_pull_coord_reference_value(pull_coord_work_t *pcrd, int coord_ind, double t)
+static void update_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double t)
{
/* With zero rate the reference value is set initially and doesn't change */
if (pcrd->params.rate != 0)
{
- double value_ref = (pcrd->params.init + pcrd->params.rate*t)*pull_conversion_factor_userinput2internal(&pcrd->params);
+ double value_ref = (pcrd->params.init + pcrd->params.rate * t)
+ * pull_conversion_factor_userinput2internal(&pcrd->params);
low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
}
}
/* Returns the dihedral angle. Updates the plane normal vectors m, n. */
-static double get_dihedral_angle_coord(PullCoordSpatialData *spatialData)
+static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
{
double phi, sign;
dvec dr32; /* store instead of dr23? */
dsvmul(-1, spatialData->dr23, dr32);
- dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
- dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
+ dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
+ dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
/* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
* we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
*/
sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
- return sign*phi;
+ return sign * phi;
}
/* Calculates pull->coord[coord_ind].value.
* This function also updates pull->coord[coord_ind].dr.
*/
-static void get_pull_coord_distance(struct pull_t *pull,
- int coord_ind,
- const t_pbc *pbc)
+static void get_pull_coord_distance(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
int m;
pcrd = &pull->coord[coord_ind];
get_pull_coord_dr(pull, coord_ind, pbc);
- PullCoordSpatialData &spatialData = pcrd->spatialData;
+ PullCoordSpatialData& spatialData = pcrd->spatialData;
switch (pcrd->params.eGeom)
{
spatialData.value = 0;
for (m = 0; m < DIM; m++)
{
- spatialData.value += spatialData.vec[m]*spatialData.dr01[m];
+ spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
}
break;
case epullgANGLE:
spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
break;
- case epullgDIHEDRAL:
- spatialData.value = get_dihedral_angle_coord(&spatialData);
- break;
+ case epullgDIHEDRAL: spatialData.value = get_dihedral_angle_coord(&spatialData); break;
case epullgANGLEAXIS:
spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
break;
- default:
- gmx_incons("Unsupported pull type in get_pull_coord_distance");
+ default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
}
}
/* Returns the deviation from the reference value.
* Updates pull->coord[coord_ind].dr, .value and .value_ref.
*/
-static double get_pull_coord_deviation(struct pull_t *pull,
- int coord_ind,
- const t_pbc *pbc,
- double t)
+static double get_pull_coord_deviation(struct pull_t* pull, int coord_ind, const t_pbc* pbc, double t)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
double dev = 0;
pcrd = &pull->coord[coord_ind];
return dev;
}
-double get_pull_coord_value(struct pull_t *pull,
- int coord_ind,
- const t_pbc *pbc)
+double get_pull_coord_value(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
{
get_pull_coord_distance(pull, coord_ind, pbc);
return pull->coord[coord_ind].spatialData.value;
}
-void clear_pull_forces(pull_t *pull)
+void clear_pull_forces(pull_t* pull)
{
/* Zeroing the forces is only required for constraint pulling.
* It can happen that multiple constraint steps need to be applied
* and therefore the constraint forces need to be accumulated.
*/
- for (pull_coord_work_t &coord : pull->coord)
+ for (pull_coord_work_t& coord : pull->coord)
{
coord.scalarForce = 0;
}
}
/* Apply constraint using SHAKE */
-static void do_constraint(struct pull_t *pull, t_pbc *pbc,
- rvec *x, rvec *v,
- gmx_bool bMaster, tensor vir,
- double dt, double t)
+static void
+do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
{
- dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
- dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
- dvec *rnew; /* current 'new' positions of the groups */
- double *dr_tot; /* the total update of the coords */
- dvec vec;
- double inpr;
- double lambda, rm, invdt = 0;
- gmx_bool bConverged_all, bConverged = FALSE;
- int niter = 0, ii, j, m, max_iter = 100;
- double a;
- dvec tmp, tmp3;
-
- snew(r_ij, pull->coord.size());
+ dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
+ dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
+ dvec* rnew; /* current 'new' positions of the groups */
+ double* dr_tot; /* the total update of the coords */
+ dvec vec;
+ double inpr;
+ double lambda, rm, invdt = 0;
+ gmx_bool bConverged_all, bConverged = FALSE;
+ int niter = 0, ii, j, m, max_iter = 100;
+ double a;
+ dvec tmp, tmp3;
+
+ snew(r_ij, pull->coord.size());
snew(dr_tot, pull->coord.size());
snew(rnew, pull->group.size());
/* Determine the constraint directions from the old positions */
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
pcrd = &pull->coord[c];
*/
get_pull_coord_distance(pull, c, pbc);
- const PullCoordSpatialData &spatialData = pcrd->spatialData;
+ const PullCoordSpatialData& spatialData = pcrd->spatialData;
if (debug)
{
- fprintf(debug, "Pull coord %zu dr %f %f %f\n",
- c, spatialData.dr01[XX], spatialData.dr01[YY], spatialData.dr01[ZZ]);
+ fprintf(debug, "Pull coord %zu dr %f %f %f\n", c, spatialData.dr01[XX],
+ spatialData.dr01[YY], spatialData.dr01[ZZ]);
}
- if (pcrd->params.eGeom == epullgDIR ||
- pcrd->params.eGeom == epullgDIRPBC)
+ if (pcrd->params.eGeom == epullgDIR || pcrd->params.eGeom == epullgDIRPBC)
{
/* Select the component along vec */
a = 0;
for (m = 0; m < DIM; m++)
{
- a += spatialData.vec[m]*spatialData.dr01[m];
+ a += spatialData.vec[m] * spatialData.dr01[m];
}
for (m = 0; m < DIM; m++)
{
- r_ij[c][m] = a*spatialData.vec[m];
+ r_ij[c][m] = a * spatialData.vec[m];
}
}
else
if (dnorm2(r_ij[c]) == 0)
{
- gmx_fatal(FARGS, "Distance for pull coordinate %zu is zero with constraint pulling, which is not allowed.", c + 1);
+ gmx_fatal(FARGS,
+ "Distance for pull coordinate %zu is zero with constraint pulling, which is "
+ "not allowed.",
+ c + 1);
}
}
/* loop over all constraints */
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
pull_group_work_t *pgrp0, *pgrp1;
dvec dr0, dr1;
- pcrd = &pull->coord[c];
+ pcrd = &pull->coord[c];
if (pcrd->params.eType != epullCONSTRAINT)
{
pgrp1 = &pull->group[pcrd->params.group[1]];
/* Get the current difference vector */
- low_get_pull_coord_dr(pull, pcrd, pbc,
- rnew[pcrd->params.group[1]],
- rnew[pcrd->params.group[0]],
- -1, unc_ij);
+ low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->params.group[1]],
+ rnew[pcrd->params.group[0]], -1, unc_ij);
if (debug)
{
fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
}
- rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
+ rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
switch (pcrd->params.eGeom)
{
case epullgDIST:
if (pcrd->value_ref <= 0)
{
- gmx_fatal(FARGS, "The pull constraint reference distance for group %zu is <= 0 (%f)", c, pcrd->value_ref);
+ gmx_fatal(
+ FARGS,
+ "The pull constraint reference distance for group %zu is <= 0 (%f)",
+ c, pcrd->value_ref);
}
{
double q, c_a, c_b, c_c;
c_a = diprod(r_ij[c], r_ij[c]);
- c_b = diprod(unc_ij, r_ij[c])*2;
+ c_b = diprod(unc_ij, r_ij[c]) * 2;
c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
if (c_b < 0)
{
- q = -0.5*(c_b - std::sqrt(c_b*c_b - 4*c_a*c_c));
- lambda = -q/c_a;
+ 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 + std::sqrt(c_b*c_b - 4*c_a*c_c));
- lambda = -c_c/q;
+ q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
+ lambda = -c_c / q;
}
if (debug)
{
- fprintf(debug,
- "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
- c_a, c_b, c_c, lambda);
+ fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b,
+ c_c, lambda);
}
}
/* The position corrections dr due to the constraints */
- dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
- dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
- dr_tot[c] += -lambda*dnorm(r_ij[c]);
+ dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
+ dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
+ dr_tot[c] += -lambda * dnorm(r_ij[c]);
break;
case epullgDIR:
case epullgDIRPBC:
for (m = 0; m < DIM; m++)
{
vec[m] = pcrd->spatialData.vec[m];
- a += unc_ij[m]*vec[m];
+ a += unc_ij[m] * vec[m];
}
/* Select only the component along the vector */
dsvmul(a, vec, unc_ij);
}
/* The position corrections dr due to the constraints */
- dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
- dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
+ dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
+ dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
dr_tot[c] += -lambda;
break;
- default:
- gmx_incons("Invalid enumeration value for eGeom");
+ default: gmx_incons("Invalid enumeration value for eGeom");
}
/* DEBUG */
g1 = pcrd->params.group[1];
low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
- fprintf(debug,
- "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
- rnew[g0][0], rnew[g0][1], rnew[g0][2],
- rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
- fprintf(debug,
- "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
- "", "", "", "", "", "", pcrd->value_ref);
- fprintf(debug,
- "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
- dr0[0], dr0[1], dr0[2],
- dr1[0], dr1[1], dr1[2],
- dnorm(tmp3));
+ fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", rnew[g0][0],
+ rnew[g0][1], rnew[g0][2], rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
+ fprintf(debug, "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", "", "", "", "", "",
+ "", pcrd->value_ref);
+ fprintf(debug, "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", dr0[0],
+ dr0[1], dr0[2], dr1[0], dr1[1], dr1[2], dnorm(tmp3));
} /* END DEBUG */
/* Update the COMs with dr */
}
/* Check if all constraints are fullfilled now */
- for (pull_coord_work_t &coord : pull->coord)
+ for (pull_coord_work_t& coord : pull->coord)
{
if (coord.params.eType != epullCONSTRAINT)
{
continue;
}
- low_get_pull_coord_dr(pull, &coord, pbc,
- rnew[coord.params.group[1]],
- rnew[coord.params.group[0]],
- -1, unc_ij);
+ low_get_pull_coord_dr(pull, &coord, pbc, rnew[coord.params.group[1]],
+ rnew[coord.params.group[0]], -1, unc_ij);
switch (coord.params.eGeom)
{
case epullgDIST:
- bConverged =
- fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
+ bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
break;
case epullgDIR:
case epullgDIRPBC:
}
inpr = diprod(unc_ij, vec);
dsvmul(inpr, vec, unc_ij);
- bConverged =
- fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
+ bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
break;
}
{
if (debug)
{
- fprintf(debug, "Pull constraint not converged: "
+ fprintf(debug,
+ "Pull constraint not converged: "
"groups %d %d,"
"d_ref = %f, current d = %f\n",
- coord.params.group[0], coord.params.group[1],
- coord.value_ref, dnorm(unc_ij));
+ coord.params.group[0], coord.params.group[1], coord.value_ref, dnorm(unc_ij));
}
bConverged_all = FALSE;
if (v)
{
- invdt = 1/dt;
+ invdt = 1 / dt;
}
/* update atoms in the groups */
for (size_t g = 0; g < pull->group.size(); g++)
{
- const pull_group_work_t *pgrp;
+ const pull_group_work_t* pgrp;
dvec dr;
pgrp = &pull->group[g];
/* update the atom positions */
auto localAtomIndices = pgrp->atomSet.localIndex();
copy_dvec(dr, tmp);
- for (gmx::index j = 0; j < localAtomIndices.size(); j++)
+ for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
{
ii = localAtomIndices[j];
if (!pgrp->localWeights.empty())
{
- dsvmul(pgrp->wscale*pgrp->localWeights[j], dr, tmp);
+ dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
}
for (m = 0; m < DIM; m++)
{
{
for (m = 0; m < DIM; m++)
{
- v[ii][m] += invdt*tmp[m];
+ v[ii][m] += invdt * tmp[m];
}
}
}
/* calculate the constraint forces, used for output and virial only */
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
pcrd = &pull->coord[c];
}
/* Accumulate the forces, in case we have multiple constraint steps */
- double force = dr_tot[c]/((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)*dt*dt);
+ double force =
+ dr_tot[c]
+ / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
+ * dt * dt);
pcrd->scalarForce += force;
if (vir != nullptr && pcrd->params.eGeom != epullgDIRPBC && bMaster)
/* Add the pull contribution to the virial */
/* We have already checked above that r_ij[c] != 0 */
- f_invr = pcrd->scalarForce/dnorm(r_ij[c]);
+ f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
for (j = 0; j < DIM; j++)
{
for (m = 0; m < DIM; m++)
{
- vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
+ vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
}
}
}
{
for (int m = 0; m < DIM; m++)
{
- vir[j][m] -= 0.5*f[j]*dr[m];
+ vir[j][m] -= 0.5 * f[j] * dr[m];
}
}
}
/* Adds the pull contribution to the virial */
-static void add_virial_coord(tensor vir,
- const pull_coord_work_t &pcrd,
- const PullCoordVectorForces &forces)
+static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
{
if (vir != nullptr && pcrd.params.eGeom != epullgDIRPBC)
{
}
}
-static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t *pcrd,
- double dev, real lambda,
- real *V, real *dVdl)
+static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
+ double dev,
+ real lambda,
+ real* V,
+ real* dVdl)
{
- real k, dkdl;
+ real k, dkdl;
- k = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
+ k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
dkdl = pcrd->params.kB - pcrd->params.k;
switch (pcrd->params.eType)
* potential is that a flat-bottom is zero above or below
the reference value.
*/
- if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0) ||
- (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
+ if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0)
+ || (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
{
dev = 0;
}
- pcrd->scalarForce = -k*dev;
- *V += 0.5* k*gmx::square(dev);
- *dVdl += 0.5*dkdl*gmx::square(dev);
+ pcrd->scalarForce = -k * dev;
+ *V += 0.5 * k * gmx::square(dev);
+ *dVdl += 0.5 * dkdl * gmx::square(dev);
break;
case epullCONST_F:
- pcrd->scalarForce = -k;
- *V += k*pcrd->spatialData.value;
- *dVdl += dkdl*pcrd->spatialData.value;
+ pcrd->scalarForce = -k;
+ *V += k * pcrd->spatialData.value;
+ *dVdl += dkdl * pcrd->spatialData.value;
break;
case epullEXTERNAL:
- gmx_incons("the scalar pull force should not be calculated internally for pull type external");
- default:
- gmx_incons("Unsupported pull type in do_pull_pot");
+ gmx_incons(
+ "the scalar pull force should not be calculated internally for pull type "
+ "external");
+ default: gmx_incons("Unsupported pull type in do_pull_pot");
}
}
-static PullCoordVectorForces
-calculateVectorForces(const pull_coord_work_t &pcrd)
+static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
{
- const t_pull_coord ¶ms = pcrd.params;
- const PullCoordSpatialData &spatialData = pcrd.spatialData;
+ const t_pull_coord& params = pcrd.params;
+ const PullCoordSpatialData& spatialData = pcrd.spatialData;
/* The geometry of the coordinate determines how the scalar force relates to the force on each group */
- PullCoordVectorForces forces;
+ PullCoordVectorForces forces;
if (params.eGeom == epullgDIST)
{
- double invdr01 = spatialData.value > 0 ? 1./spatialData.value : 0.;
+ double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
for (int m = 0; m < DIM; m++)
{
- forces.force01[m] = pcrd.scalarForce*spatialData.dr01[m]*invdr01;
+ forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
}
}
else if (params.eGeom == epullgANGLE)
* f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
*/
double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
- double b = a*cos_theta;
- double invdr01 = 1./dnorm(spatialData.dr01);
- double invdr23 = 1./dnorm(spatialData.dr23);
+ double b = a * cos_theta;
+ double invdr01 = 1. / dnorm(spatialData.dr01);
+ double invdr23 = 1. / dnorm(spatialData.dr23);
dvec normalized_dr01, normalized_dr23;
dsvmul(invdr01, spatialData.dr01, normalized_dr01);
dsvmul(invdr23, spatialData.dr23, normalized_dr23);
for (int m = 0; m < DIM; m++)
{
/* Here, f_scal is -dV/dtheta */
- forces.force01[m] = pcrd.scalarForce*invdr01*(a*normalized_dr23[m] - b*normalized_dr01[m]);
- forces.force23[m] = pcrd.scalarForce*invdr23*(a*normalized_dr01[m] - b*normalized_dr23[m]);
+ forces.force01[m] =
+ pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
+ forces.force23[m] =
+ pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
}
}
else
double invdr01;
dvec normalized_dr01;
- invdr01 = 1./dnorm(spatialData.dr01);
+ invdr01 = 1. / dnorm(spatialData.dr01);
dsvmul(invdr01, spatialData.dr01, normalized_dr01);
- a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
- b = a*cos_theta;
+ a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
+ b = a * cos_theta;
for (int m = 0; m < DIM; m++)
{
- forces.force01[m] = pcrd.scalarForce*invdr01*(a*spatialData.vec[m] - b*normalized_dr01[m]);
+ forces.force01[m] =
+ pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
}
}
else
n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
dsvmul(-1, spatialData.dr23, dr32);
sqrdist_32 = diprod(dr32, dr32);
- tol = sqrdist_32*GMX_REAL_EPS; /* Avoid tiny angles */
+ tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
if ((m2 > tol) && (n2 > tol))
{
double a_01, a_23_01, a_23_45, a_45;
double inv_dist_32, inv_sqrdist_32, dist_32;
dvec u, v;
inv_dist_32 = gmx::invsqrt(sqrdist_32);
- inv_sqrdist_32 = inv_dist_32*inv_dist_32;
- dist_32 = sqrdist_32*inv_dist_32;
+ inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
+ dist_32 = sqrdist_32 * inv_dist_32;
/* Forces on groups 0, 1 */
- a_01 = pcrd.scalarForce*dist_32/m2; /* scalarForce is -dV/dphi */
- dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
+ a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
+ dsvmul(-a_01, spatialData.planevec_m,
+ forces.force01); /* added sign to get force on group 1, not 0 */
/* Forces on groups 4, 5 */
- a_45 = -pcrd.scalarForce*dist_32/n2;
- dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
+ a_45 = -pcrd.scalarForce * dist_32 / n2;
+ dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
/* Force on groups 2, 3 (defining the axis) */
- a_23_01 = -diprod(spatialData.dr01, dr32)*inv_sqrdist_32;
- a_23_45 = -diprod(spatialData.dr45, dr32)*inv_sqrdist_32;
+ a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
+ a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
dsvmul(a_23_45, forces.force45, v);
- dvec_sub(u, v, forces.force23); /* force on group 3 */
+ dvec_sub(u, v, forces.force23); /* force on group 3 */
}
else
{
{
for (int m = 0; m < DIM; m++)
{
- forces.force01[m] = pcrd.scalarForce*spatialData.vec[m];
+ forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
}
}
using Lock = gmx::lock_guard<gmx::Mutex>;
-void register_external_pull_potential(struct pull_t *pull,
- int coord_index,
- const char *provider)
+void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
{
GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
- GMX_RELEASE_ASSERT(provider != nullptr, "register_external_pull_potential called with NULL as provider name");
+ GMX_RELEASE_ASSERT(provider != nullptr,
+ "register_external_pull_potential called with NULL as provider name");
- if (coord_index < 0 || coord_index >= static_cast<int>(pull->coord.size()))
+ if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
{
- gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d which is out of the pull coordinate range %d - %zu\n",
+ gmx_fatal(FARGS,
+ "Module '%s' attempted to register an external potential for pull coordinate %d "
+ "which is out of the pull coordinate range %d - %zu\n",
provider, coord_index + 1, 1, pull->coord.size());
}
- pull_coord_work_t *pcrd = &pull->coord[coord_index];
+ pull_coord_work_t* pcrd = &pull->coord[coord_index];
if (pcrd->params.eType != epullEXTERNAL)
{
- gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d which of type '%s', whereas external potentials are only supported with type '%s'",
- provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
+ gmx_fatal(
+ FARGS,
+ "Module '%s' attempted to register an external potential for pull coordinate %d "
+ "which of type '%s', whereas external potentials are only supported with type '%s'",
+ provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
}
- GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr, "The external potential provider string for a pull coordinate is NULL");
+ GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr,
+ "The external potential provider string for a pull coordinate is NULL");
if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider) != 0)
{
- gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d which expects the external potential to be provided by a module named '%s'",
+ gmx_fatal(FARGS,
+ "Module '%s' attempted to register an external potential for pull coordinate %d "
+ "which expects the external potential to be provided by a module named '%s'",
provider, coord_index + 1, pcrd->params.externalPotentialProvider);
}
if (pcrd->bExternalPotentialProviderHasBeenRegistered)
{
- gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d more than once",
+ gmx_fatal(FARGS,
+ "Module '%s' attempted to register an external potential for pull coordinate %d "
+ "more than once",
provider, coord_index + 1);
}
pcrd->bExternalPotentialProviderHasBeenRegistered = true;
pull->numUnregisteredExternalPotentials--;
- GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0, "Negative unregistered potentials, the pull code is inconsistent");
+ GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
+ "Negative unregistered potentials, the pull code is inconsistent");
}
-static void check_external_potential_registration(const struct pull_t *pull)
+static void check_external_potential_registration(const struct pull_t* pull)
{
if (pull->numUnregisteredExternalPotentials > 0)
{
size_t c;
for (c = 0; c < pull->coord.size(); c++)
{
- if (pull->coord[c].params.eType == epullEXTERNAL &&
- !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
+ if (pull->coord[c].params.eType == epullEXTERNAL
+ && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
{
break;
}
}
- GMX_RELEASE_ASSERT(c < pull->coord.size(), "Internal inconsistency in the pull potential provider counting");
+ GMX_RELEASE_ASSERT(c < pull->coord.size(),
+ "Internal inconsistency in the pull potential provider counting");
- gmx_fatal(FARGS, "No external provider for external pull potentials have been provided for %d pull coordinates. The first coordinate without provider is number %zu, which expects a module named '%s' to provide the external potential.",
- pull->numUnregisteredExternalPotentials,
- c + 1, pull->coord[c].params.externalPotentialProvider);
+ gmx_fatal(FARGS,
+ "No external provider for external pull potentials have been provided for %d "
+ "pull coordinates. The first coordinate without provider is number %zu, which "
+ "expects a module named '%s' to provide the external potential.",
+ pull->numUnregisteredExternalPotentials, c + 1,
+ pull->coord[c].params.externalPotentialProvider);
}
}
* potential energy is added either to the pull term or to a term
* specific to the external module.
*/
-void apply_external_pull_coord_force(struct pull_t *pull,
+void apply_external_pull_coord_force(struct pull_t* pull,
int coord_index,
double coord_force,
- const t_mdatoms *mdatoms,
- gmx::ForceWithVirial *forceWithVirial)
+ const t_mdatoms* mdatoms,
+ gmx::ForceWithVirial* forceWithVirial)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
- GMX_ASSERT(coord_index >= 0 && coord_index < static_cast<int>(pull->coord.size()), "apply_external_pull_coord_force called with coord_index out of range");
+ GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
+ "apply_external_pull_coord_force called with coord_index out of range");
if (pull->comm.bParticipate)
{
pcrd = &pull->coord[coord_index];
- GMX_RELEASE_ASSERT(pcrd->params.eType == epullEXTERNAL, "The pull force can only be set externally on pull coordinates of external type");
+ GMX_RELEASE_ASSERT(
+ pcrd->params.eType == epullEXTERNAL,
+ "The pull force can only be set externally on pull coordinates of external type");
- GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered, "apply_external_pull_coord_force called for an unregistered pull coordinate");
+ GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
+ "apply_external_pull_coord_force called for an unregistered pull coordinate");
/* Set the force */
pcrd->scalarForce = coord_force;
/* Calculate the pull potential and scalar force for a pull coordinate.
* Returns the vector forces for the pull coordinate.
*/
-static PullCoordVectorForces
-do_pull_pot_coord(struct pull_t *pull, int coord_ind, t_pbc *pbc,
- double t, real lambda,
- real *V, tensor vir, real *dVdl)
+static PullCoordVectorForces 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 = pull->coord[coord_ind];
+ pull_coord_work_t& pcrd = pull->coord[coord_ind];
assert(pcrd.params.eType != epullCONSTRAINT);
return pullCoordForces;
}
-real pull_potential(struct pull_t *pull, const t_mdatoms *md, t_pbc *pbc,
- const t_commrec *cr, double t, real lambda,
- const rvec *x, gmx::ForceWithVirial *force, real *dvdlambda)
+real pull_potential(struct pull_t* pull,
+ const t_mdatoms* md,
+ t_pbc* pbc,
+ const t_commrec* cr,
+ double t,
+ real lambda,
+ const rvec* x,
+ gmx::ForceWithVirial* force,
+ real* dvdlambda)
{
real V = 0;
pull_calc_coms(cr, pull, md, pbc, t, x, nullptr);
- rvec *f = as_rvec_array(force->force_.data());
- matrix virial = { { 0 } };
- const bool computeVirial = (force->computeVirial_ && MASTER(cr));
+ rvec* f = as_rvec_array(force->force_.data());
+ matrix virial = { { 0 } };
+ const bool computeVirial = (force->computeVirial_ && MASTER(cr));
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
pcrd = &pull->coord[c];
- /* For external potential the force is assumed to be given by an external module by a call to
- apply_pull_coord_external_force */
+ /* For external potential the force is assumed to be given by an external module by a
+ call to apply_pull_coord_external_force */
if (pcrd->params.eType == epullCONSTRAINT || pcrd->params.eType == epullEXTERNAL)
{
continue;
}
- PullCoordVectorForces pullCoordForces =
- do_pull_pot_coord(pull, c, pbc, t, lambda,
- &V,
- computeVirial ? virial : nullptr,
- &dVdl);
+ PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
+ pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
/* Distribute the force over the atoms in the pulled groups */
apply_forces_coord(pull, c, pullCoordForces, md, f);
}
}
- GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0, "Too few or too many external pull potentials have been applied the previous step");
+ GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
+ "Too few or too many external pull potentials have been applied the previous step");
/* All external pull potentials still need to be applied */
- pull->numExternalPotentialsStillToBeAppliedThisStep =
- pull->numCoordinatesWithExternalPotential;
+ pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
return (MASTER(cr) ? V : 0.0);
}
-void pull_constraint(struct pull_t *pull, const t_mdatoms *md, t_pbc *pbc,
- const t_commrec *cr, double dt, double t,
- rvec *x, rvec *xp, rvec *v, tensor vir)
+void pull_constraint(struct pull_t* pull,
+ const t_mdatoms* md,
+ t_pbc* pbc,
+ const t_commrec* cr,
+ double dt,
+ double t,
+ rvec* x,
+ rvec* xp,
+ rvec* v,
+ tensor vir)
{
assert(pull != nullptr);
}
}
-void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull)
+void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
{
- gmx_domdec_t *dd;
- pull_comm_t *comm;
- gmx_bool bMustParticipate;
+ gmx_domdec_t* dd;
+ pull_comm_t* comm;
+ gmx_bool bMustParticipate;
dd = cr->dd;
*/
bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
- for (pull_group_work_t &group : pull->group)
+ for (pull_group_work_t& group : pull->group)
{
if (!group.globalWeights.empty())
{
}
}
- GMX_ASSERT(bMustParticipate || dd != nullptr, "Either all ranks (including this rank) participate, or we use DD and need to have access to dd here");
+ GMX_ASSERT(bMustParticipate || dd != nullptr,
+ "Either all ranks (including this rank) participate, or we use DD and need to "
+ "have access to dd here");
/* We should participate if we have pull or pbc atoms */
- if (!bMustParticipate &&
- (group.atomSet.numAtomsLocal() > 0 || (group.epgrppbc == epgrppbcREFAT &&
- group.pbcAtomSet->numAtomsLocal() > 0)))
+ if (!bMustParticipate
+ && (group.atomSet.numAtomsLocal() > 0
+ || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
{
bMustParticipate = TRUE;
}
* if they needed to participate up to 20 decompositions ago.
* This avoids frequent rebuilds due to atoms jumping back and forth.
*/
- const int64_t history_count = 20;
- gmx_bool bWillParticipate;
- int count[2];
+ const int64_t history_count = 20;
+ gmx_bool bWillParticipate;
+ int count[2];
/* Increase the decomposition counter for the current call */
comm->setup_count++;
}
bWillParticipate =
- bMustParticipate ||
- (comm->bParticipate &&
- comm->must_count >= comm->setup_count - history_count);
+ bMustParticipate
+ || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
if (debug && dd != nullptr)
{
- fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
- dd->rank, gmx::boolToString(bMustParticipate), gmx::boolToString(bWillParticipate));
+ fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n", dd->rank,
+ gmx::boolToString(bMustParticipate), gmx::boolToString(bWillParticipate));
}
if (bWillParticipate)
/* If we are missing ranks or if we have 20% more ranks than needed
* we make a new sub-communicator.
*/
- if (count[1] > 0 || 6*count[0] < 5*comm->nparticipate)
+ if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
{
if (debug)
{
- fprintf(debug, "Creating new pull subcommunicator of size %d\n",
- count[0]);
+ fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
}
#if GMX_MPI
* to avoid this splitting as much as possible.
*/
assert(dd != nullptr);
- MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank,
- &comm->mpi_comm_com);
+ MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
#endif
comm->bParticipate = bWillParticipate;
/* When we use the previous COM for PBC, we need to broadcast
* the previous COM to ranks that have joined the communicator.
*/
- for (pull_group_work_t &group : pull->group)
+ for (pull_group_work_t& group : pull->group)
{
if (group.epgrppbc == epgrppbcPREVSTEPCOM)
{
GMX_ASSERT(comm->bParticipate || !MASTER(cr),
- "The master rank has to participate, as it should pass an up to date prev. COM "
+ "The master rank has to participate, as it should pass an up to "
+ "date prev. COM "
"to bcast here as well as to e.g. checkpointing");
gmx_bcast(sizeof(dvec), group.x_prev_step, cr);
pull->bSetPBCatoms = TRUE;
}
-static void init_pull_group_index(FILE *fplog, const t_commrec *cr,
- int g, pull_group_work_t *pg,
- gmx_bool bConstraint, const ivec pulldim_con,
- const gmx_mtop_t *mtop,
- const t_inputrec *ir, real lambda)
+static void init_pull_group_index(FILE* fplog,
+ const t_commrec* cr,
+ int g,
+ pull_group_work_t* pg,
+ gmx_bool bConstraint,
+ const ivec pulldim_con,
+ const gmx_mtop_t* mtop,
+ const t_inputrec* ir,
+ real lambda)
{
/* With EM and BD there are no masses in the integrator.
* But we still want to have the correct mass-weighted COMs.
* So we store the real masses in the weights.
*/
- const bool setWeights = (pg->params.nweight > 0 ||
- EI_ENERGY_MINIMIZATION(ir->eI) ||
- ir->eI == eiBD);
+ const bool setWeights = (pg->params.nweight > 0 || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
/* In parallel, store we need to extract localWeights from weights at DD time */
- std::vector<real> &weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
+ std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
- const gmx_groups_t *groups = &mtop->groups;
+ const SimulationGroups& groups = mtop->groups;
/* Count frozen dimensions and (weighted) mass */
int nfrozen = 0;
{
for (int d = 0; d < DIM; d++)
{
- if (pulldim_con[d] == 1 &&
- ir->opts.nFreeze[getGroupType(groups, egcFREEZE, ii)][d])
+ if (pulldim_con[d] == 1
+ && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
{
nfrozen++;
}
}
}
- const t_atom &atom = mtopGetAtomParameters(mtop, ii, &molb);
+ const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
real m;
if (ir->efep == efepNO)
{
}
else
{
- m = (1 - lambda)*atom.m + lambda*atom.mB;
+ m = (1 - lambda) * atom.m + lambda * atom.mB;
}
real w;
if (pg->params.nweight > 0)
if (EI_ENERGY_MINIMIZATION(ir->eI))
{
/* Move the mass to the weight */
- w *= m;
- m = 1;
+ w *= m;
+ m = 1;
}
else if (ir->eI == eiBD)
{
real mbd;
- if (ir->bd_fric != 0.0f)
+ if (ir->bd_fric != 0.0F)
{
- mbd = ir->bd_fric*ir->delta_t;
+ mbd = ir->bd_fric * ir->delta_t;
}
else
{
- if (groups->grpnr[egcTC] == nullptr)
+ if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
{
- mbd = ir->delta_t/ir->opts.tau_t[0];
+ mbd = ir->delta_t / ir->opts.tau_t[0];
}
else
{
- mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
+ mbd = ir->delta_t
+ / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
}
}
- w *= m/mbd;
- m = mbd;
+ w *= m / mbd;
+ m = mbd;
}
if (setWeights)
{
weights.push_back(w);
}
- tmass += m;
- wmass += m*w;
- wwmass += m*w*w;
+ tmass += m;
+ wmass += m * w;
+ wwmass += m * w * w;
}
if (wmass == 0)
}
if (fplog)
{
- fprintf(fplog,
- "Pull group %d: %5d atoms, mass %9.3f",
- g, pg->params.nat, tmass);
- if (pg->params.weight ||
- EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
+ fprintf(fplog, "Pull group %d: %5d atoms, mass %9.3f", g, pg->params.nat, tmass);
+ if (pg->params.weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
{
- fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
+ fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
}
if (pg->epgrppbc == epgrppbcCOS)
{
if (nfrozen == 0)
{
/* A value != 0 signals not frozen, it is updated later */
- pg->invtm = -1.0;
+ pg->invtm = -1.0;
}
else
{
int ndim = 0;
for (int d = 0; d < DIM; d++)
{
- ndim += pulldim_con[d]*pg->params.nat;
+ ndim += pulldim_con[d] * pg->params.nat;
}
if (fplog && nfrozen > 0 && nfrozen < ndim)
{
}
}
-struct pull_t *
-init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
- const gmx_mtop_t *mtop, const t_commrec *cr, gmx::LocalAtomSetManager *atomSets,
- real lambda)
+struct pull_t* init_pull(FILE* fplog,
+ const pull_params_t* pull_params,
+ const t_inputrec* ir,
+ const gmx_mtop_t* mtop,
+ const t_commrec* cr,
+ gmx::LocalAtomSetManager* atomSets,
+ real lambda)
{
- struct pull_t *pull;
- pull_comm_t *comm;
+ struct pull_t* pull;
+ pull_comm_t* comm;
pull = new pull_t;
/* Copy the pull parameters */
- pull->params = *pull_params;
+ pull->params = *pull_params;
/* Avoid pointer copies */
pull->params.group = nullptr;
pull->params.coord = nullptr;
for (int i = 0; i < pull_params->ngroup; ++i)
{
- pull->group.emplace_back(pull_params->group[i], atomSets->add({pull_params->group[i].ind, pull_params->group[i].ind+pull_params->group[i].nat}),
+ pull->group.emplace_back(pull_params->group[i],
+ atomSets->add({ pull_params->group[i].ind,
+ pull_params->group[i].ind + pull_params->group[i].nat }),
pull_params->bSetPbcRefToPrevStepCOM);
}
if (cr != nullptr && DOMAINDECOMP(cr))
{
/* Set up the global to local atom mapping for PBC atoms */
- for (pull_group_work_t &group : pull->group)
+ for (pull_group_work_t& group : pull->group)
{
if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
{
/* pbcAtomSet consists of a single atom */
- group.pbcAtomSet = gmx::compat::make_unique<gmx::LocalAtomSet>(atomSets->add({&group.params.pbcatom, &group.params.pbcatom + 1}));
+ group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
+ atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
}
}
}
pull->bXOutAverage = pull_params->bXOutAverage;
pull->bFOutAverage = pull_params->bFOutAverage;
- GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms");
+ GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0,
+ "pull group 0 is an absolute reference group and should not contain atoms");
pull->numCoordinatesWithExternalPotential = 0;
/* Construct a pull coordinate, copying all coordinate parameters */
pull->coord.emplace_back(pull_params->coord[c]);
- pull_coord_work_t *pcrd = &pull->coord.back();
+ pull_coord_work_t* pcrd = &pull->coord.back();
switch (pcrd->params.eGeom)
{
case epullgDIST:
- case epullgDIRRELATIVE: /* Direction vector is determined at each step */
+ case epullgDIRRELATIVE: /* Direction vector is determined at each step */
case epullgANGLE:
- case epullgDIHEDRAL:
- break;
+ case epullgDIHEDRAL: break;
case epullgDIR:
case epullgDIRPBC:
case epullgCYL:
* 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.",
+ 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)
{
/* Check restrictions of the constraint pull code */
- if (pcrd->params.eGeom == epullgCYL ||
- pcrd->params.eGeom == epullgDIRRELATIVE ||
- pcrd->params.eGeom == epullgANGLE ||
- pcrd->params.eGeom == epullgDIHEDRAL ||
- pcrd->params.eGeom == epullgANGLEAXIS)
+ if (pcrd->params.eGeom == epullgCYL || pcrd->params.eGeom == epullgDIRRELATIVE
+ || pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
+ || pcrd->params.eGeom == epullgANGLEAXIS)
{
- gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
- epull_names[pcrd->params.eType],
- epullg_names[pcrd->params.eGeom],
+ gmx_fatal(FARGS,
+ "Pulling of type %s can not be combined with geometry %s. Consider using "
+ "pull type %s.",
+ epull_names[pcrd->params.eType], epullg_names[pcrd->params.eGeom],
epull_names[epullUMBRELLA]);
}
{
pull->bCylinder = TRUE;
}
- else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL || pcrd->params.eGeom == epullgANGLEAXIS)
+ else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
+ || pcrd->params.eGeom == epullgANGLEAXIS)
{
pull->bAngle = TRUE;
}
for (int i = 0; i < pcrd->params.ngroup; i++)
{
int groupIndex = pcrd->params.group[i];
- if (groupIndex > 0 &&
- !(pcrd->params.eGeom == epullgCYL && i == 0))
+ if (groupIndex > 0 && !(pcrd->params.eGeom == epullgCYL && i == 0))
{
pull->group[groupIndex].needToCalcCom = true;
}
/* Initialize the constant reference value */
if (pcrd->params.eType != epullEXTERNAL)
{
- low_set_pull_coord_reference_value(pcrd, c, pcrd->params.init*pull_conversion_factor_userinput2internal(&pcrd->params));
+ low_set_pull_coord_reference_value(
+ pcrd, c,
+ pcrd->params.init * pull_conversion_factor_userinput2internal(&pcrd->params));
}
else
{
if (pcrd->params.eType == epullEXTERNAL)
{
- GMX_RELEASE_ASSERT(pcrd->params.rate == 0, "With an external potential, a pull coordinate should have rate = 0");
+ GMX_RELEASE_ASSERT(
+ pcrd->params.rate == 0,
+ "With an external potential, a pull coordinate should have rate = 0");
/* This external potential needs to be registered later */
pull->numCoordinatesWithExternalPotential++;
pcrd->bExternalPotentialProviderHasBeenRegistered = false;
}
- pull->numUnregisteredExternalPotentials =
- pull->numCoordinatesWithExternalPotential;
+ pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
pull->ePBC = ir->ePBC;
switch (pull->ePBC)
{
case epbcNONE: pull->npbcdim = 0; break;
- case epbcXY: pull->npbcdim = 2; break;
- default: pull->npbcdim = 3; break;
+ case epbcXY: pull->npbcdim = 2; break;
+ default: pull->npbcdim = 3; break;
}
if (fplog)
gmx_bool bAbs, bCos;
bAbs = FALSE;
- for (const pull_coord_work_t &coord : pull->coord)
+ for (const pull_coord_work_t& coord : pull->coord)
{
- if (pull->group[coord.params.group[0]].params.nat == 0 ||
- pull->group[coord.params.group[1]].params.nat == 0)
+ if (pull->group[coord.params.group[0]].params.nat == 0
+ || pull->group[coord.params.group[1]].params.nat == 0)
{
bAbs = TRUE;
}
}
// Don't include the reference group 0 in output, so we report ngroup-1
int numRealGroups = pull->group.size() - 1;
- GMX_RELEASE_ASSERT(numRealGroups > 0, "The reference absolute position pull group should always be present");
- fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n",
- pull->coord.size(), pull->coord.size() == 1 ? "" : "s",
- numRealGroups, numRealGroups == 1 ? "" : "s");
+ GMX_RELEASE_ASSERT(numRealGroups > 0,
+ "The reference absolute position pull group should always be present");
+ fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n", pull->coord.size(),
+ pull->coord.size() == 1 ? "" : "s", numRealGroups, numRealGroups == 1 ? "" : "s");
if (bAbs)
{
fprintf(fplog, "with an absolute reference\n");
// Don't include the reference group 0 in loop
for (size_t g = 1; g < pull->group.size(); g++)
{
- if (pull->group[g].params.nat > 1 &&
- pull->group[g].params.pbcatom < 0)
+ if (pull->group[g].params.nat > 1 && pull->group[g].params.pbcatom < 0)
{
/* We are using cosine weighting */
fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
}
}
- pull->bRefAt = FALSE;
- pull->cosdim = -1;
+ pull->bRefAt = FALSE;
+ pull->cosdim = -1;
for (size_t g = 0; g < pull->group.size(); g++)
{
- pull_group_work_t *pgrp;
+ pull_group_work_t* pgrp;
- pgrp = &pull->group[g];
+ pgrp = &pull->group[g];
if (pgrp->params.nat > 0)
{
/* There is an issue when a group is used in multiple coordinates
bConstraint = FALSE;
clear_ivec(pulldim);
clear_ivec(pulldim_con);
- for (const pull_coord_work_t &coord : pull->coord)
+ for (const pull_coord_work_t& coord : pull->coord)
{
gmx_bool bGroupUsed = FALSE;
for (int gi = 0; gi < coord.params.ngroup; gi++)
*/
switch (pgrp->epgrppbc)
{
- case epgrppbcREFAT:
- pull->bRefAt = TRUE;
- break;
+ case epgrppbcREFAT: pull->bRefAt = TRUE; break;
case epgrppbcCOS:
if (pgrp->params.weight != nullptr)
{
- gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
+ gmx_fatal(FARGS,
+ "Pull groups can not have relative weights and cosine weighting "
+ "at same time");
}
for (int m = 0; m < DIM; m++)
{
{
if (pull->cosdim >= 0 && pull->cosdim != m)
{
- gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
+ gmx_fatal(FARGS,
+ "Can only use cosine weighting with pulling in one "
+ "dimension (use mdp option pull-coord?-dim)");
}
pull->cosdim = m;
}
}
break;
- case epgrppbcNONE:
- break;
+ case epgrppbcNONE: break;
}
/* Set the indices */
- init_pull_group_index(fplog, cr, g, pgrp,
- bConstraint, pulldim_con,
- mtop, ir, lambda);
+ init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
}
else
{
/* If we use cylinder coordinates, do some initialising for them */
if (pull->bCylinder)
{
- for (const pull_coord_work_t &coord : pull->coord)
+ for (const pull_coord_work_t& coord : pull->coord)
{
if (coord.params.eGeom == epullgCYL)
{
if (pull->group[coord.params.group[0]].params.nat == 0)
{
- gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
+ gmx_fatal(FARGS,
+ "A cylinder pull group is not supported when using absolute "
+ "reference!\n");
}
}
- const auto &referenceGroup = pull->group[coord.params.group[0]];
- pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet, pull->params.bSetPbcRefToPrevStepCOM);
+ const auto& referenceGroup = pull->group[coord.params.group[0]];
+ pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet,
+ pull->params.bSetPbcRefToPrevStepCOM);
}
}
* when we have an external pull potential, since then the external
* potential provider expects each rank to have the coordinate.
*/
- comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) ||
- cr->dd->nnodes <= 32 ||
- pull->numCoordinatesWithExternalPotential > 0 ||
- getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
+ comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
+ || pull->numCoordinatesWithExternalPotential > 0
+ || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
/* This sub-commicator is not used with comm->bParticipateAll,
* so we can always initialize it to NULL.
*/
- comm->mpi_comm_com = MPI_COMM_NULL;
- comm->nparticipate = 0;
- comm->isMasterRank = (cr == nullptr || MASTER(cr));
+ comm->mpi_comm_com = MPI_COMM_NULL;
+ comm->nparticipate = 0;
+ comm->isMasterRank = (cr == nullptr || MASTER(cr));
#else
/* No MPI: 1 rank: all ranks pull */
comm->bParticipateAll = TRUE;
comm->isMasterRank = true;
#endif
- comm->bParticipate = comm->bParticipateAll;
- comm->setup_count = 0;
- comm->must_count = 0;
+ comm->bParticipate = comm->bParticipateAll;
+ comm->setup_count = 0;
+ comm->must_count = 0;
if (!comm->bParticipateAll && fplog != nullptr)
{
}
comm->pbcAtomBuffer.resize(pull->group.size());
- comm->comBuffer.resize(pull->group.size()*c_comBufferStride);
+ comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
if (pull->bCylinder)
{
- comm->cylinderBuffer.resize(pull->coord.size()*c_cylinderBufferStride);
+ comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
}
/* We still need to initialize the PBC reference coordinates */
return pull;
}
-static void destroy_pull(struct pull_t *pull)
+static void destroy_pull(struct pull_t* pull)
{
#if GMX_MPI
if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
delete pull;
}
-void preparePrevStepPullCom(const t_inputrec *ir, const t_mdatoms *md, t_state *state, const t_state *state_global, const t_commrec *cr, bool startingFromCheckpoint)
+void preparePrevStepPullCom(const t_inputrec* ir,
+ pull_t* pull_work,
+ const t_mdatoms* md,
+ t_state* state,
+ const t_state* state_global,
+ const t_commrec* cr,
+ bool startingFromCheckpoint)
{
if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
{
return;
}
- allocStatePrevStepPullCom(state, ir->pull_work);
+ allocStatePrevStepPullCom(state, pull_work);
if (startingFromCheckpoint)
{
if (MASTER(cr))
/* Only the master rank has the checkpointed COM from the previous step */
gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(), &state->pull_com_prev_step[0], cr);
}
- setPrevStepPullComFromState(ir->pull_work, state);
+ setPrevStepPullComFromState(pull_work, state);
}
else
{
t_pbc pbc;
set_pbc(&pbc, ir->ePBC, state->box);
- initPullComFromPrevStep(cr, ir->pull_work, md, &pbc, state->x.rvec_array());
- updatePrevStepPullCom(ir->pull_work, state);
+ initPullComFromPrevStep(cr, pull_work, md, &pbc, state->x.rvec_array());
+ updatePrevStepPullCom(pull_work, state);
}
}
-void finish_pull(struct pull_t *pull)
+void finish_pull(struct pull_t* pull)
{
check_external_potential_registration(pull);
destroy_pull(pull);
}
-gmx_bool pull_have_potential(const struct pull_t *pull)
+gmx_bool pull_have_potential(const struct pull_t* pull)
{
return pull->bPotential;
}
-gmx_bool pull_have_constraint(const struct pull_t *pull)
+gmx_bool pull_have_constraint(const struct pull_t* pull)
{
return pull->bConstraint;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff