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, 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.
52 #include "gromacs/commandline/filenm.h"
53 #include "gromacs/domdec/domdec_struct.h"
54 #include "gromacs/domdec/localatomset.h"
55 #include "gromacs/domdec/localatomsetmanager.h"
56 #include "gromacs/fileio/gmxfio.h"
57 #include "gromacs/gmxlib/network.h"
58 #include "gromacs/math/functions.h"
59 #include "gromacs/math/units.h"
60 #include "gromacs/math/vec.h"
61 #include "gromacs/math/vectypes.h"
62 #include "gromacs/mdlib/gmx_omp_nthreads.h"
63 #include "gromacs/mdtypes/commrec.h"
64 #include "gromacs/mdtypes/forceoutput.h"
65 #include "gromacs/mdtypes/inputrec.h"
66 #include "gromacs/mdtypes/md_enums.h"
67 #include "gromacs/mdtypes/state.h"
68 #include "gromacs/pbcutil/pbc.h"
69 #include "gromacs/topology/mtop_lookup.h"
70 #include "gromacs/topology/topology.h"
71 #include "gromacs/utility/cstringutil.h"
72 #include "gromacs/utility/exceptions.h"
73 #include "gromacs/utility/fatalerror.h"
74 #include "gromacs/utility/futil.h"
75 #include "gromacs/utility/gmxassert.h"
76 #include "gromacs/utility/mutex.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 static void low_get_pull_coord_dr(const struct pull_t* pull,
454 const pull_coord_work_t* pcrd,
461 const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
463 /* Only the first group can be an absolute reference, in that case nat=0 */
464 if (pgrp0->params.ind.empty())
466 for (int m = 0; m < DIM; m++)
468 xref[m] = pcrd->params.origin[m];
473 copy_dvec(xref, xrefr);
475 dvec dref = { 0, 0, 0 };
476 if (pcrd->params.eGeom == epullgDIRPBC)
478 for (int m = 0; m < DIM; m++)
480 dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
482 /* Add the reference position, so we use the correct periodic image */
483 dvec_inc(xrefr, dref);
486 pbc_dx_d(pbc, xg, xrefr, dr);
488 bool directional = pull_coordinate_is_directional(&pcrd->params);
490 for (int m = 0; m < DIM; m++)
492 dr[m] *= pcrd->params.dim[m];
493 if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
495 dr2 += dr[m] * dr[m];
498 /* Check if we are close to switching to another periodic image */
499 if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
501 /* Note that technically there is no issue with switching periodic
502 * image, as pbc_dx_d returns the distance to the closest periodic
503 * image. However in all cases where periodic image switches occur,
504 * the pull results will be useless in practice.
507 "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
508 "box size (%f).\n%s",
509 pcrd->params.group[0],
510 pcrd->params.group[1],
512 sqrt(0.98 * 0.98 * max_dist2),
513 pcrd->params.eGeom == epullgDIR
514 ? "You might want to consider using \"pull-geometry = "
515 "direction-periodic\" instead.\n"
519 if (pcrd->params.eGeom == epullgDIRPBC)
525 /* This function returns the distance based on the contents of the pull struct.
526 * pull->coord[coord_ind].dr, and potentially vector, are updated.
528 static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
530 pull_coord_work_t* pcrd = &pull->coord[coord_ind];
531 PullCoordSpatialData& spatialData = pcrd->spatialData;
534 /* With AWH pulling we allow for periodic pulling with geometry=direction.
535 * TODO: Store a periodicity flag instead of checking for external pull provider.
537 if (pcrd->params.eGeom == epullgDIRPBC
538 || (pcrd->params.eGeom == epullgDIR && pcrd->params.eType == epullEXTERNAL))
544 md2 = static_cast<double>(max_pull_distance2(pcrd, pbc));
547 if (pcrd->params.eGeom == epullgDIRRELATIVE)
549 /* We need to determine the pull vector */
553 const pull_group_work_t& pgrp2 = pull->group[pcrd->params.group[2]];
554 const pull_group_work_t& pgrp3 = pull->group[pcrd->params.group[3]];
556 pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
558 for (m = 0; m < DIM; m++)
560 vec[m] *= pcrd->params.dim[m];
562 spatialData.vec_len = dnorm(vec);
563 for (m = 0; m < DIM; m++)
565 spatialData.vec[m] = vec[m] / spatialData.vec_len;
570 "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
577 spatialData.vec[ZZ]);
581 pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
582 pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
584 low_get_pull_coord_dr(pull,
588 pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x,
592 if (pcrd->params.ngroup >= 4)
594 pull_group_work_t *pgrp2, *pgrp3;
595 pgrp2 = &pull->group[pcrd->params.group[2]];
596 pgrp3 = &pull->group[pcrd->params.group[3]];
598 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, md2, spatialData.dr23);
600 if (pcrd->params.ngroup >= 6)
602 pull_group_work_t *pgrp4, *pgrp5;
603 pgrp4 = &pull->group[pcrd->params.group[4]];
604 pgrp5 = &pull->group[pcrd->params.group[5]];
606 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, md2, spatialData.dr45);
610 /* Modify x so that it is periodic in [-pi, pi)
611 * It is assumed that x is in [-3pi, 3pi) so that x
612 * needs to be shifted by at most one period.
614 static void make_periodic_2pi(double* x)
626 /* This function should always be used to modify pcrd->value_ref */
627 static void low_set_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double value_ref)
629 GMX_ASSERT(pcrd->params.eType != epullEXTERNAL,
630 "The pull coord reference value should not be used with type external-potential");
632 if (pcrd->params.eGeom == epullgDIST)
637 "Pull reference distance for coordinate %d (%f) needs to be non-negative",
642 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgANGLEAXIS)
644 if (value_ref < 0 || value_ref > M_PI)
647 "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
648 "interval [0,180] deg",
650 value_ref * pull_conversion_factor_internal2userinput(&pcrd->params));
653 else if (pcrd->params.eGeom == epullgDIHEDRAL)
655 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
656 make_periodic_2pi(&value_ref);
659 pcrd->value_ref = value_ref;
662 static void update_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double t)
664 /* With zero rate the reference value is set initially and doesn't change */
665 if (pcrd->params.rate != 0)
667 double value_ref = (pcrd->params.init + pcrd->params.rate * t)
668 * pull_conversion_factor_userinput2internal(&pcrd->params);
669 low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
673 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
674 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
677 dvec dr32; /* store instead of dr23? */
679 dsvmul(-1, spatialData->dr23, dr32);
680 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
681 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
682 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
684 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
685 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
686 * Note 2: the angle between the plane normal vectors equals pi only when
687 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
688 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
690 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
694 /* Calculates pull->coord[coord_ind].value.
695 * This function also updates pull->coord[coord_ind].dr.
697 static void get_pull_coord_distance(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
699 pull_coord_work_t* pcrd;
702 pcrd = &pull->coord[coord_ind];
704 get_pull_coord_dr(pull, coord_ind, pbc);
706 PullCoordSpatialData& spatialData = pcrd->spatialData;
708 switch (pcrd->params.eGeom)
711 /* Pull along the vector between the com's */
712 spatialData.value = dnorm(spatialData.dr01);
716 case epullgDIRRELATIVE:
719 spatialData.value = 0;
720 for (m = 0; m < DIM; m++)
722 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
726 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
728 case epullgDIHEDRAL: spatialData.value = get_dihedral_angle_coord(&spatialData); break;
729 case epullgANGLEAXIS:
730 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
732 default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
736 /* Returns the deviation from the reference value.
737 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
739 static double get_pull_coord_deviation(struct pull_t* pull, int coord_ind, const t_pbc* pbc, double t)
741 pull_coord_work_t* pcrd;
744 pcrd = &pull->coord[coord_ind];
746 /* Update the reference value before computing the distance,
747 * since it is used in the distance computation with periodic pulling.
749 update_pull_coord_reference_value(pcrd, coord_ind, t);
751 get_pull_coord_distance(pull, coord_ind, pbc);
753 get_pull_coord_distance(pull, coord_ind, pbc);
755 /* Determine the deviation */
756 dev = pcrd->spatialData.value - pcrd->value_ref;
758 /* Check that values are allowed */
759 if (pcrd->params.eGeom == epullgDIST && pcrd->spatialData.value == 0)
761 /* With no vector we can not determine the direction for the force,
762 * so we set the force to zero.
766 else if (pcrd->params.eGeom == epullgDIHEDRAL)
768 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
769 Thus, the unwrapped deviation is here in (-2pi, 2pi].
770 After making it periodic, the deviation will be in [-pi, pi). */
771 make_periodic_2pi(&dev);
777 double get_pull_coord_value(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
779 get_pull_coord_distance(pull, coord_ind, pbc);
781 return pull->coord[coord_ind].spatialData.value;
784 void clear_pull_forces(pull_t* pull)
786 /* Zeroing the forces is only required for constraint pulling.
787 * It can happen that multiple constraint steps need to be applied
788 * and therefore the constraint forces need to be accumulated.
790 for (pull_coord_work_t& coord : pull->coord)
792 coord.scalarForce = 0;
796 /* Apply constraint using SHAKE */
798 do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
801 dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
802 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
803 dvec* rnew; /* current 'new' positions of the groups */
804 double* dr_tot; /* the total update of the coords */
807 double lambda, rm, invdt = 0;
808 gmx_bool bConverged_all, bConverged = FALSE;
809 int niter = 0, ii, j, m, max_iter = 100;
813 snew(r_ij, pull->coord.size());
814 snew(dr_tot, pull->coord.size());
816 snew(rnew, pull->group.size());
818 /* copy the current unconstrained positions for use in iterations. We
819 iterate until rinew[i] and rjnew[j] obey the constraints. Then
820 rinew - pull.x_unc[i] is the correction dr to group i */
821 for (size_t g = 0; g < pull->group.size(); g++)
823 copy_dvec(pull->group[g].xp, rnew[g]);
826 /* Determine the constraint directions from the old positions */
827 for (size_t c = 0; c < pull->coord.size(); c++)
829 pull_coord_work_t* pcrd;
831 pcrd = &pull->coord[c];
833 if (pcrd->params.eType != epullCONSTRAINT)
838 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
839 * We don't modify dr and value anymore, so these values are also used
842 get_pull_coord_distance(pull, c, pbc);
844 const PullCoordSpatialData& spatialData = pcrd->spatialData;
848 "Pull coord %zu dr %f %f %f\n",
850 spatialData.dr01[XX],
851 spatialData.dr01[YY],
852 spatialData.dr01[ZZ]);
855 if (pcrd->params.eGeom == epullgDIR || pcrd->params.eGeom == epullgDIRPBC)
857 /* Select the component along vec */
859 for (m = 0; m < DIM; m++)
861 a += spatialData.vec[m] * spatialData.dr01[m];
863 for (m = 0; m < DIM; m++)
865 r_ij[c][m] = a * spatialData.vec[m];
870 copy_dvec(spatialData.dr01, r_ij[c]);
873 if (dnorm2(r_ij[c]) == 0)
876 "Distance for pull coordinate %zu is zero with constraint pulling, which is "
882 bConverged_all = FALSE;
883 while (!bConverged_all && niter < max_iter)
885 bConverged_all = TRUE;
887 /* loop over all constraints */
888 for (size_t c = 0; c < pull->coord.size(); c++)
890 pull_coord_work_t* pcrd;
891 pull_group_work_t *pgrp0, *pgrp1;
894 pcrd = &pull->coord[c];
896 if (pcrd->params.eType != epullCONSTRAINT)
901 update_pull_coord_reference_value(pcrd, c, t);
903 pgrp0 = &pull->group[pcrd->params.group[0]];
904 pgrp1 = &pull->group[pcrd->params.group[1]];
906 /* Get the current difference vector */
907 low_get_pull_coord_dr(
908 pull, pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], -1, unc_ij);
912 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
915 rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
917 switch (pcrd->params.eGeom)
920 if (pcrd->value_ref <= 0)
924 "The pull constraint reference distance for group %zu is <= 0 (%f)",
930 double q, c_a, c_b, c_c;
932 c_a = diprod(r_ij[c], r_ij[c]);
933 c_b = diprod(unc_ij, r_ij[c]) * 2;
934 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
938 q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
943 q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
949 fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b, c_c, lambda);
953 /* The position corrections dr due to the constraints */
954 dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
955 dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
956 dr_tot[c] += -lambda * dnorm(r_ij[c]);
961 /* A 1-dimensional constraint along a vector */
963 for (m = 0; m < DIM; m++)
965 vec[m] = pcrd->spatialData.vec[m];
966 a += unc_ij[m] * vec[m];
968 /* Select only the component along the vector */
969 dsvmul(a, vec, unc_ij);
970 lambda = a - pcrd->value_ref;
973 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
976 /* The position corrections dr due to the constraints */
977 dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
978 dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
979 dr_tot[c] += -lambda;
981 default: gmx_incons("Invalid enumeration value for eGeom");
989 g0 = pcrd->params.group[0];
990 g1 = pcrd->params.group[1];
991 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
992 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
994 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1003 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
1012 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1022 /* Update the COMs with dr */
1023 dvec_inc(rnew[pcrd->params.group[1]], dr1);
1024 dvec_inc(rnew[pcrd->params.group[0]], dr0);
1027 /* Check if all constraints are fullfilled now */
1028 for (pull_coord_work_t& coord : pull->coord)
1030 if (coord.params.eType != epullCONSTRAINT)
1035 low_get_pull_coord_dr(
1036 pull, &coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], -1, unc_ij);
1038 switch (coord.params.eGeom)
1041 bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1046 for (m = 0; m < DIM; m++)
1048 vec[m] = coord.spatialData.vec[m];
1050 inpr = diprod(unc_ij, vec);
1051 dsvmul(inpr, vec, unc_ij);
1052 bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1061 "Pull constraint not converged: "
1063 "d_ref = %f, current d = %f\n",
1064 coord.params.group[0],
1065 coord.params.group[1],
1070 bConverged_all = FALSE;
1075 /* if after all constraints are dealt with and bConverged is still TRUE
1076 we're finished, if not we do another iteration */
1078 if (niter > max_iter)
1080 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1083 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1090 /* update atoms in the groups */
1091 for (size_t g = 0; g < pull->group.size(); g++)
1093 const pull_group_work_t* pgrp;
1096 pgrp = &pull->group[g];
1098 /* get the final constraint displacement dr for group g */
1099 dvec_sub(rnew[g], pgrp->xp, dr);
1101 if (dnorm2(dr) == 0)
1103 /* No displacement, no update necessary */
1107 /* update the atom positions */
1108 auto localAtomIndices = pgrp->atomSet.localIndex();
1110 for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1112 ii = localAtomIndices[j];
1113 if (!pgrp->localWeights.empty())
1115 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1117 for (m = 0; m < DIM; m++)
1123 for (m = 0; m < DIM; m++)
1125 v[ii][m] += invdt * tmp[m];
1131 /* calculate the constraint forces, used for output and virial only */
1132 for (size_t c = 0; c < pull->coord.size(); c++)
1134 pull_coord_work_t* pcrd;
1136 pcrd = &pull->coord[c];
1138 if (pcrd->params.eType != epullCONSTRAINT)
1143 /* Accumulate the forces, in case we have multiple constraint steps */
1146 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1148 pcrd->scalarForce += force;
1150 if (vir != nullptr && pcrd->params.eGeom != epullgDIRPBC && bMaster)
1154 /* Add the pull contribution to the virial */
1155 /* We have already checked above that r_ij[c] != 0 */
1156 f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1158 for (j = 0; j < DIM; j++)
1160 for (m = 0; m < DIM; m++)
1162 vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1168 /* finished! I hope. Give back some memory */
1174 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1176 for (int j = 0; j < DIM; j++)
1178 for (int m = 0; m < DIM; m++)
1180 vir[j][m] -= 0.5 * f[j] * dr[m];
1185 /* Adds the pull contribution to the virial */
1186 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1188 if (vir != nullptr && pcrd.params.eGeom != epullgDIRPBC)
1190 /* Add the pull contribution for each distance vector to the virial. */
1191 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1192 if (pcrd.params.ngroup >= 4)
1194 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1196 if (pcrd.params.ngroup >= 6)
1198 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1203 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1211 k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1212 dkdl = pcrd->params.kB - pcrd->params.k;
1214 switch (pcrd->params.eType)
1217 case epullFLATBOTTOM:
1218 case epullFLATBOTTOMHIGH:
1219 /* The only difference between an umbrella and a flat-bottom
1220 * potential is that a flat-bottom is zero above or below
1221 the reference value.
1223 if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0)
1224 || (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
1229 pcrd->scalarForce = -k * dev;
1230 *V += 0.5 * k * gmx::square(dev);
1231 *dVdl += 0.5 * dkdl * gmx::square(dev);
1234 pcrd->scalarForce = -k;
1235 *V += k * pcrd->spatialData.value;
1236 *dVdl += dkdl * pcrd->spatialData.value;
1240 "the scalar pull force should not be calculated internally for pull type "
1242 default: gmx_incons("Unsupported pull type in do_pull_pot");
1246 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1248 const t_pull_coord& params = pcrd.params;
1249 const PullCoordSpatialData& spatialData = pcrd.spatialData;
1251 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1252 PullCoordVectorForces forces;
1254 if (params.eGeom == epullgDIST)
1256 double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1257 for (int m = 0; m < DIM; m++)
1259 forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1262 else if (params.eGeom == epullgANGLE)
1265 double cos_theta, cos_theta2;
1267 cos_theta = cos(spatialData.value);
1268 cos_theta2 = gmx::square(cos_theta);
1270 /* The force at theta = 0, pi is undefined so we should not apply any force.
1271 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1272 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1273 * have to check for this before dividing by their norm below.
1277 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1278 * and another vector parallel to dr23:
1279 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1280 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1282 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1283 double b = a * cos_theta;
1284 double invdr01 = 1. / dnorm(spatialData.dr01);
1285 double invdr23 = 1. / dnorm(spatialData.dr23);
1286 dvec normalized_dr01, normalized_dr23;
1287 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1288 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1290 for (int m = 0; m < DIM; m++)
1292 /* Here, f_scal is -dV/dtheta */
1294 pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1296 pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1301 /* No forces to apply for ill-defined cases*/
1302 clear_dvec(forces.force01);
1303 clear_dvec(forces.force23);
1306 else if (params.eGeom == epullgANGLEAXIS)
1308 double cos_theta, cos_theta2;
1310 /* The angle-axis force is exactly the same as the angle force (above) except that in
1311 this case the second vector (dr23) is replaced by the pull vector. */
1312 cos_theta = cos(spatialData.value);
1313 cos_theta2 = gmx::square(cos_theta);
1319 dvec normalized_dr01;
1321 invdr01 = 1. / dnorm(spatialData.dr01);
1322 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1323 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1326 for (int m = 0; m < DIM; m++)
1329 pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1334 clear_dvec(forces.force01);
1337 else if (params.eGeom == epullgDIHEDRAL)
1339 double m2, n2, tol, sqrdist_32;
1341 /* Note: there is a small difference here compared to the
1342 dihedral force calculations in the bondeds (ref: Bekker 1994).
1343 There rij = ri - rj, while here dr01 = r1 - r0.
1344 However, all distance vectors occur in form of cross or inner products
1345 so that two signs cancel and we end up with the same expressions.
1346 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1348 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1349 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1350 dsvmul(-1, spatialData.dr23, dr32);
1351 sqrdist_32 = diprod(dr32, dr32);
1352 tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1353 if ((m2 > tol) && (n2 > tol))
1355 double a_01, a_23_01, a_23_45, a_45;
1356 double inv_dist_32, inv_sqrdist_32, dist_32;
1358 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1359 inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1360 dist_32 = sqrdist_32 * inv_dist_32;
1362 /* Forces on groups 0, 1 */
1363 a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
1364 dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1366 /* Forces on groups 4, 5 */
1367 a_45 = -pcrd.scalarForce * dist_32 / n2;
1368 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1370 /* Force on groups 2, 3 (defining the axis) */
1371 a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1372 a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1373 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1374 dsvmul(a_23_45, forces.force45, v);
1375 dvec_sub(u, v, forces.force23); /* force on group 3 */
1379 /* No force to apply for ill-defined cases */
1380 clear_dvec(forces.force01);
1381 clear_dvec(forces.force23);
1382 clear_dvec(forces.force45);
1387 for (int m = 0; m < DIM; m++)
1389 forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1397 /* We use a global mutex for locking access to the pull data structure
1398 * during registration of external pull potential providers.
1399 * We could use a different, local mutex for each pull object, but the overhead
1400 * is extremely small here and registration is only done during initialization.
1402 static gmx::Mutex registrationMutex;
1404 using Lock = gmx::lock_guard<gmx::Mutex>;
1406 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1408 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1409 GMX_RELEASE_ASSERT(provider != nullptr,
1410 "register_external_pull_potential called with NULL as provider name");
1412 if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1415 "Module '%s' attempted to register an external potential for pull coordinate %d "
1416 "which is out of the pull coordinate range %d - %zu\n",
1420 pull->coord.size());
1423 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1425 if (pcrd->params.eType != epullEXTERNAL)
1429 "Module '%s' attempted to register an external potential for pull coordinate %d "
1430 "which of type '%s', whereas external potentials are only supported with type '%s'",
1433 epull_names[pcrd->params.eType],
1434 epull_names[epullEXTERNAL]);
1437 GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
1438 "The external potential provider string for a pull coordinate is NULL");
1440 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
1443 "Module '%s' attempted to register an external potential for pull coordinate %d "
1444 "which expects the external potential to be provided by a module named '%s'",
1447 pcrd->params.externalPotentialProvider.c_str());
1450 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1451 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1452 * pull->numUnregisteredExternalPotentials.
1454 Lock registrationLock(registrationMutex);
1456 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1459 "Module '%s' attempted to register an external potential for pull coordinate %d "
1465 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1466 pull->numUnregisteredExternalPotentials--;
1468 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1469 "Negative unregistered potentials, the pull code is inconsistent");
1473 static void check_external_potential_registration(const struct pull_t* pull)
1475 if (pull->numUnregisteredExternalPotentials > 0)
1478 for (c = 0; c < pull->coord.size(); c++)
1480 if (pull->coord[c].params.eType == epullEXTERNAL
1481 && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1487 GMX_RELEASE_ASSERT(c < pull->coord.size(),
1488 "Internal inconsistency in the pull potential provider counting");
1491 "No external provider for external pull potentials have been provided for %d "
1492 "pull coordinates. The first coordinate without provider is number %zu, which "
1493 "expects a module named '%s' to provide the external potential.",
1494 pull->numUnregisteredExternalPotentials,
1496 pull->coord[c].params.externalPotentialProvider.c_str());
1501 /* Pull takes care of adding the forces of the external potential.
1502 * The external potential module has to make sure that the corresponding
1503 * potential energy is added either to the pull term or to a term
1504 * specific to the external module.
1506 void apply_external_pull_coord_force(struct pull_t* pull,
1510 gmx::ForceWithVirial* forceWithVirial)
1512 pull_coord_work_t* pcrd;
1514 GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1515 "apply_external_pull_coord_force called with coord_index out of range");
1517 if (pull->comm.bParticipate)
1519 pcrd = &pull->coord[coord_index];
1522 pcrd->params.eType == epullEXTERNAL,
1523 "The pull force can only be set externally on pull coordinates of external type");
1525 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1526 "apply_external_pull_coord_force called for an unregistered pull coordinate");
1529 pcrd->scalarForce = coord_force;
1531 /* Calculate the forces on the pull groups */
1532 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1534 /* Add the forces for this coordinate to the total virial and force */
1535 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1537 matrix virial = { { 0 } };
1538 add_virial_coord(virial, *pcrd, pullCoordForces);
1539 forceWithVirial->addVirialContribution(virial);
1543 pull, coord_index, pullCoordForces, masses, as_rvec_array(forceWithVirial->force_.data()));
1546 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1549 /* Calculate the pull potential and scalar force for a pull coordinate.
1550 * Returns the vector forces for the pull coordinate.
1552 static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull,
1561 pull_coord_work_t& pcrd = pull->coord[coord_ind];
1563 assert(pcrd.params.eType != epullCONSTRAINT);
1565 double dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
1567 calc_pull_coord_scalar_force_and_potential(&pcrd, dev, lambda, V, dVdl);
1569 PullCoordVectorForces pullCoordForces = calculateVectorForces(pcrd);
1571 add_virial_coord(vir, pcrd, pullCoordForces);
1573 return pullCoordForces;
1576 real pull_potential(struct pull_t* pull,
1579 const t_commrec* cr,
1583 gmx::ForceWithVirial* force,
1588 assert(pull != nullptr);
1590 /* Ideally we should check external potential registration only during
1591 * the initialization phase, but that requires another function call
1592 * that should be done exactly in the right order. So we check here.
1594 check_external_potential_registration(pull);
1596 if (pull->comm.bParticipate)
1600 pull_calc_coms(cr, pull, masses, pbc, t, x, nullptr);
1602 rvec* f = as_rvec_array(force->force_.data());
1603 matrix virial = { { 0 } };
1604 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1605 for (size_t c = 0; c < pull->coord.size(); c++)
1607 pull_coord_work_t* pcrd;
1608 pcrd = &pull->coord[c];
1610 /* For external potential the force is assumed to be given by an external module by a
1611 call to apply_pull_coord_external_force */
1612 if (pcrd->params.eType == epullCONSTRAINT || pcrd->params.eType == epullEXTERNAL)
1617 PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1618 pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1620 /* Distribute the force over the atoms in the pulled groups */
1621 apply_forces_coord(pull, c, pullCoordForces, masses, f);
1626 force->addVirialContribution(virial);
1631 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1632 "Too few or too many external pull potentials have been applied the previous step");
1633 /* All external pull potentials still need to be applied */
1634 pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1636 return (MASTER(cr) ? V : 0.0);
1639 void pull_constraint(struct pull_t* pull,
1642 const t_commrec* cr,
1650 assert(pull != nullptr);
1652 if (pull->comm.bParticipate)
1654 pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
1656 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1660 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1664 gmx_bool bMustParticipate;
1670 /* We always make the master node participate, such that it can do i/o,
1671 * add the virial and to simplify MC type extensions people might have.
1673 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1675 for (pull_group_work_t& group : pull->group)
1677 if (!group.globalWeights.empty())
1679 group.localWeights.resize(group.atomSet.numAtomsLocal());
1680 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1682 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1686 GMX_ASSERT(bMustParticipate || dd != nullptr,
1687 "Either all ranks (including this rank) participate, or we use DD and need to "
1688 "have access to dd here");
1690 /* We should participate if we have pull or pbc atoms */
1691 if (!bMustParticipate
1692 && (group.atomSet.numAtomsLocal() > 0
1693 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1695 bMustParticipate = TRUE;
1699 if (!comm->bParticipateAll)
1701 /* Keep currently not required ranks in the communicator
1702 * if they needed to participate up to 20 decompositions ago.
1703 * This avoids frequent rebuilds due to atoms jumping back and forth.
1705 const int64_t history_count = 20;
1706 gmx_bool bWillParticipate;
1709 /* Increase the decomposition counter for the current call */
1710 comm->setup_count++;
1712 if (bMustParticipate)
1714 comm->must_count = comm->setup_count;
1719 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1721 if (debug && dd != nullptr)
1724 "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
1726 gmx::boolToString(bMustParticipate),
1727 gmx::boolToString(bWillParticipate));
1730 if (bWillParticipate)
1732 /* Count the number of ranks that we want to have participating */
1734 /* Count the number of ranks that need to be added */
1735 count[1] = comm->bParticipate ? 0 : 1;
1743 /* The cost of this global operation will be less that the cost
1744 * of the extra MPI_Comm_split calls that we can avoid.
1746 gmx_sumi(2, count, cr);
1748 /* If we are missing ranks or if we have 20% more ranks than needed
1749 * we make a new sub-communicator.
1751 if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1755 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1759 if (comm->mpi_comm_com != MPI_COMM_NULL)
1761 MPI_Comm_free(&comm->mpi_comm_com);
1764 /* This might be an extremely expensive operation, so we try
1765 * to avoid this splitting as much as possible.
1767 assert(dd != nullptr);
1768 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1771 comm->bParticipate = bWillParticipate;
1772 comm->nparticipate = count[0];
1774 /* When we use the previous COM for PBC, we need to broadcast
1775 * the previous COM to ranks that have joined the communicator.
1777 for (pull_group_work_t& group : pull->group)
1779 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1781 GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1782 "The master rank has to participate, as it should pass an up to "
1784 "to bcast here as well as to e.g. checkpointing");
1786 gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
1792 /* Since the PBC of atoms might have changed, we need to update the PBC */
1793 pull->bSetPBCatoms = TRUE;
1796 static void init_pull_group_index(FILE* fplog,
1797 const t_commrec* cr,
1799 pull_group_work_t* pg,
1800 gmx_bool bConstraint,
1801 const ivec pulldim_con,
1802 const gmx_mtop_t* mtop,
1803 const t_inputrec* ir,
1806 /* With EM and BD there are no masses in the integrator.
1807 * But we still want to have the correct mass-weighted COMs.
1808 * So we store the real masses in the weights.
1810 const bool setWeights =
1811 (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
1813 /* In parallel, store we need to extract localWeights from weights at DD time */
1814 std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1816 const SimulationGroups& groups = mtop->groups;
1818 /* Count frozen dimensions and (weighted) mass */
1824 for (int i = 0; i < int(pg->params.ind.size()); i++)
1826 int ii = pg->params.ind[i];
1827 if (bConstraint && ir->opts.nFreeze)
1829 for (int d = 0; d < DIM; d++)
1831 if (pulldim_con[d] == 1
1832 && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1838 const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1840 if (ir->efep == efepNO)
1846 m = (1 - lambda) * atom.m + lambda * atom.mB;
1849 if (!pg->params.weight.empty())
1851 w = pg->params.weight[i];
1857 if (EI_ENERGY_MINIMIZATION(ir->eI))
1859 /* Move the mass to the weight */
1863 else if (ir->eI == eiBD)
1866 if (ir->bd_fric != 0.0F)
1868 mbd = ir->bd_fric * ir->delta_t;
1872 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1874 mbd = ir->delta_t / ir->opts.tau_t[0];
1879 / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1887 weights.push_back(w);
1891 wwmass += m * w * w;
1896 /* We can have single atom groups with zero mass with potential pulling
1897 * without cosine weighting.
1899 if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1901 /* With one atom the mass doesn't matter */
1907 "The total%s mass of pull group %d is zero",
1908 !pg->params.weight.empty() ? " weighted" : "",
1914 fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
1915 if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1917 fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1919 if (pg->epgrppbc == epgrppbcCOS)
1921 fprintf(fplog, ", cosine weighting will be used");
1923 fprintf(fplog, "\n");
1928 /* A value != 0 signals not frozen, it is updated later */
1934 for (int d = 0; d < DIM; d++)
1936 ndim += pulldim_con[d] * pg->params.ind.size();
1938 if (fplog && nfrozen > 0 && nfrozen < ndim)
1941 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1942 " that are subject to pulling are frozen.\n"
1943 " For constraint pulling the whole group will be frozen.\n\n",
1951 struct pull_t* init_pull(FILE* fplog,
1952 const pull_params_t* pull_params,
1953 const t_inputrec* ir,
1954 const gmx_mtop_t* mtop,
1955 const t_commrec* cr,
1956 gmx::LocalAtomSetManager* atomSets,
1959 struct pull_t* pull;
1964 /* Copy the pull parameters */
1965 pull->params = *pull_params;
1967 for (int i = 0; i < pull_params->ngroup; ++i)
1969 pull->group.emplace_back(pull_params->group[i],
1970 atomSets->add(pull_params->group[i].ind),
1971 pull_params->bSetPbcRefToPrevStepCOM);
1974 if (cr != nullptr && DOMAINDECOMP(cr))
1976 /* Set up the global to local atom mapping for PBC atoms */
1977 for (pull_group_work_t& group : pull->group)
1979 if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
1981 /* pbcAtomSet consists of a single atom */
1982 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
1983 atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
1988 pull->bPotential = FALSE;
1989 pull->bConstraint = FALSE;
1990 pull->bCylinder = FALSE;
1991 pull->bAngle = FALSE;
1992 pull->bXOutAverage = pull_params->bXOutAverage;
1993 pull->bFOutAverage = pull_params->bFOutAverage;
1995 GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
1996 "pull group 0 is an absolute reference group and should not contain atoms");
1998 pull->numCoordinatesWithExternalPotential = 0;
2000 for (int c = 0; c < pull->params.ncoord; c++)
2002 /* Construct a pull coordinate, copying all coordinate parameters */
2003 pull->coord.emplace_back(pull_params->coord[c]);
2005 pull_coord_work_t* pcrd = &pull->coord.back();
2007 switch (pcrd->params.eGeom)
2010 case epullgDIRRELATIVE: /* Direction vector is determined at each step */
2012 case epullgDIHEDRAL: break;
2016 case epullgANGLEAXIS:
2017 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2020 /* We allow reading of newer tpx files with new pull geometries,
2021 * but with the same tpx format, with old code. A new geometry
2022 * only adds a new enum value, which will be out of range for
2023 * old code. The first place we need to generate an error is
2024 * here, since the pull code can't handle this.
2025 * The enum can be used before arriving here only for printing
2026 * the string corresponding to the geometry, which will result
2027 * in printing "UNDEFINED".
2030 "Pull geometry not supported for pull coordinate %d. The geometry enum "
2031 "%d in the input is larger than that supported by the code (up to %d). "
2032 "You are probably reading a tpr file generated with a newer version of "
2033 "Gromacs with an binary from an older version of Gromacs.",
2039 if (pcrd->params.eType == epullCONSTRAINT)
2041 /* Check restrictions of the constraint pull code */
2042 if (pcrd->params.eGeom == epullgCYL || pcrd->params.eGeom == epullgDIRRELATIVE
2043 || pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
2044 || pcrd->params.eGeom == epullgANGLEAXIS)
2047 "Pulling of type %s can not be combined with geometry %s. Consider using "
2049 epull_names[pcrd->params.eType],
2050 epullg_names[pcrd->params.eGeom],
2051 epull_names[epullUMBRELLA]);
2056 "Constraint pulling can not be combined with multiple time stepping");
2058 pull->bConstraint = TRUE;
2062 pull->bPotential = TRUE;
2065 if (pcrd->params.eGeom == epullgCYL)
2067 pull->bCylinder = TRUE;
2069 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
2070 || pcrd->params.eGeom == epullgANGLEAXIS)
2072 pull->bAngle = TRUE;
2075 /* We only need to calculate the plain COM of a group
2076 * when it is not only used as a cylinder group.
2077 * Also the absolute reference group 0 needs no COM computation.
2079 for (int i = 0; i < pcrd->params.ngroup; i++)
2081 int groupIndex = pcrd->params.group[i];
2082 if (groupIndex > 0 && !(pcrd->params.eGeom == epullgCYL && i == 0))
2084 pull->group[groupIndex].needToCalcCom = true;
2088 /* With non-zero rate the reference value is set at every step */
2089 if (pcrd->params.rate == 0)
2091 /* Initialize the constant reference value */
2092 if (pcrd->params.eType != epullEXTERNAL)
2094 low_set_pull_coord_reference_value(
2095 pcrd, c, pcrd->params.init * pull_conversion_factor_userinput2internal(&pcrd->params));
2099 /* With an external pull potential, the reference value loses
2100 * it's meaning and should not be used. Setting it to zero
2101 * makes any terms dependent on it disappear.
2102 * The only issue this causes is that with cylinder pulling
2103 * no atoms of the cylinder group within the cylinder radius
2104 * should be more than half a box length away from the COM of
2105 * the pull group along the axial direction.
2107 pcrd->value_ref = 0.0;
2111 if (pcrd->params.eType == epullEXTERNAL)
2114 pcrd->params.rate == 0,
2115 "With an external potential, a pull coordinate should have rate = 0");
2117 /* This external potential needs to be registered later */
2118 pull->numCoordinatesWithExternalPotential++;
2120 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2123 pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
2124 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2126 pull->pbcType = ir->pbcType;
2127 switch (pull->pbcType)
2129 case PbcType::No: pull->npbcdim = 0; break;
2130 case PbcType::XY: pull->npbcdim = 2; break;
2131 default: pull->npbcdim = 3; break;
2136 gmx_bool bAbs, bCos;
2139 for (const pull_coord_work_t& coord : pull->coord)
2141 if (pull->group[coord.params.group[0]].params.ind.empty()
2142 || pull->group[coord.params.group[1]].params.ind.empty())
2148 fprintf(fplog, "\n");
2149 if (pull->bPotential)
2151 fprintf(fplog, "Will apply potential COM pulling\n");
2153 if (pull->bConstraint)
2155 fprintf(fplog, "Will apply constraint COM pulling\n");
2157 // Don't include the reference group 0 in output, so we report ngroup-1
2158 int numRealGroups = pull->group.size() - 1;
2159 GMX_RELEASE_ASSERT(numRealGroups > 0,
2160 "The reference absolute position pull group should always be present");
2162 "with %zu pull coordinate%s and %d group%s\n",
2164 pull->coord.size() == 1 ? "" : "s",
2166 numRealGroups == 1 ? "" : "s");
2169 fprintf(fplog, "with an absolute reference\n");
2172 // Don't include the reference group 0 in loop
2173 for (size_t g = 1; g < pull->group.size(); g++)
2175 if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
2177 /* We are using cosine weighting */
2178 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2184 please_cite(fplog, "Engin2010");
2188 pull->bRefAt = FALSE;
2190 for (size_t g = 0; g < pull->group.size(); g++)
2192 pull_group_work_t* pgrp;
2194 pgrp = &pull->group[g];
2195 if (!pgrp->params.ind.empty())
2197 /* There is an issue when a group is used in multiple coordinates
2198 * and constraints are applied in different dimensions with atoms
2199 * frozen in some, but not all dimensions.
2200 * Since there is only one mass array per group, we can't have
2201 * frozen/non-frozen atoms for different coords at the same time.
2202 * But since this is a very exotic case, we don't check for this.
2203 * A warning is printed in init_pull_group_index.
2206 gmx_bool bConstraint;
2207 ivec pulldim, pulldim_con;
2209 /* Loop over all pull coordinates to see along which dimensions
2210 * this group is pulled and if it is involved in constraints.
2212 bConstraint = FALSE;
2213 clear_ivec(pulldim);
2214 clear_ivec(pulldim_con);
2215 for (const pull_coord_work_t& coord : pull->coord)
2217 gmx_bool bGroupUsed = FALSE;
2218 for (int gi = 0; gi < coord.params.ngroup; gi++)
2220 if (coord.params.group[gi] == static_cast<int>(g))
2228 for (int m = 0; m < DIM; m++)
2230 if (coord.params.dim[m] == 1)
2234 if (coord.params.eType == epullCONSTRAINT)
2244 /* Determine if we need to take PBC into account for calculating
2245 * the COM's of the pull groups.
2247 switch (pgrp->epgrppbc)
2249 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2251 if (!pgrp->params.weight.empty())
2254 "Pull groups can not have relative weights and cosine weighting "
2257 for (int m = 0; m < DIM; m++)
2259 if (m < pull->npbcdim && pulldim[m] == 1)
2261 if (pull->cosdim >= 0 && pull->cosdim != m)
2264 "Can only use cosine weighting with pulling in one "
2265 "dimension (use mdp option pull-coord?-dim)");
2271 case epgrppbcNONE: break;
2274 /* Set the indices */
2275 init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2279 /* Absolute reference, set the inverse mass to zero.
2280 * This is only relevant (and used) with constraint pulling.
2287 /* If we use cylinder coordinates, do some initialising for them */
2288 if (pull->bCylinder)
2290 for (const pull_coord_work_t& coord : pull->coord)
2292 if (coord.params.eGeom == epullgCYL)
2294 if (pull->group[coord.params.group[0]].params.ind.empty())
2297 "A cylinder pull group is not supported when using absolute "
2301 const auto& referenceGroup = pull->group[coord.params.group[0]];
2302 pull->dyna.emplace_back(
2303 referenceGroup.params, referenceGroup.atomSet, pull->params.bSetPbcRefToPrevStepCOM);
2307 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2308 pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
2309 pull->comSums.resize(pull->nthreads);
2314 /* Use a sub-communicator when we have more than 32 ranks, but not
2315 * when we have an external pull potential, since then the external
2316 * potential provider expects each rank to have the coordinate.
2318 comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
2319 || pull->numCoordinatesWithExternalPotential > 0
2320 || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2321 /* This sub-commicator is not used with comm->bParticipateAll,
2322 * so we can always initialize it to NULL.
2324 comm->mpi_comm_com = MPI_COMM_NULL;
2325 comm->nparticipate = 0;
2326 comm->isMasterRank = (cr == nullptr || MASTER(cr));
2328 /* No MPI: 1 rank: all ranks pull */
2329 comm->bParticipateAll = TRUE;
2330 comm->isMasterRank = true;
2332 comm->bParticipate = comm->bParticipateAll;
2333 comm->setup_count = 0;
2334 comm->must_count = 0;
2336 if (!comm->bParticipateAll && fplog != nullptr)
2338 fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2341 comm->pbcAtomBuffer.resize(pull->group.size());
2342 comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
2343 if (pull->bCylinder)
2345 comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
2348 /* We still need to initialize the PBC reference coordinates */
2349 pull->bSetPBCatoms = TRUE;
2351 pull->out_x = nullptr;
2352 pull->out_f = nullptr;
2357 static void destroy_pull(struct pull_t* pull)
2360 if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2362 MPI_Comm_free(&pull->comm.mpi_comm_com);
2369 void preparePrevStepPullCom(const t_inputrec* ir,
2373 const t_state* state_global,
2374 const t_commrec* cr,
2375 bool startingFromCheckpoint)
2377 if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
2381 allocStatePrevStepPullCom(state, pull_work);
2382 if (startingFromCheckpoint)
2386 state->pull_com_prev_step = state_global->pull_com_prev_step;
2390 /* Only the master rank has the checkpointed COM from the previous step */
2391 gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(),
2392 &state->pull_com_prev_step[0],
2393 cr->mpi_comm_mygroup);
2395 setPrevStepPullComFromState(pull_work, state);
2400 set_pbc(&pbc, ir->pbcType, state->box);
2401 initPullComFromPrevStep(cr, pull_work, masses, &pbc, state->x.rvec_array());
2402 updatePrevStepPullCom(pull_work, state);
2406 void finish_pull(struct pull_t* pull)
2408 check_external_potential_registration(pull);
2412 gmx_fio_fclose(pull->out_x);
2416 gmx_fio_fclose(pull->out_f);
2422 bool pull_have_potential(const pull_t& pull)
2424 return pull.bPotential;
2427 bool pull_have_constraint(const pull_t& pull)
2429 return pull.bConstraint;
2432 bool pull_have_constraint(const pull_params_t& pullParameters)
2434 for (int c = 0; c < pullParameters.ncoord; c++)
2436 if (pullParameters.coord[c].eType == epullCONSTRAINT)