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.
40 #include "gromacs/utility/stringutil.h"
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/cstringutil.h"
75 #include "gromacs/utility/exceptions.h"
76 #include "gromacs/utility/fatalerror.h"
77 #include "gromacs/utility/futil.h"
78 #include "gromacs/utility/gmxassert.h"
79 #include "gromacs/utility/pleasecite.h"
80 #include "gromacs/utility/real.h"
81 #include "gromacs/utility/smalloc.h"
82 #include "gromacs/utility/strconvert.h"
84 #include "pull_internal.h"
88 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
91 static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
93 if (params.ind.size() <= 1)
98 else if (params.pbcatom >= 0)
100 if (setPbcRefToPrevStepCOM)
102 return epgrppbcPREVSTEPCOM;
106 return epgrppbcREFAT;
115 /* NOTE: The params initialization currently copies pointers.
116 * So the lifetime of the source, currently always inputrec,
117 * should not end before that of this object.
118 * This will be fixed when the pointers are replacted by std::vector.
120 pull_group_work_t::pull_group_work_t(const t_pull_group& params,
121 gmx::LocalAtomSet atomSet,
122 bool bSetPbcRefToPrevStepCOM) :
124 epgrppbc(groupPbcFromParams(params, bSetPbcRefToPrevStepCOM)),
125 needToCalcCom(false),
135 bool pull_coordinate_is_angletype(const t_pull_coord* pcrd)
137 return (pcrd->eGeom == PullGroupGeometry::Angle || pcrd->eGeom == PullGroupGeometry::Dihedral
138 || pcrd->eGeom == PullGroupGeometry::AngleAxis);
141 static bool pull_coordinate_is_directional(const t_pull_coord* pcrd)
143 return (pcrd->eGeom == PullGroupGeometry::Direction || pcrd->eGeom == PullGroupGeometry::DirectionPBC
144 || pcrd->eGeom == PullGroupGeometry::DirectionRelative
145 || pcrd->eGeom == PullGroupGeometry::Cylinder);
148 const char* pull_coordinate_units(const t_pull_coord* pcrd)
150 return pull_coordinate_is_angletype(pcrd) ? "deg" : "nm";
153 double pull_conversion_factor_userinput2internal(const t_pull_coord* pcrd)
155 if (pull_coordinate_is_angletype(pcrd))
157 return gmx::c_deg2Rad;
165 double pull_conversion_factor_internal2userinput(const t_pull_coord* pcrd)
167 if (pull_coordinate_is_angletype(pcrd))
169 return gmx::c_rad2Deg;
177 /* Apply forces in a mass weighted fashion for part of the pull group */
178 static void apply_forces_grp_part(const pull_group_work_t* pgrp,
186 double inv_wm = pgrp->mwscale;
188 auto localAtomIndices = pgrp->atomSet.localIndex();
189 for (int i = ind_start; i < ind_end; i++)
191 int ii = localAtomIndices[i];
192 double wmass = masses[ii];
193 if (!pgrp->localWeights.empty())
195 wmass *= pgrp->localWeights[i];
198 for (int d = 0; d < DIM; d++)
200 f[ii][d] += sign * wmass * f_pull[d] * inv_wm;
205 /* Apply forces in a mass weighted fashion */
206 static void apply_forces_grp(const pull_group_work_t* pgrp,
213 auto localAtomIndices = pgrp->atomSet.localIndex();
215 if (pgrp->params.ind.size() == 1 && pgrp->atomSet.numAtomsLocal() == 1)
217 /* Only one atom and our rank has this atom: we can skip
218 * the mass weighting, which means that this code also works
219 * for mass=0, e.g. with a virtual site.
221 for (int d = 0; d < DIM; d++)
223 f[localAtomIndices[0]][d] += sign * f_pull[d];
228 if (localAtomIndices.size() <= c_pullMaxNumLocalAtomsSingleThreaded)
230 apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), masses, f_pull, sign, f);
234 #pragma omp parallel for num_threads(nthreads) schedule(static)
235 for (int th = 0; th < nthreads; th++)
237 int ind_start = (localAtomIndices.size() * (th + 0)) / nthreads;
238 int ind_end = (localAtomIndices.size() * (th + 1)) / nthreads;
239 apply_forces_grp_part(pgrp, ind_start, ind_end, masses, f_pull, sign, f);
245 /* Apply forces in a mass weighted fashion to a cylinder group */
246 static void apply_forces_cyl_grp(const pull_group_work_t* pgrp,
247 const double dv_corr,
253 int gmx_unused nthreads)
255 double inv_wm = pgrp->mwscale;
257 auto localAtomIndices = pgrp->atomSet.localIndex();
259 /* The cylinder group is always a slab in the system, thus large.
260 * Therefore we always thread-parallelize this group.
262 int numAtomsLocal = localAtomIndices.size();
263 #pragma omp parallel for num_threads(nthreads) schedule(static)
264 for (int i = 0; i < numAtomsLocal; i++)
266 double weight = pgrp->localWeights[i];
271 int ii = localAtomIndices[i];
272 double mass = masses[ii];
273 /* The stored axial distance from the cylinder center (dv) needs
274 * to be corrected for an offset (dv_corr), which was unknown when
277 double dv_com = pgrp->dv[i] + dv_corr;
279 /* Here we not only add the pull force working along vec (f_pull),
280 * but also a radial component, due to the dependence of the weights
281 * on the radial distance.
283 for (int m = 0; m < DIM; m++)
285 f[ii][m] += sign * inv_wm * (mass * weight * f_pull[m] + pgrp->mdw[i][m] * dv_com * f_scal);
290 /* Apply torque forces in a mass weighted fashion to the groups that define
291 * the pull vector direction for pull coordinate pcrd.
293 static void apply_forces_vec_torque(const struct pull_t* pull,
294 const pull_coord_work_t* pcrd,
298 const PullCoordSpatialData& spatialData = pcrd->spatialData;
300 /* The component inpr along the pull vector is accounted for in the usual
301 * way. Here we account for the component perpendicular to vec.
304 for (int m = 0; m < DIM; m++)
306 inpr += spatialData.dr01[m] * spatialData.vec[m];
308 /* The torque force works along the component of the distance vector
309 * of between the two "usual" pull groups that is perpendicular to
310 * the pull vector. The magnitude of this force is the "usual" scale force
311 * multiplied with the ratio of the distance between the two "usual" pull
312 * groups and the distance between the two groups that define the vector.
315 for (int m = 0; m < DIM; m++)
317 f_perp[m] = (spatialData.dr01[m] - inpr * spatialData.vec[m]) / spatialData.vec_len * pcrd->scalarForce;
320 /* Apply the force to the groups defining the vector using opposite signs */
321 apply_forces_grp(&pull->group[pcrd->params.group[2]], masses, f_perp, -1, f, pull->nthreads);
322 apply_forces_grp(&pull->group[pcrd->params.group[3]], masses, f_perp, 1, f, pull->nthreads);
325 /* Apply forces in a mass weighted fashion */
326 static void apply_forces_coord(struct pull_t* pull,
328 const PullCoordVectorForces& forces,
332 /* Here it would be more efficient to use one large thread-parallel
333 * region instead of potential parallel regions within apply_forces_grp.
334 * But there could be overlap between pull groups and this would lead
338 const pull_coord_work_t& pcrd = pull->coord[coord];
340 if (pcrd.params.eGeom == PullGroupGeometry::Cylinder)
342 apply_forces_cyl_grp(&pull->dyna[coord],
343 pcrd.spatialData.cyl_dev,
351 /* Sum the force along the vector and the radial force */
353 for (int m = 0; m < DIM; m++)
355 f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
357 apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, f_tot, 1, f, pull->nthreads);
361 if (pcrd.params.eGeom == PullGroupGeometry::DirectionRelative)
363 /* We need to apply the torque forces to the pull groups
364 * that define the pull vector.
366 apply_forces_vec_torque(pull, &pcrd, masses, f);
369 if (!pull->group[pcrd.params.group[0]].params.ind.empty())
372 &pull->group[pcrd.params.group[0]], masses, forces.force01, -1, f, pull->nthreads);
374 apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, forces.force01, 1, f, pull->nthreads);
376 if (pcrd.params.ngroup >= 4)
379 &pull->group[pcrd.params.group[2]], masses, forces.force23, -1, f, pull->nthreads);
380 apply_forces_grp(&pull->group[pcrd.params.group[3]], masses, forces.force23, 1, f, pull->nthreads);
382 if (pcrd.params.ngroup >= 6)
385 &pull->group[pcrd.params.group[4]], masses, forces.force45, -1, f, pull->nthreads);
386 apply_forces_grp(&pull->group[pcrd.params.group[5]], masses, forces.force45, 1, f, pull->nthreads);
391 real max_pull_distance2(const pull_coord_work_t* pcrd, const t_pbc* pbc)
393 /* Note that this maximum distance calculation is more complex than
394 * most other cases in GROMACS, since here we have to take care of
395 * distance calculations that don't involve all three dimensions.
396 * For example, we can use distances that are larger than the
397 * box X and Y dimensions for a box that is elongated along Z.
400 real max_d2 = GMX_REAL_MAX;
402 if (pull_coordinate_is_directional(&pcrd->params))
404 /* Directional pulling along along direction pcrd->vec.
405 * Calculating the exact maximum distance is complex and bug-prone.
406 * So we take a safe approach by not allowing distances that
407 * are larger than half the distance between unit cell faces
408 * along dimensions involved in pcrd->vec.
410 for (int m = 0; m < DIM; m++)
412 if (m < pbc->ndim_ePBC && pcrd->spatialData.vec[m] != 0)
414 real imageDistance2 = gmx::square(pbc->box[m][m]);
415 for (int d = m + 1; d < DIM; d++)
417 imageDistance2 -= gmx::square(pbc->box[d][m]);
419 max_d2 = std::min(max_d2, imageDistance2);
425 /* Distance pulling along dimensions with pcrd->params.dim[d]==1.
426 * We use half the minimum box vector length of the dimensions involved.
427 * This is correct for all cases, except for corner cases with
428 * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
429 * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
430 * in such corner cases the user could get correct results,
431 * depending on the details of the setup, we avoid further
432 * code complications.
434 for (int m = 0; m < DIM; m++)
436 if (m < pbc->ndim_ePBC && pcrd->params.dim[m] != 0)
438 real imageDistance2 = gmx::square(pbc->box[m][m]);
439 for (int d = 0; d < m; d++)
441 if (pcrd->params.dim[d] != 0)
443 imageDistance2 += gmx::square(pbc->box[m][d]);
446 max_d2 = std::min(max_d2, imageDistance2);
451 return 0.25 * max_d2;
454 /* This function returns the distance based on coordinates xg and xref.
455 * Note that the pull coordinate struct pcrd is not modified.
457 * \param[in] pull The pull struct
458 * \param[in] pcrd The pull coordinate to compute a distance for
459 * \param[in] pbc The periodic boundary conditions
460 * \param[in] xg The coordinate of group 1
461 * \param[in] xref The coordinate of group 0
462 * \param[in] groupIndex0 The index of group 0 in the pcrd->params.group
463 * \param[in] groupIndex1 The index of group 1 in the pcrd->params.group
464 * \param[in] max_dist2 The maximum distance squared
465 * \param[out] dr The distance vector
467 static void low_get_pull_coord_dr(const struct pull_t* pull,
468 const pull_coord_work_t* pcrd,
472 const int groupIndex0,
473 const int groupIndex1,
474 const double max_dist2,
477 const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
479 /* Only the first group can be an absolute reference, in that case nat=0 */
480 if (pgrp0->params.ind.empty())
482 for (int m = 0; m < DIM; m++)
484 xref[m] = pcrd->params.origin[m];
489 copy_dvec(xref, xrefr);
491 dvec dref = { 0, 0, 0 };
492 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
494 for (int m = 0; m < DIM; m++)
496 dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
498 /* Add the reference position, so we use the correct periodic image */
499 dvec_inc(xrefr, dref);
502 pbc_dx_d(pbc, xg, xrefr, dr);
504 bool directional = pull_coordinate_is_directional(&pcrd->params);
506 for (int m = 0; m < DIM; m++)
508 dr[m] *= pcrd->params.dim[m];
509 if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
511 dr2 += dr[m] * dr[m];
514 /* Check if we are close to switching to another periodic image */
515 if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
517 /* Note that technically there is no issue with switching periodic
518 * image, as pbc_dx_d returns the distance to the closest periodic
519 * image. However in all cases where periodic image switches occur,
520 * the pull results will be useless in practice.
523 "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
524 "box size (%f).\n%s",
525 pcrd->params.group[groupIndex0],
526 pcrd->params.group[groupIndex1],
528 sqrt(0.98 * 0.98 * max_dist2),
529 pcrd->params.eGeom == PullGroupGeometry::Direction
530 ? "You might want to consider using \"pull-geometry = "
531 "direction-periodic\" instead.\n"
535 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
541 /* This function returns the distance based on the contents of the pull struct.
542 * pull->coord[coord_ind].dr, and potentially vector, are updated.
544 static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
546 pull_coord_work_t* pcrd = &pull->coord[coord_ind];
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 pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
599 pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
601 low_get_pull_coord_dr(
606 pcrd->params.eGeom == PullGroupGeometry::Cylinder ? pull->dyna[coord_ind].x : pgrp0->x,
612 if (pcrd->params.ngroup >= 4)
614 pull_group_work_t *pgrp2, *pgrp3;
615 pgrp2 = &pull->group[pcrd->params.group[2]];
616 pgrp3 = &pull->group[pcrd->params.group[3]];
618 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, 2, 3, md2, spatialData.dr23);
620 if (pcrd->params.ngroup >= 6)
622 pull_group_work_t *pgrp4, *pgrp5;
623 pgrp4 = &pull->group[pcrd->params.group[4]];
624 pgrp5 = &pull->group[pcrd->params.group[5]];
626 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, 4, 5, md2, spatialData.dr45);
630 /* Modify x so that it is periodic in [-pi, pi)
631 * It is assumed that x is in [-3pi, 3pi) so that x
632 * needs to be shifted by at most one period.
634 static void make_periodic_2pi(double* x)
646 /* This function should always be used to modify pcrd->value_ref */
647 static void low_set_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double value_ref)
649 GMX_ASSERT(pcrd->params.eType != PullingAlgorithm::External,
650 "The pull coord reference value should not be used with type external-potential");
652 if (pcrd->params.eGeom == PullGroupGeometry::Distance)
657 "Pull reference distance for coordinate %d (%f) needs to be non-negative",
662 else if (pcrd->params.eGeom == PullGroupGeometry::Angle
663 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
665 if (value_ref < 0 || value_ref > M_PI)
668 "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
669 "interval [0,180] deg",
671 value_ref * pull_conversion_factor_internal2userinput(&pcrd->params));
674 else if (pcrd->params.eGeom == PullGroupGeometry::Dihedral)
676 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
677 make_periodic_2pi(&value_ref);
680 pcrd->value_ref = value_ref;
683 static void update_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double t)
685 /* With zero rate the reference value is set initially and doesn't change */
686 if (pcrd->params.rate != 0)
688 double value_ref = (pcrd->params.init + pcrd->params.rate * t)
689 * pull_conversion_factor_userinput2internal(&pcrd->params);
690 low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
694 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
695 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
698 dvec dr32; /* store instead of dr23? */
700 dsvmul(-1, spatialData->dr23, dr32);
701 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
702 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
703 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
705 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
706 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
707 * Note 2: the angle between the plane normal vectors equals pi only when
708 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
709 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
711 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
715 /* Calculates pull->coord[coord_ind].value.
716 * This function also updates pull->coord[coord_ind].dr.
718 static void get_pull_coord_distance(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
720 pull_coord_work_t* pcrd;
723 pcrd = &pull->coord[coord_ind];
725 get_pull_coord_dr(pull, coord_ind, pbc);
727 PullCoordSpatialData& spatialData = pcrd->spatialData;
729 switch (pcrd->params.eGeom)
731 case PullGroupGeometry::Distance:
732 /* Pull along the vector between the com's */
733 spatialData.value = dnorm(spatialData.dr01);
735 case PullGroupGeometry::Direction:
736 case PullGroupGeometry::DirectionPBC:
737 case PullGroupGeometry::DirectionRelative:
738 case PullGroupGeometry::Cylinder:
740 spatialData.value = 0;
741 for (m = 0; m < DIM; m++)
743 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
746 case PullGroupGeometry::Angle:
747 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
749 case PullGroupGeometry::Dihedral:
750 spatialData.value = get_dihedral_angle_coord(&spatialData);
752 case PullGroupGeometry::AngleAxis:
753 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
755 default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
759 /* Returns the deviation from the reference value.
760 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
762 static double get_pull_coord_deviation(struct pull_t* pull, int coord_ind, const t_pbc* pbc, double t)
764 pull_coord_work_t* pcrd;
767 pcrd = &pull->coord[coord_ind];
769 /* Update the reference value before computing the distance,
770 * since it is used in the distance computation with periodic pulling.
772 update_pull_coord_reference_value(pcrd, coord_ind, t);
774 get_pull_coord_distance(pull, coord_ind, pbc);
776 get_pull_coord_distance(pull, coord_ind, pbc);
778 /* Determine the deviation */
779 dev = pcrd->spatialData.value - pcrd->value_ref;
781 /* Check that values are allowed */
782 if (pcrd->params.eGeom == PullGroupGeometry::Distance && pcrd->spatialData.value == 0)
784 /* With no vector we can not determine the direction for the force,
785 * so we set the force to zero.
789 else if (pcrd->params.eGeom == PullGroupGeometry::Dihedral)
791 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
792 Thus, the unwrapped deviation is here in (-2pi, 2pi].
793 After making it periodic, the deviation will be in [-pi, pi). */
794 make_periodic_2pi(&dev);
800 double get_pull_coord_value(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
802 get_pull_coord_distance(pull, coord_ind, pbc);
804 return pull->coord[coord_ind].spatialData.value;
807 void clear_pull_forces(pull_t* pull)
809 /* Zeroing the forces is only required for constraint pulling.
810 * It can happen that multiple constraint steps need to be applied
811 * and therefore the constraint forces need to be accumulated.
813 for (pull_coord_work_t& coord : pull->coord)
815 coord.scalarForce = 0;
819 /* Apply constraint using SHAKE */
821 do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
824 dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
825 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
826 dvec* rnew; /* current 'new' positions of the groups */
827 double* dr_tot; /* the total update of the coords */
830 double lambda, rm, invdt = 0;
831 gmx_bool bConverged_all, bConverged = FALSE;
832 int niter = 0, ii, j, m, max_iter = 100;
836 snew(r_ij, pull->coord.size());
837 snew(dr_tot, pull->coord.size());
839 snew(rnew, pull->group.size());
841 /* copy the current unconstrained positions for use in iterations. We
842 iterate until rinew[i] and rjnew[j] obey the constraints. Then
843 rinew - pull.x_unc[i] is the correction dr to group i */
844 for (size_t g = 0; g < pull->group.size(); g++)
846 copy_dvec(pull->group[g].xp, rnew[g]);
849 /* Determine the constraint directions from the old positions */
850 for (size_t c = 0; c < pull->coord.size(); c++)
852 pull_coord_work_t* pcrd;
854 pcrd = &pull->coord[c];
856 if (pcrd->params.eType != PullingAlgorithm::Constraint)
861 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
862 * We don't modify dr and value anymore, so these values are also used
865 get_pull_coord_distance(pull, c, pbc);
867 const PullCoordSpatialData& spatialData = pcrd->spatialData;
871 "Pull coord %zu dr %f %f %f\n",
873 spatialData.dr01[XX],
874 spatialData.dr01[YY],
875 spatialData.dr01[ZZ]);
878 if (pcrd->params.eGeom == PullGroupGeometry::Direction
879 || pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
881 /* Select the component along vec */
883 for (m = 0; m < DIM; m++)
885 a += spatialData.vec[m] * spatialData.dr01[m];
887 for (m = 0; m < DIM; m++)
889 r_ij[c][m] = a * spatialData.vec[m];
894 copy_dvec(spatialData.dr01, r_ij[c]);
897 if (dnorm2(r_ij[c]) == 0)
900 "Distance for pull coordinate %zu is zero with constraint pulling, which is "
906 bConverged_all = FALSE;
907 while (!bConverged_all && niter < max_iter)
909 bConverged_all = TRUE;
911 /* loop over all constraints */
912 for (size_t c = 0; c < pull->coord.size(); c++)
914 pull_coord_work_t* pcrd;
915 pull_group_work_t *pgrp0, *pgrp1;
918 pcrd = &pull->coord[c];
920 if (pcrd->params.eType != PullingAlgorithm::Constraint)
925 update_pull_coord_reference_value(pcrd, c, t);
927 pgrp0 = &pull->group[pcrd->params.group[0]];
928 pgrp1 = &pull->group[pcrd->params.group[1]];
930 /* Get the current difference vector */
931 low_get_pull_coord_dr(
932 pull, pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
936 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
939 rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
941 switch (pcrd->params.eGeom)
943 case PullGroupGeometry::Distance:
944 if (pcrd->value_ref <= 0)
948 "The pull constraint reference distance for group %zu is <= 0 (%f)",
954 double q, c_a, c_b, c_c;
956 c_a = diprod(r_ij[c], r_ij[c]);
957 c_b = diprod(unc_ij, r_ij[c]) * 2;
958 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
962 q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
967 q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
973 fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b, c_c, lambda);
977 /* The position corrections dr due to the constraints */
978 dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
979 dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
980 dr_tot[c] += -lambda * dnorm(r_ij[c]);
982 case PullGroupGeometry::Direction:
983 case PullGroupGeometry::DirectionPBC:
984 case PullGroupGeometry::Cylinder:
985 /* A 1-dimensional constraint along a vector */
987 for (m = 0; m < DIM; m++)
989 vec[m] = pcrd->spatialData.vec[m];
990 a += unc_ij[m] * vec[m];
992 /* Select only the component along the vector */
993 dsvmul(a, vec, unc_ij);
994 lambda = a - pcrd->value_ref;
997 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
1000 /* The position corrections dr due to the constraints */
1001 dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
1002 dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
1003 dr_tot[c] += -lambda;
1005 default: gmx_incons("Invalid enumeration value for eGeom");
1013 g0 = pcrd->params.group[0];
1014 g1 = pcrd->params.group[1];
1015 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
1016 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
1018 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1027 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
1036 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1046 /* Update the COMs with dr */
1047 dvec_inc(rnew[pcrd->params.group[1]], dr1);
1048 dvec_inc(rnew[pcrd->params.group[0]], dr0);
1051 /* Check if all constraints are fullfilled now */
1052 for (pull_coord_work_t& coord : pull->coord)
1054 if (coord.params.eType != PullingAlgorithm::Constraint)
1059 low_get_pull_coord_dr(
1060 pull, &coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
1062 switch (coord.params.eGeom)
1064 case PullGroupGeometry::Distance:
1065 bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1067 case PullGroupGeometry::Direction:
1068 case PullGroupGeometry::DirectionPBC:
1069 case PullGroupGeometry::Cylinder:
1070 for (m = 0; m < DIM; m++)
1072 vec[m] = coord.spatialData.vec[m];
1074 inpr = diprod(unc_ij, vec);
1075 dsvmul(inpr, vec, unc_ij);
1076 bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1080 gmx::formatString("Geometry %s not handled here",
1081 enumValueToString(coord.params.eGeom))
1090 "Pull constraint not converged: "
1092 "d_ref = %f, current d = %f\n",
1093 coord.params.group[0],
1094 coord.params.group[1],
1099 bConverged_all = FALSE;
1104 /* if after all constraints are dealt with and bConverged is still TRUE
1105 we're finished, if not we do another iteration */
1107 if (niter > max_iter)
1109 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1112 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1119 /* update atoms in the groups */
1120 for (size_t g = 0; g < pull->group.size(); g++)
1122 const pull_group_work_t* pgrp;
1125 pgrp = &pull->group[g];
1127 /* get the final constraint displacement dr for group g */
1128 dvec_sub(rnew[g], pgrp->xp, dr);
1130 if (dnorm2(dr) == 0)
1132 /* No displacement, no update necessary */
1136 /* update the atom positions */
1137 auto localAtomIndices = pgrp->atomSet.localIndex();
1139 for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1141 ii = localAtomIndices[j];
1142 if (!pgrp->localWeights.empty())
1144 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1146 for (m = 0; m < DIM; m++)
1152 for (m = 0; m < DIM; m++)
1154 v[ii][m] += invdt * tmp[m];
1160 /* calculate the constraint forces, used for output and virial only */
1161 for (size_t c = 0; c < pull->coord.size(); c++)
1163 pull_coord_work_t* pcrd;
1165 pcrd = &pull->coord[c];
1167 if (pcrd->params.eType != PullingAlgorithm::Constraint)
1172 /* Accumulate the forces, in case we have multiple constraint steps */
1175 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1177 pcrd->scalarForce += force;
1179 if (vir != nullptr && pcrd->params.eGeom != PullGroupGeometry::DirectionPBC && bMaster)
1183 /* Add the pull contribution to the virial */
1184 /* We have already checked above that r_ij[c] != 0 */
1185 f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1187 for (j = 0; j < DIM; j++)
1189 for (m = 0; m < DIM; m++)
1191 vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1197 /* finished! I hope. Give back some memory */
1203 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1205 for (int j = 0; j < DIM; j++)
1207 for (int m = 0; m < DIM; m++)
1209 vir[j][m] -= 0.5 * f[j] * dr[m];
1214 /* Adds the pull contribution to the virial */
1215 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1217 if (vir != nullptr && pcrd.params.eGeom != PullGroupGeometry::DirectionPBC)
1219 /* Add the pull contribution for each distance vector to the virial. */
1220 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1221 if (pcrd.params.ngroup >= 4)
1223 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1225 if (pcrd.params.ngroup >= 6)
1227 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1232 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1240 k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1241 dkdl = pcrd->params.kB - pcrd->params.k;
1243 switch (pcrd->params.eType)
1245 case PullingAlgorithm::Umbrella:
1246 case PullingAlgorithm::FlatBottom:
1247 case PullingAlgorithm::FlatBottomHigh:
1248 /* The only difference between an umbrella and a flat-bottom
1249 * potential is that a flat-bottom is zero above or below
1250 the reference value.
1252 if ((pcrd->params.eType == PullingAlgorithm::FlatBottom && dev < 0)
1253 || (pcrd->params.eType == PullingAlgorithm::FlatBottomHigh && dev > 0))
1258 pcrd->scalarForce = -k * dev;
1259 *V += 0.5 * k * gmx::square(dev);
1260 *dVdl += 0.5 * dkdl * gmx::square(dev);
1262 case PullingAlgorithm::ConstantForce:
1263 pcrd->scalarForce = -k;
1264 *V += k * pcrd->spatialData.value;
1265 *dVdl += dkdl * pcrd->spatialData.value;
1267 case PullingAlgorithm::External:
1269 "the scalar pull force should not be calculated internally for pull type "
1271 default: gmx_incons("Unsupported pull type in do_pull_pot");
1275 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1277 const t_pull_coord& params = pcrd.params;
1278 const PullCoordSpatialData& spatialData = pcrd.spatialData;
1280 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1281 PullCoordVectorForces forces;
1283 if (params.eGeom == PullGroupGeometry::Distance)
1285 double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1286 for (int m = 0; m < DIM; m++)
1288 forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1291 else if (params.eGeom == PullGroupGeometry::Angle)
1294 double cos_theta, cos_theta2;
1296 cos_theta = cos(spatialData.value);
1297 cos_theta2 = gmx::square(cos_theta);
1299 /* The force at theta = 0, pi is undefined so we should not apply any force.
1300 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1301 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1302 * have to check for this before dividing by their norm below.
1306 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1307 * and another vector parallel to dr23:
1308 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1309 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1311 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1312 double b = a * cos_theta;
1313 double invdr01 = 1. / dnorm(spatialData.dr01);
1314 double invdr23 = 1. / dnorm(spatialData.dr23);
1315 dvec normalized_dr01, normalized_dr23;
1316 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1317 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1319 for (int m = 0; m < DIM; m++)
1321 /* Here, f_scal is -dV/dtheta */
1323 pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1325 pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1330 /* No forces to apply for ill-defined cases*/
1331 clear_dvec(forces.force01);
1332 clear_dvec(forces.force23);
1335 else if (params.eGeom == PullGroupGeometry::AngleAxis)
1337 double cos_theta, cos_theta2;
1339 /* The angle-axis force is exactly the same as the angle force (above) except that in
1340 this case the second vector (dr23) is replaced by the pull vector. */
1341 cos_theta = cos(spatialData.value);
1342 cos_theta2 = gmx::square(cos_theta);
1348 dvec normalized_dr01;
1350 invdr01 = 1. / dnorm(spatialData.dr01);
1351 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1352 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1355 for (int m = 0; m < DIM; m++)
1358 pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1363 clear_dvec(forces.force01);
1366 else if (params.eGeom == PullGroupGeometry::Dihedral)
1368 double m2, n2, tol, sqrdist_32;
1370 /* Note: there is a small difference here compared to the
1371 dihedral force calculations in the bondeds (ref: Bekker 1994).
1372 There rij = ri - rj, while here dr01 = r1 - r0.
1373 However, all distance vectors occur in form of cross or inner products
1374 so that two signs cancel and we end up with the same expressions.
1375 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1377 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1378 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1379 dsvmul(-1, spatialData.dr23, dr32);
1380 sqrdist_32 = diprod(dr32, dr32);
1381 tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1382 if ((m2 > tol) && (n2 > tol))
1384 double a_01, a_23_01, a_23_45, a_45;
1385 double inv_dist_32, inv_sqrdist_32, dist_32;
1387 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1388 inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1389 dist_32 = sqrdist_32 * inv_dist_32;
1391 /* Forces on groups 0, 1 */
1392 a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
1393 dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1395 /* Forces on groups 4, 5 */
1396 a_45 = -pcrd.scalarForce * dist_32 / n2;
1397 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1399 /* Force on groups 2, 3 (defining the axis) */
1400 a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1401 a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1402 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1403 dsvmul(a_23_45, forces.force45, v);
1404 dvec_sub(u, v, forces.force23); /* force on group 3 */
1408 /* No force to apply for ill-defined cases */
1409 clear_dvec(forces.force01);
1410 clear_dvec(forces.force23);
1411 clear_dvec(forces.force45);
1416 for (int m = 0; m < DIM; m++)
1418 forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1426 /* We use a global mutex for locking access to the pull data structure
1427 * during registration of external pull potential providers.
1428 * We could use a different, local mutex for each pull object, but the overhead
1429 * is extremely small here and registration is only done during initialization.
1431 static std::mutex registrationMutex;
1433 using Lock = std::lock_guard<std::mutex>;
1435 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1437 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1438 GMX_RELEASE_ASSERT(provider != nullptr,
1439 "register_external_pull_potential called with NULL as provider name");
1441 if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1444 "Module '%s' attempted to register an external potential for pull coordinate %d "
1445 "which is out of the pull coordinate range %d - %zu\n",
1449 pull->coord.size());
1452 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1454 if (pcrd->params.eType != PullingAlgorithm::External)
1458 "Module '%s' attempted to register an external potential for pull coordinate %d "
1459 "which of type '%s', whereas external potentials are only supported with type '%s'",
1462 enumValueToString(pcrd->params.eType),
1463 enumValueToString(PullingAlgorithm::External));
1466 GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
1467 "The external potential provider string for a pull coordinate is NULL");
1469 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
1472 "Module '%s' attempted to register an external potential for pull coordinate %d "
1473 "which expects the external potential to be provided by a module named '%s'",
1476 pcrd->params.externalPotentialProvider.c_str());
1479 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1480 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1481 * pull->numUnregisteredExternalPotentials.
1483 Lock registrationLock(registrationMutex);
1485 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1488 "Module '%s' attempted to register an external potential for pull coordinate %d "
1494 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1495 pull->numUnregisteredExternalPotentials--;
1497 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1498 "Negative unregistered potentials, the pull code is inconsistent");
1502 static void check_external_potential_registration(const struct pull_t* pull)
1504 if (pull->numUnregisteredExternalPotentials > 0)
1507 for (c = 0; c < pull->coord.size(); c++)
1509 if (pull->coord[c].params.eType == PullingAlgorithm::External
1510 && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1516 GMX_RELEASE_ASSERT(c < pull->coord.size(),
1517 "Internal inconsistency in the pull potential provider counting");
1520 "No external provider for external pull potentials have been provided for %d "
1521 "pull coordinates. The first coordinate without provider is number %zu, which "
1522 "expects a module named '%s' to provide the external potential.",
1523 pull->numUnregisteredExternalPotentials,
1525 pull->coord[c].params.externalPotentialProvider.c_str());
1530 /* Pull takes care of adding the forces of the external potential.
1531 * The external potential module has to make sure that the corresponding
1532 * potential energy is added either to the pull term or to a term
1533 * specific to the external module.
1535 void apply_external_pull_coord_force(struct pull_t* pull,
1539 gmx::ForceWithVirial* forceWithVirial)
1541 pull_coord_work_t* pcrd;
1543 GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1544 "apply_external_pull_coord_force called with coord_index out of range");
1546 if (pull->comm.bParticipate)
1548 pcrd = &pull->coord[coord_index];
1551 pcrd->params.eType == PullingAlgorithm::External,
1552 "The pull force can only be set externally on pull coordinates of external type");
1554 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1555 "apply_external_pull_coord_force called for an unregistered pull coordinate");
1558 pcrd->scalarForce = coord_force;
1560 /* Calculate the forces on the pull groups */
1561 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1563 /* Add the forces for this coordinate to the total virial and force */
1564 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1566 matrix virial = { { 0 } };
1567 add_virial_coord(virial, *pcrd, pullCoordForces);
1568 forceWithVirial->addVirialContribution(virial);
1572 pull, coord_index, pullCoordForces, masses, as_rvec_array(forceWithVirial->force_.data()));
1575 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1578 /* Calculate the pull potential and scalar force for a pull coordinate.
1579 * Returns the vector forces for the pull coordinate.
1581 static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull,
1590 pull_coord_work_t& pcrd = pull->coord[coord_ind];
1592 assert(pcrd.params.eType != PullingAlgorithm::Constraint);
1594 double dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
1596 calc_pull_coord_scalar_force_and_potential(&pcrd, dev, lambda, V, dVdl);
1598 PullCoordVectorForces pullCoordForces = calculateVectorForces(pcrd);
1600 add_virial_coord(vir, pcrd, pullCoordForces);
1602 return pullCoordForces;
1605 real pull_potential(struct pull_t* pull,
1608 const t_commrec* cr,
1612 gmx::ForceWithVirial* force,
1617 assert(pull != nullptr);
1619 /* Ideally we should check external potential registration only during
1620 * the initialization phase, but that requires another function call
1621 * that should be done exactly in the right order. So we check here.
1623 check_external_potential_registration(pull);
1625 if (pull->comm.bParticipate)
1629 pull_calc_coms(cr, pull, masses, pbc, t, x, nullptr);
1631 rvec* f = as_rvec_array(force->force_.data());
1632 matrix virial = { { 0 } };
1633 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1634 for (size_t c = 0; c < pull->coord.size(); c++)
1636 pull_coord_work_t* pcrd;
1637 pcrd = &pull->coord[c];
1639 /* For external potential the force is assumed to be given by an external module by a
1640 call to apply_pull_coord_external_force */
1641 if (pcrd->params.eType == PullingAlgorithm::Constraint
1642 || pcrd->params.eType == PullingAlgorithm::External)
1647 PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1648 pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1650 /* Distribute the force over the atoms in the pulled groups */
1651 apply_forces_coord(pull, c, pullCoordForces, masses, f);
1656 force->addVirialContribution(virial);
1661 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1662 "Too few or too many external pull potentials have been applied the previous step");
1663 /* All external pull potentials still need to be applied */
1664 pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1666 return (MASTER(cr) ? V : 0.0);
1669 void pull_constraint(struct pull_t* pull,
1672 const t_commrec* cr,
1680 assert(pull != nullptr);
1682 if (pull->comm.bParticipate)
1684 pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
1686 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1690 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1694 gmx_bool bMustParticipate;
1700 /* We always make the master node participate, such that it can do i/o,
1701 * add the virial and to simplify MC type extensions people might have.
1703 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1705 for (pull_group_work_t& group : pull->group)
1707 if (!group.globalWeights.empty())
1709 group.localWeights.resize(group.atomSet.numAtomsLocal());
1710 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1712 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1716 GMX_ASSERT(bMustParticipate || dd != nullptr,
1717 "Either all ranks (including this rank) participate, or we use DD and need to "
1718 "have access to dd here");
1720 /* We should participate if we have pull or pbc atoms */
1721 if (!bMustParticipate
1722 && (group.atomSet.numAtomsLocal() > 0
1723 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1725 bMustParticipate = TRUE;
1729 if (!comm->bParticipateAll)
1731 /* Keep currently not required ranks in the communicator
1732 * if they needed to participate up to 20 decompositions ago.
1733 * This avoids frequent rebuilds due to atoms jumping back and forth.
1735 const int64_t history_count = 20;
1736 gmx_bool bWillParticipate;
1739 /* Increase the decomposition counter for the current call */
1740 comm->setup_count++;
1742 if (bMustParticipate)
1744 comm->must_count = comm->setup_count;
1749 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1751 if (debug && dd != nullptr)
1754 "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
1756 gmx::boolToString(bMustParticipate),
1757 gmx::boolToString(bWillParticipate));
1760 if (bWillParticipate)
1762 /* Count the number of ranks that we want to have participating */
1764 /* Count the number of ranks that need to be added */
1765 count[1] = comm->bParticipate ? 0 : 1;
1773 /* The cost of this global operation will be less that the cost
1774 * of the extra MPI_Comm_split calls that we can avoid.
1776 gmx_sumi(2, count, cr);
1778 /* If we are missing ranks or if we have 20% more ranks than needed
1779 * we make a new sub-communicator.
1781 if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1785 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1789 if (comm->mpi_comm_com != MPI_COMM_NULL)
1791 MPI_Comm_free(&comm->mpi_comm_com);
1794 /* This might be an extremely expensive operation, so we try
1795 * to avoid this splitting as much as possible.
1797 assert(dd != nullptr);
1798 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1801 comm->bParticipate = bWillParticipate;
1802 comm->nparticipate = count[0];
1804 /* When we use the previous COM for PBC, we need to broadcast
1805 * the previous COM to ranks that have joined the communicator.
1807 for (pull_group_work_t& group : pull->group)
1809 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1811 GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1812 "The master rank has to participate, as it should pass an up to "
1814 "to bcast here as well as to e.g. checkpointing");
1816 gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
1822 /* Since the PBC of atoms might have changed, we need to update the PBC */
1823 pull->bSetPBCatoms = TRUE;
1826 static void init_pull_group_index(FILE* fplog,
1827 const t_commrec* cr,
1829 pull_group_work_t* pg,
1830 gmx_bool bConstraint,
1831 const ivec pulldim_con,
1832 const gmx_mtop_t& mtop,
1833 const t_inputrec* ir,
1836 /* With EM and BD there are no masses in the integrator.
1837 * But we still want to have the correct mass-weighted COMs.
1838 * So we store the real masses in the weights.
1840 const bool setWeights = (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1841 || ir->eI == IntegrationAlgorithm::BD);
1843 /* In parallel, store we need to extract localWeights from weights at DD time */
1844 std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1846 const SimulationGroups& groups = mtop.groups;
1848 /* Count frozen dimensions and (weighted) mass */
1854 for (int i = 0; i < int(pg->params.ind.size()); i++)
1856 int ii = pg->params.ind[i];
1857 if (bConstraint && ir->opts.nFreeze)
1859 for (int d = 0; d < DIM; d++)
1861 if (pulldim_con[d] == 1
1862 && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1868 const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1870 if (ir->efep == FreeEnergyPerturbationType::No)
1876 m = (1 - lambda) * atom.m + lambda * atom.mB;
1879 if (!pg->params.weight.empty())
1881 w = pg->params.weight[i];
1887 if (EI_ENERGY_MINIMIZATION(ir->eI))
1889 /* Move the mass to the weight */
1893 else if (ir->eI == IntegrationAlgorithm::BD)
1896 if (ir->bd_fric != 0.0F)
1898 mbd = ir->bd_fric * ir->delta_t;
1902 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1904 mbd = ir->delta_t / ir->opts.tau_t[0];
1909 / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1917 weights.push_back(w);
1921 wwmass += m * w * w;
1926 /* We can have single atom groups with zero mass with potential pulling
1927 * without cosine weighting.
1929 if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1931 /* With one atom the mass doesn't matter */
1937 "The total%s mass of pull group %d is zero",
1938 !pg->params.weight.empty() ? " weighted" : "",
1944 fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
1945 if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1946 || ir->eI == IntegrationAlgorithm::BD)
1948 fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1950 if (pg->epgrppbc == epgrppbcCOS)
1952 fprintf(fplog, ", cosine weighting will be used");
1954 fprintf(fplog, "\n");
1959 /* A value != 0 signals not frozen, it is updated later */
1965 for (int d = 0; d < DIM; d++)
1967 ndim += pulldim_con[d] * pg->params.ind.size();
1969 if (fplog && nfrozen > 0 && nfrozen < ndim)
1972 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1973 " that are subject to pulling are frozen.\n"
1974 " For constraint pulling the whole group will be frozen.\n\n",
1982 struct pull_t* init_pull(FILE* fplog,
1983 const pull_params_t* pull_params,
1984 const t_inputrec* ir,
1985 const gmx_mtop_t& mtop,
1986 const t_commrec* cr,
1987 gmx::LocalAtomSetManager* atomSets,
1990 struct pull_t* pull;
1995 /* Copy the pull parameters */
1996 pull->params = *pull_params;
1998 for (int i = 0; i < pull_params->ngroup; ++i)
2000 pull->group.emplace_back(pull_params->group[i],
2001 atomSets->add(pull_params->group[i].ind),
2002 pull_params->bSetPbcRefToPrevStepCOM);
2005 if (cr != nullptr && DOMAINDECOMP(cr))
2007 /* Set up the global to local atom mapping for PBC atoms */
2008 for (pull_group_work_t& group : pull->group)
2010 if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
2012 /* pbcAtomSet consists of a single atom */
2013 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
2014 atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
2019 pull->bPotential = FALSE;
2020 pull->bConstraint = FALSE;
2021 pull->bCylinder = FALSE;
2022 pull->bAngle = FALSE;
2023 pull->bXOutAverage = pull_params->bXOutAverage;
2024 pull->bFOutAverage = pull_params->bFOutAverage;
2026 GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
2027 "pull group 0 is an absolute reference group and should not contain atoms");
2029 pull->numCoordinatesWithExternalPotential = 0;
2031 for (int c = 0; c < pull->params.ncoord; c++)
2033 /* Construct a pull coordinate, copying all coordinate parameters */
2034 pull->coord.emplace_back(pull_params->coord[c]);
2036 pull_coord_work_t* pcrd = &pull->coord.back();
2038 switch (pcrd->params.eGeom)
2040 case PullGroupGeometry::Distance:
2041 case PullGroupGeometry::DirectionRelative: /* Direction vector is determined at each step */
2042 case PullGroupGeometry::Angle:
2043 case PullGroupGeometry::Dihedral: break;
2044 case PullGroupGeometry::Direction:
2045 case PullGroupGeometry::DirectionPBC:
2046 case PullGroupGeometry::Cylinder:
2047 case PullGroupGeometry::AngleAxis:
2048 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2051 /* We allow reading of newer tpx files with new pull geometries,
2052 * but with the same tpx format, with old code. A new geometry
2053 * only adds a new enum value, which will be out of range for
2054 * old code. The first place we need to generate an error is
2055 * here, since the pull code can't handle this.
2056 * The enum can be used before arriving here only for printing
2057 * the string corresponding to the geometry, which will result
2058 * in printing "UNDEFINED".
2061 "Pull geometry not supported for pull coordinate %d. The geometry enum "
2062 "%s in the input is larger than that supported by the code (up to %d). "
2063 "You are probably reading a tpr file generated with a newer version of "
2064 "GROMACS with an binary from an older version of Gromacs.",
2066 enumValueToString(pcrd->params.eGeom),
2067 static_cast<int>(PullGroupGeometry::Count) - 1);
2070 if (pcrd->params.eType == PullingAlgorithm::Constraint)
2072 /* Check restrictions of the constraint pull code */
2073 if (pcrd->params.eGeom == PullGroupGeometry::Cylinder
2074 || pcrd->params.eGeom == PullGroupGeometry::DirectionRelative
2075 || pcrd->params.eGeom == PullGroupGeometry::Angle
2076 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2077 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2080 "Pulling of type %s can not be combined with geometry %s. Consider using "
2082 enumValueToString(pcrd->params.eType),
2083 enumValueToString(pcrd->params.eGeom),
2084 enumValueToString(PullingAlgorithm::Umbrella));
2089 "Constraint pulling can not be combined with multiple time stepping");
2091 pull->bConstraint = TRUE;
2095 pull->bPotential = TRUE;
2098 if (pcrd->params.eGeom == PullGroupGeometry::Cylinder)
2100 pull->bCylinder = TRUE;
2102 else if (pcrd->params.eGeom == PullGroupGeometry::Angle
2103 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2104 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2106 pull->bAngle = TRUE;
2109 /* We only need to calculate the plain COM of a group
2110 * when it is not only used as a cylinder group.
2111 * Also the absolute reference group 0 needs no COM computation.
2113 for (int i = 0; i < pcrd->params.ngroup; i++)
2115 int groupIndex = pcrd->params.group[i];
2116 if (groupIndex > 0 && !(pcrd->params.eGeom == PullGroupGeometry::Cylinder && i == 0))
2118 pull->group[groupIndex].needToCalcCom = true;
2122 /* With non-zero rate the reference value is set at every step */
2123 if (pcrd->params.rate == 0)
2125 /* Initialize the constant reference value */
2126 if (pcrd->params.eType != PullingAlgorithm::External)
2128 low_set_pull_coord_reference_value(
2129 pcrd, c, pcrd->params.init * pull_conversion_factor_userinput2internal(&pcrd->params));
2133 /* With an external pull potential, the reference value loses
2134 * it's meaning and should not be used. Setting it to zero
2135 * makes any terms dependent on it disappear.
2136 * The only issue this causes is that with cylinder pulling
2137 * no atoms of the cylinder group within the cylinder radius
2138 * should be more than half a box length away from the COM of
2139 * the pull group along the axial direction.
2141 pcrd->value_ref = 0.0;
2145 if (pcrd->params.eType == PullingAlgorithm::External)
2148 pcrd->params.rate == 0,
2149 "With an external potential, a pull coordinate should have rate = 0");
2151 /* This external potential needs to be registered later */
2152 pull->numCoordinatesWithExternalPotential++;
2154 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2157 pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
2158 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2160 pull->pbcType = ir->pbcType;
2161 switch (pull->pbcType)
2163 case PbcType::No: pull->npbcdim = 0; break;
2164 case PbcType::XY: pull->npbcdim = 2; break;
2165 default: pull->npbcdim = 3; break;
2170 gmx_bool bAbs, bCos;
2173 for (const pull_coord_work_t& coord : pull->coord)
2175 if (pull->group[coord.params.group[0]].params.ind.empty()
2176 || pull->group[coord.params.group[1]].params.ind.empty())
2182 fprintf(fplog, "\n");
2183 if (pull->bPotential)
2185 fprintf(fplog, "Will apply potential COM pulling\n");
2187 if (pull->bConstraint)
2189 fprintf(fplog, "Will apply constraint COM pulling\n");
2191 // Don't include the reference group 0 in output, so we report ngroup-1
2192 int numRealGroups = pull->group.size() - 1;
2193 GMX_RELEASE_ASSERT(numRealGroups > 0,
2194 "The reference absolute position pull group should always be present");
2196 "with %zu pull coordinate%s and %d group%s\n",
2198 pull->coord.size() == 1 ? "" : "s",
2200 numRealGroups == 1 ? "" : "s");
2203 fprintf(fplog, "with an absolute reference\n");
2206 // Don't include the reference group 0 in loop
2207 for (size_t g = 1; g < pull->group.size(); g++)
2209 if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
2211 /* We are using cosine weighting */
2212 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2218 please_cite(fplog, "Engin2010");
2222 pull->bRefAt = FALSE;
2224 for (size_t g = 0; g < pull->group.size(); g++)
2226 pull_group_work_t* pgrp;
2228 pgrp = &pull->group[g];
2229 if (!pgrp->params.ind.empty())
2231 /* There is an issue when a group is used in multiple coordinates
2232 * and constraints are applied in different dimensions with atoms
2233 * frozen in some, but not all dimensions.
2234 * Since there is only one mass array per group, we can't have
2235 * frozen/non-frozen atoms for different coords at the same time.
2236 * But since this is a very exotic case, we don't check for this.
2237 * A warning is printed in init_pull_group_index.
2240 gmx_bool bConstraint;
2241 ivec pulldim, pulldim_con;
2243 /* Loop over all pull coordinates to see along which dimensions
2244 * this group is pulled and if it is involved in constraints.
2246 bConstraint = FALSE;
2247 clear_ivec(pulldim);
2248 clear_ivec(pulldim_con);
2249 for (const pull_coord_work_t& coord : pull->coord)
2251 gmx_bool bGroupUsed = FALSE;
2252 for (int gi = 0; gi < coord.params.ngroup; gi++)
2254 if (coord.params.group[gi] == static_cast<int>(g))
2262 for (int m = 0; m < DIM; m++)
2264 if (coord.params.dim[m] == 1)
2268 if (coord.params.eType == PullingAlgorithm::Constraint)
2278 /* Determine if we need to take PBC into account for calculating
2279 * the COM's of the pull groups.
2281 switch (pgrp->epgrppbc)
2283 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2285 if (!pgrp->params.weight.empty())
2288 "Pull groups can not have relative weights and cosine weighting "
2291 for (int m = 0; m < DIM; m++)
2293 if (m < pull->npbcdim && pulldim[m] == 1)
2295 if (pull->cosdim >= 0 && pull->cosdim != m)
2298 "Can only use cosine weighting with pulling in one "
2299 "dimension (use mdp option pull-coord?-dim)");
2305 case epgrppbcNONE: break;
2308 /* Set the indices */
2309 init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2313 /* Absolute reference, set the inverse mass to zero.
2314 * This is only relevant (and used) with constraint pulling.
2321 /* If we use cylinder coordinates, do some initialising for them */
2322 if (pull->bCylinder)
2324 for (const pull_coord_work_t& coord : pull->coord)
2326 if (coord.params.eGeom == PullGroupGeometry::Cylinder)
2328 if (pull->group[coord.params.group[0]].params.ind.empty())
2331 "A cylinder pull group is not supported when using absolute "
2335 const auto& referenceGroup = pull->group[coord.params.group[0]];
2336 pull->dyna.emplace_back(
2337 referenceGroup.params, referenceGroup.atomSet, pull->params.bSetPbcRefToPrevStepCOM);
2341 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2342 pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
2343 pull->comSums.resize(pull->nthreads);
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,
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(cr, pull_work, masses, &pbc, state->x.rvec_array());
2436 updatePrevStepPullCom(pull_work, state);
2440 void finish_pull(struct pull_t* pull)
2442 check_external_potential_registration(pull);
2446 gmx_fio_fclose(pull->out_x);
2450 gmx_fio_fclose(pull->out_f);
2456 bool pull_have_potential(const pull_t& pull)
2458 return pull.bPotential;
2461 bool pull_have_constraint(const pull_t& pull)
2463 return pull.bConstraint;
2466 bool pull_have_constraint(const pull_params_t& pullParameters)
2468 for (int c = 0; c < pullParameters.ncoord; c++)
2470 if (pullParameters.coord[c].eType == PullingAlgorithm::Constraint)