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/vec.h"
63 #include "gromacs/math/vectypes.h"
64 #include "gromacs/mdlib/gmx_omp_nthreads.h"
65 #include "gromacs/mdtypes/commrec.h"
66 #include "gromacs/mdtypes/forceoutput.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/md_enums.h"
69 #include "gromacs/mdtypes/state.h"
70 #include "gromacs/pbcutil/pbc.h"
71 #include "gromacs/topology/mtop_lookup.h"
72 #include "gromacs/topology/topology.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/exceptions.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/futil.h"
77 #include "gromacs/utility/gmxassert.h"
78 #include "gromacs/utility/pleasecite.h"
79 #include "gromacs/utility/real.h"
80 #include "gromacs/utility/smalloc.h"
81 #include "gromacs/utility/strconvert.h"
83 #include "pull_internal.h"
87 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
90 static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
92 if (params.ind.size() <= 1)
97 else if (params.pbcatom >= 0)
99 if (setPbcRefToPrevStepCOM)
101 return epgrppbcPREVSTEPCOM;
105 return epgrppbcREFAT;
114 /* NOTE: The params initialization currently copies pointers.
115 * So the lifetime of the source, currently always inputrec,
116 * should not end before that of this object.
117 * This will be fixed when the pointers are replacted by std::vector.
119 pull_group_work_t::pull_group_work_t(const t_pull_group& params,
120 gmx::LocalAtomSet atomSet,
121 bool bSetPbcRefToPrevStepCOM) :
123 epgrppbc(groupPbcFromParams(params, bSetPbcRefToPrevStepCOM)),
124 needToCalcCom(false),
134 bool pull_coordinate_is_angletype(const t_pull_coord* pcrd)
136 return (pcrd->eGeom == PullGroupGeometry::Angle || pcrd->eGeom == PullGroupGeometry::Dihedral
137 || pcrd->eGeom == PullGroupGeometry::AngleAxis);
140 static bool pull_coordinate_is_directional(const t_pull_coord* pcrd)
142 return (pcrd->eGeom == PullGroupGeometry::Direction || pcrd->eGeom == PullGroupGeometry::DirectionPBC
143 || pcrd->eGeom == PullGroupGeometry::DirectionRelative
144 || pcrd->eGeom == PullGroupGeometry::Cylinder);
147 const char* pull_coordinate_units(const t_pull_coord* pcrd)
149 return pull_coordinate_is_angletype(pcrd) ? "deg" : "nm";
152 double pull_conversion_factor_userinput2internal(const t_pull_coord* pcrd)
154 if (pull_coordinate_is_angletype(pcrd))
164 double pull_conversion_factor_internal2userinput(const t_pull_coord* pcrd)
166 if (pull_coordinate_is_angletype(pcrd))
176 /* Apply forces in a mass weighted fashion for part of the pull group */
177 static void apply_forces_grp_part(const pull_group_work_t* pgrp,
185 double inv_wm = pgrp->mwscale;
187 auto localAtomIndices = pgrp->atomSet.localIndex();
188 for (int i = ind_start; i < ind_end; i++)
190 int ii = localAtomIndices[i];
191 double wmass = masses[ii];
192 if (!pgrp->localWeights.empty())
194 wmass *= pgrp->localWeights[i];
197 for (int d = 0; d < DIM; d++)
199 f[ii][d] += sign * wmass * f_pull[d] * inv_wm;
204 /* Apply forces in a mass weighted fashion */
205 static void apply_forces_grp(const pull_group_work_t* pgrp,
212 auto localAtomIndices = pgrp->atomSet.localIndex();
214 if (pgrp->params.ind.size() == 1 && pgrp->atomSet.numAtomsLocal() == 1)
216 /* Only one atom and our rank has this atom: we can skip
217 * the mass weighting, which means that this code also works
218 * for mass=0, e.g. with a virtual site.
220 for (int d = 0; d < DIM; d++)
222 f[localAtomIndices[0]][d] += sign * f_pull[d];
227 if (localAtomIndices.size() <= c_pullMaxNumLocalAtomsSingleThreaded)
229 apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), masses, f_pull, sign, f);
233 #pragma omp parallel for num_threads(nthreads) schedule(static)
234 for (int th = 0; th < nthreads; th++)
236 int ind_start = (localAtomIndices.size() * (th + 0)) / nthreads;
237 int ind_end = (localAtomIndices.size() * (th + 1)) / nthreads;
238 apply_forces_grp_part(pgrp, ind_start, ind_end, masses, f_pull, sign, f);
244 /* Apply forces in a mass weighted fashion to a cylinder group */
245 static void apply_forces_cyl_grp(const pull_group_work_t* pgrp,
246 const double dv_corr,
252 int gmx_unused nthreads)
254 double inv_wm = pgrp->mwscale;
256 auto localAtomIndices = pgrp->atomSet.localIndex();
258 /* The cylinder group is always a slab in the system, thus large.
259 * Therefore we always thread-parallelize this group.
261 int numAtomsLocal = localAtomIndices.size();
262 #pragma omp parallel for num_threads(nthreads) schedule(static)
263 for (int i = 0; i < numAtomsLocal; i++)
265 double weight = pgrp->localWeights[i];
270 int ii = localAtomIndices[i];
271 double mass = masses[ii];
272 /* The stored axial distance from the cylinder center (dv) needs
273 * to be corrected for an offset (dv_corr), which was unknown when
276 double dv_com = pgrp->dv[i] + dv_corr;
278 /* Here we not only add the pull force working along vec (f_pull),
279 * but also a radial component, due to the dependence of the weights
280 * on the radial distance.
282 for (int m = 0; m < DIM; m++)
284 f[ii][m] += sign * inv_wm * (mass * weight * f_pull[m] + pgrp->mdw[i][m] * dv_com * f_scal);
289 /* Apply torque forces in a mass weighted fashion to the groups that define
290 * the pull vector direction for pull coordinate pcrd.
292 static void apply_forces_vec_torque(const struct pull_t* pull,
293 const pull_coord_work_t* pcrd,
297 const PullCoordSpatialData& spatialData = pcrd->spatialData;
299 /* The component inpr along the pull vector is accounted for in the usual
300 * way. Here we account for the component perpendicular to vec.
303 for (int m = 0; m < DIM; m++)
305 inpr += spatialData.dr01[m] * spatialData.vec[m];
307 /* The torque force works along the component of the distance vector
308 * of between the two "usual" pull groups that is perpendicular to
309 * the pull vector. The magnitude of this force is the "usual" scale force
310 * multiplied with the ratio of the distance between the two "usual" pull
311 * groups and the distance between the two groups that define the vector.
314 for (int m = 0; m < DIM; m++)
316 f_perp[m] = (spatialData.dr01[m] - inpr * spatialData.vec[m]) / spatialData.vec_len * pcrd->scalarForce;
319 /* Apply the force to the groups defining the vector using opposite signs */
320 apply_forces_grp(&pull->group[pcrd->params.group[2]], masses, f_perp, -1, f, pull->nthreads);
321 apply_forces_grp(&pull->group[pcrd->params.group[3]], masses, f_perp, 1, f, pull->nthreads);
324 /* Apply forces in a mass weighted fashion */
325 static void apply_forces_coord(struct pull_t* pull,
327 const PullCoordVectorForces& forces,
331 /* Here it would be more efficient to use one large thread-parallel
332 * region instead of potential parallel regions within apply_forces_grp.
333 * But there could be overlap between pull groups and this would lead
337 const pull_coord_work_t& pcrd = pull->coord[coord];
339 if (pcrd.params.eGeom == PullGroupGeometry::Cylinder)
341 apply_forces_cyl_grp(&pull->dyna[coord],
342 pcrd.spatialData.cyl_dev,
350 /* Sum the force along the vector and the radial force */
352 for (int m = 0; m < DIM; m++)
354 f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
356 apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, f_tot, 1, f, pull->nthreads);
360 if (pcrd.params.eGeom == PullGroupGeometry::DirectionRelative)
362 /* We need to apply the torque forces to the pull groups
363 * that define the pull vector.
365 apply_forces_vec_torque(pull, &pcrd, masses, f);
368 if (!pull->group[pcrd.params.group[0]].params.ind.empty())
371 &pull->group[pcrd.params.group[0]], masses, forces.force01, -1, f, pull->nthreads);
373 apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, forces.force01, 1, f, pull->nthreads);
375 if (pcrd.params.ngroup >= 4)
378 &pull->group[pcrd.params.group[2]], masses, forces.force23, -1, f, pull->nthreads);
379 apply_forces_grp(&pull->group[pcrd.params.group[3]], masses, forces.force23, 1, f, pull->nthreads);
381 if (pcrd.params.ngroup >= 6)
384 &pull->group[pcrd.params.group[4]], masses, forces.force45, -1, f, pull->nthreads);
385 apply_forces_grp(&pull->group[pcrd.params.group[5]], masses, forces.force45, 1, f, pull->nthreads);
390 real max_pull_distance2(const pull_coord_work_t* pcrd, const t_pbc* pbc)
392 /* Note that this maximum distance calculation is more complex than
393 * most other cases in GROMACS, since here we have to take care of
394 * distance calculations that don't involve all three dimensions.
395 * For example, we can use distances that are larger than the
396 * box X and Y dimensions for a box that is elongated along Z.
399 real max_d2 = GMX_REAL_MAX;
401 if (pull_coordinate_is_directional(&pcrd->params))
403 /* Directional pulling along along direction pcrd->vec.
404 * Calculating the exact maximum distance is complex and bug-prone.
405 * So we take a safe approach by not allowing distances that
406 * are larger than half the distance between unit cell faces
407 * along dimensions involved in pcrd->vec.
409 for (int m = 0; m < DIM; m++)
411 if (m < pbc->ndim_ePBC && pcrd->spatialData.vec[m] != 0)
413 real imageDistance2 = gmx::square(pbc->box[m][m]);
414 for (int d = m + 1; d < DIM; d++)
416 imageDistance2 -= gmx::square(pbc->box[d][m]);
418 max_d2 = std::min(max_d2, imageDistance2);
424 /* Distance pulling along dimensions with pcrd->params.dim[d]==1.
425 * We use half the minimum box vector length of the dimensions involved.
426 * This is correct for all cases, except for corner cases with
427 * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
428 * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
429 * in such corner cases the user could get correct results,
430 * depending on the details of the setup, we avoid further
431 * code complications.
433 for (int m = 0; m < DIM; m++)
435 if (m < pbc->ndim_ePBC && pcrd->params.dim[m] != 0)
437 real imageDistance2 = gmx::square(pbc->box[m][m]);
438 for (int d = 0; d < m; d++)
440 if (pcrd->params.dim[d] != 0)
442 imageDistance2 += gmx::square(pbc->box[m][d]);
445 max_d2 = std::min(max_d2, imageDistance2);
450 return 0.25 * max_d2;
453 /* This function returns the distance based on coordinates xg and xref.
454 * Note that the pull coordinate struct pcrd is not modified.
456 * \param[in] pull The pull struct
457 * \param[in] pcrd The pull coordinate to compute a distance for
458 * \param[in] pbc The periodic boundary conditions
459 * \param[in] xg The coordinate of group 1
460 * \param[in] xref The coordinate of group 0
461 * \param[in] groupIndex0 The index of group 0 in the pcrd->params.group
462 * \param[in] groupIndex1 The index of group 1 in the pcrd->params.group
463 * \param[in] max_dist2 The maximum distance squared
464 * \param[out] dr The distance vector
466 static void low_get_pull_coord_dr(const struct pull_t* pull,
467 const pull_coord_work_t* pcrd,
471 const int groupIndex0,
472 const int groupIndex1,
473 const double max_dist2,
476 const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
478 /* Only the first group can be an absolute reference, in that case nat=0 */
479 if (pgrp0->params.ind.empty())
481 for (int m = 0; m < DIM; m++)
483 xref[m] = pcrd->params.origin[m];
488 copy_dvec(xref, xrefr);
490 dvec dref = { 0, 0, 0 };
491 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
493 for (int m = 0; m < DIM; m++)
495 dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
497 /* Add the reference position, so we use the correct periodic image */
498 dvec_inc(xrefr, dref);
501 pbc_dx_d(pbc, xg, xrefr, dr);
503 bool directional = pull_coordinate_is_directional(&pcrd->params);
505 for (int m = 0; m < DIM; m++)
507 dr[m] *= pcrd->params.dim[m];
508 if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
510 dr2 += dr[m] * dr[m];
513 /* Check if we are close to switching to another periodic image */
514 if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
516 /* Note that technically there is no issue with switching periodic
517 * image, as pbc_dx_d returns the distance to the closest periodic
518 * image. However in all cases where periodic image switches occur,
519 * the pull results will be useless in practice.
522 "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
523 "box size (%f).\n%s",
524 pcrd->params.group[groupIndex0],
525 pcrd->params.group[groupIndex1],
527 sqrt(0.98 * 0.98 * max_dist2),
528 pcrd->params.eGeom == PullGroupGeometry::Direction
529 ? "You might want to consider using \"pull-geometry = "
530 "direction-periodic\" instead.\n"
534 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
540 /* This function returns the distance based on the contents of the pull struct.
541 * pull->coord[coord_ind].dr, and potentially vector, are updated.
543 static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
545 pull_coord_work_t* pcrd = &pull->coord[coord_ind];
546 PullCoordSpatialData& spatialData = pcrd->spatialData;
549 /* With AWH pulling we allow for periodic pulling with geometry=direction.
550 * TODO: Store a periodicity flag instead of checking for external pull provider.
552 if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC
553 || (pcrd->params.eGeom == PullGroupGeometry::Direction
554 && pcrd->params.eType == PullingAlgorithm::External))
560 md2 = static_cast<double>(max_pull_distance2(pcrd, pbc));
563 if (pcrd->params.eGeom == PullGroupGeometry::DirectionRelative)
565 /* We need to determine the pull vector */
569 const pull_group_work_t& pgrp2 = pull->group[pcrd->params.group[2]];
570 const pull_group_work_t& pgrp3 = pull->group[pcrd->params.group[3]];
572 pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
574 for (m = 0; m < DIM; m++)
576 vec[m] *= pcrd->params.dim[m];
578 spatialData.vec_len = dnorm(vec);
579 for (m = 0; m < DIM; m++)
581 spatialData.vec[m] = vec[m] / spatialData.vec_len;
586 "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
593 spatialData.vec[ZZ]);
597 pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
598 pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
600 low_get_pull_coord_dr(
605 pcrd->params.eGeom == PullGroupGeometry::Cylinder ? pull->dyna[coord_ind].x : pgrp0->x,
611 if (pcrd->params.ngroup >= 4)
613 pull_group_work_t *pgrp2, *pgrp3;
614 pgrp2 = &pull->group[pcrd->params.group[2]];
615 pgrp3 = &pull->group[pcrd->params.group[3]];
617 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, 2, 3, md2, spatialData.dr23);
619 if (pcrd->params.ngroup >= 6)
621 pull_group_work_t *pgrp4, *pgrp5;
622 pgrp4 = &pull->group[pcrd->params.group[4]];
623 pgrp5 = &pull->group[pcrd->params.group[5]];
625 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, 4, 5, md2, spatialData.dr45);
629 /* Modify x so that it is periodic in [-pi, pi)
630 * It is assumed that x is in [-3pi, 3pi) so that x
631 * needs to be shifted by at most one period.
633 static void make_periodic_2pi(double* x)
645 /* This function should always be used to modify pcrd->value_ref */
646 static void low_set_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double value_ref)
648 GMX_ASSERT(pcrd->params.eType != PullingAlgorithm::External,
649 "The pull coord reference value should not be used with type external-potential");
651 if (pcrd->params.eGeom == PullGroupGeometry::Distance)
656 "Pull reference distance for coordinate %d (%f) needs to be non-negative",
661 else if (pcrd->params.eGeom == PullGroupGeometry::Angle
662 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
664 if (value_ref < 0 || value_ref > M_PI)
667 "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
668 "interval [0,180] deg",
670 value_ref * pull_conversion_factor_internal2userinput(&pcrd->params));
673 else if (pcrd->params.eGeom == PullGroupGeometry::Dihedral)
675 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
676 make_periodic_2pi(&value_ref);
679 pcrd->value_ref = value_ref;
682 static void update_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double t)
684 /* With zero rate the reference value is set initially and doesn't change */
685 if (pcrd->params.rate != 0)
687 double value_ref = (pcrd->params.init + pcrd->params.rate * t)
688 * pull_conversion_factor_userinput2internal(&pcrd->params);
689 low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
693 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
694 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
697 dvec dr32; /* store instead of dr23? */
699 dsvmul(-1, spatialData->dr23, dr32);
700 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
701 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
702 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
704 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
705 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
706 * Note 2: the angle between the plane normal vectors equals pi only when
707 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
708 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
710 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
714 /* Calculates pull->coord[coord_ind].value.
715 * This function also updates pull->coord[coord_ind].dr.
717 static void get_pull_coord_distance(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
719 pull_coord_work_t* pcrd;
722 pcrd = &pull->coord[coord_ind];
724 get_pull_coord_dr(pull, coord_ind, pbc);
726 PullCoordSpatialData& spatialData = pcrd->spatialData;
728 switch (pcrd->params.eGeom)
730 case PullGroupGeometry::Distance:
731 /* Pull along the vector between the com's */
732 spatialData.value = dnorm(spatialData.dr01);
734 case PullGroupGeometry::Direction:
735 case PullGroupGeometry::DirectionPBC:
736 case PullGroupGeometry::DirectionRelative:
737 case PullGroupGeometry::Cylinder:
739 spatialData.value = 0;
740 for (m = 0; m < DIM; m++)
742 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
745 case PullGroupGeometry::Angle:
746 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
748 case PullGroupGeometry::Dihedral:
749 spatialData.value = get_dihedral_angle_coord(&spatialData);
751 case PullGroupGeometry::AngleAxis:
752 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
754 default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
758 /* Returns the deviation from the reference value.
759 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
761 static double get_pull_coord_deviation(struct pull_t* pull, int coord_ind, const t_pbc* pbc, double t)
763 pull_coord_work_t* pcrd;
766 pcrd = &pull->coord[coord_ind];
768 /* Update the reference value before computing the distance,
769 * since it is used in the distance computation with periodic pulling.
771 update_pull_coord_reference_value(pcrd, coord_ind, t);
773 get_pull_coord_distance(pull, coord_ind, pbc);
775 get_pull_coord_distance(pull, coord_ind, pbc);
777 /* Determine the deviation */
778 dev = pcrd->spatialData.value - pcrd->value_ref;
780 /* Check that values are allowed */
781 if (pcrd->params.eGeom == PullGroupGeometry::Distance && pcrd->spatialData.value == 0)
783 /* With no vector we can not determine the direction for the force,
784 * so we set the force to zero.
788 else if (pcrd->params.eGeom == PullGroupGeometry::Dihedral)
790 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
791 Thus, the unwrapped deviation is here in (-2pi, 2pi].
792 After making it periodic, the deviation will be in [-pi, pi). */
793 make_periodic_2pi(&dev);
799 double get_pull_coord_value(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
801 get_pull_coord_distance(pull, coord_ind, pbc);
803 return pull->coord[coord_ind].spatialData.value;
806 void clear_pull_forces(pull_t* pull)
808 /* Zeroing the forces is only required for constraint pulling.
809 * It can happen that multiple constraint steps need to be applied
810 * and therefore the constraint forces need to be accumulated.
812 for (pull_coord_work_t& coord : pull->coord)
814 coord.scalarForce = 0;
818 /* Apply constraint using SHAKE */
820 do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
823 dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
824 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
825 dvec* rnew; /* current 'new' positions of the groups */
826 double* dr_tot; /* the total update of the coords */
829 double lambda, rm, invdt = 0;
830 gmx_bool bConverged_all, bConverged = FALSE;
831 int niter = 0, ii, j, m, max_iter = 100;
835 snew(r_ij, pull->coord.size());
836 snew(dr_tot, pull->coord.size());
838 snew(rnew, pull->group.size());
840 /* copy the current unconstrained positions for use in iterations. We
841 iterate until rinew[i] and rjnew[j] obey the constraints. Then
842 rinew - pull.x_unc[i] is the correction dr to group i */
843 for (size_t g = 0; g < pull->group.size(); g++)
845 copy_dvec(pull->group[g].xp, rnew[g]);
848 /* Determine the constraint directions from the old positions */
849 for (size_t c = 0; c < pull->coord.size(); c++)
851 pull_coord_work_t* pcrd;
853 pcrd = &pull->coord[c];
855 if (pcrd->params.eType != PullingAlgorithm::Constraint)
860 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
861 * We don't modify dr and value anymore, so these values are also used
864 get_pull_coord_distance(pull, c, pbc);
866 const PullCoordSpatialData& spatialData = pcrd->spatialData;
870 "Pull coord %zu dr %f %f %f\n",
872 spatialData.dr01[XX],
873 spatialData.dr01[YY],
874 spatialData.dr01[ZZ]);
877 if (pcrd->params.eGeom == PullGroupGeometry::Direction
878 || pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
880 /* Select the component along vec */
882 for (m = 0; m < DIM; m++)
884 a += spatialData.vec[m] * spatialData.dr01[m];
886 for (m = 0; m < DIM; m++)
888 r_ij[c][m] = a * spatialData.vec[m];
893 copy_dvec(spatialData.dr01, r_ij[c]);
896 if (dnorm2(r_ij[c]) == 0)
899 "Distance for pull coordinate %zu is zero with constraint pulling, which is "
905 bConverged_all = FALSE;
906 while (!bConverged_all && niter < max_iter)
908 bConverged_all = TRUE;
910 /* loop over all constraints */
911 for (size_t c = 0; c < pull->coord.size(); c++)
913 pull_coord_work_t* pcrd;
914 pull_group_work_t *pgrp0, *pgrp1;
917 pcrd = &pull->coord[c];
919 if (pcrd->params.eType != PullingAlgorithm::Constraint)
924 update_pull_coord_reference_value(pcrd, c, t);
926 pgrp0 = &pull->group[pcrd->params.group[0]];
927 pgrp1 = &pull->group[pcrd->params.group[1]];
929 /* Get the current difference vector */
930 low_get_pull_coord_dr(
931 pull, pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
935 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
938 rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
940 switch (pcrd->params.eGeom)
942 case PullGroupGeometry::Distance:
943 if (pcrd->value_ref <= 0)
947 "The pull constraint reference distance for group %zu is <= 0 (%f)",
953 double q, c_a, c_b, c_c;
955 c_a = diprod(r_ij[c], r_ij[c]);
956 c_b = diprod(unc_ij, r_ij[c]) * 2;
957 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
961 q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
966 q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
972 fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b, c_c, lambda);
976 /* The position corrections dr due to the constraints */
977 dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
978 dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
979 dr_tot[c] += -lambda * dnorm(r_ij[c]);
981 case PullGroupGeometry::Direction:
982 case PullGroupGeometry::DirectionPBC:
983 case PullGroupGeometry::Cylinder:
984 /* A 1-dimensional constraint along a vector */
986 for (m = 0; m < DIM; m++)
988 vec[m] = pcrd->spatialData.vec[m];
989 a += unc_ij[m] * vec[m];
991 /* Select only the component along the vector */
992 dsvmul(a, vec, unc_ij);
993 lambda = a - pcrd->value_ref;
996 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
999 /* The position corrections dr due to the constraints */
1000 dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
1001 dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
1002 dr_tot[c] += -lambda;
1004 default: gmx_incons("Invalid enumeration value for eGeom");
1012 g0 = pcrd->params.group[0];
1013 g1 = pcrd->params.group[1];
1014 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
1015 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
1017 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1026 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
1035 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1045 /* Update the COMs with dr */
1046 dvec_inc(rnew[pcrd->params.group[1]], dr1);
1047 dvec_inc(rnew[pcrd->params.group[0]], dr0);
1050 /* Check if all constraints are fullfilled now */
1051 for (pull_coord_work_t& coord : pull->coord)
1053 if (coord.params.eType != PullingAlgorithm::Constraint)
1058 low_get_pull_coord_dr(
1059 pull, &coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
1061 switch (coord.params.eGeom)
1063 case PullGroupGeometry::Distance:
1064 bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1066 case PullGroupGeometry::Direction:
1067 case PullGroupGeometry::DirectionPBC:
1068 case PullGroupGeometry::Cylinder:
1069 for (m = 0; m < DIM; m++)
1071 vec[m] = coord.spatialData.vec[m];
1073 inpr = diprod(unc_ij, vec);
1074 dsvmul(inpr, vec, unc_ij);
1075 bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1079 gmx::formatString("Geometry %s not handled here",
1080 enumValueToString(coord.params.eGeom))
1089 "Pull constraint not converged: "
1091 "d_ref = %f, current d = %f\n",
1092 coord.params.group[0],
1093 coord.params.group[1],
1098 bConverged_all = FALSE;
1103 /* if after all constraints are dealt with and bConverged is still TRUE
1104 we're finished, if not we do another iteration */
1106 if (niter > max_iter)
1108 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1111 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1118 /* update atoms in the groups */
1119 for (size_t g = 0; g < pull->group.size(); g++)
1121 const pull_group_work_t* pgrp;
1124 pgrp = &pull->group[g];
1126 /* get the final constraint displacement dr for group g */
1127 dvec_sub(rnew[g], pgrp->xp, dr);
1129 if (dnorm2(dr) == 0)
1131 /* No displacement, no update necessary */
1135 /* update the atom positions */
1136 auto localAtomIndices = pgrp->atomSet.localIndex();
1138 for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1140 ii = localAtomIndices[j];
1141 if (!pgrp->localWeights.empty())
1143 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1145 for (m = 0; m < DIM; m++)
1151 for (m = 0; m < DIM; m++)
1153 v[ii][m] += invdt * tmp[m];
1159 /* calculate the constraint forces, used for output and virial only */
1160 for (size_t c = 0; c < pull->coord.size(); c++)
1162 pull_coord_work_t* pcrd;
1164 pcrd = &pull->coord[c];
1166 if (pcrd->params.eType != PullingAlgorithm::Constraint)
1171 /* Accumulate the forces, in case we have multiple constraint steps */
1174 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1176 pcrd->scalarForce += force;
1178 if (vir != nullptr && pcrd->params.eGeom != PullGroupGeometry::DirectionPBC && bMaster)
1182 /* Add the pull contribution to the virial */
1183 /* We have already checked above that r_ij[c] != 0 */
1184 f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1186 for (j = 0; j < DIM; j++)
1188 for (m = 0; m < DIM; m++)
1190 vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1196 /* finished! I hope. Give back some memory */
1202 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1204 for (int j = 0; j < DIM; j++)
1206 for (int m = 0; m < DIM; m++)
1208 vir[j][m] -= 0.5 * f[j] * dr[m];
1213 /* Adds the pull contribution to the virial */
1214 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1216 if (vir != nullptr && pcrd.params.eGeom != PullGroupGeometry::DirectionPBC)
1218 /* Add the pull contribution for each distance vector to the virial. */
1219 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1220 if (pcrd.params.ngroup >= 4)
1222 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1224 if (pcrd.params.ngroup >= 6)
1226 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1231 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1239 k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1240 dkdl = pcrd->params.kB - pcrd->params.k;
1242 switch (pcrd->params.eType)
1244 case PullingAlgorithm::Umbrella:
1245 case PullingAlgorithm::FlatBottom:
1246 case PullingAlgorithm::FlatBottomHigh:
1247 /* The only difference between an umbrella and a flat-bottom
1248 * potential is that a flat-bottom is zero above or below
1249 the reference value.
1251 if ((pcrd->params.eType == PullingAlgorithm::FlatBottom && dev < 0)
1252 || (pcrd->params.eType == PullingAlgorithm::FlatBottomHigh && dev > 0))
1257 pcrd->scalarForce = -k * dev;
1258 *V += 0.5 * k * gmx::square(dev);
1259 *dVdl += 0.5 * dkdl * gmx::square(dev);
1261 case PullingAlgorithm::ConstantForce:
1262 pcrd->scalarForce = -k;
1263 *V += k * pcrd->spatialData.value;
1264 *dVdl += dkdl * pcrd->spatialData.value;
1266 case PullingAlgorithm::External:
1268 "the scalar pull force should not be calculated internally for pull type "
1270 default: gmx_incons("Unsupported pull type in do_pull_pot");
1274 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1276 const t_pull_coord& params = pcrd.params;
1277 const PullCoordSpatialData& spatialData = pcrd.spatialData;
1279 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1280 PullCoordVectorForces forces;
1282 if (params.eGeom == PullGroupGeometry::Distance)
1284 double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1285 for (int m = 0; m < DIM; m++)
1287 forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1290 else if (params.eGeom == PullGroupGeometry::Angle)
1293 double cos_theta, cos_theta2;
1295 cos_theta = cos(spatialData.value);
1296 cos_theta2 = gmx::square(cos_theta);
1298 /* The force at theta = 0, pi is undefined so we should not apply any force.
1299 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1300 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1301 * have to check for this before dividing by their norm below.
1305 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1306 * and another vector parallel to dr23:
1307 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1308 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1310 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1311 double b = a * cos_theta;
1312 double invdr01 = 1. / dnorm(spatialData.dr01);
1313 double invdr23 = 1. / dnorm(spatialData.dr23);
1314 dvec normalized_dr01, normalized_dr23;
1315 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1316 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1318 for (int m = 0; m < DIM; m++)
1320 /* Here, f_scal is -dV/dtheta */
1322 pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1324 pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1329 /* No forces to apply for ill-defined cases*/
1330 clear_dvec(forces.force01);
1331 clear_dvec(forces.force23);
1334 else if (params.eGeom == PullGroupGeometry::AngleAxis)
1336 double cos_theta, cos_theta2;
1338 /* The angle-axis force is exactly the same as the angle force (above) except that in
1339 this case the second vector (dr23) is replaced by the pull vector. */
1340 cos_theta = cos(spatialData.value);
1341 cos_theta2 = gmx::square(cos_theta);
1347 dvec normalized_dr01;
1349 invdr01 = 1. / dnorm(spatialData.dr01);
1350 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1351 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1354 for (int m = 0; m < DIM; m++)
1357 pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1362 clear_dvec(forces.force01);
1365 else if (params.eGeom == PullGroupGeometry::Dihedral)
1367 double m2, n2, tol, sqrdist_32;
1369 /* Note: there is a small difference here compared to the
1370 dihedral force calculations in the bondeds (ref: Bekker 1994).
1371 There rij = ri - rj, while here dr01 = r1 - r0.
1372 However, all distance vectors occur in form of cross or inner products
1373 so that two signs cancel and we end up with the same expressions.
1374 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1376 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1377 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1378 dsvmul(-1, spatialData.dr23, dr32);
1379 sqrdist_32 = diprod(dr32, dr32);
1380 tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1381 if ((m2 > tol) && (n2 > tol))
1383 double a_01, a_23_01, a_23_45, a_45;
1384 double inv_dist_32, inv_sqrdist_32, dist_32;
1386 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1387 inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1388 dist_32 = sqrdist_32 * inv_dist_32;
1390 /* Forces on groups 0, 1 */
1391 a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
1392 dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1394 /* Forces on groups 4, 5 */
1395 a_45 = -pcrd.scalarForce * dist_32 / n2;
1396 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1398 /* Force on groups 2, 3 (defining the axis) */
1399 a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1400 a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1401 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1402 dsvmul(a_23_45, forces.force45, v);
1403 dvec_sub(u, v, forces.force23); /* force on group 3 */
1407 /* No force to apply for ill-defined cases */
1408 clear_dvec(forces.force01);
1409 clear_dvec(forces.force23);
1410 clear_dvec(forces.force45);
1415 for (int m = 0; m < DIM; m++)
1417 forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1425 /* We use a global mutex for locking access to the pull data structure
1426 * during registration of external pull potential providers.
1427 * We could use a different, local mutex for each pull object, but the overhead
1428 * is extremely small here and registration is only done during initialization.
1430 static std::mutex registrationMutex;
1432 using Lock = std::lock_guard<std::mutex>;
1434 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1436 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1437 GMX_RELEASE_ASSERT(provider != nullptr,
1438 "register_external_pull_potential called with NULL as provider name");
1440 if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1443 "Module '%s' attempted to register an external potential for pull coordinate %d "
1444 "which is out of the pull coordinate range %d - %zu\n",
1448 pull->coord.size());
1451 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1453 if (pcrd->params.eType != PullingAlgorithm::External)
1457 "Module '%s' attempted to register an external potential for pull coordinate %d "
1458 "which of type '%s', whereas external potentials are only supported with type '%s'",
1461 enumValueToString(pcrd->params.eType),
1462 enumValueToString(PullingAlgorithm::External));
1465 GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
1466 "The external potential provider string for a pull coordinate is NULL");
1468 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
1471 "Module '%s' attempted to register an external potential for pull coordinate %d "
1472 "which expects the external potential to be provided by a module named '%s'",
1475 pcrd->params.externalPotentialProvider.c_str());
1478 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1479 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1480 * pull->numUnregisteredExternalPotentials.
1482 Lock registrationLock(registrationMutex);
1484 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1487 "Module '%s' attempted to register an external potential for pull coordinate %d "
1493 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1494 pull->numUnregisteredExternalPotentials--;
1496 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1497 "Negative unregistered potentials, the pull code is inconsistent");
1501 static void check_external_potential_registration(const struct pull_t* pull)
1503 if (pull->numUnregisteredExternalPotentials > 0)
1506 for (c = 0; c < pull->coord.size(); c++)
1508 if (pull->coord[c].params.eType == PullingAlgorithm::External
1509 && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1515 GMX_RELEASE_ASSERT(c < pull->coord.size(),
1516 "Internal inconsistency in the pull potential provider counting");
1519 "No external provider for external pull potentials have been provided for %d "
1520 "pull coordinates. The first coordinate without provider is number %zu, which "
1521 "expects a module named '%s' to provide the external potential.",
1522 pull->numUnregisteredExternalPotentials,
1524 pull->coord[c].params.externalPotentialProvider.c_str());
1529 /* Pull takes care of adding the forces of the external potential.
1530 * The external potential module has to make sure that the corresponding
1531 * potential energy is added either to the pull term or to a term
1532 * specific to the external module.
1534 void apply_external_pull_coord_force(struct pull_t* pull,
1538 gmx::ForceWithVirial* forceWithVirial)
1540 pull_coord_work_t* pcrd;
1542 GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1543 "apply_external_pull_coord_force called with coord_index out of range");
1545 if (pull->comm.bParticipate)
1547 pcrd = &pull->coord[coord_index];
1550 pcrd->params.eType == PullingAlgorithm::External,
1551 "The pull force can only be set externally on pull coordinates of external type");
1553 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1554 "apply_external_pull_coord_force called for an unregistered pull coordinate");
1557 pcrd->scalarForce = coord_force;
1559 /* Calculate the forces on the pull groups */
1560 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1562 /* Add the forces for this coordinate to the total virial and force */
1563 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1565 matrix virial = { { 0 } };
1566 add_virial_coord(virial, *pcrd, pullCoordForces);
1567 forceWithVirial->addVirialContribution(virial);
1571 pull, coord_index, pullCoordForces, masses, as_rvec_array(forceWithVirial->force_.data()));
1574 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1577 /* Calculate the pull potential and scalar force for a pull coordinate.
1578 * Returns the vector forces for the pull coordinate.
1580 static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull,
1589 pull_coord_work_t& pcrd = pull->coord[coord_ind];
1591 assert(pcrd.params.eType != PullingAlgorithm::Constraint);
1593 double dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
1595 calc_pull_coord_scalar_force_and_potential(&pcrd, dev, lambda, V, dVdl);
1597 PullCoordVectorForces pullCoordForces = calculateVectorForces(pcrd);
1599 add_virial_coord(vir, pcrd, pullCoordForces);
1601 return pullCoordForces;
1604 real pull_potential(struct pull_t* pull,
1607 const t_commrec* cr,
1611 gmx::ForceWithVirial* force,
1616 assert(pull != nullptr);
1618 /* Ideally we should check external potential registration only during
1619 * the initialization phase, but that requires another function call
1620 * that should be done exactly in the right order. So we check here.
1622 check_external_potential_registration(pull);
1624 if (pull->comm.bParticipate)
1628 pull_calc_coms(cr, pull, masses, pbc, t, x, nullptr);
1630 rvec* f = as_rvec_array(force->force_.data());
1631 matrix virial = { { 0 } };
1632 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1633 for (size_t c = 0; c < pull->coord.size(); c++)
1635 pull_coord_work_t* pcrd;
1636 pcrd = &pull->coord[c];
1638 /* For external potential the force is assumed to be given by an external module by a
1639 call to apply_pull_coord_external_force */
1640 if (pcrd->params.eType == PullingAlgorithm::Constraint
1641 || pcrd->params.eType == PullingAlgorithm::External)
1646 PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1647 pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1649 /* Distribute the force over the atoms in the pulled groups */
1650 apply_forces_coord(pull, c, pullCoordForces, masses, f);
1655 force->addVirialContribution(virial);
1660 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1661 "Too few or too many external pull potentials have been applied the previous step");
1662 /* All external pull potentials still need to be applied */
1663 pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1665 return (MASTER(cr) ? V : 0.0);
1668 void pull_constraint(struct pull_t* pull,
1671 const t_commrec* cr,
1679 assert(pull != nullptr);
1681 if (pull->comm.bParticipate)
1683 pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
1685 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1689 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1693 gmx_bool bMustParticipate;
1699 /* We always make the master node participate, such that it can do i/o,
1700 * add the virial and to simplify MC type extensions people might have.
1702 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1704 for (pull_group_work_t& group : pull->group)
1706 if (!group.globalWeights.empty())
1708 group.localWeights.resize(group.atomSet.numAtomsLocal());
1709 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1711 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1715 GMX_ASSERT(bMustParticipate || dd != nullptr,
1716 "Either all ranks (including this rank) participate, or we use DD and need to "
1717 "have access to dd here");
1719 /* We should participate if we have pull or pbc atoms */
1720 if (!bMustParticipate
1721 && (group.atomSet.numAtomsLocal() > 0
1722 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1724 bMustParticipate = TRUE;
1728 if (!comm->bParticipateAll)
1730 /* Keep currently not required ranks in the communicator
1731 * if they needed to participate up to 20 decompositions ago.
1732 * This avoids frequent rebuilds due to atoms jumping back and forth.
1734 const int64_t history_count = 20;
1735 gmx_bool bWillParticipate;
1738 /* Increase the decomposition counter for the current call */
1739 comm->setup_count++;
1741 if (bMustParticipate)
1743 comm->must_count = comm->setup_count;
1748 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1750 if (debug && dd != nullptr)
1753 "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
1755 gmx::boolToString(bMustParticipate),
1756 gmx::boolToString(bWillParticipate));
1759 if (bWillParticipate)
1761 /* Count the number of ranks that we want to have participating */
1763 /* Count the number of ranks that need to be added */
1764 count[1] = comm->bParticipate ? 0 : 1;
1772 /* The cost of this global operation will be less that the cost
1773 * of the extra MPI_Comm_split calls that we can avoid.
1775 gmx_sumi(2, count, cr);
1777 /* If we are missing ranks or if we have 20% more ranks than needed
1778 * we make a new sub-communicator.
1780 if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1784 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1788 if (comm->mpi_comm_com != MPI_COMM_NULL)
1790 MPI_Comm_free(&comm->mpi_comm_com);
1793 /* This might be an extremely expensive operation, so we try
1794 * to avoid this splitting as much as possible.
1796 assert(dd != nullptr);
1797 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1800 comm->bParticipate = bWillParticipate;
1801 comm->nparticipate = count[0];
1803 /* When we use the previous COM for PBC, we need to broadcast
1804 * the previous COM to ranks that have joined the communicator.
1806 for (pull_group_work_t& group : pull->group)
1808 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1810 GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1811 "The master rank has to participate, as it should pass an up to "
1813 "to bcast here as well as to e.g. checkpointing");
1815 gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
1821 /* Since the PBC of atoms might have changed, we need to update the PBC */
1822 pull->bSetPBCatoms = TRUE;
1825 static void init_pull_group_index(FILE* fplog,
1826 const t_commrec* cr,
1828 pull_group_work_t* pg,
1829 gmx_bool bConstraint,
1830 const ivec pulldim_con,
1831 const gmx_mtop_t* mtop,
1832 const t_inputrec* ir,
1835 /* With EM and BD there are no masses in the integrator.
1836 * But we still want to have the correct mass-weighted COMs.
1837 * So we store the real masses in the weights.
1839 const bool setWeights = (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1840 || ir->eI == IntegrationAlgorithm::BD);
1842 /* In parallel, store we need to extract localWeights from weights at DD time */
1843 std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1845 const SimulationGroups& groups = mtop->groups;
1847 /* Count frozen dimensions and (weighted) mass */
1853 for (int i = 0; i < int(pg->params.ind.size()); i++)
1855 int ii = pg->params.ind[i];
1856 if (bConstraint && ir->opts.nFreeze)
1858 for (int d = 0; d < DIM; d++)
1860 if (pulldim_con[d] == 1
1861 && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1867 const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1869 if (ir->efep == FreeEnergyPerturbationType::No)
1875 m = (1 - lambda) * atom.m + lambda * atom.mB;
1878 if (!pg->params.weight.empty())
1880 w = pg->params.weight[i];
1886 if (EI_ENERGY_MINIMIZATION(ir->eI))
1888 /* Move the mass to the weight */
1892 else if (ir->eI == IntegrationAlgorithm::BD)
1895 if (ir->bd_fric != 0.0F)
1897 mbd = ir->bd_fric * ir->delta_t;
1901 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1903 mbd = ir->delta_t / ir->opts.tau_t[0];
1908 / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1916 weights.push_back(w);
1920 wwmass += m * w * w;
1925 /* We can have single atom groups with zero mass with potential pulling
1926 * without cosine weighting.
1928 if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1930 /* With one atom the mass doesn't matter */
1936 "The total%s mass of pull group %d is zero",
1937 !pg->params.weight.empty() ? " weighted" : "",
1943 fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
1944 if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1945 || ir->eI == IntegrationAlgorithm::BD)
1947 fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1949 if (pg->epgrppbc == epgrppbcCOS)
1951 fprintf(fplog, ", cosine weighting will be used");
1953 fprintf(fplog, "\n");
1958 /* A value != 0 signals not frozen, it is updated later */
1964 for (int d = 0; d < DIM; d++)
1966 ndim += pulldim_con[d] * pg->params.ind.size();
1968 if (fplog && nfrozen > 0 && nfrozen < ndim)
1971 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1972 " that are subject to pulling are frozen.\n"
1973 " For constraint pulling the whole group will be frozen.\n\n",
1981 struct pull_t* init_pull(FILE* fplog,
1982 const pull_params_t* pull_params,
1983 const t_inputrec* ir,
1984 const gmx_mtop_t* mtop,
1985 const t_commrec* cr,
1986 gmx::LocalAtomSetManager* atomSets,
1989 struct pull_t* pull;
1994 /* Copy the pull parameters */
1995 pull->params = *pull_params;
1997 for (int i = 0; i < pull_params->ngroup; ++i)
1999 pull->group.emplace_back(pull_params->group[i],
2000 atomSets->add(pull_params->group[i].ind),
2001 pull_params->bSetPbcRefToPrevStepCOM);
2004 if (cr != nullptr && DOMAINDECOMP(cr))
2006 /* Set up the global to local atom mapping for PBC atoms */
2007 for (pull_group_work_t& group : pull->group)
2009 if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
2011 /* pbcAtomSet consists of a single atom */
2012 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
2013 atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
2018 pull->bPotential = FALSE;
2019 pull->bConstraint = FALSE;
2020 pull->bCylinder = FALSE;
2021 pull->bAngle = FALSE;
2022 pull->bXOutAverage = pull_params->bXOutAverage;
2023 pull->bFOutAverage = pull_params->bFOutAverage;
2025 GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
2026 "pull group 0 is an absolute reference group and should not contain atoms");
2028 pull->numCoordinatesWithExternalPotential = 0;
2030 for (int c = 0; c < pull->params.ncoord; c++)
2032 /* Construct a pull coordinate, copying all coordinate parameters */
2033 pull->coord.emplace_back(pull_params->coord[c]);
2035 pull_coord_work_t* pcrd = &pull->coord.back();
2037 switch (pcrd->params.eGeom)
2039 case PullGroupGeometry::Distance:
2040 case PullGroupGeometry::DirectionRelative: /* Direction vector is determined at each step */
2041 case PullGroupGeometry::Angle:
2042 case PullGroupGeometry::Dihedral: break;
2043 case PullGroupGeometry::Direction:
2044 case PullGroupGeometry::DirectionPBC:
2045 case PullGroupGeometry::Cylinder:
2046 case PullGroupGeometry::AngleAxis:
2047 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2050 /* We allow reading of newer tpx files with new pull geometries,
2051 * but with the same tpx format, with old code. A new geometry
2052 * only adds a new enum value, which will be out of range for
2053 * old code. The first place we need to generate an error is
2054 * here, since the pull code can't handle this.
2055 * The enum can be used before arriving here only for printing
2056 * the string corresponding to the geometry, which will result
2057 * in printing "UNDEFINED".
2060 "Pull geometry not supported for pull coordinate %d. The geometry enum "
2061 "%s in the input is larger than that supported by the code (up to %d). "
2062 "You are probably reading a tpr file generated with a newer version of "
2063 "GROMACS with an binary from an older version of Gromacs.",
2065 enumValueToString(pcrd->params.eGeom),
2066 static_cast<int>(PullGroupGeometry::Count) - 1);
2069 if (pcrd->params.eType == PullingAlgorithm::Constraint)
2071 /* Check restrictions of the constraint pull code */
2072 if (pcrd->params.eGeom == PullGroupGeometry::Cylinder
2073 || pcrd->params.eGeom == PullGroupGeometry::DirectionRelative
2074 || pcrd->params.eGeom == PullGroupGeometry::Angle
2075 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2076 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2079 "Pulling of type %s can not be combined with geometry %s. Consider using "
2081 enumValueToString(pcrd->params.eType),
2082 enumValueToString(pcrd->params.eGeom),
2083 enumValueToString(PullingAlgorithm::Umbrella));
2088 "Constraint pulling can not be combined with multiple time stepping");
2090 pull->bConstraint = TRUE;
2094 pull->bPotential = TRUE;
2097 if (pcrd->params.eGeom == PullGroupGeometry::Cylinder)
2099 pull->bCylinder = TRUE;
2101 else if (pcrd->params.eGeom == PullGroupGeometry::Angle
2102 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2103 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2105 pull->bAngle = TRUE;
2108 /* We only need to calculate the plain COM of a group
2109 * when it is not only used as a cylinder group.
2110 * Also the absolute reference group 0 needs no COM computation.
2112 for (int i = 0; i < pcrd->params.ngroup; i++)
2114 int groupIndex = pcrd->params.group[i];
2115 if (groupIndex > 0 && !(pcrd->params.eGeom == PullGroupGeometry::Cylinder && i == 0))
2117 pull->group[groupIndex].needToCalcCom = true;
2121 /* With non-zero rate the reference value is set at every step */
2122 if (pcrd->params.rate == 0)
2124 /* Initialize the constant reference value */
2125 if (pcrd->params.eType != PullingAlgorithm::External)
2127 low_set_pull_coord_reference_value(
2128 pcrd, c, pcrd->params.init * pull_conversion_factor_userinput2internal(&pcrd->params));
2132 /* With an external pull potential, the reference value loses
2133 * it's meaning and should not be used. Setting it to zero
2134 * makes any terms dependent on it disappear.
2135 * The only issue this causes is that with cylinder pulling
2136 * no atoms of the cylinder group within the cylinder radius
2137 * should be more than half a box length away from the COM of
2138 * the pull group along the axial direction.
2140 pcrd->value_ref = 0.0;
2144 if (pcrd->params.eType == PullingAlgorithm::External)
2147 pcrd->params.rate == 0,
2148 "With an external potential, a pull coordinate should have rate = 0");
2150 /* This external potential needs to be registered later */
2151 pull->numCoordinatesWithExternalPotential++;
2153 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2156 pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
2157 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2159 pull->pbcType = ir->pbcType;
2160 switch (pull->pbcType)
2162 case PbcType::No: pull->npbcdim = 0; break;
2163 case PbcType::XY: pull->npbcdim = 2; break;
2164 default: pull->npbcdim = 3; break;
2169 gmx_bool bAbs, bCos;
2172 for (const pull_coord_work_t& coord : pull->coord)
2174 if (pull->group[coord.params.group[0]].params.ind.empty()
2175 || pull->group[coord.params.group[1]].params.ind.empty())
2181 fprintf(fplog, "\n");
2182 if (pull->bPotential)
2184 fprintf(fplog, "Will apply potential COM pulling\n");
2186 if (pull->bConstraint)
2188 fprintf(fplog, "Will apply constraint COM pulling\n");
2190 // Don't include the reference group 0 in output, so we report ngroup-1
2191 int numRealGroups = pull->group.size() - 1;
2192 GMX_RELEASE_ASSERT(numRealGroups > 0,
2193 "The reference absolute position pull group should always be present");
2195 "with %zu pull coordinate%s and %d group%s\n",
2197 pull->coord.size() == 1 ? "" : "s",
2199 numRealGroups == 1 ? "" : "s");
2202 fprintf(fplog, "with an absolute reference\n");
2205 // Don't include the reference group 0 in loop
2206 for (size_t g = 1; g < pull->group.size(); g++)
2208 if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
2210 /* We are using cosine weighting */
2211 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2217 please_cite(fplog, "Engin2010");
2221 pull->bRefAt = FALSE;
2223 for (size_t g = 0; g < pull->group.size(); g++)
2225 pull_group_work_t* pgrp;
2227 pgrp = &pull->group[g];
2228 if (!pgrp->params.ind.empty())
2230 /* There is an issue when a group is used in multiple coordinates
2231 * and constraints are applied in different dimensions with atoms
2232 * frozen in some, but not all dimensions.
2233 * Since there is only one mass array per group, we can't have
2234 * frozen/non-frozen atoms for different coords at the same time.
2235 * But since this is a very exotic case, we don't check for this.
2236 * A warning is printed in init_pull_group_index.
2239 gmx_bool bConstraint;
2240 ivec pulldim, pulldim_con;
2242 /* Loop over all pull coordinates to see along which dimensions
2243 * this group is pulled and if it is involved in constraints.
2245 bConstraint = FALSE;
2246 clear_ivec(pulldim);
2247 clear_ivec(pulldim_con);
2248 for (const pull_coord_work_t& coord : pull->coord)
2250 gmx_bool bGroupUsed = FALSE;
2251 for (int gi = 0; gi < coord.params.ngroup; gi++)
2253 if (coord.params.group[gi] == static_cast<int>(g))
2261 for (int m = 0; m < DIM; m++)
2263 if (coord.params.dim[m] == 1)
2267 if (coord.params.eType == PullingAlgorithm::Constraint)
2277 /* Determine if we need to take PBC into account for calculating
2278 * the COM's of the pull groups.
2280 switch (pgrp->epgrppbc)
2282 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2284 if (!pgrp->params.weight.empty())
2287 "Pull groups can not have relative weights and cosine weighting "
2290 for (int m = 0; m < DIM; m++)
2292 if (m < pull->npbcdim && pulldim[m] == 1)
2294 if (pull->cosdim >= 0 && pull->cosdim != m)
2297 "Can only use cosine weighting with pulling in one "
2298 "dimension (use mdp option pull-coord?-dim)");
2304 case epgrppbcNONE: break;
2307 /* Set the indices */
2308 init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2312 /* Absolute reference, set the inverse mass to zero.
2313 * This is only relevant (and used) with constraint pulling.
2320 /* If we use cylinder coordinates, do some initialising for them */
2321 if (pull->bCylinder)
2323 for (const pull_coord_work_t& coord : pull->coord)
2325 if (coord.params.eGeom == PullGroupGeometry::Cylinder)
2327 if (pull->group[coord.params.group[0]].params.ind.empty())
2330 "A cylinder pull group is not supported when using absolute "
2334 const auto& referenceGroup = pull->group[coord.params.group[0]];
2335 pull->dyna.emplace_back(
2336 referenceGroup.params, referenceGroup.atomSet, pull->params.bSetPbcRefToPrevStepCOM);
2340 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2341 pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
2342 pull->comSums.resize(pull->nthreads);
2347 /* Use a sub-communicator when we have more than 32 ranks, but not
2348 * when we have an external pull potential, since then the external
2349 * potential provider expects each rank to have the coordinate.
2351 comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
2352 || pull->numCoordinatesWithExternalPotential > 0
2353 || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2354 /* This sub-commicator is not used with comm->bParticipateAll,
2355 * so we can always initialize it to NULL.
2357 comm->mpi_comm_com = MPI_COMM_NULL;
2358 comm->nparticipate = 0;
2359 comm->isMasterRank = (cr == nullptr || MASTER(cr));
2361 /* No MPI: 1 rank: all ranks pull */
2362 comm->bParticipateAll = TRUE;
2363 comm->isMasterRank = true;
2365 comm->bParticipate = comm->bParticipateAll;
2366 comm->setup_count = 0;
2367 comm->must_count = 0;
2369 if (!comm->bParticipateAll && fplog != nullptr)
2371 fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2374 comm->pbcAtomBuffer.resize(pull->group.size());
2375 comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
2376 if (pull->bCylinder)
2378 comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
2381 /* We still need to initialize the PBC reference coordinates */
2382 pull->bSetPBCatoms = TRUE;
2384 pull->out_x = nullptr;
2385 pull->out_f = nullptr;
2390 static void destroy_pull(struct pull_t* pull)
2393 if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2395 MPI_Comm_free(&pull->comm.mpi_comm_com);
2402 void preparePrevStepPullCom(const t_inputrec* ir,
2406 const t_state* state_global,
2407 const t_commrec* cr,
2408 bool startingFromCheckpoint)
2410 if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
2414 allocStatePrevStepPullCom(state, pull_work);
2415 if (startingFromCheckpoint)
2419 state->pull_com_prev_step = state_global->pull_com_prev_step;
2423 /* Only the master rank has the checkpointed COM from the previous step */
2424 gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(),
2425 &state->pull_com_prev_step[0],
2426 cr->mpi_comm_mygroup);
2428 setPrevStepPullComFromState(pull_work, state);
2433 set_pbc(&pbc, ir->pbcType, state->box);
2434 initPullComFromPrevStep(cr, pull_work, masses, &pbc, state->x.rvec_array());
2435 updatePrevStepPullCom(pull_work, state);
2439 void finish_pull(struct pull_t* pull)
2441 check_external_potential_registration(pull);
2445 gmx_fio_fclose(pull->out_x);
2449 gmx_fio_fclose(pull->out_f);
2455 bool pull_have_potential(const pull_t& pull)
2457 return pull.bPotential;
2460 bool pull_have_constraint(const pull_t& pull)
2462 return pull.bConstraint;
2465 bool pull_have_constraint(const pull_params_t& pullParameters)
2467 for (int c = 0; c < pullParameters.ncoord; c++)
2469 if (pullParameters.coord[c].eType == PullingAlgorithm::Constraint)