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,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
51 #include "gromacs/commandline/filenm.h"
52 #include "gromacs/domdec/domdec_struct.h"
53 #include "gromacs/domdec/localatomset.h"
54 #include "gromacs/domdec/localatomsetmanager.h"
55 #include "gromacs/fileio/gmxfio.h"
56 #include "gromacs/gmxlib/network.h"
57 #include "gromacs/math/functions.h"
58 #include "gromacs/math/units.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/math/vectypes.h"
61 #include "gromacs/mdlib/gmx_omp_nthreads.h"
62 #include "gromacs/mdtypes/commrec.h"
63 #include "gromacs/mdtypes/forceoutput.h"
64 #include "gromacs/mdtypes/inputrec.h"
65 #include "gromacs/mdtypes/md_enums.h"
66 #include "gromacs/mdtypes/mdatom.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)
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 = md->massT[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.nat == 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(), md, 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, md, 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 = md->massT[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]], md, f_perp, -1, f, pull->nthreads);
318 apply_forces_grp(&pull->group[pcrd->params.group[3]], md, 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], pcrd.spatialData.cyl_dev, md, forces.force01,
339 pcrd.scalarForce, -1, f, pull->nthreads);
341 /* Sum the force along the vector and the radial force */
343 for (int m = 0; m < DIM; m++)
345 f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
347 apply_forces_grp(&pull->group[pcrd.params.group[1]], md, f_tot, 1, f, pull->nthreads);
351 if (pcrd.params.eGeom == epullgDIRRELATIVE)
353 /* We need to apply the torque forces to the pull groups
354 * that define the pull vector.
356 apply_forces_vec_torque(pull, &pcrd, md, f);
359 if (pull->group[pcrd.params.group[0]].params.nat > 0)
361 apply_forces_grp(&pull->group[pcrd.params.group[0]], md, forces.force01, -1, f, pull->nthreads);
363 apply_forces_grp(&pull->group[pcrd.params.group[1]], md, forces.force01, 1, f, pull->nthreads);
365 if (pcrd.params.ngroup >= 4)
367 apply_forces_grp(&pull->group[pcrd.params.group[2]], md, forces.force23, -1, f, pull->nthreads);
368 apply_forces_grp(&pull->group[pcrd.params.group[3]], md, forces.force23, 1, f, pull->nthreads);
370 if (pcrd.params.ngroup >= 6)
372 apply_forces_grp(&pull->group[pcrd.params.group[4]], md, forces.force45, -1, f, pull->nthreads);
373 apply_forces_grp(&pull->group[pcrd.params.group[5]], md, forces.force45, 1, f, pull->nthreads);
378 real max_pull_distance2(const pull_coord_work_t* pcrd, const t_pbc* pbc)
380 /* Note that this maximum distance calculation is more complex than
381 * most other cases in GROMACS, since here we have to take care of
382 * distance calculations that don't involve all three dimensions.
383 * For example, we can use distances that are larger than the
384 * box X and Y dimensions for a box that is elongated along Z.
387 real max_d2 = GMX_REAL_MAX;
389 if (pull_coordinate_is_directional(&pcrd->params))
391 /* Directional pulling along along direction pcrd->vec.
392 * Calculating the exact maximum distance is complex and bug-prone.
393 * So we take a safe approach by not allowing distances that
394 * are larger than half the distance between unit cell faces
395 * along dimensions involved in pcrd->vec.
397 for (int m = 0; m < DIM; m++)
399 if (m < pbc->ndim_ePBC && pcrd->spatialData.vec[m] != 0)
401 real imageDistance2 = gmx::square(pbc->box[m][m]);
402 for (int d = m + 1; d < DIM; d++)
404 imageDistance2 -= gmx::square(pbc->box[d][m]);
406 max_d2 = std::min(max_d2, imageDistance2);
412 /* Distance pulling along dimensions with pcrd->params.dim[d]==1.
413 * We use half the minimum box vector length of the dimensions involved.
414 * This is correct for all cases, except for corner cases with
415 * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
416 * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
417 * in such corner cases the user could get correct results,
418 * depending on the details of the setup, we avoid further
419 * code complications.
421 for (int m = 0; m < DIM; m++)
423 if (m < pbc->ndim_ePBC && pcrd->params.dim[m] != 0)
425 real imageDistance2 = gmx::square(pbc->box[m][m]);
426 for (int d = 0; d < m; d++)
428 if (pcrd->params.dim[d] != 0)
430 imageDistance2 += gmx::square(pbc->box[m][d]);
433 max_d2 = std::min(max_d2, imageDistance2);
438 return 0.25 * max_d2;
441 /* This function returns the distance based on coordinates xg and xref.
442 * Note that the pull coordinate struct pcrd is not modified.
444 static void low_get_pull_coord_dr(const struct pull_t* pull,
445 const pull_coord_work_t* pcrd,
452 const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
454 /* Only the first group can be an absolute reference, in that case nat=0 */
455 if (pgrp0->params.nat == 0)
457 for (int m = 0; m < DIM; m++)
459 xref[m] = pcrd->params.origin[m];
464 copy_dvec(xref, xrefr);
466 dvec dref = { 0, 0, 0 };
467 if (pcrd->params.eGeom == epullgDIRPBC)
469 for (int m = 0; m < DIM; m++)
471 dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
473 /* Add the reference position, so we use the correct periodic image */
474 dvec_inc(xrefr, dref);
477 pbc_dx_d(pbc, xg, xrefr, dr);
479 bool directional = pull_coordinate_is_directional(&pcrd->params);
481 for (int m = 0; m < DIM; m++)
483 dr[m] *= pcrd->params.dim[m];
484 if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
486 dr2 += dr[m] * dr[m];
489 /* Check if we are close to switching to another periodic image */
490 if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
492 /* Note that technically there is no issue with switching periodic
493 * image, as pbc_dx_d returns the distance to the closest periodic
494 * image. However in all cases where periodic image switches occur,
495 * the pull results will be useless in practice.
498 "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
499 "box size (%f).\n%s",
500 pcrd->params.group[0], pcrd->params.group[1], sqrt(dr2),
501 sqrt(0.98 * 0.98 * max_dist2), pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
504 if (pcrd->params.eGeom == epullgDIRPBC)
510 /* This function returns the distance based on the contents of the pull struct.
511 * pull->coord[coord_ind].dr, and potentially vector, are updated.
513 static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
515 pull_coord_work_t* pcrd = &pull->coord[coord_ind];
516 PullCoordSpatialData& spatialData = pcrd->spatialData;
519 /* With AWH pulling we allow for periodic pulling with geometry=direction.
520 * TODO: Store a periodicity flag instead of checking for external pull provider.
522 if (pcrd->params.eGeom == epullgDIRPBC
523 || (pcrd->params.eGeom == epullgDIR && pcrd->params.eType == epullEXTERNAL))
529 md2 = static_cast<double>(max_pull_distance2(pcrd, pbc));
532 if (pcrd->params.eGeom == epullgDIRRELATIVE)
534 /* We need to determine the pull vector */
538 const pull_group_work_t& pgrp2 = pull->group[pcrd->params.group[2]];
539 const pull_group_work_t& pgrp3 = pull->group[pcrd->params.group[3]];
541 pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
543 for (m = 0; m < DIM; m++)
545 vec[m] *= pcrd->params.dim[m];
547 spatialData.vec_len = dnorm(vec);
548 for (m = 0; m < DIM; m++)
550 spatialData.vec[m] = vec[m] / spatialData.vec_len;
554 fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
555 coord_ind, vec[XX], vec[YY], vec[ZZ], spatialData.vec[XX], spatialData.vec[YY],
556 spatialData.vec[ZZ]);
560 pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
561 pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
563 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp1->x,
564 pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, md2,
567 if (pcrd->params.ngroup >= 4)
569 pull_group_work_t *pgrp2, *pgrp3;
570 pgrp2 = &pull->group[pcrd->params.group[2]];
571 pgrp3 = &pull->group[pcrd->params.group[3]];
573 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, md2, spatialData.dr23);
575 if (pcrd->params.ngroup >= 6)
577 pull_group_work_t *pgrp4, *pgrp5;
578 pgrp4 = &pull->group[pcrd->params.group[4]];
579 pgrp5 = &pull->group[pcrd->params.group[5]];
581 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, md2, spatialData.dr45);
585 /* Modify x so that it is periodic in [-pi, pi)
586 * It is assumed that x is in [-3pi, 3pi) so that x
587 * needs to be shifted by at most one period.
589 static void make_periodic_2pi(double* x)
601 /* This function should always be used to modify pcrd->value_ref */
602 static void low_set_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double value_ref)
604 GMX_ASSERT(pcrd->params.eType != epullEXTERNAL,
605 "The pull coord reference value should not be used with type external-potential");
607 if (pcrd->params.eGeom == epullgDIST)
612 "Pull reference distance for coordinate %d (%f) needs to be non-negative",
613 coord_ind + 1, value_ref);
616 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgANGLEAXIS)
618 if (value_ref < 0 || value_ref > M_PI)
621 "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
622 "interval [0,180] deg",
624 value_ref * pull_conversion_factor_internal2userinput(&pcrd->params));
627 else if (pcrd->params.eGeom == epullgDIHEDRAL)
629 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
630 make_periodic_2pi(&value_ref);
633 pcrd->value_ref = value_ref;
636 static void update_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double t)
638 /* With zero rate the reference value is set initially and doesn't change */
639 if (pcrd->params.rate != 0)
641 double value_ref = (pcrd->params.init + pcrd->params.rate * t)
642 * pull_conversion_factor_userinput2internal(&pcrd->params);
643 low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
647 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
648 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
651 dvec dr32; /* store instead of dr23? */
653 dsvmul(-1, spatialData->dr23, dr32);
654 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
655 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
656 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
658 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
659 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
660 * Note 2: the angle between the plane normal vectors equals pi only when
661 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
662 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
664 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
668 /* Calculates pull->coord[coord_ind].value.
669 * This function also updates pull->coord[coord_ind].dr.
671 static void get_pull_coord_distance(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
673 pull_coord_work_t* pcrd;
676 pcrd = &pull->coord[coord_ind];
678 get_pull_coord_dr(pull, coord_ind, pbc);
680 PullCoordSpatialData& spatialData = pcrd->spatialData;
682 switch (pcrd->params.eGeom)
685 /* Pull along the vector between the com's */
686 spatialData.value = dnorm(spatialData.dr01);
690 case epullgDIRRELATIVE:
693 spatialData.value = 0;
694 for (m = 0; m < DIM; m++)
696 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
700 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
702 case epullgDIHEDRAL: spatialData.value = get_dihedral_angle_coord(&spatialData); break;
703 case epullgANGLEAXIS:
704 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
706 default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
710 /* Returns the deviation from the reference value.
711 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
713 static double get_pull_coord_deviation(struct pull_t* pull, int coord_ind, const t_pbc* pbc, double t)
715 pull_coord_work_t* pcrd;
718 pcrd = &pull->coord[coord_ind];
720 /* Update the reference value before computing the distance,
721 * since it is used in the distance computation with periodic pulling.
723 update_pull_coord_reference_value(pcrd, coord_ind, t);
725 get_pull_coord_distance(pull, coord_ind, pbc);
727 get_pull_coord_distance(pull, coord_ind, pbc);
729 /* Determine the deviation */
730 dev = pcrd->spatialData.value - pcrd->value_ref;
732 /* Check that values are allowed */
733 if (pcrd->params.eGeom == epullgDIST && pcrd->spatialData.value == 0)
735 /* With no vector we can not determine the direction for the force,
736 * so we set the force to zero.
740 else if (pcrd->params.eGeom == epullgDIHEDRAL)
742 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
743 Thus, the unwrapped deviation is here in (-2pi, 2pi].
744 After making it periodic, the deviation will be in [-pi, pi). */
745 make_periodic_2pi(&dev);
751 double get_pull_coord_value(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
753 get_pull_coord_distance(pull, coord_ind, pbc);
755 return pull->coord[coord_ind].spatialData.value;
758 void clear_pull_forces(pull_t* pull)
760 /* Zeroing the forces is only required for constraint pulling.
761 * It can happen that multiple constraint steps need to be applied
762 * and therefore the constraint forces need to be accumulated.
764 for (pull_coord_work_t& coord : pull->coord)
766 coord.scalarForce = 0;
770 /* Apply constraint using SHAKE */
772 do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
775 dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
776 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
777 dvec* rnew; /* current 'new' positions of the groups */
778 double* dr_tot; /* the total update of the coords */
781 double lambda, rm, invdt = 0;
782 gmx_bool bConverged_all, bConverged = FALSE;
783 int niter = 0, ii, j, m, max_iter = 100;
787 snew(r_ij, pull->coord.size());
788 snew(dr_tot, pull->coord.size());
790 snew(rnew, pull->group.size());
792 /* copy the current unconstrained positions for use in iterations. We
793 iterate until rinew[i] and rjnew[j] obey the constraints. Then
794 rinew - pull.x_unc[i] is the correction dr to group i */
795 for (size_t g = 0; g < pull->group.size(); g++)
797 copy_dvec(pull->group[g].xp, rnew[g]);
800 /* Determine the constraint directions from the old positions */
801 for (size_t c = 0; c < pull->coord.size(); c++)
803 pull_coord_work_t* pcrd;
805 pcrd = &pull->coord[c];
807 if (pcrd->params.eType != epullCONSTRAINT)
812 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
813 * We don't modify dr and value anymore, so these values are also used
816 get_pull_coord_distance(pull, c, pbc);
818 const PullCoordSpatialData& spatialData = pcrd->spatialData;
821 fprintf(debug, "Pull coord %zu dr %f %f %f\n", c, spatialData.dr01[XX],
822 spatialData.dr01[YY], spatialData.dr01[ZZ]);
825 if (pcrd->params.eGeom == epullgDIR || pcrd->params.eGeom == epullgDIRPBC)
827 /* Select the component along vec */
829 for (m = 0; m < DIM; m++)
831 a += spatialData.vec[m] * spatialData.dr01[m];
833 for (m = 0; m < DIM; m++)
835 r_ij[c][m] = a * spatialData.vec[m];
840 copy_dvec(spatialData.dr01, r_ij[c]);
843 if (dnorm2(r_ij[c]) == 0)
846 "Distance for pull coordinate %zu is zero with constraint pulling, which is "
852 bConverged_all = FALSE;
853 while (!bConverged_all && niter < max_iter)
855 bConverged_all = TRUE;
857 /* loop over all constraints */
858 for (size_t c = 0; c < pull->coord.size(); c++)
860 pull_coord_work_t* pcrd;
861 pull_group_work_t *pgrp0, *pgrp1;
864 pcrd = &pull->coord[c];
866 if (pcrd->params.eType != epullCONSTRAINT)
871 update_pull_coord_reference_value(pcrd, c, t);
873 pgrp0 = &pull->group[pcrd->params.group[0]];
874 pgrp1 = &pull->group[pcrd->params.group[1]];
876 /* Get the current difference vector */
877 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->params.group[1]],
878 rnew[pcrd->params.group[0]], -1, unc_ij);
882 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
885 rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
887 switch (pcrd->params.eGeom)
890 if (pcrd->value_ref <= 0)
894 "The pull constraint reference distance for group %zu is <= 0 (%f)",
899 double q, c_a, c_b, c_c;
901 c_a = diprod(r_ij[c], r_ij[c]);
902 c_b = diprod(unc_ij, r_ij[c]) * 2;
903 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
907 q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
912 q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
918 fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b,
923 /* The position corrections dr due to the constraints */
924 dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
925 dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
926 dr_tot[c] += -lambda * dnorm(r_ij[c]);
931 /* A 1-dimensional constraint along a vector */
933 for (m = 0; m < DIM; m++)
935 vec[m] = pcrd->spatialData.vec[m];
936 a += unc_ij[m] * vec[m];
938 /* Select only the component along the vector */
939 dsvmul(a, vec, unc_ij);
940 lambda = a - pcrd->value_ref;
943 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
946 /* The position corrections dr due to the constraints */
947 dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
948 dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
949 dr_tot[c] += -lambda;
951 default: gmx_incons("Invalid enumeration value for eGeom");
959 g0 = pcrd->params.group[0];
960 g1 = pcrd->params.group[1];
961 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
962 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
963 fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", rnew[g0][0],
964 rnew[g0][1], rnew[g0][2], rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
965 fprintf(debug, "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", "", "", "", "", "",
966 "", pcrd->value_ref);
967 fprintf(debug, "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", dr0[0],
968 dr0[1], dr0[2], dr1[0], dr1[1], dr1[2], dnorm(tmp3));
971 /* Update the COMs with dr */
972 dvec_inc(rnew[pcrd->params.group[1]], dr1);
973 dvec_inc(rnew[pcrd->params.group[0]], dr0);
976 /* Check if all constraints are fullfilled now */
977 for (pull_coord_work_t& coord : pull->coord)
979 if (coord.params.eType != epullCONSTRAINT)
984 low_get_pull_coord_dr(pull, &coord, pbc, rnew[coord.params.group[1]],
985 rnew[coord.params.group[0]], -1, unc_ij);
987 switch (coord.params.eGeom)
990 bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
995 for (m = 0; m < DIM; m++)
997 vec[m] = coord.spatialData.vec[m];
999 inpr = diprod(unc_ij, vec);
1000 dsvmul(inpr, vec, unc_ij);
1001 bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1010 "Pull constraint not converged: "
1012 "d_ref = %f, current d = %f\n",
1013 coord.params.group[0], coord.params.group[1], coord.value_ref, dnorm(unc_ij));
1016 bConverged_all = FALSE;
1021 /* if after all constraints are dealt with and bConverged is still TRUE
1022 we're finished, if not we do another iteration */
1024 if (niter > max_iter)
1026 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1029 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1036 /* update atoms in the groups */
1037 for (size_t g = 0; g < pull->group.size(); g++)
1039 const pull_group_work_t* pgrp;
1042 pgrp = &pull->group[g];
1044 /* get the final constraint displacement dr for group g */
1045 dvec_sub(rnew[g], pgrp->xp, dr);
1047 if (dnorm2(dr) == 0)
1049 /* No displacement, no update necessary */
1053 /* update the atom positions */
1054 auto localAtomIndices = pgrp->atomSet.localIndex();
1056 for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1058 ii = localAtomIndices[j];
1059 if (!pgrp->localWeights.empty())
1061 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1063 for (m = 0; m < DIM; m++)
1069 for (m = 0; m < DIM; m++)
1071 v[ii][m] += invdt * tmp[m];
1077 /* calculate the constraint forces, used for output and virial only */
1078 for (size_t c = 0; c < pull->coord.size(); c++)
1080 pull_coord_work_t* pcrd;
1082 pcrd = &pull->coord[c];
1084 if (pcrd->params.eType != epullCONSTRAINT)
1089 /* Accumulate the forces, in case we have multiple constraint steps */
1092 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1094 pcrd->scalarForce += force;
1096 if (vir != nullptr && pcrd->params.eGeom != epullgDIRPBC && bMaster)
1100 /* Add the pull contribution to the virial */
1101 /* We have already checked above that r_ij[c] != 0 */
1102 f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1104 for (j = 0; j < DIM; j++)
1106 for (m = 0; m < DIM; m++)
1108 vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1114 /* finished! I hope. Give back some memory */
1120 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1122 for (int j = 0; j < DIM; j++)
1124 for (int m = 0; m < DIM; m++)
1126 vir[j][m] -= 0.5 * f[j] * dr[m];
1131 /* Adds the pull contribution to the virial */
1132 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1134 if (vir != nullptr && pcrd.params.eGeom != epullgDIRPBC)
1136 /* Add the pull contribution for each distance vector to the virial. */
1137 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1138 if (pcrd.params.ngroup >= 4)
1140 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1142 if (pcrd.params.ngroup >= 6)
1144 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1149 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1157 k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1158 dkdl = pcrd->params.kB - pcrd->params.k;
1160 switch (pcrd->params.eType)
1163 case epullFLATBOTTOM:
1164 case epullFLATBOTTOMHIGH:
1165 /* The only difference between an umbrella and a flat-bottom
1166 * potential is that a flat-bottom is zero above or below
1167 the reference value.
1169 if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0)
1170 || (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
1175 pcrd->scalarForce = -k * dev;
1176 *V += 0.5 * k * gmx::square(dev);
1177 *dVdl += 0.5 * dkdl * gmx::square(dev);
1180 pcrd->scalarForce = -k;
1181 *V += k * pcrd->spatialData.value;
1182 *dVdl += dkdl * pcrd->spatialData.value;
1186 "the scalar pull force should not be calculated internally for pull type "
1188 default: gmx_incons("Unsupported pull type in do_pull_pot");
1192 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1194 const t_pull_coord& params = pcrd.params;
1195 const PullCoordSpatialData& spatialData = pcrd.spatialData;
1197 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1198 PullCoordVectorForces forces;
1200 if (params.eGeom == epullgDIST)
1202 double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1203 for (int m = 0; m < DIM; m++)
1205 forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1208 else if (params.eGeom == epullgANGLE)
1211 double cos_theta, cos_theta2;
1213 cos_theta = cos(spatialData.value);
1214 cos_theta2 = gmx::square(cos_theta);
1216 /* The force at theta = 0, pi is undefined so we should not apply any force.
1217 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1218 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1219 * have to check for this before dividing by their norm below.
1223 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1224 * and another vector parallel to dr23:
1225 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1226 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1228 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1229 double b = a * cos_theta;
1230 double invdr01 = 1. / dnorm(spatialData.dr01);
1231 double invdr23 = 1. / dnorm(spatialData.dr23);
1232 dvec normalized_dr01, normalized_dr23;
1233 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1234 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1236 for (int m = 0; m < DIM; m++)
1238 /* Here, f_scal is -dV/dtheta */
1240 pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1242 pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1247 /* No forces to apply for ill-defined cases*/
1248 clear_dvec(forces.force01);
1249 clear_dvec(forces.force23);
1252 else if (params.eGeom == epullgANGLEAXIS)
1254 double cos_theta, cos_theta2;
1256 /* The angle-axis force is exactly the same as the angle force (above) except that in
1257 this case the second vector (dr23) is replaced by the pull vector. */
1258 cos_theta = cos(spatialData.value);
1259 cos_theta2 = gmx::square(cos_theta);
1265 dvec normalized_dr01;
1267 invdr01 = 1. / dnorm(spatialData.dr01);
1268 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1269 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1272 for (int m = 0; m < DIM; m++)
1275 pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1280 clear_dvec(forces.force01);
1283 else if (params.eGeom == epullgDIHEDRAL)
1285 double m2, n2, tol, sqrdist_32;
1287 /* Note: there is a small difference here compared to the
1288 dihedral force calculations in the bondeds (ref: Bekker 1994).
1289 There rij = ri - rj, while here dr01 = r1 - r0.
1290 However, all distance vectors occur in form of cross or inner products
1291 so that two signs cancel and we end up with the same expressions.
1292 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1294 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1295 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1296 dsvmul(-1, spatialData.dr23, dr32);
1297 sqrdist_32 = diprod(dr32, dr32);
1298 tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1299 if ((m2 > tol) && (n2 > tol))
1301 double a_01, a_23_01, a_23_45, a_45;
1302 double inv_dist_32, inv_sqrdist_32, dist_32;
1304 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1305 inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1306 dist_32 = sqrdist_32 * inv_dist_32;
1308 /* Forces on groups 0, 1 */
1309 a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
1310 dsvmul(-a_01, spatialData.planevec_m,
1311 forces.force01); /* added sign to get force on group 1, not 0 */
1313 /* Forces on groups 4, 5 */
1314 a_45 = -pcrd.scalarForce * dist_32 / n2;
1315 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1317 /* Force on groups 2, 3 (defining the axis) */
1318 a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1319 a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1320 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1321 dsvmul(a_23_45, forces.force45, v);
1322 dvec_sub(u, v, forces.force23); /* force on group 3 */
1326 /* No force to apply for ill-defined cases */
1327 clear_dvec(forces.force01);
1328 clear_dvec(forces.force23);
1329 clear_dvec(forces.force45);
1334 for (int m = 0; m < DIM; m++)
1336 forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1344 /* We use a global mutex for locking access to the pull data structure
1345 * during registration of external pull potential providers.
1346 * We could use a different, local mutex for each pull object, but the overhead
1347 * is extremely small here and registration is only done during initialization.
1349 static gmx::Mutex registrationMutex;
1351 using Lock = gmx::lock_guard<gmx::Mutex>;
1353 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1355 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1356 GMX_RELEASE_ASSERT(provider != nullptr,
1357 "register_external_pull_potential called with NULL as provider name");
1359 if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1362 "Module '%s' attempted to register an external potential for pull coordinate %d "
1363 "which is out of the pull coordinate range %d - %zu\n",
1364 provider, coord_index + 1, 1, pull->coord.size());
1367 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1369 if (pcrd->params.eType != epullEXTERNAL)
1373 "Module '%s' attempted to register an external potential for pull coordinate %d "
1374 "which of type '%s', whereas external potentials are only supported with type '%s'",
1375 provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
1378 GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr,
1379 "The external potential provider string for a pull coordinate is NULL");
1381 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider) != 0)
1384 "Module '%s' attempted to register an external potential for pull coordinate %d "
1385 "which expects the external potential to be provided by a module named '%s'",
1386 provider, coord_index + 1, pcrd->params.externalPotentialProvider);
1389 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1390 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1391 * pull->numUnregisteredExternalPotentials.
1393 Lock registrationLock(registrationMutex);
1395 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1398 "Module '%s' attempted to register an external potential for pull coordinate %d "
1400 provider, coord_index + 1);
1403 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1404 pull->numUnregisteredExternalPotentials--;
1406 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1407 "Negative unregistered potentials, the pull code is inconsistent");
1411 static void check_external_potential_registration(const struct pull_t* pull)
1413 if (pull->numUnregisteredExternalPotentials > 0)
1416 for (c = 0; c < pull->coord.size(); c++)
1418 if (pull->coord[c].params.eType == epullEXTERNAL
1419 && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1425 GMX_RELEASE_ASSERT(c < pull->coord.size(),
1426 "Internal inconsistency in the pull potential provider counting");
1429 "No external provider for external pull potentials have been provided for %d "
1430 "pull coordinates. The first coordinate without provider is number %zu, which "
1431 "expects a module named '%s' to provide the external potential.",
1432 pull->numUnregisteredExternalPotentials, c + 1,
1433 pull->coord[c].params.externalPotentialProvider);
1438 /* Pull takes care of adding the forces of the external potential.
1439 * The external potential module has to make sure that the corresponding
1440 * potential energy is added either to the pull term or to a term
1441 * specific to the external module.
1443 void apply_external_pull_coord_force(struct pull_t* pull,
1446 const t_mdatoms* mdatoms,
1447 gmx::ForceWithVirial* forceWithVirial)
1449 pull_coord_work_t* pcrd;
1451 GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1452 "apply_external_pull_coord_force called with coord_index out of range");
1454 if (pull->comm.bParticipate)
1456 pcrd = &pull->coord[coord_index];
1459 pcrd->params.eType == epullEXTERNAL,
1460 "The pull force can only be set externally on pull coordinates of external type");
1462 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1463 "apply_external_pull_coord_force called for an unregistered pull coordinate");
1466 pcrd->scalarForce = coord_force;
1468 /* Calculate the forces on the pull groups */
1469 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1471 /* Add the forces for this coordinate to the total virial and force */
1472 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1474 matrix virial = { { 0 } };
1475 add_virial_coord(virial, *pcrd, pullCoordForces);
1476 forceWithVirial->addVirialContribution(virial);
1479 apply_forces_coord(pull, coord_index, pullCoordForces, mdatoms,
1480 as_rvec_array(forceWithVirial->force_.data()));
1483 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1486 /* Calculate the pull potential and scalar force for a pull coordinate.
1487 * Returns the vector forces for the pull coordinate.
1489 static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull,
1498 pull_coord_work_t& pcrd = pull->coord[coord_ind];
1500 assert(pcrd.params.eType != epullCONSTRAINT);
1502 double dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
1504 calc_pull_coord_scalar_force_and_potential(&pcrd, dev, lambda, V, dVdl);
1506 PullCoordVectorForces pullCoordForces = calculateVectorForces(pcrd);
1508 add_virial_coord(vir, pcrd, pullCoordForces);
1510 return pullCoordForces;
1513 real pull_potential(struct pull_t* pull,
1514 const t_mdatoms* md,
1516 const t_commrec* cr,
1520 gmx::ForceWithVirial* force,
1525 assert(pull != nullptr);
1527 /* Ideally we should check external potential registration only during
1528 * the initialization phase, but that requires another function call
1529 * that should be done exactly in the right order. So we check here.
1531 check_external_potential_registration(pull);
1533 if (pull->comm.bParticipate)
1537 pull_calc_coms(cr, pull, md, pbc, t, x, nullptr);
1539 rvec* f = as_rvec_array(force->force_.data());
1540 matrix virial = { { 0 } };
1541 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1542 for (size_t c = 0; c < pull->coord.size(); c++)
1544 pull_coord_work_t* pcrd;
1545 pcrd = &pull->coord[c];
1547 /* For external potential the force is assumed to be given by an external module by a
1548 call to apply_pull_coord_external_force */
1549 if (pcrd->params.eType == epullCONSTRAINT || pcrd->params.eType == epullEXTERNAL)
1554 PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1555 pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1557 /* Distribute the force over the atoms in the pulled groups */
1558 apply_forces_coord(pull, c, pullCoordForces, md, f);
1563 force->addVirialContribution(virial);
1568 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1569 "Too few or too many external pull potentials have been applied the previous step");
1570 /* All external pull potentials still need to be applied */
1571 pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1573 return (MASTER(cr) ? V : 0.0);
1576 void pull_constraint(struct pull_t* pull,
1577 const t_mdatoms* md,
1579 const t_commrec* cr,
1587 assert(pull != nullptr);
1589 if (pull->comm.bParticipate)
1591 pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1593 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1597 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1601 gmx_bool bMustParticipate;
1607 /* We always make the master node participate, such that it can do i/o,
1608 * add the virial and to simplify MC type extensions people might have.
1610 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1612 for (pull_group_work_t& group : pull->group)
1614 if (!group.globalWeights.empty())
1616 group.localWeights.resize(group.atomSet.numAtomsLocal());
1617 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1619 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1623 GMX_ASSERT(bMustParticipate || dd != nullptr,
1624 "Either all ranks (including this rank) participate, or we use DD and need to "
1625 "have access to dd here");
1627 /* We should participate if we have pull or pbc atoms */
1628 if (!bMustParticipate
1629 && (group.atomSet.numAtomsLocal() > 0
1630 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1632 bMustParticipate = TRUE;
1636 if (!comm->bParticipateAll)
1638 /* Keep currently not required ranks in the communicator
1639 * if they needed to participate up to 20 decompositions ago.
1640 * This avoids frequent rebuilds due to atoms jumping back and forth.
1642 const int64_t history_count = 20;
1643 gmx_bool bWillParticipate;
1646 /* Increase the decomposition counter for the current call */
1647 comm->setup_count++;
1649 if (bMustParticipate)
1651 comm->must_count = comm->setup_count;
1656 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1658 if (debug && dd != nullptr)
1660 fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n", dd->rank,
1661 gmx::boolToString(bMustParticipate), gmx::boolToString(bWillParticipate));
1664 if (bWillParticipate)
1666 /* Count the number of ranks that we want to have participating */
1668 /* Count the number of ranks that need to be added */
1669 count[1] = comm->bParticipate ? 0 : 1;
1677 /* The cost of this global operation will be less that the cost
1678 * of the extra MPI_Comm_split calls that we can avoid.
1680 gmx_sumi(2, count, cr);
1682 /* If we are missing ranks or if we have 20% more ranks than needed
1683 * we make a new sub-communicator.
1685 if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1689 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1693 if (comm->mpi_comm_com != MPI_COMM_NULL)
1695 MPI_Comm_free(&comm->mpi_comm_com);
1698 /* This might be an extremely expensive operation, so we try
1699 * to avoid this splitting as much as possible.
1701 assert(dd != nullptr);
1702 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1705 comm->bParticipate = bWillParticipate;
1706 comm->nparticipate = count[0];
1708 /* When we use the previous COM for PBC, we need to broadcast
1709 * the previous COM to ranks that have joined the communicator.
1711 for (pull_group_work_t& group : pull->group)
1713 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1715 GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1716 "The master rank has to participate, as it should pass an up to "
1718 "to bcast here as well as to e.g. checkpointing");
1720 gmx_bcast(sizeof(dvec), group.x_prev_step, cr);
1726 /* Since the PBC of atoms might have changed, we need to update the PBC */
1727 pull->bSetPBCatoms = TRUE;
1730 static void init_pull_group_index(FILE* fplog,
1731 const t_commrec* cr,
1733 pull_group_work_t* pg,
1734 gmx_bool bConstraint,
1735 const ivec pulldim_con,
1736 const gmx_mtop_t* mtop,
1737 const t_inputrec* ir,
1740 /* With EM and BD there are no masses in the integrator.
1741 * But we still want to have the correct mass-weighted COMs.
1742 * So we store the real masses in the weights.
1744 const bool setWeights = (pg->params.nweight > 0 || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
1746 /* In parallel, store we need to extract localWeights from weights at DD time */
1747 std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1749 const SimulationGroups& groups = mtop->groups;
1751 /* Count frozen dimensions and (weighted) mass */
1757 for (int i = 0; i < pg->params.nat; i++)
1759 int ii = pg->params.ind[i];
1760 if (bConstraint && ir->opts.nFreeze)
1762 for (int d = 0; d < DIM; d++)
1764 if (pulldim_con[d] == 1
1765 && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1771 const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1773 if (ir->efep == efepNO)
1779 m = (1 - lambda) * atom.m + lambda * atom.mB;
1782 if (pg->params.nweight > 0)
1784 w = pg->params.weight[i];
1790 if (EI_ENERGY_MINIMIZATION(ir->eI))
1792 /* Move the mass to the weight */
1796 else if (ir->eI == eiBD)
1799 if (ir->bd_fric != 0.0F)
1801 mbd = ir->bd_fric * ir->delta_t;
1805 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1807 mbd = ir->delta_t / ir->opts.tau_t[0];
1812 / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1820 weights.push_back(w);
1824 wwmass += m * w * w;
1829 /* We can have single atom groups with zero mass with potential pulling
1830 * without cosine weighting.
1832 if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1834 /* With one atom the mass doesn't matter */
1839 gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1840 pg->params.weight ? " weighted" : "", g);
1845 fprintf(fplog, "Pull group %d: %5d atoms, mass %9.3f", g, pg->params.nat, tmass);
1846 if (pg->params.weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1848 fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1850 if (pg->epgrppbc == epgrppbcCOS)
1852 fprintf(fplog, ", cosine weighting will be used");
1854 fprintf(fplog, "\n");
1859 /* A value != 0 signals not frozen, it is updated later */
1865 for (int d = 0; d < DIM; d++)
1867 ndim += pulldim_con[d] * pg->params.nat;
1869 if (fplog && nfrozen > 0 && nfrozen < ndim)
1872 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1873 " that are subject to pulling are frozen.\n"
1874 " For constraint pulling the whole group will be frozen.\n\n",
1882 struct pull_t* init_pull(FILE* fplog,
1883 const pull_params_t* pull_params,
1884 const t_inputrec* ir,
1885 const gmx_mtop_t* mtop,
1886 const t_commrec* cr,
1887 gmx::LocalAtomSetManager* atomSets,
1890 struct pull_t* pull;
1895 /* Copy the pull parameters */
1896 pull->params = *pull_params;
1897 /* Avoid pointer copies */
1898 pull->params.group = nullptr;
1899 pull->params.coord = nullptr;
1901 for (int i = 0; i < pull_params->ngroup; ++i)
1903 pull->group.emplace_back(pull_params->group[i],
1904 atomSets->add({ pull_params->group[i].ind,
1905 pull_params->group[i].ind + pull_params->group[i].nat }),
1906 pull_params->bSetPbcRefToPrevStepCOM);
1909 if (cr != nullptr && DOMAINDECOMP(cr))
1911 /* Set up the global to local atom mapping for PBC atoms */
1912 for (pull_group_work_t& group : pull->group)
1914 if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
1916 /* pbcAtomSet consists of a single atom */
1917 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
1918 atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
1923 pull->bPotential = FALSE;
1924 pull->bConstraint = FALSE;
1925 pull->bCylinder = FALSE;
1926 pull->bAngle = FALSE;
1927 pull->bXOutAverage = pull_params->bXOutAverage;
1928 pull->bFOutAverage = pull_params->bFOutAverage;
1930 GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0,
1931 "pull group 0 is an absolute reference group and should not contain atoms");
1933 pull->numCoordinatesWithExternalPotential = 0;
1935 for (int c = 0; c < pull->params.ncoord; c++)
1937 /* Construct a pull coordinate, copying all coordinate parameters */
1938 pull->coord.emplace_back(pull_params->coord[c]);
1940 pull_coord_work_t* pcrd = &pull->coord.back();
1942 switch (pcrd->params.eGeom)
1945 case epullgDIRRELATIVE: /* Direction vector is determined at each step */
1947 case epullgDIHEDRAL: break;
1951 case epullgANGLEAXIS:
1952 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
1955 /* We allow reading of newer tpx files with new pull geometries,
1956 * but with the same tpx format, with old code. A new geometry
1957 * only adds a new enum value, which will be out of range for
1958 * old code. The first place we need to generate an error is
1959 * here, since the pull code can't handle this.
1960 * The enum can be used before arriving here only for printing
1961 * the string corresponding to the geometry, which will result
1962 * in printing "UNDEFINED".
1965 "Pull geometry not supported for pull coordinate %d. The geometry enum "
1966 "%d in the input is larger than that supported by the code (up to %d). "
1967 "You are probably reading a tpr file generated with a newer version of "
1968 "Gromacs with an binary from an older version of Gromacs.",
1969 c + 1, pcrd->params.eGeom, epullgNR - 1);
1972 if (pcrd->params.eType == epullCONSTRAINT)
1974 /* Check restrictions of the constraint pull code */
1975 if (pcrd->params.eGeom == epullgCYL || pcrd->params.eGeom == epullgDIRRELATIVE
1976 || pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
1977 || pcrd->params.eGeom == epullgANGLEAXIS)
1980 "Pulling of type %s can not be combined with geometry %s. Consider using "
1982 epull_names[pcrd->params.eType], epullg_names[pcrd->params.eGeom],
1983 epull_names[epullUMBRELLA]);
1986 pull->bConstraint = TRUE;
1990 pull->bPotential = TRUE;
1993 if (pcrd->params.eGeom == epullgCYL)
1995 pull->bCylinder = TRUE;
1997 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
1998 || pcrd->params.eGeom == epullgANGLEAXIS)
2000 pull->bAngle = TRUE;
2003 /* We only need to calculate the plain COM of a group
2004 * when it is not only used as a cylinder group.
2005 * Also the absolute reference group 0 needs no COM computation.
2007 for (int i = 0; i < pcrd->params.ngroup; i++)
2009 int groupIndex = pcrd->params.group[i];
2010 if (groupIndex > 0 && !(pcrd->params.eGeom == epullgCYL && i == 0))
2012 pull->group[groupIndex].needToCalcCom = true;
2016 /* With non-zero rate the reference value is set at every step */
2017 if (pcrd->params.rate == 0)
2019 /* Initialize the constant reference value */
2020 if (pcrd->params.eType != epullEXTERNAL)
2022 low_set_pull_coord_reference_value(
2024 pcrd->params.init * pull_conversion_factor_userinput2internal(&pcrd->params));
2028 /* With an external pull potential, the reference value loses
2029 * it's meaning and should not be used. Setting it to zero
2030 * makes any terms dependent on it disappear.
2031 * The only issue this causes is that with cylinder pulling
2032 * no atoms of the cylinder group within the cylinder radius
2033 * should be more than half a box length away from the COM of
2034 * the pull group along the axial direction.
2036 pcrd->value_ref = 0.0;
2040 if (pcrd->params.eType == epullEXTERNAL)
2043 pcrd->params.rate == 0,
2044 "With an external potential, a pull coordinate should have rate = 0");
2046 /* This external potential needs to be registered later */
2047 pull->numCoordinatesWithExternalPotential++;
2049 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2052 pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
2053 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2055 pull->ePBC = ir->ePBC;
2058 case epbcNONE: pull->npbcdim = 0; break;
2059 case epbcXY: pull->npbcdim = 2; break;
2060 default: pull->npbcdim = 3; break;
2065 gmx_bool bAbs, bCos;
2068 for (const pull_coord_work_t& coord : pull->coord)
2070 if (pull->group[coord.params.group[0]].params.nat == 0
2071 || pull->group[coord.params.group[1]].params.nat == 0)
2077 fprintf(fplog, "\n");
2078 if (pull->bPotential)
2080 fprintf(fplog, "Will apply potential COM pulling\n");
2082 if (pull->bConstraint)
2084 fprintf(fplog, "Will apply constraint COM pulling\n");
2086 // Don't include the reference group 0 in output, so we report ngroup-1
2087 int numRealGroups = pull->group.size() - 1;
2088 GMX_RELEASE_ASSERT(numRealGroups > 0,
2089 "The reference absolute position pull group should always be present");
2090 fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n", pull->coord.size(),
2091 pull->coord.size() == 1 ? "" : "s", numRealGroups, numRealGroups == 1 ? "" : "s");
2094 fprintf(fplog, "with an absolute reference\n");
2097 // Don't include the reference group 0 in loop
2098 for (size_t g = 1; g < pull->group.size(); g++)
2100 if (pull->group[g].params.nat > 1 && pull->group[g].params.pbcatom < 0)
2102 /* We are using cosine weighting */
2103 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2109 please_cite(fplog, "Engin2010");
2113 pull->bRefAt = FALSE;
2115 for (size_t g = 0; g < pull->group.size(); g++)
2117 pull_group_work_t* pgrp;
2119 pgrp = &pull->group[g];
2120 if (pgrp->params.nat > 0)
2122 /* There is an issue when a group is used in multiple coordinates
2123 * and constraints are applied in different dimensions with atoms
2124 * frozen in some, but not all dimensions.
2125 * Since there is only one mass array per group, we can't have
2126 * frozen/non-frozen atoms for different coords at the same time.
2127 * But since this is a very exotic case, we don't check for this.
2128 * A warning is printed in init_pull_group_index.
2131 gmx_bool bConstraint;
2132 ivec pulldim, pulldim_con;
2134 /* Loop over all pull coordinates to see along which dimensions
2135 * this group is pulled and if it is involved in constraints.
2137 bConstraint = FALSE;
2138 clear_ivec(pulldim);
2139 clear_ivec(pulldim_con);
2140 for (const pull_coord_work_t& coord : pull->coord)
2142 gmx_bool bGroupUsed = FALSE;
2143 for (int gi = 0; gi < coord.params.ngroup; gi++)
2145 if (coord.params.group[gi] == static_cast<int>(g))
2153 for (int m = 0; m < DIM; m++)
2155 if (coord.params.dim[m] == 1)
2159 if (coord.params.eType == epullCONSTRAINT)
2169 /* Determine if we need to take PBC into account for calculating
2170 * the COM's of the pull groups.
2172 switch (pgrp->epgrppbc)
2174 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2176 if (pgrp->params.weight != nullptr)
2179 "Pull groups can not have relative weights and cosine weighting "
2182 for (int m = 0; m < DIM; m++)
2184 if (m < pull->npbcdim && pulldim[m] == 1)
2186 if (pull->cosdim >= 0 && pull->cosdim != m)
2189 "Can only use cosine weighting with pulling in one "
2190 "dimension (use mdp option pull-coord?-dim)");
2196 case epgrppbcNONE: break;
2199 /* Set the indices */
2200 init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2204 /* Absolute reference, set the inverse mass to zero.
2205 * This is only relevant (and used) with constraint pulling.
2212 /* If we use cylinder coordinates, do some initialising for them */
2213 if (pull->bCylinder)
2215 for (const pull_coord_work_t& coord : pull->coord)
2217 if (coord.params.eGeom == epullgCYL)
2219 if (pull->group[coord.params.group[0]].params.nat == 0)
2222 "A cylinder pull group is not supported when using absolute "
2226 const auto& referenceGroup = pull->group[coord.params.group[0]];
2227 pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet,
2228 pull->params.bSetPbcRefToPrevStepCOM);
2232 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2233 pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
2234 pull->comSums.resize(pull->nthreads);
2239 /* Use a sub-communicator when we have more than 32 ranks, but not
2240 * when we have an external pull potential, since then the external
2241 * potential provider expects each rank to have the coordinate.
2243 comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
2244 || pull->numCoordinatesWithExternalPotential > 0
2245 || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2246 /* This sub-commicator is not used with comm->bParticipateAll,
2247 * so we can always initialize it to NULL.
2249 comm->mpi_comm_com = MPI_COMM_NULL;
2250 comm->nparticipate = 0;
2251 comm->isMasterRank = (cr == nullptr || MASTER(cr));
2253 /* No MPI: 1 rank: all ranks pull */
2254 comm->bParticipateAll = TRUE;
2255 comm->isMasterRank = true;
2257 comm->bParticipate = comm->bParticipateAll;
2258 comm->setup_count = 0;
2259 comm->must_count = 0;
2261 if (!comm->bParticipateAll && fplog != nullptr)
2263 fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2266 comm->pbcAtomBuffer.resize(pull->group.size());
2267 comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
2268 if (pull->bCylinder)
2270 comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
2273 /* We still need to initialize the PBC reference coordinates */
2274 pull->bSetPBCatoms = TRUE;
2276 pull->out_x = nullptr;
2277 pull->out_f = nullptr;
2282 static void destroy_pull(struct pull_t* pull)
2285 if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2287 MPI_Comm_free(&pull->comm.mpi_comm_com);
2294 void preparePrevStepPullCom(const t_inputrec* ir,
2296 const t_mdatoms* md,
2298 const t_state* state_global,
2299 const t_commrec* cr,
2300 bool startingFromCheckpoint)
2302 if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
2306 allocStatePrevStepPullCom(state, pull_work);
2307 if (startingFromCheckpoint)
2311 state->pull_com_prev_step = state_global->pull_com_prev_step;
2315 /* Only the master rank has the checkpointed COM from the previous step */
2316 gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(), &state->pull_com_prev_step[0], cr);
2318 setPrevStepPullComFromState(pull_work, state);
2323 set_pbc(&pbc, ir->ePBC, state->box);
2324 initPullComFromPrevStep(cr, pull_work, md, &pbc, state->x.rvec_array());
2325 updatePrevStepPullCom(pull_work, state);
2329 void finish_pull(struct pull_t* pull)
2331 check_external_potential_registration(pull);
2335 gmx_fio_fclose(pull->out_x);
2339 gmx_fio_fclose(pull->out_f);
2345 gmx_bool pull_have_potential(const struct pull_t* pull)
2347 return pull->bPotential;
2350 gmx_bool pull_have_constraint(const struct pull_t* pull)
2352 return pull->bConstraint;
2355 bool pull_have_constraint(const pull_params_t* pullParameters)
2357 if (pullParameters == nullptr)
2361 for (int c = 0; c < pullParameters->ncoord; c++)
2363 if (pullParameters->coord[c].eType == epullCONSTRAINT)