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.
54 #include "gromacs/commandline/filenm.h"
55 #include "gromacs/domdec/domdec_struct.h"
56 #include "gromacs/domdec/localatomset.h"
57 #include "gromacs/domdec/localatomsetmanager.h"
58 #include "gromacs/fileio/gmxfio.h"
59 #include "gromacs/gmxlib/network.h"
60 #include "gromacs/math/functions.h"
61 #include "gromacs/math/units.h"
62 #include "gromacs/math/utilities.h"
63 #include "gromacs/math/vec.h"
64 #include "gromacs/math/vectypes.h"
65 #include "gromacs/mdlib/gmx_omp_nthreads.h"
66 #include "gromacs/mdtypes/commrec.h"
67 #include "gromacs/mdtypes/forceoutput.h"
68 #include "gromacs/mdtypes/inputrec.h"
69 #include "gromacs/mdtypes/md_enums.h"
70 #include "gromacs/mdtypes/state.h"
71 #include "gromacs/pbcutil/pbc.h"
72 #include "gromacs/topology/mtop_lookup.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/arrayref.h"
75 #include "gromacs/utility/cstringutil.h"
76 #include "gromacs/utility/enumerationhelpers.h"
77 #include "gromacs/utility/exceptions.h"
78 #include "gromacs/utility/fatalerror.h"
79 #include "gromacs/utility/futil.h"
80 #include "gromacs/utility/gmxassert.h"
81 #include "gromacs/utility/pleasecite.h"
82 #include "gromacs/utility/real.h"
83 #include "gromacs/utility/smalloc.h"
84 #include "gromacs/utility/strconvert.h"
85 #include "gromacs/utility/stringutil.h"
87 #include "pull_internal.h"
88 #include "transformationcoordinate.h"
92 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
98 /*! \brief Tells whether the pull geometry is an angle type */
99 constexpr gmx::EnumerationArray<PullGroupGeometry, bool> sc_isAngleType = {
100 { false, false, false, false, false, true, true, true }
103 static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
105 if (params.ind.size() <= 1)
107 /* no PBC required */
110 else if (params.pbcatom >= 0)
112 if (setPbcRefToPrevStepCOM)
114 return epgrppbcPREVSTEPCOM;
118 return epgrppbcREFAT;
127 /* NOTE: The params initialization currently copies pointers.
128 * So the lifetime of the source, currently always inputrec,
129 * should not end before that of this object.
130 * This will be fixed when the pointers are replacted by std::vector.
132 pull_group_work_t::pull_group_work_t(const t_pull_group& params,
133 gmx::LocalAtomSet atomSet,
134 bool bSetPbcRefToPrevStepCOM,
137 epgrppbc(groupPbcFromParams(params, bSetPbcRefToPrevStepCOM)),
138 maxNumThreads(maxNumThreads),
139 needToCalcCom(false),
149 static bool pull_coordinate_is_directional(const t_pull_coord& pcrd)
151 return (pcrd.eGeom == PullGroupGeometry::Direction || pcrd.eGeom == PullGroupGeometry::DirectionPBC
152 || pcrd.eGeom == PullGroupGeometry::DirectionRelative
153 || pcrd.eGeom == PullGroupGeometry::Cylinder);
156 const char* pull_coordinate_units(const t_pull_coord& pcrd)
158 return sc_isAngleType[pcrd.eGeom] ? "deg" : "nm";
161 double pull_conversion_factor_userinput2internal(const t_pull_coord& pcrd)
163 if (sc_isAngleType[pcrd.eGeom])
165 return gmx::c_deg2Rad;
173 double pull_conversion_factor_internal2userinput(const t_pull_coord& pcrd)
175 if (sc_isAngleType[pcrd.eGeom])
177 return gmx::c_rad2Deg;
185 /* Apply forces in a mass weighted fashion for part of the pull group */
186 static void apply_forces_grp_part(const pull_group_work_t& pgrp,
189 ArrayRef<const real> masses,
194 const double inv_wm = pgrp.mwscale;
196 auto localAtomIndices = pgrp.atomSet.localIndex();
197 for (int i = ind_start; i < ind_end; i++)
199 int ii = localAtomIndices[i];
200 double wmass = masses[ii];
201 if (!pgrp.localWeights.empty())
203 wmass *= pgrp.localWeights[i];
206 for (int d = 0; d < DIM; d++)
208 f[ii][d] += sign * wmass * f_pull[d] * inv_wm;
213 /* Apply forces in a mass weighted fashion */
214 static void apply_forces_grp(const pull_group_work_t& pgrp,
215 ArrayRef<const real> masses,
220 auto localAtomIndices = pgrp.atomSet.localIndex();
222 if (pgrp.params.ind.size() == 1 && pgrp.atomSet.numAtomsLocal() == 1)
224 /* Only one atom and our rank has this atom: we can skip
225 * the mass weighting, which means that this code also works
226 * for mass=0, e.g. with a virtual site.
228 for (int d = 0; d < DIM; d++)
230 f[localAtomIndices[0]][d] += sign * f_pull[d];
235 const int numThreads = pgrp.numThreads();
238 apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), masses, f_pull, sign, f);
242 #pragma omp parallel for num_threads(numThreads) schedule(static)
243 for (int th = 0; th < numThreads; th++)
245 int ind_start = (localAtomIndices.size() * (th + 0)) / numThreads;
246 int ind_end = (localAtomIndices.size() * (th + 1)) / numThreads;
247 apply_forces_grp_part(pgrp, ind_start, ind_end, masses, f_pull, sign, f);
253 /* Apply forces in a mass weighted fashion to a cylinder group */
254 static void apply_forces_cyl_grp(const pull_group_work_t& pgrp,
255 const double dv_corr,
256 ArrayRef<const real> masses,
262 const double inv_wm = pgrp.mwscale;
264 auto localAtomIndices = pgrp.atomSet.localIndex();
266 /* The cylinder group is always a slab in the system, thus large.
267 * Therefore we always thread-parallelize this group.
269 const int numAtomsLocal = localAtomIndices.size();
270 const int gmx_unused numThreads = pgrp.numThreads();
271 #pragma omp parallel for num_threads(numThreads) schedule(static)
272 for (int i = 0; i < numAtomsLocal; i++)
274 double weight = pgrp.localWeights[i];
279 int ii = localAtomIndices[i];
280 double mass = masses[ii];
281 /* The stored axial distance from the cylinder center (dv) needs
282 * to be corrected for an offset (dv_corr), which was unknown when
285 double dv_com = pgrp.dv[i] + dv_corr;
287 /* Here we not only add the pull force working along vec (f_pull),
288 * but also a radial component, due to the dependence of the weights
289 * on the radial distance.
291 for (int m = 0; m < DIM; m++)
293 f[ii][m] += sign * inv_wm * (mass * weight * f_pull[m] + pgrp.mdw[i][m] * dv_com * f_scal);
298 /* Apply torque forces in a mass weighted fashion to the groups that define
299 * the pull vector direction for pull coordinate pcrd.
301 static void apply_forces_vec_torque(const pull_coord_work_t& pcrd,
302 ArrayRef<const pull_group_work_t> pullGroups,
303 ArrayRef<const real> masses,
306 const PullCoordSpatialData& spatialData = pcrd.spatialData;
308 /* The component inpr along the pull vector is accounted for in the usual
309 * way. Here we account for the component perpendicular to vec.
312 for (int m = 0; m < DIM; m++)
314 inpr += spatialData.dr01[m] * spatialData.vec[m];
316 /* The torque force works along the component of the distance vector
317 * of between the two "usual" pull groups that is perpendicular to
318 * the pull vector. The magnitude of this force is the "usual" scale force
319 * multiplied with the ratio of the distance between the two "usual" pull
320 * groups and the distance between the two groups that define the vector.
323 for (int m = 0; m < DIM; m++)
325 f_perp[m] = (spatialData.dr01[m] - inpr * spatialData.vec[m]) / spatialData.vec_len * pcrd.scalarForce;
328 /* Apply the force to the groups defining the vector using opposite signs */
329 apply_forces_grp(pullGroups[pcrd.params.group[2]], masses, f_perp, -1, f);
330 apply_forces_grp(pullGroups[pcrd.params.group[3]], masses, f_perp, 1, f);
333 /* Apply forces in a mass weighted fashion */
334 static void apply_forces_coord(const pull_coord_work_t& pcrd,
335 ArrayRef<const pull_group_work_t> pullGroups,
336 const PullCoordVectorForces& forces,
337 ArrayRef<const real> masses,
340 /* Here it would be more efficient to use one large thread-parallel
341 * region instead of potential parallel regions within apply_forces_grp.
342 * But there could be overlap between pull groups and this would lead
345 * This may also lead to potential issues with force redistribution
346 * for transformation pull coordinates
349 if (pcrd.params.eGeom == PullGroupGeometry::Cylinder)
351 apply_forces_cyl_grp(
352 *pcrd.dynamicGroup0, pcrd.spatialData.cyl_dev, masses, forces.force01, pcrd.scalarForce, -1, f);
354 /* Sum the force along the vector and the radial force */
356 for (int m = 0; m < DIM; m++)
358 f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
360 apply_forces_grp(pullGroups[pcrd.params.group[1]], masses, f_tot, 1, f);
362 else if (pcrd.params.eGeom == PullGroupGeometry::Transformation)
368 if (pcrd.params.eGeom == PullGroupGeometry::DirectionRelative)
370 /* We need to apply the torque forces to the pull groups
371 * that define the pull vector.
373 apply_forces_vec_torque(pcrd, pullGroups, masses, f);
376 if (!pullGroups[pcrd.params.group[0]].params.ind.empty())
378 apply_forces_grp(pullGroups[pcrd.params.group[0]], masses, forces.force01, -1, f);
380 apply_forces_grp(pullGroups[pcrd.params.group[1]], masses, forces.force01, 1, f);
382 if (pcrd.params.ngroup >= 4)
384 apply_forces_grp(pullGroups[pcrd.params.group[2]], masses, forces.force23, -1, f);
385 apply_forces_grp(pullGroups[pcrd.params.group[3]], masses, forces.force23, 1, f);
387 if (pcrd.params.ngroup >= 6)
389 apply_forces_grp(pullGroups[pcrd.params.group[4]], masses, forces.force45, -1, f);
390 apply_forces_grp(pullGroups[pcrd.params.group[5]], masses, forces.force45, 1, f);
395 real max_pull_distance2(const pull_coord_work_t& pcrd, const t_pbc& pbc)
397 /* Note that this maximum distance calculation is more complex than
398 * most other cases in GROMACS, since here we have to take care of
399 * distance calculations that don't involve all three dimensions.
400 * For example, we can use distances that are larger than the
401 * box X and Y dimensions for a box that is elongated along Z.
404 real max_d2 = GMX_REAL_MAX;
406 if (pull_coordinate_is_directional(pcrd.params))
408 /* Directional pulling along along direction pcrd.vec.
409 * Calculating the exact maximum distance is complex and bug-prone.
410 * So we take a safe approach by not allowing distances that
411 * are larger than half the distance between unit cell faces
412 * along dimensions involved in pcrd.vec.
414 for (int m = 0; m < DIM; m++)
416 if (m < pbc.ndim_ePBC && pcrd.spatialData.vec[m] != 0)
418 real imageDistance2 = gmx::square(pbc.box[m][m]);
419 for (int d = m + 1; d < DIM; d++)
421 imageDistance2 -= gmx::square(pbc.box[d][m]);
423 max_d2 = std::min(max_d2, imageDistance2);
429 /* Distance pulling along dimensions with pcrd.params.dim[d]==1.
430 * We use half the minimum box vector length of the dimensions involved.
431 * This is correct for all cases, except for corner cases with
432 * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
433 * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
434 * in such corner cases the user could get correct results,
435 * depending on the details of the setup, we avoid further
436 * code complications.
438 for (int m = 0; m < DIM; m++)
440 if (m < pbc.ndim_ePBC && pcrd.params.dim[m] != 0)
442 real imageDistance2 = gmx::square(pbc.box[m][m]);
443 for (int d = 0; d < m; d++)
445 if (pcrd.params.dim[d] != 0)
447 imageDistance2 += gmx::square(pbc.box[m][d]);
450 max_d2 = std::min(max_d2, imageDistance2);
455 return 0.25 * max_d2;
458 /* This function returns the distance based on coordinates xg and xref.
459 * Note that the pull coordinate struct pcrd is not modified.
461 * \param[in] pull The pull struct
462 * \param[in] pcrd The pull coordinate to compute a distance for
463 * \param[in] pbc The periodic boundary conditions
464 * \param[in] xg The coordinate of group 1
465 * \param[in] xref The coordinate of group 0
466 * \param[in] groupIndex0 The index of group 0 in the pcrd->params.group
467 * \param[in] groupIndex1 The index of group 1 in the pcrd->params.group
468 * \param[in] max_dist2 The maximum distance squared
469 * \param[out] dr The distance vector
471 static void low_get_pull_coord_dr(const pull_t& pull,
472 const pull_coord_work_t& pcrd,
476 const int groupIndex0,
477 const int groupIndex1,
478 const double max_dist2,
481 const pull_group_work_t& pgrp0 = pull.group[pcrd.params.group[0]];
483 // Group coordinate, to be updated with the reference position
486 /* Only the first group can be an absolute reference, in that case nat=0 */
487 if (pgrp0.params.ind.empty())
489 for (int m = 0; m < DIM; m++)
491 xrefr[m] = pcrd.params.origin[m];
496 copy_dvec(xref, xrefr);
499 dvec dref = { 0, 0, 0 };
500 if (pcrd.params.eGeom == PullGroupGeometry::DirectionPBC)
502 for (int m = 0; m < DIM; m++)
504 dref[m] = pcrd.value_ref * pcrd.spatialData.vec[m];
506 /* Add the reference position, so we use the correct periodic image */
507 dvec_inc(xrefr, dref);
510 pbc_dx_d(&pbc, xg, xrefr, dr);
512 bool directional = pull_coordinate_is_directional(pcrd.params);
514 for (int m = 0; m < DIM; m++)
516 dr[m] *= pcrd.params.dim[m];
517 if (pcrd.params.dim[m] && !(directional && pcrd.spatialData.vec[m] == 0))
519 dr2 += dr[m] * dr[m];
522 /* Check if we are close to switching to another periodic image */
523 if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
525 /* Note that technically there is no issue with switching periodic
526 * image, as pbc_dx_d returns the distance to the closest periodic
527 * image. However in all cases where periodic image switches occur,
528 * the pull results will be useless in practice.
531 "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
532 "box size (%f).\n%s",
533 pcrd.params.group[groupIndex0],
534 pcrd.params.group[groupIndex1],
536 sqrt(0.98 * 0.98 * max_dist2),
537 pcrd.params.eGeom == PullGroupGeometry::Direction
538 ? "You might want to consider using \"pull-geometry = "
539 "direction-periodic\" instead.\n"
543 if (pcrd.params.eGeom == PullGroupGeometry::DirectionPBC)
549 /* This function returns the distance based on the contents of the pull struct.
550 * pull->coord[coord_ind].dr, and potentially vector, are updated.
552 static void get_pull_coord_dr(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc& pbc)
554 const int coord_ind = pcrd->params.coordIndex;
556 PullCoordSpatialData& spatialData = pcrd->spatialData;
559 /* With AWH pulling we allow for periodic pulling with geometry=direction.
560 * TODO: Store a periodicity flag instead of checking for external pull provider.
562 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC
563 || (pcrd->params.eGeom == PullGroupGeometry::Direction
564 && pcrd->params.eType == PullingAlgorithm::External))
570 md2 = static_cast<double>(max_pull_distance2(*pcrd, pbc));
573 if (pcrd->params.eGeom == PullGroupGeometry::DirectionRelative)
575 /* We need to determine the pull vector */
579 const pull_group_work_t& pgrp2 = pull.group[pcrd->params.group[2]];
580 const pull_group_work_t& pgrp3 = pull.group[pcrd->params.group[3]];
582 pbc_dx_d(&pbc, pgrp3.x, pgrp2.x, vec);
584 for (m = 0; m < DIM; m++)
586 vec[m] *= pcrd->params.dim[m];
588 spatialData.vec_len = dnorm(vec);
589 for (m = 0; m < DIM; m++)
591 spatialData.vec[m] = vec[m] / spatialData.vec_len;
596 "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
603 spatialData.vec[ZZ]);
607 const pull_group_work_t& pgrp0 = pull.group[pcrd->params.group[0]];
608 const pull_group_work_t& pgrp1 = pull.group[pcrd->params.group[1]];
610 low_get_pull_coord_dr(
615 pcrd->params.eGeom == PullGroupGeometry::Cylinder ? pcrd->dynamicGroup0->x : pgrp0.x,
621 if (pcrd->params.ngroup >= 4)
623 const pull_group_work_t& pgrp2 = pull.group[pcrd->params.group[2]];
624 const pull_group_work_t& pgrp3 = pull.group[pcrd->params.group[3]];
626 low_get_pull_coord_dr(pull, *pcrd, pbc, pgrp3.x, pgrp2.x, 2, 3, md2, spatialData.dr23);
628 if (pcrd->params.ngroup >= 6)
630 const pull_group_work_t& pgrp4 = pull.group[pcrd->params.group[4]];
631 const pull_group_work_t& pgrp5 = pull.group[pcrd->params.group[5]];
633 low_get_pull_coord_dr(pull, *pcrd, pbc, pgrp5.x, pgrp4.x, 4, 5, md2, spatialData.dr45);
637 /* Modify x so that it is periodic in [-pi, pi)
638 * It is assumed that x is in [-3pi, 3pi) so that x
639 * needs to be shifted by at most one period.
641 static void make_periodic_2pi(double* x)
653 /* This function should always be used to get values for setting value_ref in pull_coord_work_t */
654 static double sanitizePullCoordReferenceValue(const t_pull_coord& pcrdParams, double value_ref)
656 GMX_ASSERT(pcrdParams.eType != PullingAlgorithm::External,
657 "The pull coord reference value should not be used with type external-potential");
659 if (pcrdParams.eGeom == PullGroupGeometry::Distance)
664 "Pull reference distance for coordinate %d (%f) needs to be non-negative",
665 pcrdParams.coordIndex + 1,
669 else if (pcrdParams.eGeom == PullGroupGeometry::Angle || pcrdParams.eGeom == PullGroupGeometry::AngleAxis)
671 if (value_ref < 0 || value_ref > M_PI)
674 "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
675 "interval [0,180] deg",
676 pcrdParams.coordIndex + 1,
677 value_ref * pull_conversion_factor_internal2userinput(pcrdParams));
680 else if (pcrdParams.eGeom == PullGroupGeometry::Dihedral)
682 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
683 make_periodic_2pi(&value_ref);
689 static void updatePullCoordReferenceValue(double* referenceValue, const t_pull_coord& pcrdParams, double t)
691 GMX_ASSERT(referenceValue, "Need a valid reference value object");
693 /* With zero rate the reference value is set initially and doesn't change */
694 if (pcrdParams.rate != 0)
696 const double inputValue = (pcrdParams.init + pcrdParams.rate * t)
697 * pull_conversion_factor_userinput2internal(pcrdParams);
698 *referenceValue = sanitizePullCoordReferenceValue(pcrdParams, inputValue);
702 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
703 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
706 dvec dr32; /* store instead of dr23? */
708 dsvmul(-1, spatialData->dr23, dr32);
709 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
710 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
711 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
713 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
714 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
715 * Note 2: the angle between the plane normal vectors equals pi only when
716 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
717 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
719 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
723 /* Calculates pull->coord[coord_ind].value.
724 * This function also updates pull->coord[coord_ind].dr.
726 static void get_pull_coord_distance(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc& pbc)
728 get_pull_coord_dr(pull, pcrd, pbc);
730 PullCoordSpatialData& spatialData = pcrd->spatialData;
732 switch (pcrd->params.eGeom)
734 case PullGroupGeometry::Distance:
735 /* Pull along the vector between the com's */
736 spatialData.value = dnorm(spatialData.dr01);
738 case PullGroupGeometry::Direction:
739 case PullGroupGeometry::DirectionPBC:
740 case PullGroupGeometry::DirectionRelative:
741 case PullGroupGeometry::Cylinder:
743 spatialData.value = 0;
744 for (int m = 0; m < DIM; m++)
746 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
749 case PullGroupGeometry::Angle:
750 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
752 case PullGroupGeometry::Dihedral:
753 spatialData.value = get_dihedral_angle_coord(&spatialData);
755 case PullGroupGeometry::AngleAxis:
756 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
758 case PullGroupGeometry::Transformation:
759 // Note that we would only need to pass the part of coord up to coord_ind
760 spatialData.value = gmx::getTransformationPullCoordinateValue(
761 pcrd, ArrayRef<const pull_coord_work_t>(pull.coord).subArray(0, pcrd->params.coordIndex));
763 default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
767 /* Returns the deviation from the reference value.
768 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
770 static double get_pull_coord_deviation(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc& pbc, double t)
774 /* Update the reference value before computing the distance,
775 * since it is used in the distance computation with periodic pulling.
777 updatePullCoordReferenceValue(&pcrd->value_ref, pcrd->params, t);
779 get_pull_coord_distance(pull, pcrd, pbc);
781 get_pull_coord_distance(pull, pcrd, pbc);
783 /* Determine the deviation */
784 dev = pcrd->spatialData.value - pcrd->value_ref;
786 /* Check that values are allowed */
787 if (pcrd->params.eGeom == PullGroupGeometry::Distance && pcrd->spatialData.value == 0)
789 /* With no vector we can not determine the direction for the force,
790 * so we set the force to zero.
794 else if (pcrd->params.eGeom == PullGroupGeometry::Dihedral)
796 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
797 Thus, the unwrapped deviation is here in (-2pi, 2pi].
798 After making it periodic, the deviation will be in [-pi, pi). */
799 make_periodic_2pi(&dev);
805 double get_pull_coord_value(pull_t* pull, int coordIndex, const t_pbc& pbc)
807 get_pull_coord_distance(*pull, &pull->coord[coordIndex], pbc);
809 return pull->coord[coordIndex].spatialData.value;
812 void clear_pull_forces(pull_t* pull)
814 /* Zeroing the forces is only required for constraint pulling.
815 * It can happen that multiple constraint steps need to be applied
816 * and therefore the constraint forces need to be accumulated.
818 for (pull_coord_work_t& coord : pull->coord)
820 coord.scalarForce = 0;
824 /* Apply constraint using SHAKE */
825 static void do_constraint(struct pull_t* pull,
835 dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
836 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
837 dvec* rnew; /* current 'new' positions of the groups */
838 double* dr_tot; /* the total update of the coords */
841 double lambda, rm, invdt = 0;
842 gmx_bool bConverged_all, bConverged = FALSE;
843 int niter = 0, ii, j, m, max_iter = 100;
847 snew(r_ij, pull->coord.size());
848 snew(dr_tot, pull->coord.size());
850 snew(rnew, pull->group.size());
852 /* copy the current unconstrained positions for use in iterations. We
853 iterate until rinew[i] and rjnew[j] obey the constraints. Then
854 rinew - pull.x_unc[i] is the correction dr to group i */
855 for (size_t g = 0; g < pull->group.size(); g++)
857 copy_dvec(pull->group[g].xp, rnew[g]);
860 /* Determine the constraint directions from the old positions */
861 for (size_t c = 0; c < pull->coord.size(); c++)
863 pull_coord_work_t* pcrd;
865 pcrd = &pull->coord[c];
867 if (pcrd->params.eType != PullingAlgorithm::Constraint)
872 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
873 * We don't modify dr and value anymore, so these values are also used
876 get_pull_coord_distance(*pull, pcrd, pbc);
878 const PullCoordSpatialData& spatialData = pcrd->spatialData;
882 "Pull coord %zu dr %f %f %f\n",
884 spatialData.dr01[XX],
885 spatialData.dr01[YY],
886 spatialData.dr01[ZZ]);
889 if (pcrd->params.eGeom == PullGroupGeometry::Direction
890 || pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
892 /* Select the component along vec */
894 for (m = 0; m < DIM; m++)
896 a += spatialData.vec[m] * spatialData.dr01[m];
898 for (m = 0; m < DIM; m++)
900 r_ij[c][m] = a * spatialData.vec[m];
905 copy_dvec(spatialData.dr01, r_ij[c]);
908 if (dnorm2(r_ij[c]) == 0)
911 "Distance for pull coordinate %zu is zero with constraint pulling, which is "
917 bConverged_all = FALSE;
918 while (!bConverged_all && niter < max_iter)
920 bConverged_all = TRUE;
922 /* loop over all constraints */
923 for (size_t c = 0; c < pull->coord.size(); c++)
925 pull_coord_work_t* pcrd;
926 pull_group_work_t *pgrp0, *pgrp1;
929 pcrd = &pull->coord[c];
931 if (pcrd->params.eType != PullingAlgorithm::Constraint)
936 updatePullCoordReferenceValue(&pcrd->value_ref, pcrd->params, t);
938 pgrp0 = &pull->group[pcrd->params.group[0]];
939 pgrp1 = &pull->group[pcrd->params.group[1]];
941 /* Get the current difference vector */
942 low_get_pull_coord_dr(
943 *pull, *pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
947 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
950 rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
952 switch (pcrd->params.eGeom)
954 case PullGroupGeometry::Distance:
955 if (pcrd->value_ref <= 0)
959 "The pull constraint reference distance for group %zu is <= 0 (%f)",
965 double q, c_a, c_b, c_c;
967 c_a = diprod(r_ij[c], r_ij[c]);
968 c_b = diprod(unc_ij, r_ij[c]) * 2;
969 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
973 q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
978 q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
984 fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b, c_c, lambda);
988 /* The position corrections dr due to the constraints */
989 dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
990 dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
991 dr_tot[c] += -lambda * dnorm(r_ij[c]);
993 case PullGroupGeometry::Direction:
994 case PullGroupGeometry::DirectionPBC:
995 case PullGroupGeometry::Cylinder:
996 /* A 1-dimensional constraint along a vector */
998 for (m = 0; m < DIM; m++)
1000 vec[m] = pcrd->spatialData.vec[m];
1001 a += unc_ij[m] * vec[m];
1003 /* Select only the component along the vector */
1004 dsvmul(a, vec, unc_ij);
1005 lambda = a - pcrd->value_ref;
1008 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
1011 /* The position corrections dr due to the constraints */
1012 dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
1013 dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
1014 dr_tot[c] += -lambda;
1016 case PullGroupGeometry::Transformation:
1017 GMX_RELEASE_ASSERT(false, "transformation with constraints should never occur");
1019 default: gmx_incons("Invalid enumeration value for eGeom");
1027 g0 = pcrd->params.group[0];
1028 g1 = pcrd->params.group[1];
1029 low_get_pull_coord_dr(*pull, *pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
1030 low_get_pull_coord_dr(*pull, *pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
1032 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1041 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
1050 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1060 /* Update the COMs with dr */
1061 dvec_inc(rnew[pcrd->params.group[1]], dr1);
1062 dvec_inc(rnew[pcrd->params.group[0]], dr0);
1065 /* Check if all constraints are fullfilled now */
1066 for (pull_coord_work_t& coord : pull->coord)
1068 if (coord.params.eType != PullingAlgorithm::Constraint)
1073 low_get_pull_coord_dr(
1074 *pull, coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
1076 switch (coord.params.eGeom)
1078 case PullGroupGeometry::Distance:
1079 bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1081 case PullGroupGeometry::Direction:
1082 case PullGroupGeometry::DirectionPBC:
1083 case PullGroupGeometry::Cylinder:
1084 for (m = 0; m < DIM; m++)
1086 vec[m] = coord.spatialData.vec[m];
1088 inpr = diprod(unc_ij, vec);
1089 dsvmul(inpr, vec, unc_ij);
1090 bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1094 gmx::formatString("Geometry %s not handled here",
1095 enumValueToString(coord.params.eGeom))
1104 "Pull constraint not converged: "
1106 "d_ref = %f, current d = %f\n",
1107 coord.params.group[0],
1108 coord.params.group[1],
1113 bConverged_all = FALSE;
1118 /* if after all constraints are dealt with and bConverged is still TRUE
1119 we're finished, if not we do another iteration */
1121 if (niter > max_iter)
1123 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1126 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1133 /* update atoms in the groups */
1134 for (size_t g = 0; g < pull->group.size(); g++)
1136 const pull_group_work_t* pgrp;
1139 pgrp = &pull->group[g];
1141 /* get the final constraint displacement dr for group g */
1142 dvec_sub(rnew[g], pgrp->xp, dr);
1144 if (dnorm2(dr) == 0)
1146 /* No displacement, no update necessary */
1150 /* update the atom positions */
1151 auto localAtomIndices = pgrp->atomSet.localIndex();
1153 for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1155 ii = localAtomIndices[j];
1156 if (!pgrp->localWeights.empty())
1158 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1160 for (m = 0; m < DIM; m++)
1166 for (m = 0; m < DIM; m++)
1168 v[ii][m] += invdt * tmp[m];
1174 /* calculate the constraint forces, used for output and virial only */
1175 for (size_t c = 0; c < pull->coord.size(); c++)
1177 pull_coord_work_t* pcrd;
1179 pcrd = &pull->coord[c];
1181 if (pcrd->params.eType != PullingAlgorithm::Constraint)
1186 /* Accumulate the forces, in case we have multiple constraint steps */
1189 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1191 pcrd->scalarForce += force;
1193 if (vir != nullptr && pcrd->params.eGeom != PullGroupGeometry::DirectionPBC && bMaster)
1197 /* Add the pull contribution to the virial */
1198 /* We have already checked above that r_ij[c] != 0 */
1199 f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1201 for (j = 0; j < DIM; j++)
1203 for (m = 0; m < DIM; m++)
1205 vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1211 /* finished! I hope. Give back some memory */
1217 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1219 for (int j = 0; j < DIM; j++)
1221 for (int m = 0; m < DIM; m++)
1223 vir[j][m] -= 0.5 * f[j] * dr[m];
1228 /* Adds the pull contribution to the virial */
1229 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1231 if (vir != nullptr && pcrd.params.eGeom != PullGroupGeometry::DirectionPBC)
1233 /* Add the pull contribution for each distance vector to the virial. */
1234 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1235 if (pcrd.params.ngroup >= 4)
1237 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1239 if (pcrd.params.ngroup >= 6)
1241 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1246 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1254 k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1255 dkdl = pcrd->params.kB - pcrd->params.k;
1257 switch (pcrd->params.eType)
1259 case PullingAlgorithm::Umbrella:
1260 case PullingAlgorithm::FlatBottom:
1261 case PullingAlgorithm::FlatBottomHigh:
1262 /* The only difference between an umbrella and a flat-bottom
1263 * potential is that a flat-bottom is zero above or below
1264 the reference value.
1266 if ((pcrd->params.eType == PullingAlgorithm::FlatBottom && dev < 0)
1267 || (pcrd->params.eType == PullingAlgorithm::FlatBottomHigh && dev > 0))
1272 pcrd->scalarForce = -k * dev;
1273 *V += 0.5 * k * gmx::square(dev);
1274 *dVdl += 0.5 * dkdl * gmx::square(dev);
1276 case PullingAlgorithm::ConstantForce:
1277 pcrd->scalarForce = -k;
1278 *V += k * pcrd->spatialData.value;
1279 *dVdl += dkdl * pcrd->spatialData.value;
1281 case PullingAlgorithm::External:
1283 "the scalar pull force should not be calculated internally for pull type "
1285 default: gmx_incons("Unsupported pull type in do_pull_pot");
1289 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1291 const t_pull_coord& params = pcrd.params;
1292 const PullCoordSpatialData& spatialData = pcrd.spatialData;
1294 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1295 PullCoordVectorForces forces;
1297 if (params.eGeom == PullGroupGeometry::Distance)
1299 double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1300 for (int m = 0; m < DIM; m++)
1302 forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1305 else if (params.eGeom == PullGroupGeometry::Angle)
1308 double cos_theta, cos_theta2;
1310 cos_theta = cos(spatialData.value);
1311 cos_theta2 = gmx::square(cos_theta);
1313 /* The force at theta = 0, pi is undefined so we should not apply any force.
1314 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1315 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1316 * have to check for this before dividing by their norm below.
1320 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1321 * and another vector parallel to dr23:
1322 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1323 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1325 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1326 double b = a * cos_theta;
1327 double invdr01 = 1. / dnorm(spatialData.dr01);
1328 double invdr23 = 1. / dnorm(spatialData.dr23);
1329 dvec normalized_dr01, normalized_dr23;
1330 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1331 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1333 for (int m = 0; m < DIM; m++)
1335 /* Here, f_scal is -dV/dtheta */
1337 pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1339 pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1344 /* No forces to apply for ill-defined cases*/
1345 clear_dvec(forces.force01);
1346 clear_dvec(forces.force23);
1349 else if (params.eGeom == PullGroupGeometry::AngleAxis)
1351 double cos_theta, cos_theta2;
1353 /* The angle-axis force is exactly the same as the angle force (above) except that in
1354 this case the second vector (dr23) is replaced by the pull vector. */
1355 cos_theta = cos(spatialData.value);
1356 cos_theta2 = gmx::square(cos_theta);
1362 dvec normalized_dr01;
1364 invdr01 = 1. / dnorm(spatialData.dr01);
1365 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1366 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1369 for (int m = 0; m < DIM; m++)
1372 pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1377 clear_dvec(forces.force01);
1380 else if (params.eGeom == PullGroupGeometry::Dihedral)
1382 double m2, n2, tol, sqrdist_32;
1384 /* Note: there is a small difference here compared to the
1385 dihedral force calculations in the bondeds (ref: Bekker 1994).
1386 There rij = ri - rj, while here dr01 = r1 - r0.
1387 However, all distance vectors occur in form of cross or inner products
1388 so that two signs cancel and we end up with the same expressions.
1389 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1391 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1392 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1393 dsvmul(-1, spatialData.dr23, dr32);
1394 sqrdist_32 = diprod(dr32, dr32);
1395 tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1396 if ((m2 > tol) && (n2 > tol))
1398 double a_01, a_23_01, a_23_45, a_45;
1399 double inv_dist_32, inv_sqrdist_32, dist_32;
1401 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1402 inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1403 dist_32 = sqrdist_32 * inv_dist_32;
1405 /* Forces on groups 0, 1 */
1406 a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
1407 dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1409 /* Forces on groups 4, 5 */
1410 a_45 = -pcrd.scalarForce * dist_32 / n2;
1411 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1413 /* Force on groups 2, 3 (defining the axis) */
1414 a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1415 a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1416 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1417 dsvmul(a_23_45, forces.force45, v);
1418 dvec_sub(u, v, forces.force23); /* force on group 3 */
1422 /* No force to apply for ill-defined cases */
1423 clear_dvec(forces.force01);
1424 clear_dvec(forces.force23);
1425 clear_dvec(forces.force45);
1430 for (int m = 0; m < DIM; m++)
1432 forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1440 /* We use a global mutex for locking access to the pull data structure
1441 * during registration of external pull potential providers.
1442 * We could use a different, local mutex for each pull object, but the overhead
1443 * is extremely small here and registration is only done during initialization.
1445 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
1446 static std::mutex registrationMutex;
1448 using Lock = std::lock_guard<std::mutex>;
1450 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1452 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1453 GMX_RELEASE_ASSERT(provider != nullptr,
1454 "register_external_pull_potential called with NULL as provider name");
1456 if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1459 "Module '%s' attempted to register an external potential for pull coordinate %d "
1460 "which is out of the pull coordinate range %d - %zu\n",
1464 pull->coord.size());
1467 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1469 if (pcrd->params.eType != PullingAlgorithm::External)
1473 "Module '%s' attempted to register an external potential for pull coordinate %d "
1474 "which of type '%s', whereas external potentials are only supported with type '%s'",
1477 enumValueToString(pcrd->params.eType),
1478 enumValueToString(PullingAlgorithm::External));
1481 GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
1482 "The external potential provider string for a pull coordinate is NULL");
1484 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
1487 "Module '%s' attempted to register an external potential for pull coordinate %d "
1488 "which expects the external potential to be provided by a module named '%s'",
1491 pcrd->params.externalPotentialProvider.c_str());
1494 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1495 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1496 * pull->numUnregisteredExternalPotentials.
1498 Lock registrationLock(registrationMutex);
1500 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1503 "Module '%s' attempted to register an external potential for pull coordinate %d "
1509 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1510 pull->numUnregisteredExternalPotentials--;
1512 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1513 "Negative unregistered potentials, the pull code is inconsistent");
1517 static void check_external_potential_registration(const struct pull_t* pull)
1519 if (pull->numUnregisteredExternalPotentials > 0)
1521 for (const pull_coord_work_t& pcrd : pull->coord)
1523 if (pcrd.params.eType == PullingAlgorithm::External
1524 && !pcrd.bExternalPotentialProviderHasBeenRegistered)
1527 "No external provider for external pull potentials have been provided "
1529 "pull coordinates. The first coordinate without provider is number %d, "
1531 "expects a module named '%s' to provide the external potential.",
1532 pull->numUnregisteredExternalPotentials,
1533 pcrd.params.coordIndex + 1,
1534 pcrd.params.externalPotentialProvider.c_str());
1540 /*! \brief Applies a force of any non-transformation pull coordinate
1542 * \param[in] pull The pull struct, we can't pass only the coord because of cylinder pulling
1543 * \param[in] coord_index The index of the coord to apply to force for
1544 * \param[in] coord_force The force working on this coord
1545 * \param[in] masses Atom masses
1546 * \param[in] forceWithVirial Atom force and virial object
1548 static void apply_default_pull_coord_force(pull_t* pull,
1549 const int coord_index,
1550 const double coord_force,
1551 gmx::ArrayRef<const real> masses,
1552 gmx::ForceWithVirial* forceWithVirial)
1554 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1556 pcrd->scalarForce = coord_force;
1557 /* Calculate the forces on the pull groups */
1558 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1560 /* Add the forces for this coordinate to the total virial and force */
1561 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1563 matrix virial = { { 0 } };
1564 add_virial_coord(virial, *pcrd, pullCoordForces);
1565 forceWithVirial->addVirialContribution(virial);
1567 apply_forces_coord(pull->coord[coord_index],
1571 as_rvec_array(forceWithVirial->force_.data()));
1574 /* Pull takes care of adding the forces of the external potential.
1575 * The external potential module has to make sure that the corresponding
1576 * potential energy is added either to the pull term or to a term
1577 * specific to the external module.
1579 void apply_external_pull_coord_force(struct pull_t* pull,
1582 ArrayRef<const real> masses,
1583 gmx::ForceWithVirial* forceWithVirial)
1585 pull_coord_work_t* pcrd;
1587 GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1588 "apply_external_pull_coord_force called with coord_index out of range");
1590 if (pull->comm.bParticipate)
1592 pcrd = &pull->coord[coord_index];
1595 pcrd->params.eType == PullingAlgorithm::External,
1596 "The pull force can only be set externally on pull coordinates of external type");
1598 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1599 "apply_external_pull_coord_force called for an unregistered pull coordinate");
1601 if (pcrd->params.eGeom == PullGroupGeometry::Transformation)
1603 gmx::applyTransformationPullCoordForce(
1605 gmx::ArrayRef<pull_coord_work_t>(pull->coord).subArray(0, pcrd->params.coordIndex),
1610 apply_default_pull_coord_force(pull, coord_index, coord_force, masses, forceWithVirial);
1613 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1616 /* Calculate the pull potential and scalar force for a pull coordinate.
1617 * Returns the vector forces for the pull coordinate.
1619 static PullCoordVectorForces do_pull_pot_coord(const pull_t& pull,
1620 pull_coord_work_t* pcrd,
1628 assert(pcrd->params.eType != PullingAlgorithm::Constraint);
1630 double dev = get_pull_coord_deviation(pull, pcrd, pbc, t);
1632 calc_pull_coord_scalar_force_and_potential(pcrd, dev, lambda, V, dVdl);
1634 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1636 add_virial_coord(vir, *pcrd, pullCoordForces);
1638 return pullCoordForces;
1641 real pull_potential(struct pull_t* pull,
1642 ArrayRef<const real> masses,
1644 const t_commrec* cr,
1647 ArrayRef<const RVec> x,
1648 gmx::ForceWithVirial* force,
1653 assert(pull != nullptr);
1655 /* Ideally we should check external potential registration only during
1656 * the initialization phase, but that requires another function call
1657 * that should be done exactly in the right order. So we check here.
1659 check_external_potential_registration(pull);
1661 if (pull->comm.bParticipate)
1665 pull_calc_coms(cr, pull, masses, pbc, t, x, {});
1667 rvec* f = as_rvec_array(force->force_.data());
1668 matrix virial = { { 0 } };
1669 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1670 for (pull_coord_work_t& pcrd : pull->coord)
1672 /* For external potential the force is assumed to be given by an external module by a
1673 call to apply_pull_coord_external_force */
1674 if (pcrd.params.eType == PullingAlgorithm::Constraint
1675 || pcrd.params.eType == PullingAlgorithm::External)
1679 PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1680 *pull, &pcrd, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1682 if (pcrd.params.eGeom == PullGroupGeometry::Transformation)
1684 gmx::applyTransformationPullCoordForce(
1686 gmx::ArrayRef<pull_coord_work_t>(pull->coord).subArray(0, pcrd.params.coordIndex),
1691 /* Distribute the force over the atoms in the pulled groups */
1692 apply_forces_coord(pcrd, pull->group, pullCoordForces, masses, f);
1698 force->addVirialContribution(virial);
1703 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1704 "Too few or too many external pull potentials have been applied the previous step");
1705 /* All external pull potentials still need to be applied */
1706 pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1708 return (MASTER(cr) ? V : 0.0);
1711 void pull_constraint(struct pull_t* pull,
1712 ArrayRef<const real> masses,
1714 const t_commrec* cr,
1722 assert(pull != nullptr);
1724 if (pull->comm.bParticipate)
1726 pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
1728 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1732 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1736 gmx_bool bMustParticipate;
1742 /* We always make the master node participate, such that it can do i/o,
1743 * add the virial and to simplify MC type extensions people might have.
1745 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1747 for (pull_group_work_t& group : pull->group)
1749 if (!group.globalWeights.empty())
1751 group.localWeights.resize(group.atomSet.numAtomsLocal());
1752 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1754 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1758 GMX_ASSERT(bMustParticipate || dd != nullptr,
1759 "Either all ranks (including this rank) participate, or we use DD and need to "
1760 "have access to dd here");
1762 /* We should participate if we have pull or pbc atoms */
1763 if (!bMustParticipate
1764 && (group.atomSet.numAtomsLocal() > 0
1765 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1767 bMustParticipate = TRUE;
1771 if (!comm->bParticipateAll)
1773 /* Keep currently not required ranks in the communicator
1774 * if they needed to participate up to 20 decompositions ago.
1775 * This avoids frequent rebuilds due to atoms jumping back and forth.
1777 const int64_t history_count = 20;
1778 gmx_bool bWillParticipate;
1781 /* Increase the decomposition counter for the current call */
1782 comm->setup_count++;
1784 if (bMustParticipate)
1786 comm->must_count = comm->setup_count;
1791 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1793 if (debug && dd != nullptr)
1796 "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
1798 gmx::boolToString(bMustParticipate),
1799 gmx::boolToString(bWillParticipate));
1802 if (bWillParticipate)
1804 /* Count the number of ranks that we want to have participating */
1806 /* Count the number of ranks that need to be added */
1807 count[1] = comm->bParticipate ? 0 : 1;
1815 /* The cost of this global operation will be less that the cost
1816 * of the extra MPI_Comm_split calls that we can avoid.
1818 gmx_sumi(2, count, cr);
1820 /* If we are missing ranks or if we have 20% more ranks than needed
1821 * we make a new sub-communicator.
1823 if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1827 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1831 if (comm->mpi_comm_com != MPI_COMM_NULL)
1833 MPI_Comm_free(&comm->mpi_comm_com);
1836 /* This might be an extremely expensive operation, so we try
1837 * to avoid this splitting as much as possible.
1839 assert(dd != nullptr);
1840 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1843 comm->bParticipate = bWillParticipate;
1844 comm->nparticipate = count[0];
1846 /* When we use the previous COM for PBC, we need to broadcast
1847 * the previous COM to ranks that have joined the communicator.
1849 for (pull_group_work_t& group : pull->group)
1851 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1853 GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1854 "The master rank has to participate, as it should pass an up to "
1856 "to bcast here as well as to e.g. checkpointing");
1858 gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
1864 /* Since the PBC of atoms might have changed, we need to update the PBC */
1865 pull->bSetPBCatoms = TRUE;
1868 static void init_pull_group_index(FILE* fplog,
1869 const t_commrec* cr,
1871 pull_group_work_t* pg,
1872 gmx_bool bConstraint,
1873 const ivec pulldim_con,
1874 const gmx_mtop_t& mtop,
1875 const t_inputrec* ir,
1878 /* With EM and BD there are no masses in the integrator.
1879 * But we still want to have the correct mass-weighted COMs.
1880 * So we store the real masses in the weights.
1882 const bool setWeights = (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1883 || ir->eI == IntegrationAlgorithm::BD);
1885 /* In parallel, store we need to extract localWeights from weights at DD time */
1886 std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1888 const SimulationGroups& groups = mtop.groups;
1890 /* Count frozen dimensions and (weighted) mass */
1896 for (int i = 0; i < int(pg->params.ind.size()); i++)
1898 int ii = pg->params.ind[i];
1899 if (bConstraint && ir->opts.nFreeze)
1901 for (int d = 0; d < DIM; d++)
1903 if (pulldim_con[d] == 1
1904 && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1910 const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1912 if (ir->efep == FreeEnergyPerturbationType::No)
1918 m = (1 - lambda) * atom.m + lambda * atom.mB;
1921 if (!pg->params.weight.empty())
1923 w = pg->params.weight[i];
1929 if (EI_ENERGY_MINIMIZATION(ir->eI))
1931 /* Move the mass to the weight */
1935 else if (ir->eI == IntegrationAlgorithm::BD)
1938 if (ir->bd_fric != 0.0F)
1940 mbd = ir->bd_fric * ir->delta_t;
1944 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1946 mbd = ir->delta_t / ir->opts.tau_t[0];
1951 / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1959 weights.push_back(w);
1963 wwmass += m * w * w;
1968 /* We can have single atom groups with zero mass with potential pulling
1969 * without cosine weighting.
1971 if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1973 /* With one atom the mass doesn't matter */
1979 "The total%s mass of pull group %d is zero",
1980 !pg->params.weight.empty() ? " weighted" : "",
1986 fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
1987 if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1988 || ir->eI == IntegrationAlgorithm::BD)
1990 fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1992 if (pg->epgrppbc == epgrppbcCOS)
1994 fprintf(fplog, ", cosine weighting will be used");
1996 fprintf(fplog, "\n");
2001 /* A value != 0 signals not frozen, it is updated later */
2007 for (int d = 0; d < DIM; d++)
2009 ndim += pulldim_con[d] * pg->params.ind.size();
2011 if (fplog && nfrozen > 0 && nfrozen < ndim)
2014 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
2015 " that are subject to pulling are frozen.\n"
2016 " For constraint pulling the whole group will be frozen.\n\n",
2024 struct pull_t* init_pull(FILE* fplog,
2025 const pull_params_t* pull_params,
2026 const t_inputrec* ir,
2027 const gmx_mtop_t& mtop,
2028 const t_commrec* cr,
2029 gmx::LocalAtomSetManager* atomSets,
2032 struct pull_t* pull;
2037 /* Copy the pull parameters */
2038 pull->params = *pull_params;
2040 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2041 const int maxNumThreads = std::max(1, gmx_omp_nthreads_get(ModuleMultiThread::Default));
2043 for (int i = 0; i < pull_params->ngroup; ++i)
2045 pull->group.emplace_back(pull_params->group[i],
2046 atomSets->add(pull_params->group[i].ind),
2047 pull_params->bSetPbcRefToPrevStepCOM,
2051 if (cr != nullptr && DOMAINDECOMP(cr))
2053 /* Set up the global to local atom mapping for PBC atoms */
2054 for (pull_group_work_t& group : pull->group)
2056 if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
2058 /* pbcAtomSet consists of a single atom */
2059 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
2060 atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
2065 pull->bPotential = FALSE;
2066 pull->bConstraint = FALSE;
2067 pull->bCylinder = FALSE;
2068 pull->bAngle = FALSE;
2069 pull->bXOutAverage = pull_params->bXOutAverage;
2070 pull->bFOutAverage = pull_params->bFOutAverage;
2072 GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
2073 "pull group 0 is an absolute reference group and should not contain atoms");
2075 pull->numCoordinatesWithExternalPotential = 0;
2077 for (int c = 0; c < pull->params.ncoord; c++)
2079 GMX_RELEASE_ASSERT(pull_params->coord[c].coordIndex == c,
2080 "The stored index should match the position in the vector");
2082 /* Construct a pull coordinate, copying all coordinate parameters */
2083 pull->coord.emplace_back(pull_params->coord[c]);
2085 pull_coord_work_t* pcrd = &pull->coord.back();
2087 switch (pcrd->params.eGeom)
2089 case PullGroupGeometry::Distance:
2090 case PullGroupGeometry::DirectionRelative: /* Direction vector is determined at each step */
2091 case PullGroupGeometry::Angle:
2092 case PullGroupGeometry::Dihedral: break;
2093 case PullGroupGeometry::Direction:
2094 case PullGroupGeometry::DirectionPBC:
2095 case PullGroupGeometry::Cylinder:
2096 case PullGroupGeometry::Transformation:
2097 case PullGroupGeometry::AngleAxis:
2098 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2101 /* We allow reading of newer tpx files with new pull geometries,
2102 * but with the same tpx format, with old code. A new geometry
2103 * only adds a new enum value, which will be out of range for
2104 * old code. The first place we need to generate an error is
2105 * here, since the pull code can't handle this.
2106 * The enum can be used before arriving here only for printing
2107 * the string corresponding to the geometry, which will result
2108 * in printing "UNDEFINED".
2111 "Pull geometry not supported for pull coordinate %d. The geometry enum "
2112 "%s in the input is larger than that supported by the code (up to %d). "
2113 "You are probably reading a tpr file generated with a newer version of "
2114 "GROMACS with an binary from an older version of Gromacs.",
2116 enumValueToString(pcrd->params.eGeom),
2117 static_cast<int>(PullGroupGeometry::Count) - 1);
2120 if (pcrd->params.eType == PullingAlgorithm::Constraint)
2122 /* Check restrictions of the constraint pull code */
2123 if (pcrd->params.eGeom == PullGroupGeometry::Cylinder
2124 || pcrd->params.eGeom == PullGroupGeometry::DirectionRelative
2125 || pcrd->params.eGeom == PullGroupGeometry::Angle
2126 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2127 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2130 "Pulling of type %s can not be combined with geometry %s. Consider using "
2132 enumValueToString(pcrd->params.eType),
2133 enumValueToString(pcrd->params.eGeom),
2134 enumValueToString(PullingAlgorithm::Umbrella));
2139 "Constraint pulling can not be combined with multiple time stepping");
2141 pull->bConstraint = TRUE;
2145 pull->bPotential = TRUE;
2148 if (pcrd->params.eGeom == PullGroupGeometry::Cylinder)
2150 pull->bCylinder = TRUE;
2152 else if (pcrd->params.eGeom == PullGroupGeometry::Angle
2153 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2154 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2156 pull->bAngle = TRUE;
2159 /* We only need to calculate the plain COM of a group
2160 * when it is not only used as a cylinder group.
2161 * Also the absolute reference group 0 needs no COM computation.
2163 for (int i = 0; i < pcrd->params.ngroup; i++)
2165 int groupIndex = pcrd->params.group[i];
2166 if (groupIndex > 0 && !(pcrd->params.eGeom == PullGroupGeometry::Cylinder && i == 0))
2168 pull->group[groupIndex].needToCalcCom = true;
2172 /* With non-zero rate the reference value is set at every step */
2173 if (pcrd->params.rate == 0)
2175 /* Initialize the constant reference value */
2176 if (pcrd->params.eType != PullingAlgorithm::External)
2178 pcrd->value_ref = sanitizePullCoordReferenceValue(
2180 pcrd->params.init * pull_conversion_factor_userinput2internal(pcrd->params));
2184 /* With an external pull potential, the reference value loses
2185 * it's meaning and should not be used. Setting it to zero
2186 * makes any terms dependent on it disappear.
2187 * The only issue this causes is that with cylinder pulling
2188 * no atoms of the cylinder group within the cylinder radius
2189 * should be more than half a box length away from the COM of
2190 * the pull group along the axial direction.
2192 pcrd->value_ref = 0.0;
2196 if (pcrd->params.eType == PullingAlgorithm::External)
2199 pcrd->params.rate == 0,
2200 "With an external potential, a pull coordinate should have rate = 0");
2202 /* This external potential needs to be registered later */
2203 pull->numCoordinatesWithExternalPotential++;
2205 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2208 pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
2209 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2211 pull->pbcType = ir->pbcType;
2212 switch (pull->pbcType)
2214 case PbcType::No: pull->npbcdim = 0; break;
2215 case PbcType::XY: pull->npbcdim = 2; break;
2216 default: pull->npbcdim = 3; break;
2221 gmx_bool bAbs, bCos;
2224 for (const pull_coord_work_t& coord : pull->coord)
2226 if (pull->group[coord.params.group[0]].params.ind.empty()
2227 || pull->group[coord.params.group[1]].params.ind.empty())
2233 fprintf(fplog, "\n");
2234 if (pull->bPotential)
2236 fprintf(fplog, "Will apply potential COM pulling\n");
2238 if (pull->bConstraint)
2240 fprintf(fplog, "Will apply constraint COM pulling\n");
2242 // Don't include the reference group 0 in output, so we report ngroup-1
2243 int numRealGroups = pull->group.size() - 1;
2244 GMX_RELEASE_ASSERT(numRealGroups > 0,
2245 "The reference absolute position pull group should always be present");
2247 "with %zu pull coordinate%s and %d group%s\n",
2249 pull->coord.size() == 1 ? "" : "s",
2251 numRealGroups == 1 ? "" : "s");
2254 fprintf(fplog, "with an absolute reference\n");
2257 // Don't include the reference group 0 in loop
2258 for (size_t g = 1; g < pull->group.size(); g++)
2260 if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
2262 /* We are using cosine weighting */
2263 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2269 please_cite(fplog, "Engin2010");
2273 pull->bRefAt = FALSE;
2275 for (size_t g = 0; g < pull->group.size(); g++)
2277 pull_group_work_t* pgrp;
2279 pgrp = &pull->group[g];
2280 if (!pgrp->params.ind.empty())
2282 /* There is an issue when a group is used in multiple coordinates
2283 * and constraints are applied in different dimensions with atoms
2284 * frozen in some, but not all dimensions.
2285 * Since there is only one mass array per group, we can't have
2286 * frozen/non-frozen atoms for different coords at the same time.
2287 * But since this is a very exotic case, we don't check for this.
2288 * A warning is printed in init_pull_group_index.
2291 gmx_bool bConstraint;
2292 ivec pulldim, pulldim_con;
2294 /* Loop over all pull coordinates to see along which dimensions
2295 * this group is pulled and if it is involved in constraints.
2297 bConstraint = FALSE;
2298 clear_ivec(pulldim);
2299 clear_ivec(pulldim_con);
2300 for (const pull_coord_work_t& coord : pull->coord)
2302 gmx_bool bGroupUsed = FALSE;
2303 for (int gi = 0; gi < coord.params.ngroup; gi++)
2305 if (coord.params.group[gi] == static_cast<int>(g))
2313 for (int m = 0; m < DIM; m++)
2315 if (coord.params.dim[m] == 1)
2319 if (coord.params.eType == PullingAlgorithm::Constraint)
2329 /* Determine if we need to take PBC into account for calculating
2330 * the COM's of the pull groups.
2332 switch (pgrp->epgrppbc)
2334 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2336 if (!pgrp->params.weight.empty())
2339 "Pull groups can not have relative weights and cosine weighting "
2342 for (int m = 0; m < DIM; m++)
2344 if (m < pull->npbcdim && pulldim[m] == 1)
2346 if (pull->cosdim >= 0 && pull->cosdim != m)
2349 "Can only use cosine weighting with pulling in one "
2350 "dimension (use mdp option pull-coord?-dim)");
2356 case epgrppbcNONE: break;
2359 /* Set the indices */
2360 init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2364 /* Absolute reference, set the inverse mass to zero.
2365 * This is only relevant (and used) with constraint pulling.
2372 /* If we use cylinder coordinates, do some initialising for them */
2373 if (pull->bCylinder)
2375 for (pull_coord_work_t& coord : pull->coord)
2377 if (coord.params.eGeom == PullGroupGeometry::Cylinder)
2379 if (pull->group[coord.params.group[0]].params.ind.empty())
2382 "A cylinder pull group is not supported when using absolute "
2386 const auto& group0 = pull->group[coord.params.group[0]];
2387 coord.dynamicGroup0 = std::make_unique<pull_group_work_t>(
2388 group0.params, group0.atomSet, pull->params.bSetPbcRefToPrevStepCOM, maxNumThreads);
2392 pull->comSums.resize(maxNumThreads);
2397 /* Use a sub-communicator when we have more than 32 ranks, but not
2398 * when we have an external pull potential, since then the external
2399 * potential provider expects each rank to have the coordinate.
2401 comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
2402 || pull->numCoordinatesWithExternalPotential > 0
2403 || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2404 /* This sub-commicator is not used with comm->bParticipateAll,
2405 * so we can always initialize it to NULL.
2407 comm->mpi_comm_com = MPI_COMM_NULL;
2408 comm->nparticipate = 0;
2409 comm->isMasterRank = (cr == nullptr || MASTER(cr));
2411 /* No MPI: 1 rank: all ranks pull */
2412 comm->bParticipateAll = TRUE;
2413 comm->isMasterRank = true;
2415 comm->bParticipate = comm->bParticipateAll;
2416 comm->setup_count = 0;
2417 comm->must_count = 0;
2419 if (!comm->bParticipateAll && fplog != nullptr)
2421 fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2424 comm->pbcAtomBuffer.resize(pull->group.size());
2425 comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
2426 if (pull->bCylinder)
2428 comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
2431 /* We still need to initialize the PBC reference coordinates */
2432 pull->bSetPBCatoms = TRUE;
2434 pull->out_x = nullptr;
2435 pull->out_f = nullptr;
2440 static void destroy_pull(struct pull_t* pull)
2443 if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2445 MPI_Comm_free(&pull->comm.mpi_comm_com);
2452 void preparePrevStepPullCom(const t_inputrec* ir,
2454 ArrayRef<const real> masses,
2456 const t_state* state_global,
2457 const t_commrec* cr,
2458 bool startingFromCheckpoint)
2460 if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
2464 allocStatePrevStepPullCom(state, pull_work);
2465 if (startingFromCheckpoint)
2469 state->pull_com_prev_step = state_global->pull_com_prev_step;
2473 /* Only the master rank has the checkpointed COM from the previous step */
2474 gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(),
2475 &state->pull_com_prev_step[0],
2476 cr->mpi_comm_mygroup);
2478 setPrevStepPullComFromState(pull_work, state);
2483 set_pbc(&pbc, ir->pbcType, state->box);
2484 initPullComFromPrevStep(
2485 cr, pull_work, masses, pbc, state->x.arrayRefWithPadding().unpaddedArrayRef());
2486 updatePrevStepPullCom(pull_work, state);
2490 void finish_pull(struct pull_t* pull)
2492 check_external_potential_registration(pull);
2496 gmx_fio_fclose(pull->out_x);
2500 gmx_fio_fclose(pull->out_f);
2506 bool pull_have_potential(const pull_t& pull)
2508 return pull.bPotential;
2511 bool pull_have_constraint(const pull_t& pull)
2513 return pull.bConstraint;
2516 bool pull_have_constraint(const pull_params_t& pullParameters)
2518 for (int c = 0; c < pullParameters.ncoord; c++)
2520 if (pullParameters.coord[c].eType == PullingAlgorithm::Constraint)