2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
53 #include "gromacs/commandline/filenm.h"
54 #include "gromacs/domdec/domdec_struct.h"
55 #include "gromacs/domdec/localatomset.h"
56 #include "gromacs/domdec/localatomsetmanager.h"
57 #include "gromacs/fileio/gmxfio.h"
58 #include "gromacs/gmxlib/network.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/utilities.h"
62 #include "gromacs/math/vec.h"
63 #include "gromacs/math/vectypes.h"
64 #include "gromacs/mdlib/gmx_omp_nthreads.h"
65 #include "gromacs/mdtypes/commrec.h"
66 #include "gromacs/mdtypes/forceoutput.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/md_enums.h"
69 #include "gromacs/mdtypes/state.h"
70 #include "gromacs/pbcutil/pbc.h"
71 #include "gromacs/topology/mtop_lookup.h"
72 #include "gromacs/topology/topology.h"
73 #include "gromacs/utility/arrayref.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/enumerationhelpers.h"
76 #include "gromacs/utility/exceptions.h"
77 #include "gromacs/utility/fatalerror.h"
78 #include "gromacs/utility/futil.h"
79 #include "gromacs/utility/gmxassert.h"
80 #include "gromacs/utility/pleasecite.h"
81 #include "gromacs/utility/real.h"
82 #include "gromacs/utility/smalloc.h"
83 #include "gromacs/utility/strconvert.h"
84 #include "gromacs/utility/stringutil.h"
86 #include "pull_internal.h"
90 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
96 /*! \brief Tells whether the pull geometry is an angle type */
97 constexpr gmx::EnumerationArray<PullGroupGeometry, bool> sc_isAngleType = {
98 { false, false, false, false, false, true, true, true }
101 static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
103 if (params.ind.size() <= 1)
105 /* no PBC required */
108 else if (params.pbcatom >= 0)
110 if (setPbcRefToPrevStepCOM)
112 return epgrppbcPREVSTEPCOM;
116 return epgrppbcREFAT;
125 /* NOTE: The params initialization currently copies pointers.
126 * So the lifetime of the source, currently always inputrec,
127 * should not end before that of this object.
128 * This will be fixed when the pointers are replacted by std::vector.
130 pull_group_work_t::pull_group_work_t(const t_pull_group& params,
131 gmx::LocalAtomSet atomSet,
132 bool bSetPbcRefToPrevStepCOM,
135 epgrppbc(groupPbcFromParams(params, bSetPbcRefToPrevStepCOM)),
136 maxNumThreads(maxNumThreads),
137 needToCalcCom(false),
147 static bool pull_coordinate_is_directional(const t_pull_coord* pcrd)
149 return (pcrd->eGeom == PullGroupGeometry::Direction || pcrd->eGeom == PullGroupGeometry::DirectionPBC
150 || pcrd->eGeom == PullGroupGeometry::DirectionRelative
151 || pcrd->eGeom == PullGroupGeometry::Cylinder);
154 const char* pull_coordinate_units(const t_pull_coord& pcrd)
156 return sc_isAngleType[pcrd.eGeom] ? "deg" : "nm";
159 double pull_conversion_factor_userinput2internal(const t_pull_coord& pcrd)
161 if (sc_isAngleType[pcrd.eGeom])
163 return gmx::c_deg2Rad;
171 double pull_conversion_factor_internal2userinput(const t_pull_coord& pcrd)
173 if (sc_isAngleType[pcrd.eGeom])
175 return gmx::c_rad2Deg;
183 /* Apply forces in a mass weighted fashion for part of the pull group */
184 static void apply_forces_grp_part(const pull_group_work_t& pgrp,
187 ArrayRef<const real> masses,
192 const double inv_wm = pgrp.mwscale;
194 auto localAtomIndices = pgrp.atomSet.localIndex();
195 for (int i = ind_start; i < ind_end; i++)
197 int ii = localAtomIndices[i];
198 double wmass = masses[ii];
199 if (!pgrp.localWeights.empty())
201 wmass *= pgrp.localWeights[i];
204 for (int d = 0; d < DIM; d++)
206 f[ii][d] += sign * wmass * f_pull[d] * inv_wm;
211 /* Apply forces in a mass weighted fashion */
212 static void apply_forces_grp(const pull_group_work_t& pgrp,
213 ArrayRef<const real> masses,
218 auto localAtomIndices = pgrp.atomSet.localIndex();
220 if (pgrp.params.ind.size() == 1 && pgrp.atomSet.numAtomsLocal() == 1)
222 /* Only one atom and our rank has this atom: we can skip
223 * the mass weighting, which means that this code also works
224 * for mass=0, e.g. with a virtual site.
226 for (int d = 0; d < DIM; d++)
228 f[localAtomIndices[0]][d] += sign * f_pull[d];
233 const int numThreads = pgrp.numThreads();
236 apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), masses, f_pull, sign, f);
240 #pragma omp parallel for num_threads(numThreads) schedule(static)
241 for (int th = 0; th < numThreads; th++)
243 int ind_start = (localAtomIndices.size() * (th + 0)) / numThreads;
244 int ind_end = (localAtomIndices.size() * (th + 1)) / numThreads;
245 apply_forces_grp_part(pgrp, ind_start, ind_end, masses, f_pull, sign, f);
251 /* Apply forces in a mass weighted fashion to a cylinder group */
252 static void apply_forces_cyl_grp(const pull_group_work_t& pgrp,
253 const double dv_corr,
254 ArrayRef<const real> masses,
260 const double inv_wm = pgrp.mwscale;
262 auto localAtomIndices = pgrp.atomSet.localIndex();
264 /* The cylinder group is always a slab in the system, thus large.
265 * Therefore we always thread-parallelize this group.
267 const int numAtomsLocal = localAtomIndices.size();
268 const int gmx_unused numThreads = pgrp.numThreads();
269 #pragma omp parallel for num_threads(numThreads) schedule(static)
270 for (int i = 0; i < numAtomsLocal; i++)
272 double weight = pgrp.localWeights[i];
277 int ii = localAtomIndices[i];
278 double mass = masses[ii];
279 /* The stored axial distance from the cylinder center (dv) needs
280 * to be corrected for an offset (dv_corr), which was unknown when
283 double dv_com = pgrp.dv[i] + dv_corr;
285 /* Here we not only add the pull force working along vec (f_pull),
286 * but also a radial component, due to the dependence of the weights
287 * on the radial distance.
289 for (int m = 0; m < DIM; m++)
291 f[ii][m] += sign * inv_wm * (mass * weight * f_pull[m] + pgrp.mdw[i][m] * dv_com * f_scal);
296 /* Apply torque forces in a mass weighted fashion to the groups that define
297 * the pull vector direction for pull coordinate pcrd.
299 static void apply_forces_vec_torque(const pull_coord_work_t& pcrd,
300 ArrayRef<const pull_group_work_t> pullGroups,
301 ArrayRef<const real> masses,
304 const PullCoordSpatialData& spatialData = pcrd.spatialData;
306 /* The component inpr along the pull vector is accounted for in the usual
307 * way. Here we account for the component perpendicular to vec.
310 for (int m = 0; m < DIM; m++)
312 inpr += spatialData.dr01[m] * spatialData.vec[m];
314 /* The torque force works along the component of the distance vector
315 * of between the two "usual" pull groups that is perpendicular to
316 * the pull vector. The magnitude of this force is the "usual" scale force
317 * multiplied with the ratio of the distance between the two "usual" pull
318 * groups and the distance between the two groups that define the vector.
321 for (int m = 0; m < DIM; m++)
323 f_perp[m] = (spatialData.dr01[m] - inpr * spatialData.vec[m]) / spatialData.vec_len * pcrd.scalarForce;
326 /* Apply the force to the groups defining the vector using opposite signs */
327 apply_forces_grp(pullGroups[pcrd.params.group[2]], masses, f_perp, -1, f);
328 apply_forces_grp(pullGroups[pcrd.params.group[3]], masses, f_perp, 1, f);
331 /* Apply forces in a mass weighted fashion */
332 static void apply_forces_coord(const pull_coord_work_t& pcrd,
333 ArrayRef<const pull_group_work_t> pullGroups,
334 const PullCoordVectorForces& forces,
335 ArrayRef<const real> masses,
338 /* Here it would be more efficient to use one large thread-parallel
339 * region instead of potential parallel regions within apply_forces_grp.
340 * But there could be overlap between pull groups and this would lead
344 if (pcrd.params.eGeom == PullGroupGeometry::Cylinder)
346 apply_forces_cyl_grp(
347 *pcrd.dynamicGroup0, pcrd.spatialData.cyl_dev, masses, forces.force01, pcrd.scalarForce, -1, f);
349 /* Sum the force along the vector and the radial force */
351 for (int m = 0; m < DIM; m++)
353 f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
355 apply_forces_grp(pullGroups[pcrd.params.group[1]], masses, f_tot, 1, f);
359 if (pcrd.params.eGeom == PullGroupGeometry::DirectionRelative)
361 /* We need to apply the torque forces to the pull groups
362 * that define the pull vector.
364 apply_forces_vec_torque(pcrd, pullGroups, masses, f);
367 if (!pullGroups[pcrd.params.group[0]].params.ind.empty())
369 apply_forces_grp(pullGroups[pcrd.params.group[0]], masses, forces.force01, -1, f);
371 apply_forces_grp(pullGroups[pcrd.params.group[1]], masses, forces.force01, 1, f);
373 if (pcrd.params.ngroup >= 4)
375 apply_forces_grp(pullGroups[pcrd.params.group[2]], masses, forces.force23, -1, f);
376 apply_forces_grp(pullGroups[pcrd.params.group[3]], masses, forces.force23, 1, f);
378 if (pcrd.params.ngroup >= 6)
380 apply_forces_grp(pullGroups[pcrd.params.group[4]], masses, forces.force45, -1, f);
381 apply_forces_grp(pullGroups[pcrd.params.group[5]], masses, forces.force45, 1, f);
386 real max_pull_distance2(const pull_coord_work_t* pcrd, const t_pbc* pbc)
388 /* Note that this maximum distance calculation is more complex than
389 * most other cases in GROMACS, since here we have to take care of
390 * distance calculations that don't involve all three dimensions.
391 * For example, we can use distances that are larger than the
392 * box X and Y dimensions for a box that is elongated along Z.
395 real max_d2 = GMX_REAL_MAX;
397 if (pull_coordinate_is_directional(&pcrd->params))
399 /* Directional pulling along along direction pcrd->vec.
400 * Calculating the exact maximum distance is complex and bug-prone.
401 * So we take a safe approach by not allowing distances that
402 * are larger than half the distance between unit cell faces
403 * along dimensions involved in pcrd->vec.
405 for (int m = 0; m < DIM; m++)
407 if (m < pbc->ndim_ePBC && pcrd->spatialData.vec[m] != 0)
409 real imageDistance2 = gmx::square(pbc->box[m][m]);
410 for (int d = m + 1; d < DIM; d++)
412 imageDistance2 -= gmx::square(pbc->box[d][m]);
414 max_d2 = std::min(max_d2, imageDistance2);
420 /* Distance pulling along dimensions with pcrd->params.dim[d]==1.
421 * We use half the minimum box vector length of the dimensions involved.
422 * This is correct for all cases, except for corner cases with
423 * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
424 * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
425 * in such corner cases the user could get correct results,
426 * depending on the details of the setup, we avoid further
427 * code complications.
429 for (int m = 0; m < DIM; m++)
431 if (m < pbc->ndim_ePBC && pcrd->params.dim[m] != 0)
433 real imageDistance2 = gmx::square(pbc->box[m][m]);
434 for (int d = 0; d < m; d++)
436 if (pcrd->params.dim[d] != 0)
438 imageDistance2 += gmx::square(pbc->box[m][d]);
441 max_d2 = std::min(max_d2, imageDistance2);
446 return 0.25 * max_d2;
449 /* This function returns the distance based on coordinates xg and xref.
450 * Note that the pull coordinate struct pcrd is not modified.
452 * \param[in] pull The pull struct
453 * \param[in] pcrd The pull coordinate to compute a distance for
454 * \param[in] pbc The periodic boundary conditions
455 * \param[in] xg The coordinate of group 1
456 * \param[in] xref The coordinate of group 0
457 * \param[in] groupIndex0 The index of group 0 in the pcrd->params.group
458 * \param[in] groupIndex1 The index of group 1 in the pcrd->params.group
459 * \param[in] max_dist2 The maximum distance squared
460 * \param[out] dr The distance vector
462 static void low_get_pull_coord_dr(const pull_t& pull,
463 const pull_coord_work_t* pcrd,
467 const int groupIndex0,
468 const int groupIndex1,
469 const double max_dist2,
472 const pull_group_work_t& pgrp0 = pull.group[pcrd->params.group[0]];
474 // Group coordinate, to be updated with the reference position
477 /* Only the first group can be an absolute reference, in that case nat=0 */
478 if (pgrp0.params.ind.empty())
480 for (int m = 0; m < DIM; m++)
482 xrefr[m] = pcrd->params.origin[m];
487 copy_dvec(xref, xrefr);
490 dvec dref = { 0, 0, 0 };
491 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
493 for (int m = 0; m < DIM; m++)
495 dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
497 /* Add the reference position, so we use the correct periodic image */
498 dvec_inc(xrefr, dref);
501 pbc_dx_d(pbc, xg, xrefr, dr);
503 bool directional = pull_coordinate_is_directional(&pcrd->params);
505 for (int m = 0; m < DIM; m++)
507 dr[m] *= pcrd->params.dim[m];
508 if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
510 dr2 += dr[m] * dr[m];
513 /* Check if we are close to switching to another periodic image */
514 if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
516 /* Note that technically there is no issue with switching periodic
517 * image, as pbc_dx_d returns the distance to the closest periodic
518 * image. However in all cases where periodic image switches occur,
519 * the pull results will be useless in practice.
522 "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
523 "box size (%f).\n%s",
524 pcrd->params.group[groupIndex0],
525 pcrd->params.group[groupIndex1],
527 sqrt(0.98 * 0.98 * max_dist2),
528 pcrd->params.eGeom == PullGroupGeometry::Direction
529 ? "You might want to consider using \"pull-geometry = "
530 "direction-periodic\" instead.\n"
534 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
540 /* This function returns the distance based on the contents of the pull struct.
541 * pull->coord[coord_ind].dr, and potentially vector, are updated.
543 static void get_pull_coord_dr(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc* pbc)
545 const int coord_ind = pcrd->params.coordIndex;
547 PullCoordSpatialData& spatialData = pcrd->spatialData;
550 /* With AWH pulling we allow for periodic pulling with geometry=direction.
551 * TODO: Store a periodicity flag instead of checking for external pull provider.
553 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC
554 || (pcrd->params.eGeom == PullGroupGeometry::Direction
555 && pcrd->params.eType == PullingAlgorithm::External))
561 md2 = static_cast<double>(max_pull_distance2(pcrd, pbc));
564 if (pcrd->params.eGeom == PullGroupGeometry::DirectionRelative)
566 /* We need to determine the pull vector */
570 const pull_group_work_t& pgrp2 = pull.group[pcrd->params.group[2]];
571 const pull_group_work_t& pgrp3 = pull.group[pcrd->params.group[3]];
573 pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
575 for (m = 0; m < DIM; m++)
577 vec[m] *= pcrd->params.dim[m];
579 spatialData.vec_len = dnorm(vec);
580 for (m = 0; m < DIM; m++)
582 spatialData.vec[m] = vec[m] / spatialData.vec_len;
587 "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
594 spatialData.vec[ZZ]);
598 const pull_group_work_t& pgrp0 = pull.group[pcrd->params.group[0]];
599 const pull_group_work_t& pgrp1 = pull.group[pcrd->params.group[1]];
601 low_get_pull_coord_dr(
606 pcrd->params.eGeom == PullGroupGeometry::Cylinder ? pcrd->dynamicGroup0->x : pgrp0.x,
612 if (pcrd->params.ngroup >= 4)
614 const pull_group_work_t& pgrp2 = pull.group[pcrd->params.group[2]];
615 const pull_group_work_t& pgrp3 = pull.group[pcrd->params.group[3]];
617 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3.x, pgrp2.x, 2, 3, md2, spatialData.dr23);
619 if (pcrd->params.ngroup >= 6)
621 const pull_group_work_t& pgrp4 = pull.group[pcrd->params.group[4]];
622 const pull_group_work_t& pgrp5 = pull.group[pcrd->params.group[5]];
624 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5.x, pgrp4.x, 4, 5, md2, spatialData.dr45);
628 /* Modify x so that it is periodic in [-pi, pi)
629 * It is assumed that x is in [-3pi, 3pi) so that x
630 * needs to be shifted by at most one period.
632 static void make_periodic_2pi(double* x)
644 /* This function should always be used to get values for setting value_ref in pull_coord_work_t */
645 static double sanitizePullCoordReferenceValue(const t_pull_coord& pcrdParams, double value_ref)
647 GMX_ASSERT(pcrdParams.eType != PullingAlgorithm::External,
648 "The pull coord reference value should not be used with type external-potential");
650 if (pcrdParams.eGeom == PullGroupGeometry::Distance)
655 "Pull reference distance for coordinate %d (%f) needs to be non-negative",
656 pcrdParams.coordIndex + 1,
660 else if (pcrdParams.eGeom == PullGroupGeometry::Angle || pcrdParams.eGeom == PullGroupGeometry::AngleAxis)
662 if (value_ref < 0 || value_ref > M_PI)
665 "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
666 "interval [0,180] deg",
667 pcrdParams.coordIndex + 1,
668 value_ref * pull_conversion_factor_internal2userinput(pcrdParams));
671 else if (pcrdParams.eGeom == PullGroupGeometry::Dihedral)
673 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
674 make_periodic_2pi(&value_ref);
680 static void updatePullCoordReferenceValue(double* referenceValue, const t_pull_coord& pcrdParams, double t)
682 GMX_ASSERT(referenceValue, "Need a valid reference value object");
684 /* With zero rate the reference value is set initially and doesn't change */
685 if (pcrdParams.rate != 0)
687 const double inputValue = (pcrdParams.init + pcrdParams.rate * t)
688 * pull_conversion_factor_userinput2internal(pcrdParams);
689 *referenceValue = sanitizePullCoordReferenceValue(pcrdParams, inputValue);
693 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
694 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
697 dvec dr32; /* store instead of dr23? */
699 dsvmul(-1, spatialData->dr23, dr32);
700 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
701 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
702 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
704 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
705 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
706 * Note 2: the angle between the plane normal vectors equals pi only when
707 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
708 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
710 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
714 /* Calculates pull->coord[coord_ind].value.
715 * This function also updates pull->coord[coord_ind].dr.
717 static void get_pull_coord_distance(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc* pbc)
719 get_pull_coord_dr(pull, pcrd, pbc);
721 PullCoordSpatialData& spatialData = pcrd->spatialData;
723 switch (pcrd->params.eGeom)
725 case PullGroupGeometry::Distance:
726 /* Pull along the vector between the com's */
727 spatialData.value = dnorm(spatialData.dr01);
729 case PullGroupGeometry::Direction:
730 case PullGroupGeometry::DirectionPBC:
731 case PullGroupGeometry::DirectionRelative:
732 case PullGroupGeometry::Cylinder:
734 spatialData.value = 0;
735 for (int m = 0; m < DIM; m++)
737 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
740 case PullGroupGeometry::Angle:
741 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
743 case PullGroupGeometry::Dihedral:
744 spatialData.value = get_dihedral_angle_coord(&spatialData);
746 case PullGroupGeometry::AngleAxis:
747 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
749 default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
753 /* Returns the deviation from the reference value.
754 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
756 static double get_pull_coord_deviation(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc* pbc, double t)
760 /* Update the reference value before computing the distance,
761 * since it is used in the distance computation with periodic pulling.
763 updatePullCoordReferenceValue(&pcrd->value_ref, pcrd->params, t);
765 get_pull_coord_distance(pull, pcrd, pbc);
767 get_pull_coord_distance(pull, pcrd, pbc);
769 /* Determine the deviation */
770 dev = pcrd->spatialData.value - pcrd->value_ref;
772 /* Check that values are allowed */
773 if (pcrd->params.eGeom == PullGroupGeometry::Distance && pcrd->spatialData.value == 0)
775 /* With no vector we can not determine the direction for the force,
776 * so we set the force to zero.
780 else if (pcrd->params.eGeom == PullGroupGeometry::Dihedral)
782 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
783 Thus, the unwrapped deviation is here in (-2pi, 2pi].
784 After making it periodic, the deviation will be in [-pi, pi). */
785 make_periodic_2pi(&dev);
791 double get_pull_coord_value(pull_t* pull, int coord_ind, const t_pbc* pbc)
793 get_pull_coord_distance(*pull, &pull->coord[coord_ind], pbc);
795 return pull->coord[coord_ind].spatialData.value;
798 void clear_pull_forces(pull_t* pull)
800 /* Zeroing the forces is only required for constraint pulling.
801 * It can happen that multiple constraint steps need to be applied
802 * and therefore the constraint forces need to be accumulated.
804 for (pull_coord_work_t& coord : pull->coord)
806 coord.scalarForce = 0;
810 /* Apply constraint using SHAKE */
811 static void do_constraint(struct pull_t* pull,
821 dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
822 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
823 dvec* rnew; /* current 'new' positions of the groups */
824 double* dr_tot; /* the total update of the coords */
827 double lambda, rm, invdt = 0;
828 gmx_bool bConverged_all, bConverged = FALSE;
829 int niter = 0, ii, j, m, max_iter = 100;
833 snew(r_ij, pull->coord.size());
834 snew(dr_tot, pull->coord.size());
836 snew(rnew, pull->group.size());
838 /* copy the current unconstrained positions for use in iterations. We
839 iterate until rinew[i] and rjnew[j] obey the constraints. Then
840 rinew - pull.x_unc[i] is the correction dr to group i */
841 for (size_t g = 0; g < pull->group.size(); g++)
843 copy_dvec(pull->group[g].xp, rnew[g]);
846 /* Determine the constraint directions from the old positions */
847 for (size_t c = 0; c < pull->coord.size(); c++)
849 pull_coord_work_t* pcrd;
851 pcrd = &pull->coord[c];
853 if (pcrd->params.eType != PullingAlgorithm::Constraint)
858 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
859 * We don't modify dr and value anymore, so these values are also used
862 get_pull_coord_distance(*pull, pcrd, pbc);
864 const PullCoordSpatialData& spatialData = pcrd->spatialData;
868 "Pull coord %zu dr %f %f %f\n",
870 spatialData.dr01[XX],
871 spatialData.dr01[YY],
872 spatialData.dr01[ZZ]);
875 if (pcrd->params.eGeom == PullGroupGeometry::Direction
876 || pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
878 /* Select the component along vec */
880 for (m = 0; m < DIM; m++)
882 a += spatialData.vec[m] * spatialData.dr01[m];
884 for (m = 0; m < DIM; m++)
886 r_ij[c][m] = a * spatialData.vec[m];
891 copy_dvec(spatialData.dr01, r_ij[c]);
894 if (dnorm2(r_ij[c]) == 0)
897 "Distance for pull coordinate %zu is zero with constraint pulling, which is "
903 bConverged_all = FALSE;
904 while (!bConverged_all && niter < max_iter)
906 bConverged_all = TRUE;
908 /* loop over all constraints */
909 for (size_t c = 0; c < pull->coord.size(); c++)
911 pull_coord_work_t* pcrd;
912 pull_group_work_t *pgrp0, *pgrp1;
915 pcrd = &pull->coord[c];
917 if (pcrd->params.eType != PullingAlgorithm::Constraint)
922 updatePullCoordReferenceValue(&pcrd->value_ref, pcrd->params, t);
924 pgrp0 = &pull->group[pcrd->params.group[0]];
925 pgrp1 = &pull->group[pcrd->params.group[1]];
927 /* Get the current difference vector */
928 low_get_pull_coord_dr(
929 *pull, pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
933 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
936 rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
938 switch (pcrd->params.eGeom)
940 case PullGroupGeometry::Distance:
941 if (pcrd->value_ref <= 0)
945 "The pull constraint reference distance for group %zu is <= 0 (%f)",
951 double q, c_a, c_b, c_c;
953 c_a = diprod(r_ij[c], r_ij[c]);
954 c_b = diprod(unc_ij, r_ij[c]) * 2;
955 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
959 q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
964 q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
970 fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b, c_c, lambda);
974 /* The position corrections dr due to the constraints */
975 dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
976 dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
977 dr_tot[c] += -lambda * dnorm(r_ij[c]);
979 case PullGroupGeometry::Direction:
980 case PullGroupGeometry::DirectionPBC:
981 case PullGroupGeometry::Cylinder:
982 /* A 1-dimensional constraint along a vector */
984 for (m = 0; m < DIM; m++)
986 vec[m] = pcrd->spatialData.vec[m];
987 a += unc_ij[m] * vec[m];
989 /* Select only the component along the vector */
990 dsvmul(a, vec, unc_ij);
991 lambda = a - pcrd->value_ref;
994 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
997 /* The position corrections dr due to the constraints */
998 dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
999 dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
1000 dr_tot[c] += -lambda;
1002 default: gmx_incons("Invalid enumeration value for eGeom");
1010 g0 = pcrd->params.group[0];
1011 g1 = pcrd->params.group[1];
1012 low_get_pull_coord_dr(*pull, pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
1013 low_get_pull_coord_dr(*pull, pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
1015 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1024 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
1033 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1043 /* Update the COMs with dr */
1044 dvec_inc(rnew[pcrd->params.group[1]], dr1);
1045 dvec_inc(rnew[pcrd->params.group[0]], dr0);
1048 /* Check if all constraints are fullfilled now */
1049 for (pull_coord_work_t& coord : pull->coord)
1051 if (coord.params.eType != PullingAlgorithm::Constraint)
1056 low_get_pull_coord_dr(
1057 *pull, &coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
1059 switch (coord.params.eGeom)
1061 case PullGroupGeometry::Distance:
1062 bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1064 case PullGroupGeometry::Direction:
1065 case PullGroupGeometry::DirectionPBC:
1066 case PullGroupGeometry::Cylinder:
1067 for (m = 0; m < DIM; m++)
1069 vec[m] = coord.spatialData.vec[m];
1071 inpr = diprod(unc_ij, vec);
1072 dsvmul(inpr, vec, unc_ij);
1073 bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1077 gmx::formatString("Geometry %s not handled here",
1078 enumValueToString(coord.params.eGeom))
1087 "Pull constraint not converged: "
1089 "d_ref = %f, current d = %f\n",
1090 coord.params.group[0],
1091 coord.params.group[1],
1096 bConverged_all = FALSE;
1101 /* if after all constraints are dealt with and bConverged is still TRUE
1102 we're finished, if not we do another iteration */
1104 if (niter > max_iter)
1106 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1109 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1116 /* update atoms in the groups */
1117 for (size_t g = 0; g < pull->group.size(); g++)
1119 const pull_group_work_t* pgrp;
1122 pgrp = &pull->group[g];
1124 /* get the final constraint displacement dr for group g */
1125 dvec_sub(rnew[g], pgrp->xp, dr);
1127 if (dnorm2(dr) == 0)
1129 /* No displacement, no update necessary */
1133 /* update the atom positions */
1134 auto localAtomIndices = pgrp->atomSet.localIndex();
1136 for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1138 ii = localAtomIndices[j];
1139 if (!pgrp->localWeights.empty())
1141 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1143 for (m = 0; m < DIM; m++)
1149 for (m = 0; m < DIM; m++)
1151 v[ii][m] += invdt * tmp[m];
1157 /* calculate the constraint forces, used for output and virial only */
1158 for (size_t c = 0; c < pull->coord.size(); c++)
1160 pull_coord_work_t* pcrd;
1162 pcrd = &pull->coord[c];
1164 if (pcrd->params.eType != PullingAlgorithm::Constraint)
1169 /* Accumulate the forces, in case we have multiple constraint steps */
1172 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1174 pcrd->scalarForce += force;
1176 if (vir != nullptr && pcrd->params.eGeom != PullGroupGeometry::DirectionPBC && bMaster)
1180 /* Add the pull contribution to the virial */
1181 /* We have already checked above that r_ij[c] != 0 */
1182 f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1184 for (j = 0; j < DIM; j++)
1186 for (m = 0; m < DIM; m++)
1188 vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1194 /* finished! I hope. Give back some memory */
1200 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1202 for (int j = 0; j < DIM; j++)
1204 for (int m = 0; m < DIM; m++)
1206 vir[j][m] -= 0.5 * f[j] * dr[m];
1211 /* Adds the pull contribution to the virial */
1212 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1214 if (vir != nullptr && pcrd.params.eGeom != PullGroupGeometry::DirectionPBC)
1216 /* Add the pull contribution for each distance vector to the virial. */
1217 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1218 if (pcrd.params.ngroup >= 4)
1220 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1222 if (pcrd.params.ngroup >= 6)
1224 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1229 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1237 k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1238 dkdl = pcrd->params.kB - pcrd->params.k;
1240 switch (pcrd->params.eType)
1242 case PullingAlgorithm::Umbrella:
1243 case PullingAlgorithm::FlatBottom:
1244 case PullingAlgorithm::FlatBottomHigh:
1245 /* The only difference between an umbrella and a flat-bottom
1246 * potential is that a flat-bottom is zero above or below
1247 the reference value.
1249 if ((pcrd->params.eType == PullingAlgorithm::FlatBottom && dev < 0)
1250 || (pcrd->params.eType == PullingAlgorithm::FlatBottomHigh && dev > 0))
1255 pcrd->scalarForce = -k * dev;
1256 *V += 0.5 * k * gmx::square(dev);
1257 *dVdl += 0.5 * dkdl * gmx::square(dev);
1259 case PullingAlgorithm::ConstantForce:
1260 pcrd->scalarForce = -k;
1261 *V += k * pcrd->spatialData.value;
1262 *dVdl += dkdl * pcrd->spatialData.value;
1264 case PullingAlgorithm::External:
1266 "the scalar pull force should not be calculated internally for pull type "
1268 default: gmx_incons("Unsupported pull type in do_pull_pot");
1272 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1274 const t_pull_coord& params = pcrd.params;
1275 const PullCoordSpatialData& spatialData = pcrd.spatialData;
1277 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1278 PullCoordVectorForces forces;
1280 if (params.eGeom == PullGroupGeometry::Distance)
1282 double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1283 for (int m = 0; m < DIM; m++)
1285 forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1288 else if (params.eGeom == PullGroupGeometry::Angle)
1291 double cos_theta, cos_theta2;
1293 cos_theta = cos(spatialData.value);
1294 cos_theta2 = gmx::square(cos_theta);
1296 /* The force at theta = 0, pi is undefined so we should not apply any force.
1297 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1298 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1299 * have to check for this before dividing by their norm below.
1303 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1304 * and another vector parallel to dr23:
1305 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1306 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1308 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1309 double b = a * cos_theta;
1310 double invdr01 = 1. / dnorm(spatialData.dr01);
1311 double invdr23 = 1. / dnorm(spatialData.dr23);
1312 dvec normalized_dr01, normalized_dr23;
1313 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1314 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1316 for (int m = 0; m < DIM; m++)
1318 /* Here, f_scal is -dV/dtheta */
1320 pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1322 pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1327 /* No forces to apply for ill-defined cases*/
1328 clear_dvec(forces.force01);
1329 clear_dvec(forces.force23);
1332 else if (params.eGeom == PullGroupGeometry::AngleAxis)
1334 double cos_theta, cos_theta2;
1336 /* The angle-axis force is exactly the same as the angle force (above) except that in
1337 this case the second vector (dr23) is replaced by the pull vector. */
1338 cos_theta = cos(spatialData.value);
1339 cos_theta2 = gmx::square(cos_theta);
1345 dvec normalized_dr01;
1347 invdr01 = 1. / dnorm(spatialData.dr01);
1348 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1349 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1352 for (int m = 0; m < DIM; m++)
1355 pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1360 clear_dvec(forces.force01);
1363 else if (params.eGeom == PullGroupGeometry::Dihedral)
1365 double m2, n2, tol, sqrdist_32;
1367 /* Note: there is a small difference here compared to the
1368 dihedral force calculations in the bondeds (ref: Bekker 1994).
1369 There rij = ri - rj, while here dr01 = r1 - r0.
1370 However, all distance vectors occur in form of cross or inner products
1371 so that two signs cancel and we end up with the same expressions.
1372 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1374 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1375 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1376 dsvmul(-1, spatialData.dr23, dr32);
1377 sqrdist_32 = diprod(dr32, dr32);
1378 tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1379 if ((m2 > tol) && (n2 > tol))
1381 double a_01, a_23_01, a_23_45, a_45;
1382 double inv_dist_32, inv_sqrdist_32, dist_32;
1384 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1385 inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1386 dist_32 = sqrdist_32 * inv_dist_32;
1388 /* Forces on groups 0, 1 */
1389 a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
1390 dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1392 /* Forces on groups 4, 5 */
1393 a_45 = -pcrd.scalarForce * dist_32 / n2;
1394 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1396 /* Force on groups 2, 3 (defining the axis) */
1397 a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1398 a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1399 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1400 dsvmul(a_23_45, forces.force45, v);
1401 dvec_sub(u, v, forces.force23); /* force on group 3 */
1405 /* No force to apply for ill-defined cases */
1406 clear_dvec(forces.force01);
1407 clear_dvec(forces.force23);
1408 clear_dvec(forces.force45);
1413 for (int m = 0; m < DIM; m++)
1415 forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1423 /* We use a global mutex for locking access to the pull data structure
1424 * during registration of external pull potential providers.
1425 * We could use a different, local mutex for each pull object, but the overhead
1426 * is extremely small here and registration is only done during initialization.
1428 static std::mutex registrationMutex;
1430 using Lock = std::lock_guard<std::mutex>;
1432 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1434 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1435 GMX_RELEASE_ASSERT(provider != nullptr,
1436 "register_external_pull_potential called with NULL as provider name");
1438 if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1441 "Module '%s' attempted to register an external potential for pull coordinate %d "
1442 "which is out of the pull coordinate range %d - %zu\n",
1446 pull->coord.size());
1449 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1451 if (pcrd->params.eType != PullingAlgorithm::External)
1455 "Module '%s' attempted to register an external potential for pull coordinate %d "
1456 "which of type '%s', whereas external potentials are only supported with type '%s'",
1459 enumValueToString(pcrd->params.eType),
1460 enumValueToString(PullingAlgorithm::External));
1463 GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
1464 "The external potential provider string for a pull coordinate is NULL");
1466 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
1469 "Module '%s' attempted to register an external potential for pull coordinate %d "
1470 "which expects the external potential to be provided by a module named '%s'",
1473 pcrd->params.externalPotentialProvider.c_str());
1476 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1477 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1478 * pull->numUnregisteredExternalPotentials.
1480 Lock registrationLock(registrationMutex);
1482 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1485 "Module '%s' attempted to register an external potential for pull coordinate %d "
1491 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1492 pull->numUnregisteredExternalPotentials--;
1494 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1495 "Negative unregistered potentials, the pull code is inconsistent");
1499 static void check_external_potential_registration(const struct pull_t* pull)
1501 if (pull->numUnregisteredExternalPotentials > 0)
1504 for (c = 0; c < pull->coord.size(); c++)
1506 if (pull->coord[c].params.eType == PullingAlgorithm::External
1507 && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1513 GMX_RELEASE_ASSERT(c < pull->coord.size(),
1514 "Internal inconsistency in the pull potential provider counting");
1517 "No external provider for external pull potentials have been provided for %d "
1518 "pull coordinates. The first coordinate without provider is number %zu, which "
1519 "expects a module named '%s' to provide the external potential.",
1520 pull->numUnregisteredExternalPotentials,
1522 pull->coord[c].params.externalPotentialProvider.c_str());
1526 /* Pull takes care of adding the forces of the external potential.
1527 * The external potential module has to make sure that the corresponding
1528 * potential energy is added either to the pull term or to a term
1529 * specific to the external module.
1531 void apply_external_pull_coord_force(struct pull_t* pull,
1534 ArrayRef<const real> masses,
1535 gmx::ForceWithVirial* forceWithVirial)
1537 pull_coord_work_t* pcrd;
1539 GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1540 "apply_external_pull_coord_force called with coord_index out of range");
1542 if (pull->comm.bParticipate)
1544 pcrd = &pull->coord[coord_index];
1547 pcrd->params.eType == PullingAlgorithm::External,
1548 "The pull force can only be set externally on pull coordinates of external type");
1550 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1551 "apply_external_pull_coord_force called for an unregistered pull coordinate");
1554 pcrd->scalarForce = coord_force;
1556 /* Calculate the forces on the pull groups */
1557 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1559 /* Add the forces for this coordinate to the total virial and force */
1560 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1562 matrix virial = { { 0 } };
1563 add_virial_coord(virial, *pcrd, pullCoordForces);
1564 forceWithVirial->addVirialContribution(virial);
1568 *pcrd, pull->group, pullCoordForces, masses, as_rvec_array(forceWithVirial->force_.data()));
1571 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1574 /* Calculate the pull potential and scalar force for a pull coordinate.
1575 * Returns the vector forces for the pull coordinate.
1577 static PullCoordVectorForces do_pull_pot_coord(const pull_t& pull,
1578 pull_coord_work_t* pcrd,
1586 assert(pcrd->params.eType != PullingAlgorithm::Constraint);
1588 double dev = get_pull_coord_deviation(pull, pcrd, pbc, t);
1590 calc_pull_coord_scalar_force_and_potential(pcrd, dev, lambda, V, dVdl);
1592 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1594 add_virial_coord(vir, *pcrd, pullCoordForces);
1596 return pullCoordForces;
1599 real pull_potential(struct pull_t* pull,
1600 ArrayRef<const real> masses,
1602 const t_commrec* cr,
1605 ArrayRef<const RVec> x,
1606 gmx::ForceWithVirial* force,
1611 assert(pull != nullptr);
1613 /* Ideally we should check external potential registration only during
1614 * the initialization phase, but that requires another function call
1615 * that should be done exactly in the right order. So we check here.
1617 check_external_potential_registration(pull);
1619 if (pull->comm.bParticipate)
1623 pull_calc_coms(cr, pull, masses, pbc, t, x, {});
1625 rvec* f = as_rvec_array(force->force_.data());
1626 matrix virial = { { 0 } };
1627 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1628 for (size_t c = 0; c < pull->coord.size(); c++)
1630 pull_coord_work_t* pcrd;
1631 pcrd = &pull->coord[c];
1633 /* For external potential the force is assumed to be given by an external module by a
1634 call to apply_pull_coord_external_force */
1635 if (pcrd->params.eType == PullingAlgorithm::Constraint
1636 || pcrd->params.eType == PullingAlgorithm::External)
1641 PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1642 *pull, pcrd, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1644 /* Distribute the force over the atoms in the pulled groups */
1645 apply_forces_coord(*pcrd, pull->group, pullCoordForces, masses, f);
1650 force->addVirialContribution(virial);
1655 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1656 "Too few or too many external pull potentials have been applied the previous step");
1657 /* All external pull potentials still need to be applied */
1658 pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1660 return (MASTER(cr) ? V : 0.0);
1663 void pull_constraint(struct pull_t* pull,
1664 ArrayRef<const real> masses,
1666 const t_commrec* cr,
1674 assert(pull != nullptr);
1676 if (pull->comm.bParticipate)
1678 pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
1680 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1684 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1688 gmx_bool bMustParticipate;
1694 /* We always make the master node participate, such that it can do i/o,
1695 * add the virial and to simplify MC type extensions people might have.
1697 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1699 for (pull_group_work_t& group : pull->group)
1701 if (!group.globalWeights.empty())
1703 group.localWeights.resize(group.atomSet.numAtomsLocal());
1704 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1706 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1710 GMX_ASSERT(bMustParticipate || dd != nullptr,
1711 "Either all ranks (including this rank) participate, or we use DD and need to "
1712 "have access to dd here");
1714 /* We should participate if we have pull or pbc atoms */
1715 if (!bMustParticipate
1716 && (group.atomSet.numAtomsLocal() > 0
1717 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1719 bMustParticipate = TRUE;
1723 if (!comm->bParticipateAll)
1725 /* Keep currently not required ranks in the communicator
1726 * if they needed to participate up to 20 decompositions ago.
1727 * This avoids frequent rebuilds due to atoms jumping back and forth.
1729 const int64_t history_count = 20;
1730 gmx_bool bWillParticipate;
1733 /* Increase the decomposition counter for the current call */
1734 comm->setup_count++;
1736 if (bMustParticipate)
1738 comm->must_count = comm->setup_count;
1743 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1745 if (debug && dd != nullptr)
1748 "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
1750 gmx::boolToString(bMustParticipate),
1751 gmx::boolToString(bWillParticipate));
1754 if (bWillParticipate)
1756 /* Count the number of ranks that we want to have participating */
1758 /* Count the number of ranks that need to be added */
1759 count[1] = comm->bParticipate ? 0 : 1;
1767 /* The cost of this global operation will be less that the cost
1768 * of the extra MPI_Comm_split calls that we can avoid.
1770 gmx_sumi(2, count, cr);
1772 /* If we are missing ranks or if we have 20% more ranks than needed
1773 * we make a new sub-communicator.
1775 if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1779 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1783 if (comm->mpi_comm_com != MPI_COMM_NULL)
1785 MPI_Comm_free(&comm->mpi_comm_com);
1788 /* This might be an extremely expensive operation, so we try
1789 * to avoid this splitting as much as possible.
1791 assert(dd != nullptr);
1792 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1795 comm->bParticipate = bWillParticipate;
1796 comm->nparticipate = count[0];
1798 /* When we use the previous COM for PBC, we need to broadcast
1799 * the previous COM to ranks that have joined the communicator.
1801 for (pull_group_work_t& group : pull->group)
1803 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1805 GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1806 "The master rank has to participate, as it should pass an up to "
1808 "to bcast here as well as to e.g. checkpointing");
1810 gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
1816 /* Since the PBC of atoms might have changed, we need to update the PBC */
1817 pull->bSetPBCatoms = TRUE;
1820 static void init_pull_group_index(FILE* fplog,
1821 const t_commrec* cr,
1823 pull_group_work_t* pg,
1824 gmx_bool bConstraint,
1825 const ivec pulldim_con,
1826 const gmx_mtop_t& mtop,
1827 const t_inputrec* ir,
1830 /* With EM and BD there are no masses in the integrator.
1831 * But we still want to have the correct mass-weighted COMs.
1832 * So we store the real masses in the weights.
1834 const bool setWeights = (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1835 || ir->eI == IntegrationAlgorithm::BD);
1837 /* In parallel, store we need to extract localWeights from weights at DD time */
1838 std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1840 const SimulationGroups& groups = mtop.groups;
1842 /* Count frozen dimensions and (weighted) mass */
1848 for (int i = 0; i < int(pg->params.ind.size()); i++)
1850 int ii = pg->params.ind[i];
1851 if (bConstraint && ir->opts.nFreeze)
1853 for (int d = 0; d < DIM; d++)
1855 if (pulldim_con[d] == 1
1856 && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1862 const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1864 if (ir->efep == FreeEnergyPerturbationType::No)
1870 m = (1 - lambda) * atom.m + lambda * atom.mB;
1873 if (!pg->params.weight.empty())
1875 w = pg->params.weight[i];
1881 if (EI_ENERGY_MINIMIZATION(ir->eI))
1883 /* Move the mass to the weight */
1887 else if (ir->eI == IntegrationAlgorithm::BD)
1890 if (ir->bd_fric != 0.0F)
1892 mbd = ir->bd_fric * ir->delta_t;
1896 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1898 mbd = ir->delta_t / ir->opts.tau_t[0];
1903 / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1911 weights.push_back(w);
1915 wwmass += m * w * w;
1920 /* We can have single atom groups with zero mass with potential pulling
1921 * without cosine weighting.
1923 if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1925 /* With one atom the mass doesn't matter */
1931 "The total%s mass of pull group %d is zero",
1932 !pg->params.weight.empty() ? " weighted" : "",
1938 fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
1939 if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1940 || ir->eI == IntegrationAlgorithm::BD)
1942 fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1944 if (pg->epgrppbc == epgrppbcCOS)
1946 fprintf(fplog, ", cosine weighting will be used");
1948 fprintf(fplog, "\n");
1953 /* A value != 0 signals not frozen, it is updated later */
1959 for (int d = 0; d < DIM; d++)
1961 ndim += pulldim_con[d] * pg->params.ind.size();
1963 if (fplog && nfrozen > 0 && nfrozen < ndim)
1966 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1967 " that are subject to pulling are frozen.\n"
1968 " For constraint pulling the whole group will be frozen.\n\n",
1976 struct pull_t* init_pull(FILE* fplog,
1977 const pull_params_t* pull_params,
1978 const t_inputrec* ir,
1979 const gmx_mtop_t& mtop,
1980 const t_commrec* cr,
1981 gmx::LocalAtomSetManager* atomSets,
1984 struct pull_t* pull;
1989 /* Copy the pull parameters */
1990 pull->params = *pull_params;
1992 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
1993 const int maxNumThreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
1995 for (int i = 0; i < pull_params->ngroup; ++i)
1997 pull->group.emplace_back(pull_params->group[i],
1998 atomSets->add(pull_params->group[i].ind),
1999 pull_params->bSetPbcRefToPrevStepCOM,
2003 if (cr != nullptr && DOMAINDECOMP(cr))
2005 /* Set up the global to local atom mapping for PBC atoms */
2006 for (pull_group_work_t& group : pull->group)
2008 if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
2010 /* pbcAtomSet consists of a single atom */
2011 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
2012 atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
2017 pull->bPotential = FALSE;
2018 pull->bConstraint = FALSE;
2019 pull->bCylinder = FALSE;
2020 pull->bAngle = FALSE;
2021 pull->bXOutAverage = pull_params->bXOutAverage;
2022 pull->bFOutAverage = pull_params->bFOutAverage;
2024 GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
2025 "pull group 0 is an absolute reference group and should not contain atoms");
2027 pull->numCoordinatesWithExternalPotential = 0;
2029 for (int c = 0; c < pull->params.ncoord; c++)
2031 GMX_RELEASE_ASSERT(pull_params->coord[c].coordIndex == c,
2032 "The stored index should match the position in the vector");
2034 /* Construct a pull coordinate, copying all coordinate parameters */
2035 pull->coord.emplace_back(pull_params->coord[c]);
2037 pull_coord_work_t* pcrd = &pull->coord.back();
2039 switch (pcrd->params.eGeom)
2041 case PullGroupGeometry::Distance:
2042 case PullGroupGeometry::DirectionRelative: /* Direction vector is determined at each step */
2043 case PullGroupGeometry::Angle:
2044 case PullGroupGeometry::Dihedral: break;
2045 case PullGroupGeometry::Direction:
2046 case PullGroupGeometry::DirectionPBC:
2047 case PullGroupGeometry::Cylinder:
2048 case PullGroupGeometry::AngleAxis:
2049 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2052 /* We allow reading of newer tpx files with new pull geometries,
2053 * but with the same tpx format, with old code. A new geometry
2054 * only adds a new enum value, which will be out of range for
2055 * old code. The first place we need to generate an error is
2056 * here, since the pull code can't handle this.
2057 * The enum can be used before arriving here only for printing
2058 * the string corresponding to the geometry, which will result
2059 * in printing "UNDEFINED".
2062 "Pull geometry not supported for pull coordinate %d. The geometry enum "
2063 "%s in the input is larger than that supported by the code (up to %d). "
2064 "You are probably reading a tpr file generated with a newer version of "
2065 "GROMACS with an binary from an older version of Gromacs.",
2067 enumValueToString(pcrd->params.eGeom),
2068 static_cast<int>(PullGroupGeometry::Count) - 1);
2071 if (pcrd->params.eType == PullingAlgorithm::Constraint)
2073 /* Check restrictions of the constraint pull code */
2074 if (pcrd->params.eGeom == PullGroupGeometry::Cylinder
2075 || pcrd->params.eGeom == PullGroupGeometry::DirectionRelative
2076 || pcrd->params.eGeom == PullGroupGeometry::Angle
2077 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2078 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2081 "Pulling of type %s can not be combined with geometry %s. Consider using "
2083 enumValueToString(pcrd->params.eType),
2084 enumValueToString(pcrd->params.eGeom),
2085 enumValueToString(PullingAlgorithm::Umbrella));
2090 "Constraint pulling can not be combined with multiple time stepping");
2092 pull->bConstraint = TRUE;
2096 pull->bPotential = TRUE;
2099 if (pcrd->params.eGeom == PullGroupGeometry::Cylinder)
2101 pull->bCylinder = TRUE;
2103 else if (pcrd->params.eGeom == PullGroupGeometry::Angle
2104 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2105 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2107 pull->bAngle = TRUE;
2110 /* We only need to calculate the plain COM of a group
2111 * when it is not only used as a cylinder group.
2112 * Also the absolute reference group 0 needs no COM computation.
2114 for (int i = 0; i < pcrd->params.ngroup; i++)
2116 int groupIndex = pcrd->params.group[i];
2117 if (groupIndex > 0 && !(pcrd->params.eGeom == PullGroupGeometry::Cylinder && i == 0))
2119 pull->group[groupIndex].needToCalcCom = true;
2123 /* With non-zero rate the reference value is set at every step */
2124 if (pcrd->params.rate == 0)
2126 /* Initialize the constant reference value */
2127 if (pcrd->params.eType != PullingAlgorithm::External)
2129 pcrd->value_ref = sanitizePullCoordReferenceValue(
2131 pcrd->params.init * pull_conversion_factor_userinput2internal(pcrd->params));
2135 /* With an external pull potential, the reference value loses
2136 * it's meaning and should not be used. Setting it to zero
2137 * makes any terms dependent on it disappear.
2138 * The only issue this causes is that with cylinder pulling
2139 * no atoms of the cylinder group within the cylinder radius
2140 * should be more than half a box length away from the COM of
2141 * the pull group along the axial direction.
2143 pcrd->value_ref = 0.0;
2147 if (pcrd->params.eType == PullingAlgorithm::External)
2150 pcrd->params.rate == 0,
2151 "With an external potential, a pull coordinate should have rate = 0");
2153 /* This external potential needs to be registered later */
2154 pull->numCoordinatesWithExternalPotential++;
2156 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2159 pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
2160 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2162 pull->pbcType = ir->pbcType;
2163 switch (pull->pbcType)
2165 case PbcType::No: pull->npbcdim = 0; break;
2166 case PbcType::XY: pull->npbcdim = 2; break;
2167 default: pull->npbcdim = 3; break;
2172 gmx_bool bAbs, bCos;
2175 for (const pull_coord_work_t& coord : pull->coord)
2177 if (pull->group[coord.params.group[0]].params.ind.empty()
2178 || pull->group[coord.params.group[1]].params.ind.empty())
2184 fprintf(fplog, "\n");
2185 if (pull->bPotential)
2187 fprintf(fplog, "Will apply potential COM pulling\n");
2189 if (pull->bConstraint)
2191 fprintf(fplog, "Will apply constraint COM pulling\n");
2193 // Don't include the reference group 0 in output, so we report ngroup-1
2194 int numRealGroups = pull->group.size() - 1;
2195 GMX_RELEASE_ASSERT(numRealGroups > 0,
2196 "The reference absolute position pull group should always be present");
2198 "with %zu pull coordinate%s and %d group%s\n",
2200 pull->coord.size() == 1 ? "" : "s",
2202 numRealGroups == 1 ? "" : "s");
2205 fprintf(fplog, "with an absolute reference\n");
2208 // Don't include the reference group 0 in loop
2209 for (size_t g = 1; g < pull->group.size(); g++)
2211 if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
2213 /* We are using cosine weighting */
2214 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2220 please_cite(fplog, "Engin2010");
2224 pull->bRefAt = FALSE;
2226 for (size_t g = 0; g < pull->group.size(); g++)
2228 pull_group_work_t* pgrp;
2230 pgrp = &pull->group[g];
2231 if (!pgrp->params.ind.empty())
2233 /* There is an issue when a group is used in multiple coordinates
2234 * and constraints are applied in different dimensions with atoms
2235 * frozen in some, but not all dimensions.
2236 * Since there is only one mass array per group, we can't have
2237 * frozen/non-frozen atoms for different coords at the same time.
2238 * But since this is a very exotic case, we don't check for this.
2239 * A warning is printed in init_pull_group_index.
2242 gmx_bool bConstraint;
2243 ivec pulldim, pulldim_con;
2245 /* Loop over all pull coordinates to see along which dimensions
2246 * this group is pulled and if it is involved in constraints.
2248 bConstraint = FALSE;
2249 clear_ivec(pulldim);
2250 clear_ivec(pulldim_con);
2251 for (const pull_coord_work_t& coord : pull->coord)
2253 gmx_bool bGroupUsed = FALSE;
2254 for (int gi = 0; gi < coord.params.ngroup; gi++)
2256 if (coord.params.group[gi] == static_cast<int>(g))
2264 for (int m = 0; m < DIM; m++)
2266 if (coord.params.dim[m] == 1)
2270 if (coord.params.eType == PullingAlgorithm::Constraint)
2280 /* Determine if we need to take PBC into account for calculating
2281 * the COM's of the pull groups.
2283 switch (pgrp->epgrppbc)
2285 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2287 if (!pgrp->params.weight.empty())
2290 "Pull groups can not have relative weights and cosine weighting "
2293 for (int m = 0; m < DIM; m++)
2295 if (m < pull->npbcdim && pulldim[m] == 1)
2297 if (pull->cosdim >= 0 && pull->cosdim != m)
2300 "Can only use cosine weighting with pulling in one "
2301 "dimension (use mdp option pull-coord?-dim)");
2307 case epgrppbcNONE: break;
2310 /* Set the indices */
2311 init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2315 /* Absolute reference, set the inverse mass to zero.
2316 * This is only relevant (and used) with constraint pulling.
2323 /* If we use cylinder coordinates, do some initialising for them */
2324 if (pull->bCylinder)
2326 for (pull_coord_work_t& coord : pull->coord)
2328 if (coord.params.eGeom == PullGroupGeometry::Cylinder)
2330 if (pull->group[coord.params.group[0]].params.ind.empty())
2333 "A cylinder pull group is not supported when using absolute "
2337 const auto& group0 = pull->group[coord.params.group[0]];
2338 coord.dynamicGroup0 = std::make_unique<pull_group_work_t>(
2339 group0.params, group0.atomSet, pull->params.bSetPbcRefToPrevStepCOM, maxNumThreads);
2343 pull->comSums.resize(maxNumThreads);
2348 /* Use a sub-communicator when we have more than 32 ranks, but not
2349 * when we have an external pull potential, since then the external
2350 * potential provider expects each rank to have the coordinate.
2352 comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
2353 || pull->numCoordinatesWithExternalPotential > 0
2354 || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2355 /* This sub-commicator is not used with comm->bParticipateAll,
2356 * so we can always initialize it to NULL.
2358 comm->mpi_comm_com = MPI_COMM_NULL;
2359 comm->nparticipate = 0;
2360 comm->isMasterRank = (cr == nullptr || MASTER(cr));
2362 /* No MPI: 1 rank: all ranks pull */
2363 comm->bParticipateAll = TRUE;
2364 comm->isMasterRank = true;
2366 comm->bParticipate = comm->bParticipateAll;
2367 comm->setup_count = 0;
2368 comm->must_count = 0;
2370 if (!comm->bParticipateAll && fplog != nullptr)
2372 fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2375 comm->pbcAtomBuffer.resize(pull->group.size());
2376 comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
2377 if (pull->bCylinder)
2379 comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
2382 /* We still need to initialize the PBC reference coordinates */
2383 pull->bSetPBCatoms = TRUE;
2385 pull->out_x = nullptr;
2386 pull->out_f = nullptr;
2391 static void destroy_pull(struct pull_t* pull)
2394 if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2396 MPI_Comm_free(&pull->comm.mpi_comm_com);
2403 void preparePrevStepPullCom(const t_inputrec* ir,
2405 ArrayRef<const real> masses,
2407 const t_state* state_global,
2408 const t_commrec* cr,
2409 bool startingFromCheckpoint)
2411 if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
2415 allocStatePrevStepPullCom(state, pull_work);
2416 if (startingFromCheckpoint)
2420 state->pull_com_prev_step = state_global->pull_com_prev_step;
2424 /* Only the master rank has the checkpointed COM from the previous step */
2425 gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(),
2426 &state->pull_com_prev_step[0],
2427 cr->mpi_comm_mygroup);
2429 setPrevStepPullComFromState(pull_work, state);
2434 set_pbc(&pbc, ir->pbcType, state->box);
2435 initPullComFromPrevStep(
2436 cr, pull_work, masses, &pbc, state->x.arrayRefWithPadding().unpaddedArrayRef());
2437 updatePrevStepPullCom(pull_work, state);
2441 void finish_pull(struct pull_t* pull)
2443 check_external_potential_registration(pull);
2447 gmx_fio_fclose(pull->out_x);
2451 gmx_fio_fclose(pull->out_f);
2457 bool pull_have_potential(const pull_t& pull)
2459 return pull.bPotential;
2462 bool pull_have_constraint(const pull_t& pull)
2464 return pull.bConstraint;
2467 bool pull_have_constraint(const pull_params_t& pullParameters)
2469 for (int c = 0; c < pullParameters.ncoord; c++)
2471 if (pullParameters.coord[c].eType == PullingAlgorithm::Constraint)