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)
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.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(), 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], pcrd.spatialData.cyl_dev, masses, 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]], masses, 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, masses, f);
359 if (pull->group[pcrd.params.group[0]].params.nat > 0)
361 apply_forces_grp(&pull->group[pcrd.params.group[0]], masses, forces.force01, -1, f,
364 apply_forces_grp(&pull->group[pcrd.params.group[1]], masses, forces.force01, 1, f, pull->nthreads);
366 if (pcrd.params.ngroup >= 4)
368 apply_forces_grp(&pull->group[pcrd.params.group[2]], masses, forces.force23, -1, f,
370 apply_forces_grp(&pull->group[pcrd.params.group[3]], masses, forces.force23, 1, f,
373 if (pcrd.params.ngroup >= 6)
375 apply_forces_grp(&pull->group[pcrd.params.group[4]], masses, forces.force45, -1, f,
377 apply_forces_grp(&pull->group[pcrd.params.group[5]], masses, forces.force45, 1, f,
383 real max_pull_distance2(const pull_coord_work_t* pcrd, const t_pbc* pbc)
385 /* Note that this maximum distance calculation is more complex than
386 * most other cases in GROMACS, since here we have to take care of
387 * distance calculations that don't involve all three dimensions.
388 * For example, we can use distances that are larger than the
389 * box X and Y dimensions for a box that is elongated along Z.
392 real max_d2 = GMX_REAL_MAX;
394 if (pull_coordinate_is_directional(&pcrd->params))
396 /* Directional pulling along along direction pcrd->vec.
397 * Calculating the exact maximum distance is complex and bug-prone.
398 * So we take a safe approach by not allowing distances that
399 * are larger than half the distance between unit cell faces
400 * along dimensions involved in pcrd->vec.
402 for (int m = 0; m < DIM; m++)
404 if (m < pbc->ndim_ePBC && pcrd->spatialData.vec[m] != 0)
406 real imageDistance2 = gmx::square(pbc->box[m][m]);
407 for (int d = m + 1; d < DIM; d++)
409 imageDistance2 -= gmx::square(pbc->box[d][m]);
411 max_d2 = std::min(max_d2, imageDistance2);
417 /* Distance pulling along dimensions with pcrd->params.dim[d]==1.
418 * We use half the minimum box vector length of the dimensions involved.
419 * This is correct for all cases, except for corner cases with
420 * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
421 * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
422 * in such corner cases the user could get correct results,
423 * depending on the details of the setup, we avoid further
424 * code complications.
426 for (int m = 0; m < DIM; m++)
428 if (m < pbc->ndim_ePBC && pcrd->params.dim[m] != 0)
430 real imageDistance2 = gmx::square(pbc->box[m][m]);
431 for (int d = 0; d < m; d++)
433 if (pcrd->params.dim[d] != 0)
435 imageDistance2 += gmx::square(pbc->box[m][d]);
438 max_d2 = std::min(max_d2, imageDistance2);
443 return 0.25 * max_d2;
446 /* This function returns the distance based on coordinates xg and xref.
447 * Note that the pull coordinate struct pcrd is not modified.
449 static void low_get_pull_coord_dr(const struct pull_t* pull,
450 const pull_coord_work_t* pcrd,
457 const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
459 /* Only the first group can be an absolute reference, in that case nat=0 */
460 if (pgrp0->params.nat == 0)
462 for (int m = 0; m < DIM; m++)
464 xref[m] = pcrd->params.origin[m];
469 copy_dvec(xref, xrefr);
471 dvec dref = { 0, 0, 0 };
472 if (pcrd->params.eGeom == epullgDIRPBC)
474 for (int m = 0; m < DIM; m++)
476 dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
478 /* Add the reference position, so we use the correct periodic image */
479 dvec_inc(xrefr, dref);
482 pbc_dx_d(pbc, xg, xrefr, dr);
484 bool directional = pull_coordinate_is_directional(&pcrd->params);
486 for (int m = 0; m < DIM; m++)
488 dr[m] *= pcrd->params.dim[m];
489 if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
491 dr2 += dr[m] * dr[m];
494 /* Check if we are close to switching to another periodic image */
495 if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
497 /* Note that technically there is no issue with switching periodic
498 * image, as pbc_dx_d returns the distance to the closest periodic
499 * image. However in all cases where periodic image switches occur,
500 * the pull results will be useless in practice.
503 "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
504 "box size (%f).\n%s",
505 pcrd->params.group[0], pcrd->params.group[1], sqrt(dr2),
506 sqrt(0.98 * 0.98 * max_dist2), pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
509 if (pcrd->params.eGeom == epullgDIRPBC)
515 /* This function returns the distance based on the contents of the pull struct.
516 * pull->coord[coord_ind].dr, and potentially vector, are updated.
518 static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
520 pull_coord_work_t* pcrd = &pull->coord[coord_ind];
521 PullCoordSpatialData& spatialData = pcrd->spatialData;
524 /* With AWH pulling we allow for periodic pulling with geometry=direction.
525 * TODO: Store a periodicity flag instead of checking for external pull provider.
527 if (pcrd->params.eGeom == epullgDIRPBC
528 || (pcrd->params.eGeom == epullgDIR && pcrd->params.eType == epullEXTERNAL))
534 md2 = static_cast<double>(max_pull_distance2(pcrd, pbc));
537 if (pcrd->params.eGeom == epullgDIRRELATIVE)
539 /* We need to determine the pull vector */
543 const pull_group_work_t& pgrp2 = pull->group[pcrd->params.group[2]];
544 const pull_group_work_t& pgrp3 = pull->group[pcrd->params.group[3]];
546 pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
548 for (m = 0; m < DIM; m++)
550 vec[m] *= pcrd->params.dim[m];
552 spatialData.vec_len = dnorm(vec);
553 for (m = 0; m < DIM; m++)
555 spatialData.vec[m] = vec[m] / spatialData.vec_len;
559 fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
560 coord_ind, vec[XX], vec[YY], vec[ZZ], spatialData.vec[XX], spatialData.vec[YY],
561 spatialData.vec[ZZ]);
565 pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
566 pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
568 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp1->x,
569 pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, md2,
572 if (pcrd->params.ngroup >= 4)
574 pull_group_work_t *pgrp2, *pgrp3;
575 pgrp2 = &pull->group[pcrd->params.group[2]];
576 pgrp3 = &pull->group[pcrd->params.group[3]];
578 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, md2, spatialData.dr23);
580 if (pcrd->params.ngroup >= 6)
582 pull_group_work_t *pgrp4, *pgrp5;
583 pgrp4 = &pull->group[pcrd->params.group[4]];
584 pgrp5 = &pull->group[pcrd->params.group[5]];
586 low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, md2, spatialData.dr45);
590 /* Modify x so that it is periodic in [-pi, pi)
591 * It is assumed that x is in [-3pi, 3pi) so that x
592 * needs to be shifted by at most one period.
594 static void make_periodic_2pi(double* x)
606 /* This function should always be used to modify pcrd->value_ref */
607 static void low_set_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double value_ref)
609 GMX_ASSERT(pcrd->params.eType != epullEXTERNAL,
610 "The pull coord reference value should not be used with type external-potential");
612 if (pcrd->params.eGeom == epullgDIST)
617 "Pull reference distance for coordinate %d (%f) needs to be non-negative",
618 coord_ind + 1, value_ref);
621 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgANGLEAXIS)
623 if (value_ref < 0 || value_ref > M_PI)
626 "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
627 "interval [0,180] deg",
629 value_ref * pull_conversion_factor_internal2userinput(&pcrd->params));
632 else if (pcrd->params.eGeom == epullgDIHEDRAL)
634 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
635 make_periodic_2pi(&value_ref);
638 pcrd->value_ref = value_ref;
641 static void update_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double t)
643 /* With zero rate the reference value is set initially and doesn't change */
644 if (pcrd->params.rate != 0)
646 double value_ref = (pcrd->params.init + pcrd->params.rate * t)
647 * pull_conversion_factor_userinput2internal(&pcrd->params);
648 low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
652 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
653 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
656 dvec dr32; /* store instead of dr23? */
658 dsvmul(-1, spatialData->dr23, dr32);
659 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
660 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
661 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
663 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
664 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
665 * Note 2: the angle between the plane normal vectors equals pi only when
666 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
667 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
669 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
673 /* Calculates pull->coord[coord_ind].value.
674 * This function also updates pull->coord[coord_ind].dr.
676 static void get_pull_coord_distance(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
678 pull_coord_work_t* pcrd;
681 pcrd = &pull->coord[coord_ind];
683 get_pull_coord_dr(pull, coord_ind, pbc);
685 PullCoordSpatialData& spatialData = pcrd->spatialData;
687 switch (pcrd->params.eGeom)
690 /* Pull along the vector between the com's */
691 spatialData.value = dnorm(spatialData.dr01);
695 case epullgDIRRELATIVE:
698 spatialData.value = 0;
699 for (m = 0; m < DIM; m++)
701 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
705 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
707 case epullgDIHEDRAL: spatialData.value = get_dihedral_angle_coord(&spatialData); break;
708 case epullgANGLEAXIS:
709 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
711 default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
715 /* Returns the deviation from the reference value.
716 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
718 static double get_pull_coord_deviation(struct pull_t* pull, int coord_ind, const t_pbc* pbc, double t)
720 pull_coord_work_t* pcrd;
723 pcrd = &pull->coord[coord_ind];
725 /* Update the reference value before computing the distance,
726 * since it is used in the distance computation with periodic pulling.
728 update_pull_coord_reference_value(pcrd, coord_ind, t);
730 get_pull_coord_distance(pull, coord_ind, pbc);
732 get_pull_coord_distance(pull, coord_ind, pbc);
734 /* Determine the deviation */
735 dev = pcrd->spatialData.value - pcrd->value_ref;
737 /* Check that values are allowed */
738 if (pcrd->params.eGeom == epullgDIST && pcrd->spatialData.value == 0)
740 /* With no vector we can not determine the direction for the force,
741 * so we set the force to zero.
745 else if (pcrd->params.eGeom == epullgDIHEDRAL)
747 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
748 Thus, the unwrapped deviation is here in (-2pi, 2pi].
749 After making it periodic, the deviation will be in [-pi, pi). */
750 make_periodic_2pi(&dev);
756 double get_pull_coord_value(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
758 get_pull_coord_distance(pull, coord_ind, pbc);
760 return pull->coord[coord_ind].spatialData.value;
763 void clear_pull_forces(pull_t* pull)
765 /* Zeroing the forces is only required for constraint pulling.
766 * It can happen that multiple constraint steps need to be applied
767 * and therefore the constraint forces need to be accumulated.
769 for (pull_coord_work_t& coord : pull->coord)
771 coord.scalarForce = 0;
775 /* Apply constraint using SHAKE */
777 do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
780 dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
781 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
782 dvec* rnew; /* current 'new' positions of the groups */
783 double* dr_tot; /* the total update of the coords */
786 double lambda, rm, invdt = 0;
787 gmx_bool bConverged_all, bConverged = FALSE;
788 int niter = 0, ii, j, m, max_iter = 100;
792 snew(r_ij, pull->coord.size());
793 snew(dr_tot, pull->coord.size());
795 snew(rnew, pull->group.size());
797 /* copy the current unconstrained positions for use in iterations. We
798 iterate until rinew[i] and rjnew[j] obey the constraints. Then
799 rinew - pull.x_unc[i] is the correction dr to group i */
800 for (size_t g = 0; g < pull->group.size(); g++)
802 copy_dvec(pull->group[g].xp, rnew[g]);
805 /* Determine the constraint directions from the old positions */
806 for (size_t c = 0; c < pull->coord.size(); c++)
808 pull_coord_work_t* pcrd;
810 pcrd = &pull->coord[c];
812 if (pcrd->params.eType != epullCONSTRAINT)
817 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
818 * We don't modify dr and value anymore, so these values are also used
821 get_pull_coord_distance(pull, c, pbc);
823 const PullCoordSpatialData& spatialData = pcrd->spatialData;
826 fprintf(debug, "Pull coord %zu dr %f %f %f\n", c, spatialData.dr01[XX],
827 spatialData.dr01[YY], spatialData.dr01[ZZ]);
830 if (pcrd->params.eGeom == epullgDIR || pcrd->params.eGeom == epullgDIRPBC)
832 /* Select the component along vec */
834 for (m = 0; m < DIM; m++)
836 a += spatialData.vec[m] * spatialData.dr01[m];
838 for (m = 0; m < DIM; m++)
840 r_ij[c][m] = a * spatialData.vec[m];
845 copy_dvec(spatialData.dr01, r_ij[c]);
848 if (dnorm2(r_ij[c]) == 0)
851 "Distance for pull coordinate %zu is zero with constraint pulling, which is "
857 bConverged_all = FALSE;
858 while (!bConverged_all && niter < max_iter)
860 bConverged_all = TRUE;
862 /* loop over all constraints */
863 for (size_t c = 0; c < pull->coord.size(); c++)
865 pull_coord_work_t* pcrd;
866 pull_group_work_t *pgrp0, *pgrp1;
869 pcrd = &pull->coord[c];
871 if (pcrd->params.eType != epullCONSTRAINT)
876 update_pull_coord_reference_value(pcrd, c, t);
878 pgrp0 = &pull->group[pcrd->params.group[0]];
879 pgrp1 = &pull->group[pcrd->params.group[1]];
881 /* Get the current difference vector */
882 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->params.group[1]],
883 rnew[pcrd->params.group[0]], -1, unc_ij);
887 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
890 rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
892 switch (pcrd->params.eGeom)
895 if (pcrd->value_ref <= 0)
899 "The pull constraint reference distance for group %zu is <= 0 (%f)",
904 double q, c_a, c_b, c_c;
906 c_a = diprod(r_ij[c], r_ij[c]);
907 c_b = diprod(unc_ij, r_ij[c]) * 2;
908 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
912 q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
917 q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
923 fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b,
928 /* The position corrections dr due to the constraints */
929 dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
930 dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
931 dr_tot[c] += -lambda * dnorm(r_ij[c]);
936 /* A 1-dimensional constraint along a vector */
938 for (m = 0; m < DIM; m++)
940 vec[m] = pcrd->spatialData.vec[m];
941 a += unc_ij[m] * vec[m];
943 /* Select only the component along the vector */
944 dsvmul(a, vec, unc_ij);
945 lambda = a - pcrd->value_ref;
948 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
951 /* The position corrections dr due to the constraints */
952 dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
953 dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
954 dr_tot[c] += -lambda;
956 default: gmx_incons("Invalid enumeration value for eGeom");
964 g0 = pcrd->params.group[0];
965 g1 = pcrd->params.group[1];
966 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
967 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
968 fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", rnew[g0][0],
969 rnew[g0][1], rnew[g0][2], rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
970 fprintf(debug, "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", "", "", "", "", "",
971 "", pcrd->value_ref);
972 fprintf(debug, "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", dr0[0],
973 dr0[1], dr0[2], dr1[0], dr1[1], dr1[2], dnorm(tmp3));
976 /* Update the COMs with dr */
977 dvec_inc(rnew[pcrd->params.group[1]], dr1);
978 dvec_inc(rnew[pcrd->params.group[0]], dr0);
981 /* Check if all constraints are fullfilled now */
982 for (pull_coord_work_t& coord : pull->coord)
984 if (coord.params.eType != epullCONSTRAINT)
989 low_get_pull_coord_dr(pull, &coord, pbc, rnew[coord.params.group[1]],
990 rnew[coord.params.group[0]], -1, unc_ij);
992 switch (coord.params.eGeom)
995 bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1000 for (m = 0; m < DIM; m++)
1002 vec[m] = coord.spatialData.vec[m];
1004 inpr = diprod(unc_ij, vec);
1005 dsvmul(inpr, vec, unc_ij);
1006 bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1015 "Pull constraint not converged: "
1017 "d_ref = %f, current d = %f\n",
1018 coord.params.group[0], coord.params.group[1], coord.value_ref, dnorm(unc_ij));
1021 bConverged_all = FALSE;
1026 /* if after all constraints are dealt with and bConverged is still TRUE
1027 we're finished, if not we do another iteration */
1029 if (niter > max_iter)
1031 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1034 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1041 /* update atoms in the groups */
1042 for (size_t g = 0; g < pull->group.size(); g++)
1044 const pull_group_work_t* pgrp;
1047 pgrp = &pull->group[g];
1049 /* get the final constraint displacement dr for group g */
1050 dvec_sub(rnew[g], pgrp->xp, dr);
1052 if (dnorm2(dr) == 0)
1054 /* No displacement, no update necessary */
1058 /* update the atom positions */
1059 auto localAtomIndices = pgrp->atomSet.localIndex();
1061 for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1063 ii = localAtomIndices[j];
1064 if (!pgrp->localWeights.empty())
1066 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1068 for (m = 0; m < DIM; m++)
1074 for (m = 0; m < DIM; m++)
1076 v[ii][m] += invdt * tmp[m];
1082 /* calculate the constraint forces, used for output and virial only */
1083 for (size_t c = 0; c < pull->coord.size(); c++)
1085 pull_coord_work_t* pcrd;
1087 pcrd = &pull->coord[c];
1089 if (pcrd->params.eType != epullCONSTRAINT)
1094 /* Accumulate the forces, in case we have multiple constraint steps */
1097 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1099 pcrd->scalarForce += force;
1101 if (vir != nullptr && pcrd->params.eGeom != epullgDIRPBC && bMaster)
1105 /* Add the pull contribution to the virial */
1106 /* We have already checked above that r_ij[c] != 0 */
1107 f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1109 for (j = 0; j < DIM; j++)
1111 for (m = 0; m < DIM; m++)
1113 vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1119 /* finished! I hope. Give back some memory */
1125 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1127 for (int j = 0; j < DIM; j++)
1129 for (int m = 0; m < DIM; m++)
1131 vir[j][m] -= 0.5 * f[j] * dr[m];
1136 /* Adds the pull contribution to the virial */
1137 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1139 if (vir != nullptr && pcrd.params.eGeom != epullgDIRPBC)
1141 /* Add the pull contribution for each distance vector to the virial. */
1142 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1143 if (pcrd.params.ngroup >= 4)
1145 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1147 if (pcrd.params.ngroup >= 6)
1149 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1154 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1162 k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1163 dkdl = pcrd->params.kB - pcrd->params.k;
1165 switch (pcrd->params.eType)
1168 case epullFLATBOTTOM:
1169 case epullFLATBOTTOMHIGH:
1170 /* The only difference between an umbrella and a flat-bottom
1171 * potential is that a flat-bottom is zero above or below
1172 the reference value.
1174 if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0)
1175 || (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
1180 pcrd->scalarForce = -k * dev;
1181 *V += 0.5 * k * gmx::square(dev);
1182 *dVdl += 0.5 * dkdl * gmx::square(dev);
1185 pcrd->scalarForce = -k;
1186 *V += k * pcrd->spatialData.value;
1187 *dVdl += dkdl * pcrd->spatialData.value;
1191 "the scalar pull force should not be calculated internally for pull type "
1193 default: gmx_incons("Unsupported pull type in do_pull_pot");
1197 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1199 const t_pull_coord& params = pcrd.params;
1200 const PullCoordSpatialData& spatialData = pcrd.spatialData;
1202 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1203 PullCoordVectorForces forces;
1205 if (params.eGeom == epullgDIST)
1207 double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1208 for (int m = 0; m < DIM; m++)
1210 forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1213 else if (params.eGeom == epullgANGLE)
1216 double cos_theta, cos_theta2;
1218 cos_theta = cos(spatialData.value);
1219 cos_theta2 = gmx::square(cos_theta);
1221 /* The force at theta = 0, pi is undefined so we should not apply any force.
1222 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1223 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1224 * have to check for this before dividing by their norm below.
1228 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1229 * and another vector parallel to dr23:
1230 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1231 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1233 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1234 double b = a * cos_theta;
1235 double invdr01 = 1. / dnorm(spatialData.dr01);
1236 double invdr23 = 1. / dnorm(spatialData.dr23);
1237 dvec normalized_dr01, normalized_dr23;
1238 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1239 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1241 for (int m = 0; m < DIM; m++)
1243 /* Here, f_scal is -dV/dtheta */
1245 pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1247 pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1252 /* No forces to apply for ill-defined cases*/
1253 clear_dvec(forces.force01);
1254 clear_dvec(forces.force23);
1257 else if (params.eGeom == epullgANGLEAXIS)
1259 double cos_theta, cos_theta2;
1261 /* The angle-axis force is exactly the same as the angle force (above) except that in
1262 this case the second vector (dr23) is replaced by the pull vector. */
1263 cos_theta = cos(spatialData.value);
1264 cos_theta2 = gmx::square(cos_theta);
1270 dvec normalized_dr01;
1272 invdr01 = 1. / dnorm(spatialData.dr01);
1273 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1274 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1277 for (int m = 0; m < DIM; m++)
1280 pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1285 clear_dvec(forces.force01);
1288 else if (params.eGeom == epullgDIHEDRAL)
1290 double m2, n2, tol, sqrdist_32;
1292 /* Note: there is a small difference here compared to the
1293 dihedral force calculations in the bondeds (ref: Bekker 1994).
1294 There rij = ri - rj, while here dr01 = r1 - r0.
1295 However, all distance vectors occur in form of cross or inner products
1296 so that two signs cancel and we end up with the same expressions.
1297 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1299 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1300 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1301 dsvmul(-1, spatialData.dr23, dr32);
1302 sqrdist_32 = diprod(dr32, dr32);
1303 tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1304 if ((m2 > tol) && (n2 > tol))
1306 double a_01, a_23_01, a_23_45, a_45;
1307 double inv_dist_32, inv_sqrdist_32, dist_32;
1309 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1310 inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1311 dist_32 = sqrdist_32 * inv_dist_32;
1313 /* Forces on groups 0, 1 */
1314 a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
1315 dsvmul(-a_01, spatialData.planevec_m,
1316 forces.force01); /* added sign to get force on group 1, not 0 */
1318 /* Forces on groups 4, 5 */
1319 a_45 = -pcrd.scalarForce * dist_32 / n2;
1320 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1322 /* Force on groups 2, 3 (defining the axis) */
1323 a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1324 a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1325 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1326 dsvmul(a_23_45, forces.force45, v);
1327 dvec_sub(u, v, forces.force23); /* force on group 3 */
1331 /* No force to apply for ill-defined cases */
1332 clear_dvec(forces.force01);
1333 clear_dvec(forces.force23);
1334 clear_dvec(forces.force45);
1339 for (int m = 0; m < DIM; m++)
1341 forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1349 /* We use a global mutex for locking access to the pull data structure
1350 * during registration of external pull potential providers.
1351 * We could use a different, local mutex for each pull object, but the overhead
1352 * is extremely small here and registration is only done during initialization.
1354 static gmx::Mutex registrationMutex;
1356 using Lock = gmx::lock_guard<gmx::Mutex>;
1358 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1360 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1361 GMX_RELEASE_ASSERT(provider != nullptr,
1362 "register_external_pull_potential called with NULL as provider name");
1364 if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1367 "Module '%s' attempted to register an external potential for pull coordinate %d "
1368 "which is out of the pull coordinate range %d - %zu\n",
1369 provider, coord_index + 1, 1, pull->coord.size());
1372 pull_coord_work_t* pcrd = &pull->coord[coord_index];
1374 if (pcrd->params.eType != epullEXTERNAL)
1378 "Module '%s' attempted to register an external potential for pull coordinate %d "
1379 "which of type '%s', whereas external potentials are only supported with type '%s'",
1380 provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
1383 GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr,
1384 "The external potential provider string for a pull coordinate is NULL");
1386 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider) != 0)
1389 "Module '%s' attempted to register an external potential for pull coordinate %d "
1390 "which expects the external potential to be provided by a module named '%s'",
1391 provider, coord_index + 1, pcrd->params.externalPotentialProvider);
1394 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1395 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1396 * pull->numUnregisteredExternalPotentials.
1398 Lock registrationLock(registrationMutex);
1400 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1403 "Module '%s' attempted to register an external potential for pull coordinate %d "
1405 provider, coord_index + 1);
1408 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1409 pull->numUnregisteredExternalPotentials--;
1411 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1412 "Negative unregistered potentials, the pull code is inconsistent");
1416 static void check_external_potential_registration(const struct pull_t* pull)
1418 if (pull->numUnregisteredExternalPotentials > 0)
1421 for (c = 0; c < pull->coord.size(); c++)
1423 if (pull->coord[c].params.eType == epullEXTERNAL
1424 && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1430 GMX_RELEASE_ASSERT(c < pull->coord.size(),
1431 "Internal inconsistency in the pull potential provider counting");
1434 "No external provider for external pull potentials have been provided for %d "
1435 "pull coordinates. The first coordinate without provider is number %zu, which "
1436 "expects a module named '%s' to provide the external potential.",
1437 pull->numUnregisteredExternalPotentials, c + 1,
1438 pull->coord[c].params.externalPotentialProvider);
1443 /* Pull takes care of adding the forces of the external potential.
1444 * The external potential module has to make sure that the corresponding
1445 * potential energy is added either to the pull term or to a term
1446 * specific to the external module.
1448 void apply_external_pull_coord_force(struct pull_t* pull,
1452 gmx::ForceWithVirial* forceWithVirial)
1454 pull_coord_work_t* pcrd;
1456 GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1457 "apply_external_pull_coord_force called with coord_index out of range");
1459 if (pull->comm.bParticipate)
1461 pcrd = &pull->coord[coord_index];
1464 pcrd->params.eType == epullEXTERNAL,
1465 "The pull force can only be set externally on pull coordinates of external type");
1467 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1468 "apply_external_pull_coord_force called for an unregistered pull coordinate");
1471 pcrd->scalarForce = coord_force;
1473 /* Calculate the forces on the pull groups */
1474 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1476 /* Add the forces for this coordinate to the total virial and force */
1477 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1479 matrix virial = { { 0 } };
1480 add_virial_coord(virial, *pcrd, pullCoordForces);
1481 forceWithVirial->addVirialContribution(virial);
1484 apply_forces_coord(pull, coord_index, pullCoordForces, masses,
1485 as_rvec_array(forceWithVirial->force_.data()));
1488 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1491 /* Calculate the pull potential and scalar force for a pull coordinate.
1492 * Returns the vector forces for the pull coordinate.
1494 static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull,
1503 pull_coord_work_t& pcrd = pull->coord[coord_ind];
1505 assert(pcrd.params.eType != epullCONSTRAINT);
1507 double dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
1509 calc_pull_coord_scalar_force_and_potential(&pcrd, dev, lambda, V, dVdl);
1511 PullCoordVectorForces pullCoordForces = calculateVectorForces(pcrd);
1513 add_virial_coord(vir, pcrd, pullCoordForces);
1515 return pullCoordForces;
1518 real pull_potential(struct pull_t* pull,
1521 const t_commrec* cr,
1525 gmx::ForceWithVirial* force,
1530 assert(pull != nullptr);
1532 /* Ideally we should check external potential registration only during
1533 * the initialization phase, but that requires another function call
1534 * that should be done exactly in the right order. So we check here.
1536 check_external_potential_registration(pull);
1538 if (pull->comm.bParticipate)
1542 pull_calc_coms(cr, pull, masses, pbc, t, x, nullptr);
1544 rvec* f = as_rvec_array(force->force_.data());
1545 matrix virial = { { 0 } };
1546 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1547 for (size_t c = 0; c < pull->coord.size(); c++)
1549 pull_coord_work_t* pcrd;
1550 pcrd = &pull->coord[c];
1552 /* For external potential the force is assumed to be given by an external module by a
1553 call to apply_pull_coord_external_force */
1554 if (pcrd->params.eType == epullCONSTRAINT || pcrd->params.eType == epullEXTERNAL)
1559 PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1560 pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1562 /* Distribute the force over the atoms in the pulled groups */
1563 apply_forces_coord(pull, c, pullCoordForces, masses, f);
1568 force->addVirialContribution(virial);
1573 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1574 "Too few or too many external pull potentials have been applied the previous step");
1575 /* All external pull potentials still need to be applied */
1576 pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1578 return (MASTER(cr) ? V : 0.0);
1581 void pull_constraint(struct pull_t* pull,
1584 const t_commrec* cr,
1592 assert(pull != nullptr);
1594 if (pull->comm.bParticipate)
1596 pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
1598 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1602 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1606 gmx_bool bMustParticipate;
1612 /* We always make the master node participate, such that it can do i/o,
1613 * add the virial and to simplify MC type extensions people might have.
1615 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1617 for (pull_group_work_t& group : pull->group)
1619 if (!group.globalWeights.empty())
1621 group.localWeights.resize(group.atomSet.numAtomsLocal());
1622 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1624 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1628 GMX_ASSERT(bMustParticipate || dd != nullptr,
1629 "Either all ranks (including this rank) participate, or we use DD and need to "
1630 "have access to dd here");
1632 /* We should participate if we have pull or pbc atoms */
1633 if (!bMustParticipate
1634 && (group.atomSet.numAtomsLocal() > 0
1635 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1637 bMustParticipate = TRUE;
1641 if (!comm->bParticipateAll)
1643 /* Keep currently not required ranks in the communicator
1644 * if they needed to participate up to 20 decompositions ago.
1645 * This avoids frequent rebuilds due to atoms jumping back and forth.
1647 const int64_t history_count = 20;
1648 gmx_bool bWillParticipate;
1651 /* Increase the decomposition counter for the current call */
1652 comm->setup_count++;
1654 if (bMustParticipate)
1656 comm->must_count = comm->setup_count;
1661 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1663 if (debug && dd != nullptr)
1665 fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n", dd->rank,
1666 gmx::boolToString(bMustParticipate), gmx::boolToString(bWillParticipate));
1669 if (bWillParticipate)
1671 /* Count the number of ranks that we want to have participating */
1673 /* Count the number of ranks that need to be added */
1674 count[1] = comm->bParticipate ? 0 : 1;
1682 /* The cost of this global operation will be less that the cost
1683 * of the extra MPI_Comm_split calls that we can avoid.
1685 gmx_sumi(2, count, cr);
1687 /* If we are missing ranks or if we have 20% more ranks than needed
1688 * we make a new sub-communicator.
1690 if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1694 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1698 if (comm->mpi_comm_com != MPI_COMM_NULL)
1700 MPI_Comm_free(&comm->mpi_comm_com);
1703 /* This might be an extremely expensive operation, so we try
1704 * to avoid this splitting as much as possible.
1706 assert(dd != nullptr);
1707 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1710 comm->bParticipate = bWillParticipate;
1711 comm->nparticipate = count[0];
1713 /* When we use the previous COM for PBC, we need to broadcast
1714 * the previous COM to ranks that have joined the communicator.
1716 for (pull_group_work_t& group : pull->group)
1718 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1720 GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1721 "The master rank has to participate, as it should pass an up to "
1723 "to bcast here as well as to e.g. checkpointing");
1725 gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
1731 /* Since the PBC of atoms might have changed, we need to update the PBC */
1732 pull->bSetPBCatoms = TRUE;
1735 static void init_pull_group_index(FILE* fplog,
1736 const t_commrec* cr,
1738 pull_group_work_t* pg,
1739 gmx_bool bConstraint,
1740 const ivec pulldim_con,
1741 const gmx_mtop_t* mtop,
1742 const t_inputrec* ir,
1745 /* With EM and BD there are no masses in the integrator.
1746 * But we still want to have the correct mass-weighted COMs.
1747 * So we store the real masses in the weights.
1749 const bool setWeights = (pg->params.nweight > 0 || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
1751 /* In parallel, store we need to extract localWeights from weights at DD time */
1752 std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1754 const SimulationGroups& groups = mtop->groups;
1756 /* Count frozen dimensions and (weighted) mass */
1762 for (int i = 0; i < pg->params.nat; i++)
1764 int ii = pg->params.ind[i];
1765 if (bConstraint && ir->opts.nFreeze)
1767 for (int d = 0; d < DIM; d++)
1769 if (pulldim_con[d] == 1
1770 && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1776 const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1778 if (ir->efep == efepNO)
1784 m = (1 - lambda) * atom.m + lambda * atom.mB;
1787 if (pg->params.nweight > 0)
1789 w = pg->params.weight[i];
1795 if (EI_ENERGY_MINIMIZATION(ir->eI))
1797 /* Move the mass to the weight */
1801 else if (ir->eI == eiBD)
1804 if (ir->bd_fric != 0.0F)
1806 mbd = ir->bd_fric * ir->delta_t;
1810 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1812 mbd = ir->delta_t / ir->opts.tau_t[0];
1817 / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1825 weights.push_back(w);
1829 wwmass += m * w * w;
1834 /* We can have single atom groups with zero mass with potential pulling
1835 * without cosine weighting.
1837 if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1839 /* With one atom the mass doesn't matter */
1844 gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1845 pg->params.weight ? " weighted" : "", g);
1850 fprintf(fplog, "Pull group %d: %5d atoms, mass %9.3f", g, pg->params.nat, tmass);
1851 if (pg->params.weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1853 fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1855 if (pg->epgrppbc == epgrppbcCOS)
1857 fprintf(fplog, ", cosine weighting will be used");
1859 fprintf(fplog, "\n");
1864 /* A value != 0 signals not frozen, it is updated later */
1870 for (int d = 0; d < DIM; d++)
1872 ndim += pulldim_con[d] * pg->params.nat;
1874 if (fplog && nfrozen > 0 && nfrozen < ndim)
1877 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1878 " that are subject to pulling are frozen.\n"
1879 " For constraint pulling the whole group will be frozen.\n\n",
1887 struct pull_t* init_pull(FILE* fplog,
1888 const pull_params_t* pull_params,
1889 const t_inputrec* ir,
1890 const gmx_mtop_t* mtop,
1891 const t_commrec* cr,
1892 gmx::LocalAtomSetManager* atomSets,
1895 struct pull_t* pull;
1900 /* Copy the pull parameters */
1901 pull->params = *pull_params;
1902 /* Avoid pointer copies */
1903 pull->params.group = nullptr;
1904 pull->params.coord = nullptr;
1906 for (int i = 0; i < pull_params->ngroup; ++i)
1908 pull->group.emplace_back(pull_params->group[i],
1909 atomSets->add({ pull_params->group[i].ind,
1910 pull_params->group[i].ind + pull_params->group[i].nat }),
1911 pull_params->bSetPbcRefToPrevStepCOM);
1914 if (cr != nullptr && DOMAINDECOMP(cr))
1916 /* Set up the global to local atom mapping for PBC atoms */
1917 for (pull_group_work_t& group : pull->group)
1919 if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
1921 /* pbcAtomSet consists of a single atom */
1922 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
1923 atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
1928 pull->bPotential = FALSE;
1929 pull->bConstraint = FALSE;
1930 pull->bCylinder = FALSE;
1931 pull->bAngle = FALSE;
1932 pull->bXOutAverage = pull_params->bXOutAverage;
1933 pull->bFOutAverage = pull_params->bFOutAverage;
1935 GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0,
1936 "pull group 0 is an absolute reference group and should not contain atoms");
1938 pull->numCoordinatesWithExternalPotential = 0;
1940 for (int c = 0; c < pull->params.ncoord; c++)
1942 /* Construct a pull coordinate, copying all coordinate parameters */
1943 pull->coord.emplace_back(pull_params->coord[c]);
1945 pull_coord_work_t* pcrd = &pull->coord.back();
1947 switch (pcrd->params.eGeom)
1950 case epullgDIRRELATIVE: /* Direction vector is determined at each step */
1952 case epullgDIHEDRAL: break;
1956 case epullgANGLEAXIS:
1957 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
1960 /* We allow reading of newer tpx files with new pull geometries,
1961 * but with the same tpx format, with old code. A new geometry
1962 * only adds a new enum value, which will be out of range for
1963 * old code. The first place we need to generate an error is
1964 * here, since the pull code can't handle this.
1965 * The enum can be used before arriving here only for printing
1966 * the string corresponding to the geometry, which will result
1967 * in printing "UNDEFINED".
1970 "Pull geometry not supported for pull coordinate %d. The geometry enum "
1971 "%d in the input is larger than that supported by the code (up to %d). "
1972 "You are probably reading a tpr file generated with a newer version of "
1973 "Gromacs with an binary from an older version of Gromacs.",
1974 c + 1, pcrd->params.eGeom, epullgNR - 1);
1977 if (pcrd->params.eType == epullCONSTRAINT)
1979 /* Check restrictions of the constraint pull code */
1980 if (pcrd->params.eGeom == epullgCYL || pcrd->params.eGeom == epullgDIRRELATIVE
1981 || pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
1982 || pcrd->params.eGeom == epullgANGLEAXIS)
1985 "Pulling of type %s can not be combined with geometry %s. Consider using "
1987 epull_names[pcrd->params.eType], epullg_names[pcrd->params.eGeom],
1988 epull_names[epullUMBRELLA]);
1991 pull->bConstraint = TRUE;
1995 pull->bPotential = TRUE;
1998 if (pcrd->params.eGeom == epullgCYL)
2000 pull->bCylinder = TRUE;
2002 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
2003 || pcrd->params.eGeom == epullgANGLEAXIS)
2005 pull->bAngle = TRUE;
2008 /* We only need to calculate the plain COM of a group
2009 * when it is not only used as a cylinder group.
2010 * Also the absolute reference group 0 needs no COM computation.
2012 for (int i = 0; i < pcrd->params.ngroup; i++)
2014 int groupIndex = pcrd->params.group[i];
2015 if (groupIndex > 0 && !(pcrd->params.eGeom == epullgCYL && i == 0))
2017 pull->group[groupIndex].needToCalcCom = true;
2021 /* With non-zero rate the reference value is set at every step */
2022 if (pcrd->params.rate == 0)
2024 /* Initialize the constant reference value */
2025 if (pcrd->params.eType != epullEXTERNAL)
2027 low_set_pull_coord_reference_value(
2029 pcrd->params.init * pull_conversion_factor_userinput2internal(&pcrd->params));
2033 /* With an external pull potential, the reference value loses
2034 * it's meaning and should not be used. Setting it to zero
2035 * makes any terms dependent on it disappear.
2036 * The only issue this causes is that with cylinder pulling
2037 * no atoms of the cylinder group within the cylinder radius
2038 * should be more than half a box length away from the COM of
2039 * the pull group along the axial direction.
2041 pcrd->value_ref = 0.0;
2045 if (pcrd->params.eType == epullEXTERNAL)
2048 pcrd->params.rate == 0,
2049 "With an external potential, a pull coordinate should have rate = 0");
2051 /* This external potential needs to be registered later */
2052 pull->numCoordinatesWithExternalPotential++;
2054 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2057 pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
2058 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2060 pull->pbcType = ir->pbcType;
2061 switch (pull->pbcType)
2063 case PbcType::No: pull->npbcdim = 0; break;
2064 case PbcType::XY: pull->npbcdim = 2; break;
2065 default: pull->npbcdim = 3; break;
2070 gmx_bool bAbs, bCos;
2073 for (const pull_coord_work_t& coord : pull->coord)
2075 if (pull->group[coord.params.group[0]].params.nat == 0
2076 || pull->group[coord.params.group[1]].params.nat == 0)
2082 fprintf(fplog, "\n");
2083 if (pull->bPotential)
2085 fprintf(fplog, "Will apply potential COM pulling\n");
2087 if (pull->bConstraint)
2089 fprintf(fplog, "Will apply constraint COM pulling\n");
2091 // Don't include the reference group 0 in output, so we report ngroup-1
2092 int numRealGroups = pull->group.size() - 1;
2093 GMX_RELEASE_ASSERT(numRealGroups > 0,
2094 "The reference absolute position pull group should always be present");
2095 fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n", pull->coord.size(),
2096 pull->coord.size() == 1 ? "" : "s", numRealGroups, numRealGroups == 1 ? "" : "s");
2099 fprintf(fplog, "with an absolute reference\n");
2102 // Don't include the reference group 0 in loop
2103 for (size_t g = 1; g < pull->group.size(); g++)
2105 if (pull->group[g].params.nat > 1 && pull->group[g].params.pbcatom < 0)
2107 /* We are using cosine weighting */
2108 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2114 please_cite(fplog, "Engin2010");
2118 pull->bRefAt = FALSE;
2120 for (size_t g = 0; g < pull->group.size(); g++)
2122 pull_group_work_t* pgrp;
2124 pgrp = &pull->group[g];
2125 if (pgrp->params.nat > 0)
2127 /* There is an issue when a group is used in multiple coordinates
2128 * and constraints are applied in different dimensions with atoms
2129 * frozen in some, but not all dimensions.
2130 * Since there is only one mass array per group, we can't have
2131 * frozen/non-frozen atoms for different coords at the same time.
2132 * But since this is a very exotic case, we don't check for this.
2133 * A warning is printed in init_pull_group_index.
2136 gmx_bool bConstraint;
2137 ivec pulldim, pulldim_con;
2139 /* Loop over all pull coordinates to see along which dimensions
2140 * this group is pulled and if it is involved in constraints.
2142 bConstraint = FALSE;
2143 clear_ivec(pulldim);
2144 clear_ivec(pulldim_con);
2145 for (const pull_coord_work_t& coord : pull->coord)
2147 gmx_bool bGroupUsed = FALSE;
2148 for (int gi = 0; gi < coord.params.ngroup; gi++)
2150 if (coord.params.group[gi] == static_cast<int>(g))
2158 for (int m = 0; m < DIM; m++)
2160 if (coord.params.dim[m] == 1)
2164 if (coord.params.eType == epullCONSTRAINT)
2174 /* Determine if we need to take PBC into account for calculating
2175 * the COM's of the pull groups.
2177 switch (pgrp->epgrppbc)
2179 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2181 if (pgrp->params.weight != nullptr)
2184 "Pull groups can not have relative weights and cosine weighting "
2187 for (int m = 0; m < DIM; m++)
2189 if (m < pull->npbcdim && pulldim[m] == 1)
2191 if (pull->cosdim >= 0 && pull->cosdim != m)
2194 "Can only use cosine weighting with pulling in one "
2195 "dimension (use mdp option pull-coord?-dim)");
2201 case epgrppbcNONE: break;
2204 /* Set the indices */
2205 init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2209 /* Absolute reference, set the inverse mass to zero.
2210 * This is only relevant (and used) with constraint pulling.
2217 /* If we use cylinder coordinates, do some initialising for them */
2218 if (pull->bCylinder)
2220 for (const pull_coord_work_t& coord : pull->coord)
2222 if (coord.params.eGeom == epullgCYL)
2224 if (pull->group[coord.params.group[0]].params.nat == 0)
2227 "A cylinder pull group is not supported when using absolute "
2231 const auto& referenceGroup = pull->group[coord.params.group[0]];
2232 pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet,
2233 pull->params.bSetPbcRefToPrevStepCOM);
2237 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2238 pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
2239 pull->comSums.resize(pull->nthreads);
2244 /* Use a sub-communicator when we have more than 32 ranks, but not
2245 * when we have an external pull potential, since then the external
2246 * potential provider expects each rank to have the coordinate.
2248 comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
2249 || pull->numCoordinatesWithExternalPotential > 0
2250 || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2251 /* This sub-commicator is not used with comm->bParticipateAll,
2252 * so we can always initialize it to NULL.
2254 comm->mpi_comm_com = MPI_COMM_NULL;
2255 comm->nparticipate = 0;
2256 comm->isMasterRank = (cr == nullptr || MASTER(cr));
2258 /* No MPI: 1 rank: all ranks pull */
2259 comm->bParticipateAll = TRUE;
2260 comm->isMasterRank = true;
2262 comm->bParticipate = comm->bParticipateAll;
2263 comm->setup_count = 0;
2264 comm->must_count = 0;
2266 if (!comm->bParticipateAll && fplog != nullptr)
2268 fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2271 comm->pbcAtomBuffer.resize(pull->group.size());
2272 comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
2273 if (pull->bCylinder)
2275 comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
2278 /* We still need to initialize the PBC reference coordinates */
2279 pull->bSetPBCatoms = TRUE;
2281 pull->out_x = nullptr;
2282 pull->out_f = nullptr;
2287 static void destroy_pull(struct pull_t* pull)
2290 if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2292 MPI_Comm_free(&pull->comm.mpi_comm_com);
2299 void preparePrevStepPullCom(const t_inputrec* ir,
2303 const t_state* state_global,
2304 const t_commrec* cr,
2305 bool startingFromCheckpoint)
2307 if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
2311 allocStatePrevStepPullCom(state, pull_work);
2312 if (startingFromCheckpoint)
2316 state->pull_com_prev_step = state_global->pull_com_prev_step;
2320 /* Only the master rank has the checkpointed COM from the previous step */
2321 gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(),
2322 &state->pull_com_prev_step[0], cr->mpi_comm_mygroup);
2324 setPrevStepPullComFromState(pull_work, state);
2329 set_pbc(&pbc, ir->pbcType, state->box);
2330 initPullComFromPrevStep(cr, pull_work, masses, &pbc, state->x.rvec_array());
2331 updatePrevStepPullCom(pull_work, state);
2335 void finish_pull(struct pull_t* pull)
2337 check_external_potential_registration(pull);
2341 gmx_fio_fclose(pull->out_x);
2345 gmx_fio_fclose(pull->out_f);
2351 gmx_bool pull_have_potential(const struct pull_t* pull)
2353 return pull->bPotential;
2356 gmx_bool pull_have_constraint(const struct pull_t* pull)
2358 return pull->bConstraint;
2361 bool pull_have_constraint(const pull_params_t* pullParameters)
2363 if (pullParameters == nullptr)
2367 for (int c = 0; c < pullParameters->ncoord; c++)
2369 if (pullParameters->coord[c].eType == epullCONSTRAINT)