2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
53 #include "gromacs/commandline/filenm.h"
54 #include "gromacs/domdec/domdec_struct.h"
55 #include "gromacs/domdec/localatomset.h"
56 #include "gromacs/domdec/localatomsetmanager.h"
57 #include "gromacs/fileio/gmxfio.h"
58 #include "gromacs/gmxlib/network.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/math/vectypes.h"
63 #include "gromacs/mdlib/gmx_omp_nthreads.h"
64 #include "gromacs/mdtypes/commrec.h"
65 #include "gromacs/mdtypes/forceoutput.h"
66 #include "gromacs/mdtypes/inputrec.h"
67 #include "gromacs/mdtypes/md_enums.h"
68 #include "gromacs/mdtypes/state.h"
69 #include "gromacs/pbcutil/pbc.h"
70 #include "gromacs/topology/mtop_lookup.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/utility/cstringutil.h"
73 #include "gromacs/utility/exceptions.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/futil.h"
76 #include "gromacs/utility/gmxassert.h"
77 #include "gromacs/utility/pleasecite.h"
78 #include "gromacs/utility/real.h"
79 #include "gromacs/utility/smalloc.h"
80 #include "gromacs/utility/strconvert.h"
82 #include "pull_internal.h"
86 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
89 static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
91 if (params.ind.size() <= 1)
96 else if (params.pbcatom >= 0)
98 if (setPbcRefToPrevStepCOM)
100 return epgrppbcPREVSTEPCOM;
104 return epgrppbcREFAT;
113 /* NOTE: The params initialization currently copies pointers.
114 * So the lifetime of the source, currently always inputrec,
115 * should not end before that of this object.
116 * This will be fixed when the pointers are replacted by std::vector.
118 pull_group_work_t::pull_group_work_t(const t_pull_group& params,
119 gmx::LocalAtomSet atomSet,
120 bool bSetPbcRefToPrevStepCOM) :
122 epgrppbc(groupPbcFromParams(params, bSetPbcRefToPrevStepCOM)),
123 needToCalcCom(false),
133 bool pull_coordinate_is_angletype(const t_pull_coord* pcrd)
135 return (pcrd->eGeom == epullgANGLE || pcrd->eGeom == epullgDIHEDRAL || pcrd->eGeom == epullgANGLEAXIS);
138 static bool pull_coordinate_is_directional(const t_pull_coord* pcrd)
140 return (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC
141 || pcrd->eGeom == epullgDIRRELATIVE || pcrd->eGeom == epullgCYL);
144 const char* pull_coordinate_units(const t_pull_coord* pcrd)
146 return pull_coordinate_is_angletype(pcrd) ? "deg" : "nm";
149 double pull_conversion_factor_userinput2internal(const t_pull_coord* pcrd)
151 if (pull_coordinate_is_angletype(pcrd))
161 double pull_conversion_factor_internal2userinput(const t_pull_coord* pcrd)
163 if (pull_coordinate_is_angletype(pcrd))
173 /* Apply forces in a mass weighted fashion for part of the pull group */
174 static void apply_forces_grp_part(const pull_group_work_t* pgrp,
182 double inv_wm = pgrp->mwscale;
184 auto localAtomIndices = pgrp->atomSet.localIndex();
185 for (int i = ind_start; i < ind_end; i++)
187 int ii = localAtomIndices[i];
188 double wmass = masses[ii];
189 if (!pgrp->localWeights.empty())
191 wmass *= pgrp->localWeights[i];
194 for (int d = 0; d < DIM; d++)
196 f[ii][d] += sign * wmass * f_pull[d] * inv_wm;
201 /* Apply forces in a mass weighted fashion */
202 static void apply_forces_grp(const pull_group_work_t* pgrp,
209 auto localAtomIndices = pgrp->atomSet.localIndex();
211 if (pgrp->params.ind.size() == 1 && pgrp->atomSet.numAtomsLocal() == 1)
213 /* Only one atom and our rank has this atom: we can skip
214 * the mass weighting, which means that this code also works
215 * for mass=0, e.g. with a virtual site.
217 for (int d = 0; d < DIM; d++)
219 f[localAtomIndices[0]][d] += sign * f_pull[d];
224 if (localAtomIndices.size() <= c_pullMaxNumLocalAtomsSingleThreaded)
226 apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), masses, f_pull, sign, f);
230 #pragma omp parallel for num_threads(nthreads) schedule(static)
231 for (int th = 0; th < nthreads; th++)
233 int ind_start = (localAtomIndices.size() * (th + 0)) / nthreads;
234 int ind_end = (localAtomIndices.size() * (th + 1)) / nthreads;
235 apply_forces_grp_part(pgrp, ind_start, ind_end, masses, f_pull, sign, f);
241 /* Apply forces in a mass weighted fashion to a cylinder group */
242 static void apply_forces_cyl_grp(const pull_group_work_t* pgrp,
243 const double dv_corr,
249 int gmx_unused nthreads)
251 double inv_wm = pgrp->mwscale;
253 auto localAtomIndices = pgrp->atomSet.localIndex();
255 /* The cylinder group is always a slab in the system, thus large.
256 * Therefore we always thread-parallelize this group.
258 int numAtomsLocal = localAtomIndices.size();
259 #pragma omp parallel for num_threads(nthreads) schedule(static)
260 for (int i = 0; i < numAtomsLocal; i++)
262 double weight = pgrp->localWeights[i];
267 int ii = localAtomIndices[i];
268 double mass = masses[ii];
269 /* The stored axial distance from the cylinder center (dv) needs
270 * to be corrected for an offset (dv_corr), which was unknown when
273 double dv_com = pgrp->dv[i] + dv_corr;
275 /* Here we not only add the pull force working along vec (f_pull),
276 * but also a radial component, due to the dependence of the weights
277 * on the radial distance.
279 for (int m = 0; m < DIM; m++)
281 f[ii][m] += sign * inv_wm * (mass * weight * f_pull[m] + pgrp->mdw[i][m] * dv_com * f_scal);
286 /* Apply torque forces in a mass weighted fashion to the groups that define
287 * the pull vector direction for pull coordinate pcrd.
289 static void apply_forces_vec_torque(const struct pull_t* pull,
290 const pull_coord_work_t* pcrd,
294 const PullCoordSpatialData& spatialData = pcrd->spatialData;
296 /* The component inpr along the pull vector is accounted for in the usual
297 * way. Here we account for the component perpendicular to vec.
300 for (int m = 0; m < DIM; m++)
302 inpr += spatialData.dr01[m] * spatialData.vec[m];
304 /* The torque force works along the component of the distance vector
305 * of between the two "usual" pull groups that is perpendicular to
306 * the pull vector. The magnitude of this force is the "usual" scale force
307 * multiplied with the ratio of the distance between the two "usual" pull
308 * groups and the distance between the two groups that define the vector.
311 for (int m = 0; m < DIM; m++)
313 f_perp[m] = (spatialData.dr01[m] - inpr * spatialData.vec[m]) / spatialData.vec_len * pcrd->scalarForce;
316 /* Apply the force to the groups defining the vector using opposite signs */
317 apply_forces_grp(&pull->group[pcrd->params.group[2]], masses, f_perp, -1, f, pull->nthreads);
318 apply_forces_grp(&pull->group[pcrd->params.group[3]], masses, f_perp, 1, f, pull->nthreads);
321 /* Apply forces in a mass weighted fashion */
322 static void apply_forces_coord(struct pull_t* pull,
324 const PullCoordVectorForces& forces,
328 /* Here it would be more efficient to use one large thread-parallel
329 * region instead of potential parallel regions within apply_forces_grp.
330 * But there could be overlap between pull groups and this would lead
334 const pull_coord_work_t& pcrd = pull->coord[coord];
336 if (pcrd.params.eGeom == epullgCYL)
338 apply_forces_cyl_grp(&pull->dyna[coord],
339 pcrd.spatialData.cyl_dev,
347 /* Sum the force along the vector and the radial force */
349 for (int m = 0; m < DIM; m++)
351 f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
353 apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, f_tot, 1, f, pull->nthreads);
357 if (pcrd.params.eGeom == epullgDIRRELATIVE)
359 /* We need to apply the torque forces to the pull groups
360 * that define the pull vector.
362 apply_forces_vec_torque(pull, &pcrd, masses, f);
365 if (!pull->group[pcrd.params.group[0]].params.ind.empty())
368 &pull->group[pcrd.params.group[0]], masses, forces.force01, -1, f, pull->nthreads);
370 apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, forces.force01, 1, f, pull->nthreads);
372 if (pcrd.params.ngroup >= 4)
375 &pull->group[pcrd.params.group[2]], masses, forces.force23, -1, f, pull->nthreads);
376 apply_forces_grp(&pull->group[pcrd.params.group[3]], masses, forces.force23, 1, f, pull->nthreads);
378 if (pcrd.params.ngroup >= 6)
381 &pull->group[pcrd.params.group[4]], masses, forces.force45, -1, f, pull->nthreads);
382 apply_forces_grp(&pull->group[pcrd.params.group[5]], masses, forces.force45, 1, f, pull->nthreads);
387 real max_pull_distance2(const pull_coord_work_t* pcrd, const t_pbc* pbc)
389 /* Note that this maximum distance calculation is more complex than
390 * most other cases in GROMACS, since here we have to take care of
391 * distance calculations that don't involve all three dimensions.
392 * For example, we can use distances that are larger than the
393 * box X and Y dimensions for a box that is elongated along Z.
396 real max_d2 = GMX_REAL_MAX;
398 if (pull_coordinate_is_directional(&pcrd->params))
400 /* Directional pulling along along direction pcrd->vec.
401 * Calculating the exact maximum distance is complex and bug-prone.
402 * So we take a safe approach by not allowing distances that
403 * are larger than half the distance between unit cell faces
404 * along dimensions involved in pcrd->vec.
406 for (int m = 0; m < DIM; m++)
408 if (m < pbc->ndim_ePBC && pcrd->spatialData.vec[m] != 0)
410 real imageDistance2 = gmx::square(pbc->box[m][m]);
411 for (int d = m + 1; d < DIM; d++)
413 imageDistance2 -= gmx::square(pbc->box[d][m]);
415 max_d2 = std::min(max_d2, imageDistance2);
421 /* Distance pulling along dimensions with pcrd->params.dim[d]==1.
422 * We use half the minimum box vector length of the dimensions involved.
423 * This is correct for all cases, except for corner cases with
424 * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
425 * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
426 * in such corner cases the user could get correct results,
427 * depending on the details of the setup, we avoid further
428 * code complications.
430 for (int m = 0; m < DIM; m++)
432 if (m < pbc->ndim_ePBC && pcrd->params.dim[m] != 0)
434 real imageDistance2 = gmx::square(pbc->box[m][m]);
435 for (int d = 0; d < m; d++)
437 if (pcrd->params.dim[d] != 0)
439 imageDistance2 += gmx::square(pbc->box[m][d]);
442 max_d2 = std::min(max_d2, imageDistance2);
447 return 0.25 * max_d2;
450 /* This function returns the distance based on coordinates xg and xref.
451 * Note that the pull coordinate struct pcrd is not modified.
453 * \param[in] pull The pull struct
454 * \param[in] pcrd The pull coordinate to compute a distance for
455 * \param[in] pbc The periodic boundary conditions
456 * \param[in] xg The coordinate of group 1
457 * \param[in] xref The coordinate of group 0
458 * \param[in] groupIndex0 The index of group 0 in the pcrd->params.group
459 * \param[in] groupIndex1 The index of group 1 in the pcrd->params.group
460 * \param[in] max_dist2 The maximum distance squared
461 * \param[out] dr The distance vector
463 static void low_get_pull_coord_dr(const struct pull_t* pull,
464 const pull_coord_work_t* pcrd,
468 const int groupIndex0,
469 const int groupIndex1,
470 const double max_dist2,
473 const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
475 /* Only the first group can be an absolute reference, in that case nat=0 */
476 if (pgrp0->params.ind.empty())
478 for (int m = 0; m < DIM; m++)
480 xref[m] = pcrd->params.origin[m];
485 copy_dvec(xref, xrefr);
487 dvec dref = { 0, 0, 0 };
488 if (pcrd->params.eGeom == epullgDIRPBC)
490 for (int m = 0; m < DIM; m++)
492 dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
494 /* Add the reference position, so we use the correct periodic image */
495 dvec_inc(xrefr, dref);
498 pbc_dx_d(pbc, xg, xrefr, dr);
500 bool directional = pull_coordinate_is_directional(&pcrd->params);
502 for (int m = 0; m < DIM; m++)
504 dr[m] *= pcrd->params.dim[m];
505 if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
507 dr2 += dr[m] * dr[m];
510 /* Check if we are close to switching to another periodic image */
511 if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
513 /* Note that technically there is no issue with switching periodic
514 * image, as pbc_dx_d returns the distance to the closest periodic
515 * image. However in all cases where periodic image switches occur,
516 * the pull results will be useless in practice.
519 "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
520 "box size (%f).\n%s",
521 pcrd->params.group[groupIndex0],
522 pcrd->params.group[groupIndex1],
524 sqrt(0.98 * 0.98 * max_dist2),
525 pcrd->params.eGeom == epullgDIR
526 ? "You might want to consider using \"pull-geometry = "
527 "direction-periodic\" instead.\n"
531 if (pcrd->params.eGeom == epullgDIRPBC)
537 /* This function returns the distance based on the contents of the pull struct.
538 * pull->coord[coord_ind].dr, and potentially vector, are updated.
540 static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
542 pull_coord_work_t* pcrd = &pull->coord[coord_ind];
543 PullCoordSpatialData& spatialData = pcrd->spatialData;
546 /* With AWH pulling we allow for periodic pulling with geometry=direction.
547 * TODO: Store a periodicity flag instead of checking for external pull provider.
549 if (pcrd->params.eGeom == epullgDIRPBC
550 || (pcrd->params.eGeom == epullgDIR && pcrd->params.eType == epullEXTERNAL))
556 md2 = static_cast<double>(max_pull_distance2(pcrd, pbc));
559 if (pcrd->params.eGeom == epullgDIRRELATIVE)
561 /* We need to determine the pull vector */
565 const pull_group_work_t& pgrp2 = pull->group[pcrd->params.group[2]];
566 const pull_group_work_t& pgrp3 = pull->group[pcrd->params.group[3]];
568 pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
570 for (m = 0; m < DIM; m++)
572 vec[m] *= pcrd->params.dim[m];
574 spatialData.vec_len = dnorm(vec);
575 for (m = 0; m < DIM; m++)
577 spatialData.vec[m] = vec[m] / spatialData.vec_len;
582 "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
589 spatialData.vec[ZZ]);
593 pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
594 pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
596 low_get_pull_coord_dr(pull,
600 pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x,
606 if (pcrd->params.ngroup >= 4)
608 pull_group_work_t *pgrp2, *pgrp3;
609 pgrp2 = &pull->group[pcrd->params.group[2]];
610 pgrp3 = &pull->group[pcrd->params.group[3]];
612 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, 2, 3, md2, spatialData.dr23);
614 if (pcrd->params.ngroup >= 6)
616 pull_group_work_t *pgrp4, *pgrp5;
617 pgrp4 = &pull->group[pcrd->params.group[4]];
618 pgrp5 = &pull->group[pcrd->params.group[5]];
620 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, 4, 5, md2, spatialData.dr45);
624 /* Modify x so that it is periodic in [-pi, pi)
625 * It is assumed that x is in [-3pi, 3pi) so that x
626 * needs to be shifted by at most one period.
628 static void make_periodic_2pi(double* x)
640 /* This function should always be used to modify pcrd->value_ref */
641 static void low_set_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double value_ref)
643 GMX_ASSERT(pcrd->params.eType != epullEXTERNAL,
644 "The pull coord reference value should not be used with type external-potential");
646 if (pcrd->params.eGeom == epullgDIST)
651 "Pull reference distance for coordinate %d (%f) needs to be non-negative",
656 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgANGLEAXIS)
658 if (value_ref < 0 || value_ref > M_PI)
661 "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
662 "interval [0,180] deg",
664 value_ref * pull_conversion_factor_internal2userinput(&pcrd->params));
667 else if (pcrd->params.eGeom == epullgDIHEDRAL)
669 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
670 make_periodic_2pi(&value_ref);
673 pcrd->value_ref = value_ref;
676 static void update_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double t)
678 /* With zero rate the reference value is set initially and doesn't change */
679 if (pcrd->params.rate != 0)
681 double value_ref = (pcrd->params.init + pcrd->params.rate * t)
682 * pull_conversion_factor_userinput2internal(&pcrd->params);
683 low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
687 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
688 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
691 dvec dr32; /* store instead of dr23? */
693 dsvmul(-1, spatialData->dr23, dr32);
694 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
695 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
696 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
698 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
699 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
700 * Note 2: the angle between the plane normal vectors equals pi only when
701 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
702 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
704 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
708 /* Calculates pull->coord[coord_ind].value.
709 * This function also updates pull->coord[coord_ind].dr.
711 static void get_pull_coord_distance(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
713 pull_coord_work_t* pcrd;
716 pcrd = &pull->coord[coord_ind];
718 get_pull_coord_dr(pull, coord_ind, pbc);
720 PullCoordSpatialData& spatialData = pcrd->spatialData;
722 switch (pcrd->params.eGeom)
725 /* Pull along the vector between the com's */
726 spatialData.value = dnorm(spatialData.dr01);
730 case epullgDIRRELATIVE:
733 spatialData.value = 0;
734 for (m = 0; m < DIM; m++)
736 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
740 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
742 case epullgDIHEDRAL: spatialData.value = get_dihedral_angle_coord(&spatialData); break;
743 case epullgANGLEAXIS:
744 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
746 default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
750 /* Returns the deviation from the reference value.
751 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
753 static double get_pull_coord_deviation(struct pull_t* pull, int coord_ind, const t_pbc* pbc, double t)
755 pull_coord_work_t* pcrd;
758 pcrd = &pull->coord[coord_ind];
760 /* Update the reference value before computing the distance,
761 * since it is used in the distance computation with periodic pulling.
763 update_pull_coord_reference_value(pcrd, coord_ind, t);
765 get_pull_coord_distance(pull, coord_ind, pbc);
767 get_pull_coord_distance(pull, coord_ind, pbc);
769 /* Determine the deviation */
770 dev = pcrd->spatialData.value - pcrd->value_ref;
772 /* Check that values are allowed */
773 if (pcrd->params.eGeom == epullgDIST && pcrd->spatialData.value == 0)
775 /* With no vector we can not determine the direction for the force,
776 * so we set the force to zero.
780 else if (pcrd->params.eGeom == epullgDIHEDRAL)
782 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
783 Thus, the unwrapped deviation is here in (-2pi, 2pi].
784 After making it periodic, the deviation will be in [-pi, pi). */
785 make_periodic_2pi(&dev);
791 double get_pull_coord_value(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
793 get_pull_coord_distance(pull, coord_ind, pbc);
795 return pull->coord[coord_ind].spatialData.value;
798 void clear_pull_forces(pull_t* pull)
800 /* Zeroing the forces is only required for constraint pulling.
801 * It can happen that multiple constraint steps need to be applied
802 * and therefore the constraint forces need to be accumulated.
804 for (pull_coord_work_t& coord : pull->coord)
806 coord.scalarForce = 0;
810 /* Apply constraint using SHAKE */
812 do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
815 dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
816 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
817 dvec* rnew; /* current 'new' positions of the groups */
818 double* dr_tot; /* the total update of the coords */
821 double lambda, rm, invdt = 0;
822 gmx_bool bConverged_all, bConverged = FALSE;
823 int niter = 0, ii, j, m, max_iter = 100;
827 snew(r_ij, pull->coord.size());
828 snew(dr_tot, pull->coord.size());
830 snew(rnew, pull->group.size());
832 /* copy the current unconstrained positions for use in iterations. We
833 iterate until rinew[i] and rjnew[j] obey the constraints. Then
834 rinew - pull.x_unc[i] is the correction dr to group i */
835 for (size_t g = 0; g < pull->group.size(); g++)
837 copy_dvec(pull->group[g].xp, rnew[g]);
840 /* Determine the constraint directions from the old positions */
841 for (size_t c = 0; c < pull->coord.size(); c++)
843 pull_coord_work_t* pcrd;
845 pcrd = &pull->coord[c];
847 if (pcrd->params.eType != epullCONSTRAINT)
852 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
853 * We don't modify dr and value anymore, so these values are also used
856 get_pull_coord_distance(pull, c, pbc);
858 const PullCoordSpatialData& spatialData = pcrd->spatialData;
862 "Pull coord %zu dr %f %f %f\n",
864 spatialData.dr01[XX],
865 spatialData.dr01[YY],
866 spatialData.dr01[ZZ]);
869 if (pcrd->params.eGeom == epullgDIR || pcrd->params.eGeom == epullgDIRPBC)
871 /* Select the component along vec */
873 for (m = 0; m < DIM; m++)
875 a += spatialData.vec[m] * spatialData.dr01[m];
877 for (m = 0; m < DIM; m++)
879 r_ij[c][m] = a * spatialData.vec[m];
884 copy_dvec(spatialData.dr01, r_ij[c]);
887 if (dnorm2(r_ij[c]) == 0)
890 "Distance for pull coordinate %zu is zero with constraint pulling, which is "
896 bConverged_all = FALSE;
897 while (!bConverged_all && niter < max_iter)
899 bConverged_all = TRUE;
901 /* loop over all constraints */
902 for (size_t c = 0; c < pull->coord.size(); c++)
904 pull_coord_work_t* pcrd;
905 pull_group_work_t *pgrp0, *pgrp1;
908 pcrd = &pull->coord[c];
910 if (pcrd->params.eType != epullCONSTRAINT)
915 update_pull_coord_reference_value(pcrd, c, t);
917 pgrp0 = &pull->group[pcrd->params.group[0]];
918 pgrp1 = &pull->group[pcrd->params.group[1]];
920 /* Get the current difference vector */
921 low_get_pull_coord_dr(
922 pull, pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
926 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
929 rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
931 switch (pcrd->params.eGeom)
934 if (pcrd->value_ref <= 0)
938 "The pull constraint reference distance for group %zu is <= 0 (%f)",
944 double q, c_a, c_b, c_c;
946 c_a = diprod(r_ij[c], r_ij[c]);
947 c_b = diprod(unc_ij, r_ij[c]) * 2;
948 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
952 q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
957 q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
963 fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b, c_c, lambda);
967 /* The position corrections dr due to the constraints */
968 dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
969 dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
970 dr_tot[c] += -lambda * dnorm(r_ij[c]);
975 /* A 1-dimensional constraint along a vector */
977 for (m = 0; m < DIM; m++)
979 vec[m] = pcrd->spatialData.vec[m];
980 a += unc_ij[m] * vec[m];
982 /* Select only the component along the vector */
983 dsvmul(a, vec, unc_ij);
984 lambda = a - pcrd->value_ref;
987 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
990 /* The position corrections dr due to the constraints */
991 dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
992 dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
993 dr_tot[c] += -lambda;
995 default: gmx_incons("Invalid enumeration value for eGeom");
1003 g0 = pcrd->params.group[0];
1004 g1 = pcrd->params.group[1];
1005 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
1006 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
1008 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1017 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
1026 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1036 /* Update the COMs with dr */
1037 dvec_inc(rnew[pcrd->params.group[1]], dr1);
1038 dvec_inc(rnew[pcrd->params.group[0]], dr0);
1041 /* Check if all constraints are fullfilled now */
1042 for (pull_coord_work_t& coord : pull->coord)
1044 if (coord.params.eType != epullCONSTRAINT)
1049 low_get_pull_coord_dr(
1050 pull, &coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
1052 switch (coord.params.eGeom)
1055 bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1060 for (m = 0; m < DIM; m++)
1062 vec[m] = coord.spatialData.vec[m];
1064 inpr = diprod(unc_ij, vec);
1065 dsvmul(inpr, vec, unc_ij);
1066 bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1075 "Pull constraint not converged: "
1077 "d_ref = %f, current d = %f\n",
1078 coord.params.group[0],
1079 coord.params.group[1],
1084 bConverged_all = FALSE;
1089 /* if after all constraints are dealt with and bConverged is still TRUE
1090 we're finished, if not we do another iteration */
1092 if (niter > max_iter)
1094 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1097 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1104 /* update atoms in the groups */
1105 for (size_t g = 0; g < pull->group.size(); g++)
1107 const pull_group_work_t* pgrp;
1110 pgrp = &pull->group[g];
1112 /* get the final constraint displacement dr for group g */
1113 dvec_sub(rnew[g], pgrp->xp, dr);
1115 if (dnorm2(dr) == 0)
1117 /* No displacement, no update necessary */
1121 /* update the atom positions */
1122 auto localAtomIndices = pgrp->atomSet.localIndex();
1124 for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1126 ii = localAtomIndices[j];
1127 if (!pgrp->localWeights.empty())
1129 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1131 for (m = 0; m < DIM; m++)
1137 for (m = 0; m < DIM; m++)
1139 v[ii][m] += invdt * tmp[m];
1145 /* calculate the constraint forces, used for output and virial only */
1146 for (size_t c = 0; c < pull->coord.size(); c++)
1148 pull_coord_work_t* pcrd;
1150 pcrd = &pull->coord[c];
1152 if (pcrd->params.eType != epullCONSTRAINT)
1157 /* Accumulate the forces, in case we have multiple constraint steps */
1160 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1162 pcrd->scalarForce += force;
1164 if (vir != nullptr && pcrd->params.eGeom != epullgDIRPBC && bMaster)
1168 /* Add the pull contribution to the virial */
1169 /* We have already checked above that r_ij[c] != 0 */
1170 f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1172 for (j = 0; j < DIM; j++)
1174 for (m = 0; m < DIM; m++)
1176 vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1182 /* finished! I hope. Give back some memory */
1188 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1190 for (int j = 0; j < DIM; j++)
1192 for (int m = 0; m < DIM; m++)
1194 vir[j][m] -= 0.5 * f[j] * dr[m];
1199 /* Adds the pull contribution to the virial */
1200 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1202 if (vir != nullptr && pcrd.params.eGeom != epullgDIRPBC)
1204 /* Add the pull contribution for each distance vector to the virial. */
1205 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1206 if (pcrd.params.ngroup >= 4)
1208 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1210 if (pcrd.params.ngroup >= 6)
1212 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1217 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1225 k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1226 dkdl = pcrd->params.kB - pcrd->params.k;
1228 switch (pcrd->params.eType)
1231 case epullFLATBOTTOM:
1232 case epullFLATBOTTOMHIGH:
1233 /* The only difference between an umbrella and a flat-bottom
1234 * potential is that a flat-bottom is zero above or below
1235 the reference value.
1237 if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0)
1238 || (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
1243 pcrd->scalarForce = -k * dev;
1244 *V += 0.5 * k * gmx::square(dev);
1245 *dVdl += 0.5 * dkdl * gmx::square(dev);
1248 pcrd->scalarForce = -k;
1249 *V += k * pcrd->spatialData.value;
1250 *dVdl += dkdl * pcrd->spatialData.value;
1254 "the scalar pull force should not be calculated internally for pull type "
1256 default: gmx_incons("Unsupported pull type in do_pull_pot");
1260 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1262 const t_pull_coord& params = pcrd.params;
1263 const PullCoordSpatialData& spatialData = pcrd.spatialData;
1265 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1266 PullCoordVectorForces forces;
1268 if (params.eGeom == epullgDIST)
1270 double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1271 for (int m = 0; m < DIM; m++)
1273 forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1276 else if (params.eGeom == epullgANGLE)
1279 double cos_theta, cos_theta2;
1281 cos_theta = cos(spatialData.value);
1282 cos_theta2 = gmx::square(cos_theta);
1284 /* The force at theta = 0, pi is undefined so we should not apply any force.
1285 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1286 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1287 * have to check for this before dividing by their norm below.
1291 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1292 * and another vector parallel to dr23:
1293 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1294 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1296 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1297 double b = a * cos_theta;
1298 double invdr01 = 1. / dnorm(spatialData.dr01);
1299 double invdr23 = 1. / dnorm(spatialData.dr23);
1300 dvec normalized_dr01, normalized_dr23;
1301 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1302 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1304 for (int m = 0; m < DIM; m++)
1306 /* Here, f_scal is -dV/dtheta */
1308 pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1310 pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1315 /* No forces to apply for ill-defined cases*/
1316 clear_dvec(forces.force01);
1317 clear_dvec(forces.force23);
1320 else if (params.eGeom == epullgANGLEAXIS)
1322 double cos_theta, cos_theta2;
1324 /* The angle-axis force is exactly the same as the angle force (above) except that in
1325 this case the second vector (dr23) is replaced by the pull vector. */
1326 cos_theta = cos(spatialData.value);
1327 cos_theta2 = gmx::square(cos_theta);
1333 dvec normalized_dr01;
1335 invdr01 = 1. / dnorm(spatialData.dr01);
1336 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1337 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1340 for (int m = 0; m < DIM; m++)
1343 pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1348 clear_dvec(forces.force01);
1351 else if (params.eGeom == epullgDIHEDRAL)
1353 double m2, n2, tol, sqrdist_32;
1355 /* Note: there is a small difference here compared to the
1356 dihedral force calculations in the bondeds (ref: Bekker 1994).
1357 There rij = ri - rj, while here dr01 = r1 - r0.
1358 However, all distance vectors occur in form of cross or inner products
1359 so that two signs cancel and we end up with the same expressions.
1360 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1362 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1363 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1364 dsvmul(-1, spatialData.dr23, dr32);
1365 sqrdist_32 = diprod(dr32, dr32);
1366 tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1367 if ((m2 > tol) && (n2 > tol))
1369 double a_01, a_23_01, a_23_45, a_45;
1370 double inv_dist_32, inv_sqrdist_32, dist_32;
1372 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1373 inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1374 dist_32 = sqrdist_32 * inv_dist_32;
1376 /* Forces on groups 0, 1 */
1377 a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
1378 dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1380 /* Forces on groups 4, 5 */
1381 a_45 = -pcrd.scalarForce * dist_32 / n2;
1382 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1384 /* Force on groups 2, 3 (defining the axis) */
1385 a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1386 a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1387 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1388 dsvmul(a_23_45, forces.force45, v);
1389 dvec_sub(u, v, forces.force23); /* force on group 3 */
1393 /* No force to apply for ill-defined cases */
1394 clear_dvec(forces.force01);
1395 clear_dvec(forces.force23);
1396 clear_dvec(forces.force45);
1401 for (int m = 0; m < DIM; m++)
1403 forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1411 /* We use a global mutex for locking access to the pull data structure
1412 * during registration of external pull potential providers.
1413 * We could use a different, local mutex for each pull object, but the overhead
1414 * is extremely small here and registration is only done during initialization.
1416 static std::mutex registrationMutex;
1418 using Lock = std::lock_guard<std::mutex>;
1420 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1422 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1423 GMX_RELEASE_ASSERT(provider != nullptr,
1424 "register_external_pull_potential called with NULL as provider name");
1426 if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1429 "Module '%s' attempted to register an external potential for pull coordinate %d "
1430 "which is out of the pull coordinate range %d - %zu\n",
1434 pull->coord.size());
1437 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1439 if (pcrd->params.eType != epullEXTERNAL)
1443 "Module '%s' attempted to register an external potential for pull coordinate %d "
1444 "which of type '%s', whereas external potentials are only supported with type '%s'",
1447 epull_names[pcrd->params.eType],
1448 epull_names[epullEXTERNAL]);
1451 GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
1452 "The external potential provider string for a pull coordinate is NULL");
1454 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
1457 "Module '%s' attempted to register an external potential for pull coordinate %d "
1458 "which expects the external potential to be provided by a module named '%s'",
1461 pcrd->params.externalPotentialProvider.c_str());
1464 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1465 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1466 * pull->numUnregisteredExternalPotentials.
1468 Lock registrationLock(registrationMutex);
1470 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1473 "Module '%s' attempted to register an external potential for pull coordinate %d "
1479 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1480 pull->numUnregisteredExternalPotentials--;
1482 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1483 "Negative unregistered potentials, the pull code is inconsistent");
1487 static void check_external_potential_registration(const struct pull_t* pull)
1489 if (pull->numUnregisteredExternalPotentials > 0)
1492 for (c = 0; c < pull->coord.size(); c++)
1494 if (pull->coord[c].params.eType == epullEXTERNAL
1495 && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1501 GMX_RELEASE_ASSERT(c < pull->coord.size(),
1502 "Internal inconsistency in the pull potential provider counting");
1505 "No external provider for external pull potentials have been provided for %d "
1506 "pull coordinates. The first coordinate without provider is number %zu, which "
1507 "expects a module named '%s' to provide the external potential.",
1508 pull->numUnregisteredExternalPotentials,
1510 pull->coord[c].params.externalPotentialProvider.c_str());
1515 /* Pull takes care of adding the forces of the external potential.
1516 * The external potential module has to make sure that the corresponding
1517 * potential energy is added either to the pull term or to a term
1518 * specific to the external module.
1520 void apply_external_pull_coord_force(struct pull_t* pull,
1524 gmx::ForceWithVirial* forceWithVirial)
1526 pull_coord_work_t* pcrd;
1528 GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1529 "apply_external_pull_coord_force called with coord_index out of range");
1531 if (pull->comm.bParticipate)
1533 pcrd = &pull->coord[coord_index];
1536 pcrd->params.eType == epullEXTERNAL,
1537 "The pull force can only be set externally on pull coordinates of external type");
1539 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1540 "apply_external_pull_coord_force called for an unregistered pull coordinate");
1543 pcrd->scalarForce = coord_force;
1545 /* Calculate the forces on the pull groups */
1546 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1548 /* Add the forces for this coordinate to the total virial and force */
1549 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1551 matrix virial = { { 0 } };
1552 add_virial_coord(virial, *pcrd, pullCoordForces);
1553 forceWithVirial->addVirialContribution(virial);
1557 pull, coord_index, pullCoordForces, masses, as_rvec_array(forceWithVirial->force_.data()));
1560 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1563 /* Calculate the pull potential and scalar force for a pull coordinate.
1564 * Returns the vector forces for the pull coordinate.
1566 static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull,
1575 pull_coord_work_t& pcrd = pull->coord[coord_ind];
1577 assert(pcrd.params.eType != epullCONSTRAINT);
1579 double dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
1581 calc_pull_coord_scalar_force_and_potential(&pcrd, dev, lambda, V, dVdl);
1583 PullCoordVectorForces pullCoordForces = calculateVectorForces(pcrd);
1585 add_virial_coord(vir, pcrd, pullCoordForces);
1587 return pullCoordForces;
1590 real pull_potential(struct pull_t* pull,
1593 const t_commrec* cr,
1597 gmx::ForceWithVirial* force,
1602 assert(pull != nullptr);
1604 /* Ideally we should check external potential registration only during
1605 * the initialization phase, but that requires another function call
1606 * that should be done exactly in the right order. So we check here.
1608 check_external_potential_registration(pull);
1610 if (pull->comm.bParticipate)
1614 pull_calc_coms(cr, pull, masses, pbc, t, x, nullptr);
1616 rvec* f = as_rvec_array(force->force_.data());
1617 matrix virial = { { 0 } };
1618 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1619 for (size_t c = 0; c < pull->coord.size(); c++)
1621 pull_coord_work_t* pcrd;
1622 pcrd = &pull->coord[c];
1624 /* For external potential the force is assumed to be given by an external module by a
1625 call to apply_pull_coord_external_force */
1626 if (pcrd->params.eType == epullCONSTRAINT || pcrd->params.eType == epullEXTERNAL)
1631 PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1632 pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1634 /* Distribute the force over the atoms in the pulled groups */
1635 apply_forces_coord(pull, c, pullCoordForces, masses, f);
1640 force->addVirialContribution(virial);
1645 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1646 "Too few or too many external pull potentials have been applied the previous step");
1647 /* All external pull potentials still need to be applied */
1648 pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1650 return (MASTER(cr) ? V : 0.0);
1653 void pull_constraint(struct pull_t* pull,
1656 const t_commrec* cr,
1664 assert(pull != nullptr);
1666 if (pull->comm.bParticipate)
1668 pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
1670 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1674 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1678 gmx_bool bMustParticipate;
1684 /* We always make the master node participate, such that it can do i/o,
1685 * add the virial and to simplify MC type extensions people might have.
1687 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1689 for (pull_group_work_t& group : pull->group)
1691 if (!group.globalWeights.empty())
1693 group.localWeights.resize(group.atomSet.numAtomsLocal());
1694 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1696 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1700 GMX_ASSERT(bMustParticipate || dd != nullptr,
1701 "Either all ranks (including this rank) participate, or we use DD and need to "
1702 "have access to dd here");
1704 /* We should participate if we have pull or pbc atoms */
1705 if (!bMustParticipate
1706 && (group.atomSet.numAtomsLocal() > 0
1707 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1709 bMustParticipate = TRUE;
1713 if (!comm->bParticipateAll)
1715 /* Keep currently not required ranks in the communicator
1716 * if they needed to participate up to 20 decompositions ago.
1717 * This avoids frequent rebuilds due to atoms jumping back and forth.
1719 const int64_t history_count = 20;
1720 gmx_bool bWillParticipate;
1723 /* Increase the decomposition counter for the current call */
1724 comm->setup_count++;
1726 if (bMustParticipate)
1728 comm->must_count = comm->setup_count;
1733 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1735 if (debug && dd != nullptr)
1738 "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
1740 gmx::boolToString(bMustParticipate),
1741 gmx::boolToString(bWillParticipate));
1744 if (bWillParticipate)
1746 /* Count the number of ranks that we want to have participating */
1748 /* Count the number of ranks that need to be added */
1749 count[1] = comm->bParticipate ? 0 : 1;
1757 /* The cost of this global operation will be less that the cost
1758 * of the extra MPI_Comm_split calls that we can avoid.
1760 gmx_sumi(2, count, cr);
1762 /* If we are missing ranks or if we have 20% more ranks than needed
1763 * we make a new sub-communicator.
1765 if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1769 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1773 if (comm->mpi_comm_com != MPI_COMM_NULL)
1775 MPI_Comm_free(&comm->mpi_comm_com);
1778 /* This might be an extremely expensive operation, so we try
1779 * to avoid this splitting as much as possible.
1781 assert(dd != nullptr);
1782 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1785 comm->bParticipate = bWillParticipate;
1786 comm->nparticipate = count[0];
1788 /* When we use the previous COM for PBC, we need to broadcast
1789 * the previous COM to ranks that have joined the communicator.
1791 for (pull_group_work_t& group : pull->group)
1793 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1795 GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1796 "The master rank has to participate, as it should pass an up to "
1798 "to bcast here as well as to e.g. checkpointing");
1800 gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
1806 /* Since the PBC of atoms might have changed, we need to update the PBC */
1807 pull->bSetPBCatoms = TRUE;
1810 static void init_pull_group_index(FILE* fplog,
1811 const t_commrec* cr,
1813 pull_group_work_t* pg,
1814 gmx_bool bConstraint,
1815 const ivec pulldim_con,
1816 const gmx_mtop_t* mtop,
1817 const t_inputrec* ir,
1820 /* With EM and BD there are no masses in the integrator.
1821 * But we still want to have the correct mass-weighted COMs.
1822 * So we store the real masses in the weights.
1824 const bool setWeights =
1825 (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
1827 /* In parallel, store we need to extract localWeights from weights at DD time */
1828 std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1830 const SimulationGroups& groups = mtop->groups;
1832 /* Count frozen dimensions and (weighted) mass */
1838 for (int i = 0; i < int(pg->params.ind.size()); i++)
1840 int ii = pg->params.ind[i];
1841 if (bConstraint && ir->opts.nFreeze)
1843 for (int d = 0; d < DIM; d++)
1845 if (pulldim_con[d] == 1
1846 && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1852 const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1854 if (ir->efep == efepNO)
1860 m = (1 - lambda) * atom.m + lambda * atom.mB;
1863 if (!pg->params.weight.empty())
1865 w = pg->params.weight[i];
1871 if (EI_ENERGY_MINIMIZATION(ir->eI))
1873 /* Move the mass to the weight */
1877 else if (ir->eI == eiBD)
1880 if (ir->bd_fric != 0.0F)
1882 mbd = ir->bd_fric * ir->delta_t;
1886 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1888 mbd = ir->delta_t / ir->opts.tau_t[0];
1893 / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1901 weights.push_back(w);
1905 wwmass += m * w * w;
1910 /* We can have single atom groups with zero mass with potential pulling
1911 * without cosine weighting.
1913 if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1915 /* With one atom the mass doesn't matter */
1921 "The total%s mass of pull group %d is zero",
1922 !pg->params.weight.empty() ? " weighted" : "",
1928 fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
1929 if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1931 fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1933 if (pg->epgrppbc == epgrppbcCOS)
1935 fprintf(fplog, ", cosine weighting will be used");
1937 fprintf(fplog, "\n");
1942 /* A value != 0 signals not frozen, it is updated later */
1948 for (int d = 0; d < DIM; d++)
1950 ndim += pulldim_con[d] * pg->params.ind.size();
1952 if (fplog && nfrozen > 0 && nfrozen < ndim)
1955 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1956 " that are subject to pulling are frozen.\n"
1957 " For constraint pulling the whole group will be frozen.\n\n",
1965 struct pull_t* init_pull(FILE* fplog,
1966 const pull_params_t* pull_params,
1967 const t_inputrec* ir,
1968 const gmx_mtop_t* mtop,
1969 const t_commrec* cr,
1970 gmx::LocalAtomSetManager* atomSets,
1973 struct pull_t* pull;
1978 /* Copy the pull parameters */
1979 pull->params = *pull_params;
1981 for (int i = 0; i < pull_params->ngroup; ++i)
1983 pull->group.emplace_back(pull_params->group[i],
1984 atomSets->add(pull_params->group[i].ind),
1985 pull_params->bSetPbcRefToPrevStepCOM);
1988 if (cr != nullptr && DOMAINDECOMP(cr))
1990 /* Set up the global to local atom mapping for PBC atoms */
1991 for (pull_group_work_t& group : pull->group)
1993 if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
1995 /* pbcAtomSet consists of a single atom */
1996 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
1997 atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
2002 pull->bPotential = FALSE;
2003 pull->bConstraint = FALSE;
2004 pull->bCylinder = FALSE;
2005 pull->bAngle = FALSE;
2006 pull->bXOutAverage = pull_params->bXOutAverage;
2007 pull->bFOutAverage = pull_params->bFOutAverage;
2009 GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
2010 "pull group 0 is an absolute reference group and should not contain atoms");
2012 pull->numCoordinatesWithExternalPotential = 0;
2014 for (int c = 0; c < pull->params.ncoord; c++)
2016 /* Construct a pull coordinate, copying all coordinate parameters */
2017 pull->coord.emplace_back(pull_params->coord[c]);
2019 pull_coord_work_t* pcrd = &pull->coord.back();
2021 switch (pcrd->params.eGeom)
2024 case epullgDIRRELATIVE: /* Direction vector is determined at each step */
2026 case epullgDIHEDRAL: break;
2030 case epullgANGLEAXIS:
2031 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2034 /* We allow reading of newer tpx files with new pull geometries,
2035 * but with the same tpx format, with old code. A new geometry
2036 * only adds a new enum value, which will be out of range for
2037 * old code. The first place we need to generate an error is
2038 * here, since the pull code can't handle this.
2039 * The enum can be used before arriving here only for printing
2040 * the string corresponding to the geometry, which will result
2041 * in printing "UNDEFINED".
2044 "Pull geometry not supported for pull coordinate %d. The geometry enum "
2045 "%d in the input is larger than that supported by the code (up to %d). "
2046 "You are probably reading a tpr file generated with a newer version of "
2047 "Gromacs with an binary from an older version of Gromacs.",
2053 if (pcrd->params.eType == epullCONSTRAINT)
2055 /* Check restrictions of the constraint pull code */
2056 if (pcrd->params.eGeom == epullgCYL || pcrd->params.eGeom == epullgDIRRELATIVE
2057 || pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
2058 || pcrd->params.eGeom == epullgANGLEAXIS)
2061 "Pulling of type %s can not be combined with geometry %s. Consider using "
2063 epull_names[pcrd->params.eType],
2064 epullg_names[pcrd->params.eGeom],
2065 epull_names[epullUMBRELLA]);
2070 "Constraint pulling can not be combined with multiple time stepping");
2072 pull->bConstraint = TRUE;
2076 pull->bPotential = TRUE;
2079 if (pcrd->params.eGeom == epullgCYL)
2081 pull->bCylinder = TRUE;
2083 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
2084 || pcrd->params.eGeom == epullgANGLEAXIS)
2086 pull->bAngle = TRUE;
2089 /* We only need to calculate the plain COM of a group
2090 * when it is not only used as a cylinder group.
2091 * Also the absolute reference group 0 needs no COM computation.
2093 for (int i = 0; i < pcrd->params.ngroup; i++)
2095 int groupIndex = pcrd->params.group[i];
2096 if (groupIndex > 0 && !(pcrd->params.eGeom == epullgCYL && i == 0))
2098 pull->group[groupIndex].needToCalcCom = true;
2102 /* With non-zero rate the reference value is set at every step */
2103 if (pcrd->params.rate == 0)
2105 /* Initialize the constant reference value */
2106 if (pcrd->params.eType != epullEXTERNAL)
2108 low_set_pull_coord_reference_value(
2109 pcrd, c, pcrd->params.init * pull_conversion_factor_userinput2internal(&pcrd->params));
2113 /* With an external pull potential, the reference value loses
2114 * it's meaning and should not be used. Setting it to zero
2115 * makes any terms dependent on it disappear.
2116 * The only issue this causes is that with cylinder pulling
2117 * no atoms of the cylinder group within the cylinder radius
2118 * should be more than half a box length away from the COM of
2119 * the pull group along the axial direction.
2121 pcrd->value_ref = 0.0;
2125 if (pcrd->params.eType == epullEXTERNAL)
2128 pcrd->params.rate == 0,
2129 "With an external potential, a pull coordinate should have rate = 0");
2131 /* This external potential needs to be registered later */
2132 pull->numCoordinatesWithExternalPotential++;
2134 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2137 pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
2138 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2140 pull->pbcType = ir->pbcType;
2141 switch (pull->pbcType)
2143 case PbcType::No: pull->npbcdim = 0; break;
2144 case PbcType::XY: pull->npbcdim = 2; break;
2145 default: pull->npbcdim = 3; break;
2150 gmx_bool bAbs, bCos;
2153 for (const pull_coord_work_t& coord : pull->coord)
2155 if (pull->group[coord.params.group[0]].params.ind.empty()
2156 || pull->group[coord.params.group[1]].params.ind.empty())
2162 fprintf(fplog, "\n");
2163 if (pull->bPotential)
2165 fprintf(fplog, "Will apply potential COM pulling\n");
2167 if (pull->bConstraint)
2169 fprintf(fplog, "Will apply constraint COM pulling\n");
2171 // Don't include the reference group 0 in output, so we report ngroup-1
2172 int numRealGroups = pull->group.size() - 1;
2173 GMX_RELEASE_ASSERT(numRealGroups > 0,
2174 "The reference absolute position pull group should always be present");
2176 "with %zu pull coordinate%s and %d group%s\n",
2178 pull->coord.size() == 1 ? "" : "s",
2180 numRealGroups == 1 ? "" : "s");
2183 fprintf(fplog, "with an absolute reference\n");
2186 // Don't include the reference group 0 in loop
2187 for (size_t g = 1; g < pull->group.size(); g++)
2189 if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
2191 /* We are using cosine weighting */
2192 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2198 please_cite(fplog, "Engin2010");
2202 pull->bRefAt = FALSE;
2204 for (size_t g = 0; g < pull->group.size(); g++)
2206 pull_group_work_t* pgrp;
2208 pgrp = &pull->group[g];
2209 if (!pgrp->params.ind.empty())
2211 /* There is an issue when a group is used in multiple coordinates
2212 * and constraints are applied in different dimensions with atoms
2213 * frozen in some, but not all dimensions.
2214 * Since there is only one mass array per group, we can't have
2215 * frozen/non-frozen atoms for different coords at the same time.
2216 * But since this is a very exotic case, we don't check for this.
2217 * A warning is printed in init_pull_group_index.
2220 gmx_bool bConstraint;
2221 ivec pulldim, pulldim_con;
2223 /* Loop over all pull coordinates to see along which dimensions
2224 * this group is pulled and if it is involved in constraints.
2226 bConstraint = FALSE;
2227 clear_ivec(pulldim);
2228 clear_ivec(pulldim_con);
2229 for (const pull_coord_work_t& coord : pull->coord)
2231 gmx_bool bGroupUsed = FALSE;
2232 for (int gi = 0; gi < coord.params.ngroup; gi++)
2234 if (coord.params.group[gi] == static_cast<int>(g))
2242 for (int m = 0; m < DIM; m++)
2244 if (coord.params.dim[m] == 1)
2248 if (coord.params.eType == epullCONSTRAINT)
2258 /* Determine if we need to take PBC into account for calculating
2259 * the COM's of the pull groups.
2261 switch (pgrp->epgrppbc)
2263 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2265 if (!pgrp->params.weight.empty())
2268 "Pull groups can not have relative weights and cosine weighting "
2271 for (int m = 0; m < DIM; m++)
2273 if (m < pull->npbcdim && pulldim[m] == 1)
2275 if (pull->cosdim >= 0 && pull->cosdim != m)
2278 "Can only use cosine weighting with pulling in one "
2279 "dimension (use mdp option pull-coord?-dim)");
2285 case epgrppbcNONE: break;
2288 /* Set the indices */
2289 init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2293 /* Absolute reference, set the inverse mass to zero.
2294 * This is only relevant (and used) with constraint pulling.
2301 /* If we use cylinder coordinates, do some initialising for them */
2302 if (pull->bCylinder)
2304 for (const pull_coord_work_t& coord : pull->coord)
2306 if (coord.params.eGeom == epullgCYL)
2308 if (pull->group[coord.params.group[0]].params.ind.empty())
2311 "A cylinder pull group is not supported when using absolute "
2315 const auto& referenceGroup = pull->group[coord.params.group[0]];
2316 pull->dyna.emplace_back(
2317 referenceGroup.params, referenceGroup.atomSet, pull->params.bSetPbcRefToPrevStepCOM);
2321 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2322 pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
2323 pull->comSums.resize(pull->nthreads);
2328 /* Use a sub-communicator when we have more than 32 ranks, but not
2329 * when we have an external pull potential, since then the external
2330 * potential provider expects each rank to have the coordinate.
2332 comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
2333 || pull->numCoordinatesWithExternalPotential > 0
2334 || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2335 /* This sub-commicator is not used with comm->bParticipateAll,
2336 * so we can always initialize it to NULL.
2338 comm->mpi_comm_com = MPI_COMM_NULL;
2339 comm->nparticipate = 0;
2340 comm->isMasterRank = (cr == nullptr || MASTER(cr));
2342 /* No MPI: 1 rank: all ranks pull */
2343 comm->bParticipateAll = TRUE;
2344 comm->isMasterRank = true;
2346 comm->bParticipate = comm->bParticipateAll;
2347 comm->setup_count = 0;
2348 comm->must_count = 0;
2350 if (!comm->bParticipateAll && fplog != nullptr)
2352 fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2355 comm->pbcAtomBuffer.resize(pull->group.size());
2356 comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
2357 if (pull->bCylinder)
2359 comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
2362 /* We still need to initialize the PBC reference coordinates */
2363 pull->bSetPBCatoms = TRUE;
2365 pull->out_x = nullptr;
2366 pull->out_f = nullptr;
2371 static void destroy_pull(struct pull_t* pull)
2374 if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2376 MPI_Comm_free(&pull->comm.mpi_comm_com);
2383 void preparePrevStepPullCom(const t_inputrec* ir,
2387 const t_state* state_global,
2388 const t_commrec* cr,
2389 bool startingFromCheckpoint)
2391 if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
2395 allocStatePrevStepPullCom(state, pull_work);
2396 if (startingFromCheckpoint)
2400 state->pull_com_prev_step = state_global->pull_com_prev_step;
2404 /* Only the master rank has the checkpointed COM from the previous step */
2405 gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(),
2406 &state->pull_com_prev_step[0],
2407 cr->mpi_comm_mygroup);
2409 setPrevStepPullComFromState(pull_work, state);
2414 set_pbc(&pbc, ir->pbcType, state->box);
2415 initPullComFromPrevStep(cr, pull_work, masses, &pbc, state->x.rvec_array());
2416 updatePrevStepPullCom(pull_work, state);
2420 void finish_pull(struct pull_t* pull)
2422 check_external_potential_registration(pull);
2426 gmx_fio_fclose(pull->out_x);
2430 gmx_fio_fclose(pull->out_f);
2436 bool pull_have_potential(const pull_t& pull)
2438 return pull.bPotential;
2441 bool pull_have_constraint(const pull_t& pull)
2443 return pull.bConstraint;
2446 bool pull_have_constraint(const pull_params_t& pullParameters)
2448 for (int c = 0; c < pullParameters.ncoord; c++)
2450 if (pullParameters.coord[c].eType == epullCONSTRAINT)