src/gromacs/pulling/pull.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "pull.h"
  41
  42 #include "config.h"
  43
  44 #include <cassert>
  45 #include <cinttypes>
  46 #include <cmath>
  47 #include <cstdlib>
  48
  49 #include <algorithm>
  50 #include <memory>
  51 #include <mutex>
  52 #include <string>
  53
  54 #include "gromacs/commandline/filenm.h"
  55 #include "gromacs/domdec/domdec_struct.h"
  56 #include "gromacs/domdec/localatomset.h"
  57 #include "gromacs/domdec/localatomsetmanager.h"
  58 #include "gromacs/fileio/gmxfio.h"
  59 #include "gromacs/gmxlib/network.h"
  60 #include "gromacs/math/functions.h"
  61 #include "gromacs/math/units.h"
  62 #include "gromacs/math/utilities.h"
  63 #include "gromacs/math/vec.h"
  64 #include "gromacs/math/vectypes.h"
  65 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  66 #include "gromacs/mdtypes/commrec.h"
  67 #include "gromacs/mdtypes/forceoutput.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/md_enums.h"
  70 #include "gromacs/mdtypes/state.h"
  71 #include "gromacs/pbcutil/pbc.h"
  72 #include "gromacs/topology/mtop_lookup.h"
  73 #include "gromacs/topology/topology.h"
  74 #include "gromacs/utility/arrayref.h"
  75 #include "gromacs/utility/cstringutil.h"
  76 #include "gromacs/utility/enumerationhelpers.h"
  77 #include "gromacs/utility/exceptions.h"
  78 #include "gromacs/utility/fatalerror.h"
  79 #include "gromacs/utility/futil.h"
  80 #include "gromacs/utility/gmxassert.h"
  81 #include "gromacs/utility/pleasecite.h"
  82 #include "gromacs/utility/real.h"
  83 #include "gromacs/utility/smalloc.h"
  84 #include "gromacs/utility/strconvert.h"
  85 #include "gromacs/utility/stringutil.h"
  86
  87 #include "pull_internal.h"
  88 #include "transformationcoordinate.h"
  89
  90 namespace gmx
  91 {
  92 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
  93 } // namespace gmx
  94
  95 using gmx::ArrayRef;
  96 using gmx::RVec;
  97
  98 /*! \brief Tells whether the pull geometry is an angle type */
  99 constexpr gmx::EnumerationArray<PullGroupGeometry, bool> sc_isAngleType = {
 100     { false, false, false, false, false, true, true, true }
 101 };
 102
 103 static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
 104 {
 105     if (params.ind.size() <= 1)
 106     {
 107         /* no PBC required */
 108         return epgrppbcNONE;
 109     }
 110     else if (params.pbcatom >= 0)
 111     {
 112         if (setPbcRefToPrevStepCOM)
 113         {
 114             return epgrppbcPREVSTEPCOM;
 115         }
 116         else
 117         {
 118             return epgrppbcREFAT;
 119         }
 120     }
 121     else
 122     {
 123         return epgrppbcCOS;
 124     }
 125 }
 126
 127 /* NOTE: The params initialization currently copies pointers.
 128  *       So the lifetime of the source, currently always inputrec,
 129  *       should not end before that of this object.
 130  *       This will be fixed when the pointers are replacted by std::vector.
 131  */
 132 pull_group_work_t::pull_group_work_t(const t_pull_group& params,
 133                                      gmx::LocalAtomSet   atomSet,
 134                                      bool                bSetPbcRefToPrevStepCOM,
 135                                      int                 maxNumThreads) :
 136     params(params),
 137     epgrppbc(groupPbcFromParams(params, bSetPbcRefToPrevStepCOM)),
 138     maxNumThreads(maxNumThreads),
 139     needToCalcCom(false),
 140     atomSet(atomSet),
 141     mwscale(0),
 142     wscale(0),
 143     invtm(0)
 144 {
 145     clear_dvec(x);
 146     clear_dvec(xp);
 147 };
 148
 149 static bool pull_coordinate_is_directional(const t_pull_coord& pcrd)
 150 {
 151     return (pcrd.eGeom == PullGroupGeometry::Direction || pcrd.eGeom == PullGroupGeometry::DirectionPBC
 152             || pcrd.eGeom == PullGroupGeometry::DirectionRelative
 153             || pcrd.eGeom == PullGroupGeometry::Cylinder);
 154 }
 155
 156 const char* pull_coordinate_units(const t_pull_coord& pcrd)
 157 {
 158     return sc_isAngleType[pcrd.eGeom] ? "deg" : "nm";
 159 }
 160
 161 double pull_conversion_factor_userinput2internal(const t_pull_coord& pcrd)
 162 {
 163     if (sc_isAngleType[pcrd.eGeom])
 164     {
 165         return gmx::c_deg2Rad;
 166     }
 167     else
 168     {
 169         return 1.0;
 170     }
 171 }
 172
 173 double pull_conversion_factor_internal2userinput(const t_pull_coord& pcrd)
 174 {
 175     if (sc_isAngleType[pcrd.eGeom])
 176     {
 177         return gmx::c_rad2Deg;
 178     }
 179     else
 180     {
 181         return 1.0;
 182     }
 183 }
 184
 185 /* Apply forces in a mass weighted fashion for part of the pull group */
 186 static void apply_forces_grp_part(const pull_group_work_t& pgrp,
 187                                   int                      ind_start,
 188                                   int                      ind_end,
 189                                   ArrayRef<const real>     masses,
 190                                   const dvec               f_pull,
 191                                   const int                sign,
 192                                   rvec*                    f)
 193 {
 194     const double inv_wm = pgrp.mwscale;
 195
 196     auto localAtomIndices = pgrp.atomSet.localIndex();
 197     for (int i = ind_start; i < ind_end; i++)
 198     {
 199         int    ii    = localAtomIndices[i];
 200         double wmass = masses[ii];
 201         if (!pgrp.localWeights.empty())
 202         {
 203             wmass *= pgrp.localWeights[i];
 204         }
 205
 206         for (int d = 0; d < DIM; d++)
 207         {
 208             f[ii][d] += sign * wmass * f_pull[d] * inv_wm;
 209         }
 210     }
 211 }
 212
 213 /* Apply forces in a mass weighted fashion */
 214 static void apply_forces_grp(const pull_group_work_t& pgrp,
 215                              ArrayRef<const real>     masses,
 216                              const dvec               f_pull,
 217                              const int                sign,
 218                              rvec*                    f)
 219 {
 220     auto localAtomIndices = pgrp.atomSet.localIndex();
 221
 222     if (pgrp.params.ind.size() == 1 && pgrp.atomSet.numAtomsLocal() == 1)
 223     {
 224         /* Only one atom and our rank has this atom: we can skip
 225          * the mass weighting, which means that this code also works
 226          * for mass=0, e.g. with a virtual site.
 227          */
 228         for (int d = 0; d < DIM; d++)
 229         {
 230             f[localAtomIndices[0]][d] += sign * f_pull[d];
 231         }
 232     }
 233     else
 234     {
 235         const int numThreads = pgrp.numThreads();
 236         if (numThreads == 1)
 237         {
 238             apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), masses, f_pull, sign, f);
 239         }
 240         else
 241         {
 242 #pragma omp parallel for num_threads(numThreads) schedule(static)
 243             for (int th = 0; th < numThreads; th++)
 244             {
 245                 int ind_start = (localAtomIndices.size() * (th + 0)) / numThreads;
 246                 int ind_end   = (localAtomIndices.size() * (th + 1)) / numThreads;
 247                 apply_forces_grp_part(pgrp, ind_start, ind_end, masses, f_pull, sign, f);
 248             }
 249         }
 250     }
 251 }
 252
 253 /* Apply forces in a mass weighted fashion to a cylinder group */
 254 static void apply_forces_cyl_grp(const pull_group_work_t& pgrp,
 255                                  const double             dv_corr,
 256                                  ArrayRef<const real>     masses,
 257                                  const dvec               f_pull,
 258                                  double                   f_scal,
 259                                  const int                sign,
 260                                  rvec*                    f)
 261 {
 262     const double inv_wm = pgrp.mwscale;
 263
 264     auto localAtomIndices = pgrp.atomSet.localIndex();
 265
 266     /* The cylinder group is always a slab in the system, thus large.
 267      * Therefore we always thread-parallelize this group.
 268      */
 269     const int            numAtomsLocal = localAtomIndices.size();
 270     const int gmx_unused numThreads    = pgrp.numThreads();
 271 #pragma omp parallel for num_threads(numThreads) schedule(static)
 272     for (int i = 0; i < numAtomsLocal; i++)
 273     {
 274         double weight = pgrp.localWeights[i];
 275         if (weight == 0)
 276         {
 277             continue;
 278         }
 279         int    ii   = localAtomIndices[i];
 280         double mass = masses[ii];
 281         /* The stored axial distance from the cylinder center (dv) needs
 282          * to be corrected for an offset (dv_corr), which was unknown when
 283          * we calculated dv.
 284          */
 285         double dv_com = pgrp.dv[i] + dv_corr;
 286
 287         /* Here we not only add the pull force working along vec (f_pull),
 288          * but also a radial component, due to the dependence of the weights
 289          * on the radial distance.
 290          */
 291         for (int m = 0; m < DIM; m++)
 292         {
 293             f[ii][m] += sign * inv_wm * (mass * weight * f_pull[m] + pgrp.mdw[i][m] * dv_com * f_scal);
 294         }
 295     }
 296 }
 297
 298 /* Apply torque forces in a mass weighted fashion to the groups that define
 299  * the pull vector direction for pull coordinate pcrd.
 300  */
 301 static void apply_forces_vec_torque(const pull_coord_work_t&          pcrd,
 302                                     ArrayRef<const pull_group_work_t> pullGroups,
 303                                     ArrayRef<const real>              masses,
 304                                     rvec*                             f)
 305 {
 306     const PullCoordSpatialData& spatialData = pcrd.spatialData;
 307
 308     /* The component inpr along the pull vector is accounted for in the usual
 309      * way. Here we account for the component perpendicular to vec.
 310      */
 311     double inpr = 0;
 312     for (int m = 0; m < DIM; m++)
 313     {
 314         inpr += spatialData.dr01[m] * spatialData.vec[m];
 315     }
 316     /* The torque force works along the component of the distance vector
 317      * of between the two "usual" pull groups that is perpendicular to
 318      * the pull vector. The magnitude of this force is the "usual" scale force
 319      * multiplied with the ratio of the distance between the two "usual" pull
 320      * groups and the distance between the two groups that define the vector.
 321      */
 322     dvec f_perp;
 323     for (int m = 0; m < DIM; m++)
 324     {
 325         f_perp[m] = (spatialData.dr01[m] - inpr * spatialData.vec[m]) / spatialData.vec_len * pcrd.scalarForce;
 326     }
 327
 328     /* Apply the force to the groups defining the vector using opposite signs */
 329     apply_forces_grp(pullGroups[pcrd.params.group[2]], masses, f_perp, -1, f);
 330     apply_forces_grp(pullGroups[pcrd.params.group[3]], masses, f_perp, 1, f);
 331 }
 332
 333 /* Apply forces in a mass weighted fashion */
 334 static void apply_forces_coord(const pull_coord_work_t&          pcrd,
 335                                ArrayRef<const pull_group_work_t> pullGroups,
 336                                const PullCoordVectorForces&      forces,
 337                                ArrayRef<const real>              masses,
 338                                rvec*                             f)
 339 {
 340     /* Here it would be more efficient to use one large thread-parallel
 341      * region instead of potential parallel regions within apply_forces_grp.
 342      * But there could be overlap between pull groups and this would lead
 343      * to data races.
 344      *
 345      * This may also lead to potential issues with force redistribution
 346      * for transformation pull coordinates
 347      */
 348
 349     if (pcrd.params.eGeom == PullGroupGeometry::Cylinder)
 350     {
 351         apply_forces_cyl_grp(
 352                 *pcrd.dynamicGroup0, pcrd.spatialData.cyl_dev, masses, forces.force01, pcrd.scalarForce, -1, f);
 353
 354         /* Sum the force along the vector and the radial force */
 355         dvec f_tot;
 356         for (int m = 0; m < DIM; m++)
 357         {
 358             f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
 359         }
 360         apply_forces_grp(pullGroups[pcrd.params.group[1]], masses, f_tot, 1, f);
 361     }
 362     else if (pcrd.params.eGeom == PullGroupGeometry::Transformation)
 363     {
 364         return;
 365     }
 366     else
 367     {
 368         if (pcrd.params.eGeom == PullGroupGeometry::DirectionRelative)
 369         {
 370             /* We need to apply the torque forces to the pull groups
 371              * that define the pull vector.
 372              */
 373             apply_forces_vec_torque(pcrd, pullGroups, masses, f);
 374         }
 375
 376         if (!pullGroups[pcrd.params.group[0]].params.ind.empty())
 377         {
 378             apply_forces_grp(pullGroups[pcrd.params.group[0]], masses, forces.force01, -1, f);
 379         }
 380         apply_forces_grp(pullGroups[pcrd.params.group[1]], masses, forces.force01, 1, f);
 381
 382         if (pcrd.params.ngroup >= 4)
 383         {
 384             apply_forces_grp(pullGroups[pcrd.params.group[2]], masses, forces.force23, -1, f);
 385             apply_forces_grp(pullGroups[pcrd.params.group[3]], masses, forces.force23, 1, f);
 386         }
 387         if (pcrd.params.ngroup >= 6)
 388         {
 389             apply_forces_grp(pullGroups[pcrd.params.group[4]], masses, forces.force45, -1, f);
 390             apply_forces_grp(pullGroups[pcrd.params.group[5]], masses, forces.force45, 1, f);
 391         }
 392     }
 393 }
 394
 395 real max_pull_distance2(const pull_coord_work_t& pcrd, const t_pbc& pbc)
 396 {
 397     /* Note that this maximum distance calculation is more complex than
 398      * most other cases in GROMACS, since here we have to take care of
 399      * distance calculations that don't involve all three dimensions.
 400      * For example, we can use distances that are larger than the
 401      * box X and Y dimensions for a box that is elongated along Z.
 402      */
 403
 404     real max_d2 = GMX_REAL_MAX;
 405
 406     if (pull_coordinate_is_directional(pcrd.params))
 407     {
 408         /* Directional pulling along along direction pcrd.vec.
 409          * Calculating the exact maximum distance is complex and bug-prone.
 410          * So we take a safe approach by not allowing distances that
 411          * are larger than half the distance between unit cell faces
 412          * along dimensions involved in pcrd.vec.
 413          */
 414         for (int m = 0; m < DIM; m++)
 415         {
 416             if (m < pbc.ndim_ePBC && pcrd.spatialData.vec[m] != 0)
 417             {
 418                 real imageDistance2 = gmx::square(pbc.box[m][m]);
 419                 for (int d = m + 1; d < DIM; d++)
 420                 {
 421                     imageDistance2 -= gmx::square(pbc.box[d][m]);
 422                 }
 423                 max_d2 = std::min(max_d2, imageDistance2);
 424             }
 425         }
 426     }
 427     else
 428     {
 429         /* Distance pulling along dimensions with pcrd.params.dim[d]==1.
 430          * We use half the minimum box vector length of the dimensions involved.
 431          * This is correct for all cases, except for corner cases with
 432          * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
 433          * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
 434          * in such corner cases the user could get correct results,
 435          * depending on the details of the setup, we avoid further
 436          * code complications.
 437          */
 438         for (int m = 0; m < DIM; m++)
 439         {
 440             if (m < pbc.ndim_ePBC && pcrd.params.dim[m] != 0)
 441             {
 442                 real imageDistance2 = gmx::square(pbc.box[m][m]);
 443                 for (int d = 0; d < m; d++)
 444                 {
 445                     if (pcrd.params.dim[d] != 0)
 446                     {
 447                         imageDistance2 += gmx::square(pbc.box[m][d]);
 448                     }
 449                 }
 450                 max_d2 = std::min(max_d2, imageDistance2);
 451             }
 452         }
 453     }
 454
 455     return 0.25 * max_d2;
 456 }
 457
 458 /* This function returns the distance based on coordinates xg and xref.
 459  * Note that the pull coordinate struct pcrd is not modified.
 460  *
 461  * \param[in]  pull  The pull struct
 462  * \param[in]  pcrd  The pull coordinate to compute a distance for
 463  * \param[in]  pbc   The periodic boundary conditions
 464  * \param[in]  xg    The coordinate of group 1
 465  * \param[in]  xref  The coordinate of group 0
 466  * \param[in]  groupIndex0  The index of group 0 in the pcrd->params.group
 467  * \param[in]  groupIndex1  The index of group 1 in the pcrd->params.group
 468  * \param[in]  max_dist2    The maximum distance squared
 469  * \param[out] dr           The distance vector
 470  */
 471 static void low_get_pull_coord_dr(const pull_t&            pull,
 472                                   const pull_coord_work_t& pcrd,
 473                                   const t_pbc&             pbc,
 474                                   const dvec               xg,
 475                                   const dvec               xref,
 476                                   const int                groupIndex0,
 477                                   const int                groupIndex1,
 478                                   const double             max_dist2,
 479                                   dvec                     dr)
 480 {
 481     const pull_group_work_t& pgrp0 = pull.group[pcrd.params.group[0]];
 482
 483     // Group coordinate, to be updated with the reference position
 484     dvec xrefr;
 485
 486     /* Only the first group can be an absolute reference, in that case nat=0 */
 487     if (pgrp0.params.ind.empty())
 488     {
 489         for (int m = 0; m < DIM; m++)
 490         {
 491             xrefr[m] = pcrd.params.origin[m];
 492         }
 493     }
 494     else
 495     {
 496         copy_dvec(xref, xrefr);
 497     }
 498
 499     dvec dref = { 0, 0, 0 };
 500     if (pcrd.params.eGeom == PullGroupGeometry::DirectionPBC)
 501     {
 502         for (int m = 0; m < DIM; m++)
 503         {
 504             dref[m] = pcrd.value_ref * pcrd.spatialData.vec[m];
 505         }
 506         /* Add the reference position, so we use the correct periodic image */
 507         dvec_inc(xrefr, dref);
 508     }
 509
 510     pbc_dx_d(&pbc, xg, xrefr, dr);
 511
 512     bool   directional = pull_coordinate_is_directional(pcrd.params);
 513     double dr2         = 0;
 514     for (int m = 0; m < DIM; m++)
 515     {
 516         dr[m] *= pcrd.params.dim[m];
 517         if (pcrd.params.dim[m] && !(directional && pcrd.spatialData.vec[m] == 0))
 518         {
 519             dr2 += dr[m] * dr[m];
 520         }
 521     }
 522     /* Check if we are close to switching to another periodic image */
 523     if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
 524     {
 525         /* Note that technically there is no issue with switching periodic
 526          * image, as pbc_dx_d returns the distance to the closest periodic
 527          * image. However in all cases where periodic image switches occur,
 528          * the pull results will be useless in practice.
 529          */
 530         gmx_fatal(FARGS,
 531                   "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
 532                   "box size (%f).\n%s",
 533                   pcrd.params.group[groupIndex0],
 534                   pcrd.params.group[groupIndex1],
 535                   sqrt(dr2),
 536                   sqrt(0.98 * 0.98 * max_dist2),
 537                   pcrd.params.eGeom == PullGroupGeometry::Direction
 538                           ? "You might want to consider using \"pull-geometry = "
 539                             "direction-periodic\" instead.\n"
 540                           : "");
 541     }
 542
 543     if (pcrd.params.eGeom == PullGroupGeometry::DirectionPBC)
 544     {
 545         dvec_inc(dr, dref);
 546     }
 547 }
 548
 549 /* This function returns the distance based on the contents of the pull struct.
 550  * pull->coord[coord_ind].dr, and potentially vector, are updated.
 551  */
 552 static void get_pull_coord_dr(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc& pbc)
 553 {
 554     const int coord_ind = pcrd->params.coordIndex;
 555
 556     PullCoordSpatialData& spatialData = pcrd->spatialData;
 557
 558     double md2;
 559     /* With AWH pulling we allow for periodic pulling with geometry=direction.
 560      * TODO: Store a periodicity flag instead of checking for external pull provider.
 561      */
 562     if (pcrd->params.eGeom == PullGroupGeometry::DirectionPBC
 563         || (pcrd->params.eGeom == PullGroupGeometry::Direction
 564             && pcrd->params.eType == PullingAlgorithm::External))
 565     {
 566         md2 = -1;
 567     }
 568     else
 569     {
 570         md2 = static_cast<double>(max_pull_distance2(*pcrd, pbc));
 571     }
 572
 573     if (pcrd->params.eGeom == PullGroupGeometry::DirectionRelative)
 574     {
 575         /* We need to determine the pull vector */
 576         dvec vec;
 577         int  m;
 578
 579         const pull_group_work_t& pgrp2 = pull.group[pcrd->params.group[2]];
 580         const pull_group_work_t& pgrp3 = pull.group[pcrd->params.group[3]];
 581
 582         pbc_dx_d(&pbc, pgrp3.x, pgrp2.x, vec);
 583
 584         for (m = 0; m < DIM; m++)
 585         {
 586             vec[m] *= pcrd->params.dim[m];
 587         }
 588         spatialData.vec_len = dnorm(vec);
 589         for (m = 0; m < DIM; m++)
 590         {
 591             spatialData.vec[m] = vec[m] / spatialData.vec_len;
 592         }
 593         if (debug)
 594         {
 595             fprintf(debug,
 596                     "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
 597                     coord_ind,
 598                     vec[XX],
 599                     vec[YY],
 600                     vec[ZZ],
 601                     spatialData.vec[XX],
 602                     spatialData.vec[YY],
 603                     spatialData.vec[ZZ]);
 604         }
 605     }
 606
 607     const pull_group_work_t& pgrp0 = pull.group[pcrd->params.group[0]];
 608     const pull_group_work_t& pgrp1 = pull.group[pcrd->params.group[1]];
 609
 610     low_get_pull_coord_dr(
 611             pull,
 612             *pcrd,
 613             pbc,
 614             pgrp1.x,
 615             pcrd->params.eGeom == PullGroupGeometry::Cylinder ? pcrd->dynamicGroup0->x : pgrp0.x,
 616             0,
 617             1,
 618             md2,
 619             spatialData.dr01);
 620
 621     if (pcrd->params.ngroup >= 4)
 622     {
 623         const pull_group_work_t& pgrp2 = pull.group[pcrd->params.group[2]];
 624         const pull_group_work_t& pgrp3 = pull.group[pcrd->params.group[3]];
 625
 626         low_get_pull_coord_dr(pull, *pcrd, pbc, pgrp3.x, pgrp2.x, 2, 3, md2, spatialData.dr23);
 627     }
 628     if (pcrd->params.ngroup >= 6)
 629     {
 630         const pull_group_work_t& pgrp4 = pull.group[pcrd->params.group[4]];
 631         const pull_group_work_t& pgrp5 = pull.group[pcrd->params.group[5]];
 632
 633         low_get_pull_coord_dr(pull, *pcrd, pbc, pgrp5.x, pgrp4.x, 4, 5, md2, spatialData.dr45);
 634     }
 635 }
 636
 637 /* Modify x so that it is periodic in [-pi, pi)
 638  * It is assumed that x is in [-3pi, 3pi) so that x
 639  * needs to be shifted by at most one period.
 640  */
 641 static void make_periodic_2pi(double* x)
 642 {
 643     if (*x >= M_PI)
 644     {
 645         *x -= M_2PI;
 646     }
 647     else if (*x < -M_PI)
 648     {
 649         *x += M_2PI;
 650     }
 651 }
 652
 653 /* This function should always be used to get values for setting value_ref in pull_coord_work_t */
 654 static double sanitizePullCoordReferenceValue(const t_pull_coord& pcrdParams, double value_ref)
 655 {
 656     GMX_ASSERT(pcrdParams.eType != PullingAlgorithm::External,
 657                "The pull coord reference value should not be used with type external-potential");
 658
 659     if (pcrdParams.eGeom == PullGroupGeometry::Distance)
 660     {
 661         if (value_ref < 0)
 662         {
 663             gmx_fatal(FARGS,
 664                       "Pull reference distance for coordinate %d (%f) needs to be non-negative",
 665                       pcrdParams.coordIndex + 1,
 666                       value_ref);
 667         }
 668     }
 669     else if (pcrdParams.eGeom == PullGroupGeometry::Angle || pcrdParams.eGeom == PullGroupGeometry::AngleAxis)
 670     {
 671         if (value_ref < 0 || value_ref > M_PI)
 672         {
 673             gmx_fatal(FARGS,
 674                       "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
 675                       "interval [0,180] deg",
 676                       pcrdParams.coordIndex + 1,
 677                       value_ref * pull_conversion_factor_internal2userinput(pcrdParams));
 678         }
 679     }
 680     else if (pcrdParams.eGeom == PullGroupGeometry::Dihedral)
 681     {
 682         /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
 683         make_periodic_2pi(&value_ref);
 684     }
 685
 686     return value_ref;
 687 }
 688
 689 static void updatePullCoordReferenceValue(double* referenceValue, const t_pull_coord& pcrdParams, double t)
 690 {
 691     GMX_ASSERT(referenceValue, "Need a valid reference value object");
 692
 693     /* With zero rate the reference value is set initially and doesn't change */
 694     if (pcrdParams.rate != 0)
 695     {
 696         const double inputValue = (pcrdParams.init + pcrdParams.rate * t)
 697                                   * pull_conversion_factor_userinput2internal(pcrdParams);
 698         *referenceValue = sanitizePullCoordReferenceValue(pcrdParams, inputValue);
 699     }
 700 }
 701
 702 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
 703 static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
 704 {
 705     double phi, sign;
 706     dvec   dr32; /* store instead of dr23? */
 707
 708     dsvmul(-1, spatialData->dr23, dr32);
 709     dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
 710     dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
 711     phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
 712
 713     /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
 714      * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
 715      * Note 2: the angle between the plane normal vectors equals pi only when
 716      * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
 717      * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
 718      */
 719     sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
 720     return sign * phi;
 721 }
 722
 723 /* Calculates pull->coord[coord_ind].value.
 724  * This function also updates pull->coord[coord_ind].dr.
 725  */
 726 static void get_pull_coord_distance(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc& pbc)
 727 {
 728     get_pull_coord_dr(pull, pcrd, pbc);
 729
 730     PullCoordSpatialData& spatialData = pcrd->spatialData;
 731
 732     switch (pcrd->params.eGeom)
 733     {
 734         case PullGroupGeometry::Distance:
 735             /* Pull along the vector between the com's */
 736             spatialData.value = dnorm(spatialData.dr01);
 737             break;
 738         case PullGroupGeometry::Direction:
 739         case PullGroupGeometry::DirectionPBC:
 740         case PullGroupGeometry::DirectionRelative:
 741         case PullGroupGeometry::Cylinder:
 742             /* Pull along vec */
 743             spatialData.value = 0;
 744             for (int m = 0; m < DIM; m++)
 745             {
 746                 spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
 747             }
 748             break;
 749         case PullGroupGeometry::Angle:
 750             spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
 751             break;
 752         case PullGroupGeometry::Dihedral:
 753             spatialData.value = get_dihedral_angle_coord(&spatialData);
 754             break;
 755         case PullGroupGeometry::AngleAxis:
 756             spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
 757             break;
 758         case PullGroupGeometry::Transformation:
 759             // Note that we would only need to pass the part of coord up to coord_ind
 760             spatialData.value = gmx::getTransformationPullCoordinateValue(
 761                     pcrd, ArrayRef<const pull_coord_work_t>(pull.coord).subArray(0, pcrd->params.coordIndex));
 762             break;
 763         default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
 764     }
 765 }
 766
 767 /* Returns the deviation from the reference value.
 768  * Updates pull->coord[coord_ind].dr, .value and .value_ref.
 769  */
 770 static double get_pull_coord_deviation(const pull_t& pull, pull_coord_work_t* pcrd, const t_pbc& pbc, double t)
 771 {
 772     double dev = 0;
 773
 774     /* Update the reference value before computing the distance,
 775      * since it is used in the distance computation with periodic pulling.
 776      */
 777     updatePullCoordReferenceValue(&pcrd->value_ref, pcrd->params, t);
 778
 779     get_pull_coord_distance(pull, pcrd, pbc);
 780
 781     get_pull_coord_distance(pull, pcrd, pbc);
 782
 783     /* Determine the deviation */
 784     dev = pcrd->spatialData.value - pcrd->value_ref;
 785
 786     /* Check that values are allowed */
 787     if (pcrd->params.eGeom == PullGroupGeometry::Distance && pcrd->spatialData.value == 0)
 788     {
 789         /* With no vector we can not determine the direction for the force,
 790          * so we set the force to zero.
 791          */
 792         dev = 0;
 793     }
 794     else if (pcrd->params.eGeom == PullGroupGeometry::Dihedral)
 795     {
 796         /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
 797            Thus, the unwrapped deviation is here in (-2pi, 2pi].
 798            After making it periodic, the deviation will be in [-pi, pi). */
 799         make_periodic_2pi(&dev);
 800     }
 801
 802     return dev;
 803 }
 804
 805 double get_pull_coord_value(pull_t* pull, int coordIndex, const t_pbc& pbc)
 806 {
 807     get_pull_coord_distance(*pull, &pull->coord[coordIndex], pbc);
 808
 809     return pull->coord[coordIndex].spatialData.value;
 810 }
 811
 812 void clear_pull_forces(pull_t* pull)
 813 {
 814     /* Zeroing the forces is only required for constraint pulling.
 815      * It can happen that multiple constraint steps need to be applied
 816      * and therefore the constraint forces need to be accumulated.
 817      */
 818     for (pull_coord_work_t& coord : pull->coord)
 819     {
 820         coord.scalarForce = 0;
 821     }
 822 }
 823
 824 /* Apply constraint using SHAKE */
 825 static void do_constraint(struct pull_t* pull,
 826                           const t_pbc&   pbc,
 827                           ArrayRef<RVec> x,
 828                           ArrayRef<RVec> v,
 829                           gmx_bool       bMaster,
 830                           tensor         vir,
 831                           double         dt,
 832                           double         t)
 833 {
 834
 835     dvec*    r_ij;   /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
 836     dvec     unc_ij; /* xp[i] com of i this step, before constr.   -> unc_ij  */
 837     dvec*    rnew;   /* current 'new' positions of the groups */
 838     double*  dr_tot; /* the total update of the coords */
 839     dvec     vec;
 840     double   inpr;
 841     double   lambda, rm, invdt = 0;
 842     gmx_bool bConverged_all, bConverged = FALSE;
 843     int      niter = 0, ii, j, m, max_iter = 100;
 844     double   a;
 845     dvec     tmp, tmp3;
 846
 847     snew(r_ij, pull->coord.size());
 848     snew(dr_tot, pull->coord.size());
 849
 850     snew(rnew, pull->group.size());
 851
 852     /* copy the current unconstrained positions for use in iterations. We
 853        iterate until rinew[i] and rjnew[j] obey the constraints. Then
 854        rinew - pull.x_unc[i] is the correction dr to group i */
 855     for (size_t g = 0; g < pull->group.size(); g++)
 856     {
 857         copy_dvec(pull->group[g].xp, rnew[g]);
 858     }
 859
 860     /* Determine the constraint directions from the old positions */
 861     for (size_t c = 0; c < pull->coord.size(); c++)
 862     {
 863         pull_coord_work_t* pcrd;
 864
 865         pcrd = &pull->coord[c];
 866
 867         if (pcrd->params.eType != PullingAlgorithm::Constraint)
 868         {
 869             continue;
 870         }
 871
 872         /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
 873          * We don't modify dr and value anymore, so these values are also used
 874          * for printing.
 875          */
 876         get_pull_coord_distance(*pull, pcrd, pbc);
 877
 878         const PullCoordSpatialData& spatialData = pcrd->spatialData;
 879         if (debug)
 880         {
 881             fprintf(debug,
 882                     "Pull coord %zu dr %f %f %f\n",
 883                     c,
 884                     spatialData.dr01[XX],
 885                     spatialData.dr01[YY],
 886                     spatialData.dr01[ZZ]);
 887         }
 888
 889         if (pcrd->params.eGeom == PullGroupGeometry::Direction
 890             || pcrd->params.eGeom == PullGroupGeometry::DirectionPBC)
 891         {
 892             /* Select the component along vec */
 893             a = 0;
 894             for (m = 0; m < DIM; m++)
 895             {
 896                 a += spatialData.vec[m] * spatialData.dr01[m];
 897             }
 898             for (m = 0; m < DIM; m++)
 899             {
 900                 r_ij[c][m] = a * spatialData.vec[m];
 901             }
 902         }
 903         else
 904         {
 905             copy_dvec(spatialData.dr01, r_ij[c]);
 906         }
 907
 908         if (dnorm2(r_ij[c]) == 0)
 909         {
 910             gmx_fatal(FARGS,
 911                       "Distance for pull coordinate %zu is zero with constraint pulling, which is "
 912                       "not allowed.",
 913                       c + 1);
 914         }
 915     }
 916
 917     bConverged_all = FALSE;
 918     while (!bConverged_all && niter < max_iter)
 919     {
 920         bConverged_all = TRUE;
 921
 922         /* loop over all constraints */
 923         for (size_t c = 0; c < pull->coord.size(); c++)
 924         {
 925             pull_coord_work_t* pcrd;
 926             pull_group_work_t *pgrp0, *pgrp1;
 927             dvec               dr0, dr1;
 928
 929             pcrd = &pull->coord[c];
 930
 931             if (pcrd->params.eType != PullingAlgorithm::Constraint)
 932             {
 933                 continue;
 934             }
 935
 936             updatePullCoordReferenceValue(&pcrd->value_ref, pcrd->params, t);
 937
 938             pgrp0 = &pull->group[pcrd->params.group[0]];
 939             pgrp1 = &pull->group[pcrd->params.group[1]];
 940
 941             /* Get the current difference vector */
 942             low_get_pull_coord_dr(
 943                     *pull, *pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
 944
 945             if (debug)
 946             {
 947                 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
 948             }
 949
 950             rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
 951
 952             switch (pcrd->params.eGeom)
 953             {
 954                 case PullGroupGeometry::Distance:
 955                     if (pcrd->value_ref <= 0)
 956                     {
 957                         gmx_fatal(
 958                                 FARGS,
 959                                 "The pull constraint reference distance for group %zu is <= 0 (%f)",
 960                                 c,
 961                                 pcrd->value_ref);
 962                     }
 963
 964                     {
 965                         double q, c_a, c_b, c_c;
 966
 967                         c_a = diprod(r_ij[c], r_ij[c]);
 968                         c_b = diprod(unc_ij, r_ij[c]) * 2;
 969                         c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
 970
 971                         if (c_b < 0)
 972                         {
 973                             q      = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
 974                             lambda = -q / c_a;
 975                         }
 976                         else
 977                         {
 978                             q      = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
 979                             lambda = -c_c / q;
 980                         }
 981
 982                         if (debug)
 983                         {
 984                             fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b, c_c, lambda);
 985                         }
 986                     }
 987
 988                     /* The position corrections dr due to the constraints */
 989                     dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
 990                     dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
 991                     dr_tot[c] += -lambda * dnorm(r_ij[c]);
 992                     break;
 993                 case PullGroupGeometry::Direction:
 994                 case PullGroupGeometry::DirectionPBC:
 995                 case PullGroupGeometry::Cylinder:
 996                     /* A 1-dimensional constraint along a vector */
 997                     a = 0;
 998                     for (m = 0; m < DIM; m++)
 999                     {
1000                         vec[m] = pcrd->spatialData.vec[m];
1001                         a += unc_ij[m] * vec[m];
1002                     }
1003                     /* Select only the component along the vector */
1004                     dsvmul(a, vec, unc_ij);
1005                     lambda = a - pcrd->value_ref;
1006                     if (debug)
1007                     {
1008                         fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
1009                     }
1010
1011                     /* The position corrections dr due to the constraints */
1012                     dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
1013                     dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
1014                     dr_tot[c] += -lambda;
1015                     break;
1016                 case PullGroupGeometry::Transformation:
1017                     GMX_RELEASE_ASSERT(false, "transformation with constraints should never occur");
1018                     break;
1019                 default: gmx_incons("Invalid enumeration value for eGeom");
1020             }
1021
1022             /* DEBUG */
1023             if (debug)
1024             {
1025                 int g0, g1;
1026
1027                 g0 = pcrd->params.group[0];
1028                 g1 = pcrd->params.group[1];
1029                 low_get_pull_coord_dr(*pull, *pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
1030                 low_get_pull_coord_dr(*pull, *pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
1031                 fprintf(debug,
1032                         "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1033                         rnew[g0][0],
1034                         rnew[g0][1],
1035                         rnew[g0][2],
1036                         rnew[g1][0],
1037                         rnew[g1][1],
1038                         rnew[g1][2],
1039                         dnorm(tmp));
1040                 fprintf(debug,
1041                         "Pull ref %8s %8s %8s   %8s %8s %8s d: %8.5f\n",
1042                         "",
1043                         "",
1044                         "",
1045                         "",
1046                         "",
1047                         "",
1048                         pcrd->value_ref);
1049                 fprintf(debug,
1050                         "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1051                         dr0[0],
1052                         dr0[1],
1053                         dr0[2],
1054                         dr1[0],
1055                         dr1[1],
1056                         dr1[2],
1057                         dnorm(tmp3));
1058             } /* END DEBUG */
1059
1060             /* Update the COMs with dr */
1061             dvec_inc(rnew[pcrd->params.group[1]], dr1);
1062             dvec_inc(rnew[pcrd->params.group[0]], dr0);
1063         }
1064
1065         /* Check if all constraints are fullfilled now */
1066         for (pull_coord_work_t& coord : pull->coord)
1067         {
1068             if (coord.params.eType != PullingAlgorithm::Constraint)
1069             {
1070                 continue;
1071             }
1072
1073             low_get_pull_coord_dr(
1074                     *pull, coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
1075
1076             switch (coord.params.eGeom)
1077             {
1078                 case PullGroupGeometry::Distance:
1079                     bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1080                     break;
1081                 case PullGroupGeometry::Direction:
1082                 case PullGroupGeometry::DirectionPBC:
1083                 case PullGroupGeometry::Cylinder:
1084                     for (m = 0; m < DIM; m++)
1085                     {
1086                         vec[m] = coord.spatialData.vec[m];
1087                     }
1088                     inpr = diprod(unc_ij, vec);
1089                     dsvmul(inpr, vec, unc_ij);
1090                     bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1091                     break;
1092                 default:
1093                     GMX_ASSERT(false,
1094                                gmx::formatString("Geometry %s not handled here",
1095                                                  enumValueToString(coord.params.eGeom))
1096                                        .c_str());
1097             }
1098
1099             if (!bConverged)
1100             {
1101                 if (debug)
1102                 {
1103                     fprintf(debug,
1104                             "Pull constraint not converged: "
1105                             "groups %d %d,"
1106                             "d_ref = %f, current d = %f\n",
1107                             coord.params.group[0],
1108                             coord.params.group[1],
1109                             coord.value_ref,
1110                             dnorm(unc_ij));
1111                 }
1112
1113                 bConverged_all = FALSE;
1114             }
1115         }
1116
1117         niter++;
1118         /* if after all constraints are dealt with and bConverged is still TRUE
1119            we're finished, if not we do another iteration */
1120     }
1121     if (niter > max_iter)
1122     {
1123         gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1124     }
1125
1126     /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1127
1128     if (!v.empty())
1129     {
1130         invdt = 1 / dt;
1131     }
1132
1133     /* update atoms in the groups */
1134     for (size_t g = 0; g < pull->group.size(); g++)
1135     {
1136         const pull_group_work_t* pgrp;
1137         dvec                     dr;
1138
1139         pgrp = &pull->group[g];
1140
1141         /* get the final constraint displacement dr for group g */
1142         dvec_sub(rnew[g], pgrp->xp, dr);
1143
1144         if (dnorm2(dr) == 0)
1145         {
1146             /* No displacement, no update necessary */
1147             continue;
1148         }
1149
1150         /* update the atom positions */
1151         auto localAtomIndices = pgrp->atomSet.localIndex();
1152         copy_dvec(dr, tmp);
1153         for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
1154         {
1155             ii = localAtomIndices[j];
1156             if (!pgrp->localWeights.empty())
1157             {
1158                 dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
1159             }
1160             for (m = 0; m < DIM; m++)
1161             {
1162                 x[ii][m] += tmp[m];
1163             }
1164             if (!v.empty())
1165             {
1166                 for (m = 0; m < DIM; m++)
1167                 {
1168                     v[ii][m] += invdt * tmp[m];
1169                 }
1170             }
1171         }
1172     }
1173
1174     /* calculate the constraint forces, used for output and virial only */
1175     for (size_t c = 0; c < pull->coord.size(); c++)
1176     {
1177         pull_coord_work_t* pcrd;
1178
1179         pcrd = &pull->coord[c];
1180
1181         if (pcrd->params.eType != PullingAlgorithm::Constraint)
1182         {
1183             continue;
1184         }
1185
1186         /* Accumulate the forces, in case we have multiple constraint steps */
1187         double force =
1188                 dr_tot[c]
1189                 / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
1190                    * dt * dt);
1191         pcrd->scalarForce += force;
1192
1193         if (vir != nullptr && pcrd->params.eGeom != PullGroupGeometry::DirectionPBC && bMaster)
1194         {
1195             double f_invr;
1196
1197             /* Add the pull contribution to the virial */
1198             /* We have already checked above that r_ij[c] != 0 */
1199             f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
1200
1201             for (j = 0; j < DIM; j++)
1202             {
1203                 for (m = 0; m < DIM; m++)
1204                 {
1205                     vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
1206                 }
1207             }
1208         }
1209     }
1210
1211     /* finished! I hope. Give back some memory */
1212     sfree(r_ij);
1213     sfree(dr_tot);
1214     sfree(rnew);
1215 }
1216
1217 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1218 {
1219     for (int j = 0; j < DIM; j++)
1220     {
1221         for (int m = 0; m < DIM; m++)
1222         {
1223             vir[j][m] -= 0.5 * f[j] * dr[m];
1224         }
1225     }
1226 }
1227
1228 /* Adds the pull contribution to the virial */
1229 static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
1230 {
1231     if (vir != nullptr && pcrd.params.eGeom != PullGroupGeometry::DirectionPBC)
1232     {
1233         /* Add the pull contribution for each distance vector to the virial. */
1234         add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1235         if (pcrd.params.ngroup >= 4)
1236         {
1237             add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1238         }
1239         if (pcrd.params.ngroup >= 6)
1240         {
1241             add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1242         }
1243     }
1244 }
1245
1246 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
1247                                                        double             dev,
1248                                                        real               lambda,
1249                                                        real*              V,
1250                                                        real*              dVdl)
1251 {
1252     real k, dkdl;
1253
1254     k    = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
1255     dkdl = pcrd->params.kB - pcrd->params.k;
1256
1257     switch (pcrd->params.eType)
1258     {
1259         case PullingAlgorithm::Umbrella:
1260         case PullingAlgorithm::FlatBottom:
1261         case PullingAlgorithm::FlatBottomHigh:
1262             /* The only difference between an umbrella and a flat-bottom
1263              * potential is that a flat-bottom is zero above or below
1264                the reference value.
1265              */
1266             if ((pcrd->params.eType == PullingAlgorithm::FlatBottom && dev < 0)
1267                 || (pcrd->params.eType == PullingAlgorithm::FlatBottomHigh && dev > 0))
1268             {
1269                 dev = 0;
1270             }
1271
1272             pcrd->scalarForce = -k * dev;
1273             *V += 0.5 * k * gmx::square(dev);
1274             *dVdl += 0.5 * dkdl * gmx::square(dev);
1275             break;
1276         case PullingAlgorithm::ConstantForce:
1277             pcrd->scalarForce = -k;
1278             *V += k * pcrd->spatialData.value;
1279             *dVdl += dkdl * pcrd->spatialData.value;
1280             break;
1281         case PullingAlgorithm::External:
1282             gmx_incons(
1283                     "the scalar pull force should not be calculated internally for pull type "
1284                     "external");
1285         default: gmx_incons("Unsupported pull type in do_pull_pot");
1286     }
1287 }
1288
1289 static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
1290 {
1291     const t_pull_coord&         params      = pcrd.params;
1292     const PullCoordSpatialData& spatialData = pcrd.spatialData;
1293
1294     /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1295     PullCoordVectorForces forces;
1296
1297     if (params.eGeom == PullGroupGeometry::Distance)
1298     {
1299         double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
1300         for (int m = 0; m < DIM; m++)
1301         {
1302             forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
1303         }
1304     }
1305     else if (params.eGeom == PullGroupGeometry::Angle)
1306     {
1307
1308         double cos_theta, cos_theta2;
1309
1310         cos_theta  = cos(spatialData.value);
1311         cos_theta2 = gmx::square(cos_theta);
1312
1313         /* The force at theta = 0, pi is undefined so we should not apply any force.
1314          * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1315          * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1316          * have to check for this before dividing by their norm below.
1317          */
1318         if (cos_theta2 < 1)
1319         {
1320             /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1321              * and another vector parallel to dr23:
1322              * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1323              * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1324              */
1325             double a       = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1326             double b       = a * cos_theta;
1327             double invdr01 = 1. / dnorm(spatialData.dr01);
1328             double invdr23 = 1. / dnorm(spatialData.dr23);
1329             dvec   normalized_dr01, normalized_dr23;
1330             dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1331             dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1332
1333             for (int m = 0; m < DIM; m++)
1334             {
1335                 /* Here, f_scal is -dV/dtheta */
1336                 forces.force01[m] =
1337                         pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
1338                 forces.force23[m] =
1339                         pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
1340             }
1341         }
1342         else
1343         {
1344             /* No forces to apply for ill-defined cases*/
1345             clear_dvec(forces.force01);
1346             clear_dvec(forces.force23);
1347         }
1348     }
1349     else if (params.eGeom == PullGroupGeometry::AngleAxis)
1350     {
1351         double cos_theta, cos_theta2;
1352
1353         /* The angle-axis force is exactly the same as the angle force (above) except that in
1354            this case the second vector (dr23) is replaced by the pull vector. */
1355         cos_theta  = cos(spatialData.value);
1356         cos_theta2 = gmx::square(cos_theta);
1357
1358         if (cos_theta2 < 1)
1359         {
1360             double a, b;
1361             double invdr01;
1362             dvec   normalized_dr01;
1363
1364             invdr01 = 1. / dnorm(spatialData.dr01);
1365             dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1366             a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1367             b = a * cos_theta;
1368
1369             for (int m = 0; m < DIM; m++)
1370             {
1371                 forces.force01[m] =
1372                         pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
1373             }
1374         }
1375         else
1376         {
1377             clear_dvec(forces.force01);
1378         }
1379     }
1380     else if (params.eGeom == PullGroupGeometry::Dihedral)
1381     {
1382         double m2, n2, tol, sqrdist_32;
1383         dvec   dr32;
1384         /* Note: there is a small difference here compared to the
1385            dihedral force calculations in the bondeds (ref: Bekker 1994).
1386            There rij = ri - rj, while here dr01 = r1 - r0.
1387            However, all distance vectors occur in form of cross or inner products
1388            so that two signs cancel and we end up with the same expressions.
1389            Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1390          */
1391         m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1392         n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1393         dsvmul(-1, spatialData.dr23, dr32);
1394         sqrdist_32 = diprod(dr32, dr32);
1395         tol        = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
1396         if ((m2 > tol) && (n2 > tol))
1397         {
1398             double a_01, a_23_01, a_23_45, a_45;
1399             double inv_dist_32, inv_sqrdist_32, dist_32;
1400             dvec   u, v;
1401             inv_dist_32    = gmx::invsqrt(sqrdist_32);
1402             inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
1403             dist_32        = sqrdist_32 * inv_dist_32;
1404
1405             /* Forces on groups 0, 1 */
1406             a_01 = pcrd.scalarForce * dist_32 / m2;                /* scalarForce is -dV/dphi */
1407             dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1408
1409             /* Forces on groups 4, 5 */
1410             a_45 = -pcrd.scalarForce * dist_32 / n2;
1411             dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1412
1413             /* Force on groups 2, 3 (defining the axis) */
1414             a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
1415             a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
1416             dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1417             dsvmul(a_23_45, forces.force45, v);
1418             dvec_sub(u, v, forces.force23); /* force on group 3 */
1419         }
1420         else
1421         {
1422             /* No force to apply for ill-defined cases */
1423             clear_dvec(forces.force01);
1424             clear_dvec(forces.force23);
1425             clear_dvec(forces.force45);
1426         }
1427     }
1428     else
1429     {
1430         for (int m = 0; m < DIM; m++)
1431         {
1432             forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
1433         }
1434     }
1435
1436     return forces;
1437 }
1438
1439
1440 /* We use a global mutex for locking access to the pull data structure
1441  * during registration of external pull potential providers.
1442  * We could use a different, local mutex for each pull object, but the overhead
1443  * is extremely small here and registration is only done during initialization.
1444  */
1445 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
1446 static std::mutex registrationMutex;
1447
1448 using Lock = std::lock_guard<std::mutex>;
1449
1450 void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
1451 {
1452     GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1453     GMX_RELEASE_ASSERT(provider != nullptr,
1454                        "register_external_pull_potential called with NULL as provider name");
1455
1456     if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
1457     {
1458         gmx_fatal(FARGS,
1459                   "Module '%s' attempted to register an external potential for pull coordinate %d "
1460                   "which is out of the pull coordinate range %d - %zu\n",
1461                   provider,
1462                   coord_index + 1,
1463                   1,
1464                   pull->coord.size());
1465     }
1466
1467     pull_coord_work_t* pcrd = &pull->coord[coord_index];
1468
1469     if (pcrd->params.eType != PullingAlgorithm::External)
1470     {
1471         gmx_fatal(
1472                 FARGS,
1473                 "Module '%s' attempted to register an external potential for pull coordinate %d "
1474                 "which of type '%s', whereas external potentials are only supported with type '%s'",
1475                 provider,
1476                 coord_index + 1,
1477                 enumValueToString(pcrd->params.eType),
1478                 enumValueToString(PullingAlgorithm::External));
1479     }
1480
1481     GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
1482                        "The external potential provider string for a pull coordinate is NULL");
1483
1484     if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
1485     {
1486         gmx_fatal(FARGS,
1487                   "Module '%s' attempted to register an external potential for pull coordinate %d "
1488                   "which expects the external potential to be provided by a module named '%s'",
1489                   provider,
1490                   coord_index + 1,
1491                   pcrd->params.externalPotentialProvider.c_str());
1492     }
1493
1494     /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1495      * pcrd->bExternalPotentialProviderHasBeenRegistered and
1496      * pull->numUnregisteredExternalPotentials.
1497      */
1498     Lock registrationLock(registrationMutex);
1499
1500     if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1501     {
1502         gmx_fatal(FARGS,
1503                   "Module '%s' attempted to register an external potential for pull coordinate %d "
1504                   "more than once",
1505                   provider,
1506                   coord_index + 1);
1507     }
1508
1509     pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1510     pull->numUnregisteredExternalPotentials--;
1511
1512     GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
1513                        "Negative unregistered potentials, the pull code is inconsistent");
1514 }
1515
1516
1517 static void check_external_potential_registration(const struct pull_t* pull)
1518 {
1519     if (pull->numUnregisteredExternalPotentials > 0)
1520     {
1521         for (const pull_coord_work_t& pcrd : pull->coord)
1522         {
1523             if (pcrd.params.eType == PullingAlgorithm::External
1524                 && !pcrd.bExternalPotentialProviderHasBeenRegistered)
1525             {
1526                 gmx_fatal(FARGS,
1527                           "No external provider for external pull potentials have been provided "
1528                           "for %d "
1529                           "pull coordinates. The first coordinate without provider is number %d, "
1530                           "which "
1531                           "expects a module named '%s' to provide the external potential.",
1532                           pull->numUnregisteredExternalPotentials,
1533                           pcrd.params.coordIndex + 1,
1534                           pcrd.params.externalPotentialProvider.c_str());
1535             }
1536         }
1537     }
1538 }
1539
1540 /*! \brief Applies a force of any non-transformation pull coordinate
1541  *
1542  * \param[in] pull         The pull struct, we can't pass only the coord because of cylinder pulling
1543  * \param[in] coord_index  The index of the coord to apply to force for
1544  * \param[in] coord_force  The force working on this coord
1545  * \param[in] masses       Atom masses
1546  * \param[in] forceWithVirial  Atom force and virial object
1547  */
1548 static void apply_default_pull_coord_force(pull_t*                   pull,
1549                                            const int                 coord_index,
1550                                            const double              coord_force,
1551                                            gmx::ArrayRef<const real> masses,
1552                                            gmx::ForceWithVirial*     forceWithVirial)
1553 {
1554     pull_coord_work_t* pcrd = &pull->coord[coord_index];
1555     /* Set the force */
1556     pcrd->scalarForce = coord_force;
1557     /* Calculate the forces on the pull groups */
1558     PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1559
1560     /* Add the forces for this coordinate to the total virial and force */
1561     if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1562     {
1563         matrix virial = { { 0 } };
1564         add_virial_coord(virial, *pcrd, pullCoordForces);
1565         forceWithVirial->addVirialContribution(virial);
1566     }
1567     apply_forces_coord(pull->coord[coord_index],
1568                        pull->group,
1569                        pullCoordForces,
1570                        masses,
1571                        as_rvec_array(forceWithVirial->force_.data()));
1572 }
1573
1574 /* Pull takes care of adding the forces of the external potential.
1575  * The external potential module  has to make sure that the corresponding
1576  * potential energy is added either to the pull term or to a term
1577  * specific to the external module.
1578  */
1579 void apply_external_pull_coord_force(struct pull_t*        pull,
1580                                      int                   coord_index,
1581                                      double                coord_force,
1582                                      ArrayRef<const real>  masses,
1583                                      gmx::ForceWithVirial* forceWithVirial)
1584 {
1585     pull_coord_work_t* pcrd;
1586
1587     GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
1588                "apply_external_pull_coord_force called with coord_index out of range");
1589
1590     if (pull->comm.bParticipate)
1591     {
1592         pcrd = &pull->coord[coord_index];
1593
1594         GMX_RELEASE_ASSERT(
1595                 pcrd->params.eType == PullingAlgorithm::External,
1596                 "The pull force can only be set externally on pull coordinates of external type");
1597
1598         GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
1599                    "apply_external_pull_coord_force called for an unregistered pull coordinate");
1600
1601         if (pcrd->params.eGeom == PullGroupGeometry::Transformation)
1602         {
1603             gmx::applyTransformationPullCoordForce(
1604                     pcrd,
1605                     gmx::ArrayRef<pull_coord_work_t>(pull->coord).subArray(0, pcrd->params.coordIndex),
1606                     coord_force);
1607         }
1608         else
1609         {
1610             apply_default_pull_coord_force(pull, coord_index, coord_force, masses, forceWithVirial);
1611         }
1612     }
1613     pull->numExternalPotentialsStillToBeAppliedThisStep--;
1614 }
1615
1616 /* Calculate the pull potential and scalar force for a pull coordinate.
1617  * Returns the vector forces for the pull coordinate.
1618  */
1619 static PullCoordVectorForces do_pull_pot_coord(const pull_t&      pull,
1620                                                pull_coord_work_t* pcrd,
1621                                                const t_pbc&       pbc,
1622                                                double             t,
1623                                                real               lambda,
1624                                                real*              V,
1625                                                tensor             vir,
1626                                                real*              dVdl)
1627 {
1628     assert(pcrd->params.eType != PullingAlgorithm::Constraint);
1629
1630     double dev = get_pull_coord_deviation(pull, pcrd, pbc, t);
1631
1632     calc_pull_coord_scalar_force_and_potential(pcrd, dev, lambda, V, dVdl);
1633
1634     PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1635
1636     add_virial_coord(vir, *pcrd, pullCoordForces);
1637
1638     return pullCoordForces;
1639 }
1640
1641 real pull_potential(struct pull_t*        pull,
1642                     ArrayRef<const real>  masses,
1643                     const t_pbc&          pbc,
1644                     const t_commrec*      cr,
1645                     const double          t,
1646                     const real            lambda,
1647                     ArrayRef<const RVec>  x,
1648                     gmx::ForceWithVirial* force,
1649                     real*                 dvdlambda)
1650 {
1651     real V = 0;
1652
1653     assert(pull != nullptr);
1654
1655     /* Ideally we should check external potential registration only during
1656      * the initialization phase, but that requires another function call
1657      * that should be done exactly in the right order. So we check here.
1658      */
1659     check_external_potential_registration(pull);
1660
1661     if (pull->comm.bParticipate)
1662     {
1663         real dVdl = 0;
1664
1665         pull_calc_coms(cr, pull, masses, pbc, t, x, {});
1666
1667         rvec*      f             = as_rvec_array(force->force_.data());
1668         matrix     virial        = { { 0 } };
1669         const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1670         for (pull_coord_work_t& pcrd : pull->coord)
1671         {
1672             /* For external potential the force is assumed to be given by an external module by a
1673                call to apply_pull_coord_external_force */
1674             if (pcrd.params.eType == PullingAlgorithm::Constraint
1675                 || pcrd.params.eType == PullingAlgorithm::External)
1676             {
1677                 continue;
1678             }
1679             PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
1680                     *pull, &pcrd, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
1681
1682             if (pcrd.params.eGeom == PullGroupGeometry::Transformation)
1683             {
1684                 gmx::applyTransformationPullCoordForce(
1685                         &pcrd,
1686                         gmx::ArrayRef<pull_coord_work_t>(pull->coord).subArray(0, pcrd.params.coordIndex),
1687                         pcrd.scalarForce);
1688             }
1689             else
1690             {
1691                 /* Distribute the force over the atoms in the pulled groups */
1692                 apply_forces_coord(pcrd, pull->group, pullCoordForces, masses, f);
1693             }
1694         }
1695
1696         if (MASTER(cr))
1697         {
1698             force->addVirialContribution(virial);
1699             *dvdlambda += dVdl;
1700         }
1701     }
1702
1703     GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
1704                "Too few or too many external pull potentials have been applied the previous step");
1705     /* All external pull potentials still need to be applied */
1706     pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
1707
1708     return (MASTER(cr) ? V : 0.0);
1709 }
1710
1711 void pull_constraint(struct pull_t*       pull,
1712                      ArrayRef<const real> masses,
1713                      const t_pbc&         pbc,
1714                      const t_commrec*     cr,
1715                      double               dt,
1716                      double               t,
1717                      ArrayRef<RVec>       x,
1718                      ArrayRef<RVec>       xp,
1719                      ArrayRef<RVec>       v,
1720                      tensor               vir)
1721 {
1722     assert(pull != nullptr);
1723
1724     if (pull->comm.bParticipate)
1725     {
1726         pull_calc_coms(cr, pull, masses, pbc, t, x, xp);
1727
1728         do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1729     }
1730 }
1731
1732 void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
1733 {
1734     gmx_domdec_t* dd;
1735     pull_comm_t*  comm;
1736     gmx_bool      bMustParticipate;
1737
1738     dd = cr->dd;
1739
1740     comm = &pull->comm;
1741
1742     /* We always make the master node participate, such that it can do i/o,
1743      * add the virial and to simplify MC type extensions people might have.
1744      */
1745     bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1746
1747     for (pull_group_work_t& group : pull->group)
1748     {
1749         if (!group.globalWeights.empty())
1750         {
1751             group.localWeights.resize(group.atomSet.numAtomsLocal());
1752             for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1753             {
1754                 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1755             }
1756         }
1757
1758         GMX_ASSERT(bMustParticipate || dd != nullptr,
1759                    "Either all ranks (including this rank) participate, or we use DD and need to "
1760                    "have access to dd here");
1761
1762         /* We should participate if we have pull or pbc atoms */
1763         if (!bMustParticipate
1764             && (group.atomSet.numAtomsLocal() > 0
1765                 || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
1766         {
1767             bMustParticipate = TRUE;
1768         }
1769     }
1770
1771     if (!comm->bParticipateAll)
1772     {
1773         /* Keep currently not required ranks in the communicator
1774          * if they needed to participate up to 20 decompositions ago.
1775          * This avoids frequent rebuilds due to atoms jumping back and forth.
1776          */
1777         const int64_t history_count = 20;
1778         gmx_bool      bWillParticipate;
1779         int           count[2];
1780
1781         /* Increase the decomposition counter for the current call */
1782         comm->setup_count++;
1783
1784         if (bMustParticipate)
1785         {
1786             comm->must_count = comm->setup_count;
1787         }
1788
1789         bWillParticipate =
1790                 bMustParticipate
1791                 || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
1792
1793         if (debug && dd != nullptr)
1794         {
1795             fprintf(debug,
1796                     "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
1797                     dd->rank,
1798                     gmx::boolToString(bMustParticipate),
1799                     gmx::boolToString(bWillParticipate));
1800         }
1801
1802         if (bWillParticipate)
1803         {
1804             /* Count the number of ranks that we want to have participating */
1805             count[0] = 1;
1806             /* Count the number of ranks that need to be added */
1807             count[1] = comm->bParticipate ? 0 : 1;
1808         }
1809         else
1810         {
1811             count[0] = 0;
1812             count[1] = 0;
1813         }
1814
1815         /* The cost of this global operation will be less that the cost
1816          * of the extra MPI_Comm_split calls that we can avoid.
1817          */
1818         gmx_sumi(2, count, cr);
1819
1820         /* If we are missing ranks or if we have 20% more ranks than needed
1821          * we make a new sub-communicator.
1822          */
1823         if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
1824         {
1825             if (debug)
1826             {
1827                 fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
1828             }
1829
1830 #if GMX_MPI
1831             if (comm->mpi_comm_com != MPI_COMM_NULL)
1832             {
1833                 MPI_Comm_free(&comm->mpi_comm_com);
1834             }
1835
1836             /* This might be an extremely expensive operation, so we try
1837              * to avoid this splitting as much as possible.
1838              */
1839             assert(dd != nullptr);
1840             MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
1841 #endif
1842
1843             comm->bParticipate = bWillParticipate;
1844             comm->nparticipate = count[0];
1845
1846             /* When we use the previous COM for PBC, we need to broadcast
1847              * the previous COM to ranks that have joined the communicator.
1848              */
1849             for (pull_group_work_t& group : pull->group)
1850             {
1851                 if (group.epgrppbc == epgrppbcPREVSTEPCOM)
1852                 {
1853                     GMX_ASSERT(comm->bParticipate || !MASTER(cr),
1854                                "The master rank has to participate, as it should pass an up to "
1855                                "date prev. COM "
1856                                "to bcast here as well as to e.g. checkpointing");
1857
1858                     gmx_bcast(sizeof(dvec), group.x_prev_step, cr->mpi_comm_mygroup);
1859                 }
1860             }
1861         }
1862     }
1863
1864     /* Since the PBC of atoms might have changed, we need to update the PBC */
1865     pull->bSetPBCatoms = TRUE;
1866 }
1867
1868 static void init_pull_group_index(FILE*              fplog,
1869                                   const t_commrec*   cr,
1870                                   int                g,
1871                                   pull_group_work_t* pg,
1872                                   gmx_bool           bConstraint,
1873                                   const ivec         pulldim_con,
1874                                   const gmx_mtop_t&  mtop,
1875                                   const t_inputrec*  ir,
1876                                   real               lambda)
1877 {
1878     /* With EM and BD there are no masses in the integrator.
1879      * But we still want to have the correct mass-weighted COMs.
1880      * So we store the real masses in the weights.
1881      */
1882     const bool setWeights = (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1883                              || ir->eI == IntegrationAlgorithm::BD);
1884
1885     /* In parallel, store we need to extract localWeights from weights at DD time */
1886     std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1887
1888     const SimulationGroups& groups = mtop.groups;
1889
1890     /* Count frozen dimensions and (weighted) mass */
1891     int    nfrozen = 0;
1892     double tmass   = 0;
1893     double wmass   = 0;
1894     double wwmass  = 0;
1895     int    molb    = 0;
1896     for (int i = 0; i < int(pg->params.ind.size()); i++)
1897     {
1898         int ii = pg->params.ind[i];
1899         if (bConstraint && ir->opts.nFreeze)
1900         {
1901             for (int d = 0; d < DIM; d++)
1902             {
1903                 if (pulldim_con[d] == 1
1904                     && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
1905                 {
1906                     nfrozen++;
1907                 }
1908             }
1909         }
1910         const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
1911         real          m;
1912         if (ir->efep == FreeEnergyPerturbationType::No)
1913         {
1914             m = atom.m;
1915         }
1916         else
1917         {
1918             m = (1 - lambda) * atom.m + lambda * atom.mB;
1919         }
1920         real w;
1921         if (!pg->params.weight.empty())
1922         {
1923             w = pg->params.weight[i];
1924         }
1925         else
1926         {
1927             w = 1;
1928         }
1929         if (EI_ENERGY_MINIMIZATION(ir->eI))
1930         {
1931             /* Move the mass to the weight */
1932             w *= m;
1933             m = 1;
1934         }
1935         else if (ir->eI == IntegrationAlgorithm::BD)
1936         {
1937             real mbd;
1938             if (ir->bd_fric != 0.0F)
1939             {
1940                 mbd = ir->bd_fric * ir->delta_t;
1941             }
1942             else
1943             {
1944                 if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
1945                 {
1946                     mbd = ir->delta_t / ir->opts.tau_t[0];
1947                 }
1948                 else
1949                 {
1950                     mbd = ir->delta_t
1951                           / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
1952                 }
1953             }
1954             w *= m / mbd;
1955             m = mbd;
1956         }
1957         if (setWeights)
1958         {
1959             weights.push_back(w);
1960         }
1961         tmass += m;
1962         wmass += m * w;
1963         wwmass += m * w * w;
1964     }
1965
1966     if (wmass == 0)
1967     {
1968         /* We can have single atom groups with zero mass with potential pulling
1969          * without cosine weighting.
1970          */
1971         if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1972         {
1973             /* With one atom the mass doesn't matter */
1974             wwmass = 1;
1975         }
1976         else
1977         {
1978             gmx_fatal(FARGS,
1979                       "The total%s mass of pull group %d is zero",
1980                       !pg->params.weight.empty() ? " weighted" : "",
1981                       g);
1982         }
1983     }
1984     if (fplog)
1985     {
1986         fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
1987         if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI)
1988             || ir->eI == IntegrationAlgorithm::BD)
1989         {
1990             fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
1991         }
1992         if (pg->epgrppbc == epgrppbcCOS)
1993         {
1994             fprintf(fplog, ", cosine weighting will be used");
1995         }
1996         fprintf(fplog, "\n");
1997     }
1998
1999     if (nfrozen == 0)
2000     {
2001         /* A value != 0 signals not frozen, it is updated later */
2002         pg->invtm = -1.0;
2003     }
2004     else
2005     {
2006         int ndim = 0;
2007         for (int d = 0; d < DIM; d++)
2008         {
2009             ndim += pulldim_con[d] * pg->params.ind.size();
2010         }
2011         if (fplog && nfrozen > 0 && nfrozen < ndim)
2012         {
2013             fprintf(fplog,
2014                     "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
2015                     "         that are subject to pulling are frozen.\n"
2016                     "         For constraint pulling the whole group will be frozen.\n\n",
2017                     g);
2018         }
2019         pg->invtm  = 0.0;
2020         pg->wscale = 1.0;
2021     }
2022 }
2023
2024 struct pull_t* init_pull(FILE*                     fplog,
2025                          const pull_params_t*      pull_params,
2026                          const t_inputrec*         ir,
2027                          const gmx_mtop_t&         mtop,
2028                          const t_commrec*          cr,
2029                          gmx::LocalAtomSetManager* atomSets,
2030                          real                      lambda)
2031 {
2032     struct pull_t* pull;
2033     pull_comm_t*   comm;
2034
2035     pull = new pull_t;
2036
2037     /* Copy the pull parameters */
2038     pull->params = *pull_params;
2039
2040     /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2041     const int maxNumThreads = std::max(1, gmx_omp_nthreads_get(ModuleMultiThread::Default));
2042
2043     for (int i = 0; i < pull_params->ngroup; ++i)
2044     {
2045         pull->group.emplace_back(pull_params->group[i],
2046                                  atomSets->add(pull_params->group[i].ind),
2047                                  pull_params->bSetPbcRefToPrevStepCOM,
2048                                  maxNumThreads);
2049     }
2050
2051     if (cr != nullptr && haveDDAtomOrdering(*cr))
2052     {
2053         /* Set up the global to local atom mapping for PBC atoms */
2054         for (pull_group_work_t& group : pull->group)
2055         {
2056             if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
2057             {
2058                 /* pbcAtomSet consists of a single atom */
2059                 group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
2060                         atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
2061             }
2062         }
2063     }
2064
2065     pull->bPotential   = FALSE;
2066     pull->bConstraint  = FALSE;
2067     pull->bCylinder    = FALSE;
2068     pull->bAngle       = FALSE;
2069     pull->bXOutAverage = pull_params->bXOutAverage;
2070     pull->bFOutAverage = pull_params->bFOutAverage;
2071
2072     GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
2073                        "pull group 0 is an absolute reference group and should not contain atoms");
2074
2075     pull->numCoordinatesWithExternalPotential = 0;
2076
2077     for (int c = 0; c < pull->params.ncoord; c++)
2078     {
2079         GMX_RELEASE_ASSERT(pull_params->coord[c].coordIndex == c,
2080                            "The stored index should match the position in the vector");
2081
2082         /* Construct a pull coordinate, copying all coordinate parameters */
2083         pull->coord.emplace_back(pull_params->coord[c]);
2084
2085         pull_coord_work_t* pcrd = &pull->coord.back();
2086
2087         switch (pcrd->params.eGeom)
2088         {
2089             case PullGroupGeometry::Distance:
2090             case PullGroupGeometry::DirectionRelative: /* Direction vector is determined at each step */
2091             case PullGroupGeometry::Angle:
2092             case PullGroupGeometry::Dihedral: break;
2093             case PullGroupGeometry::Direction:
2094             case PullGroupGeometry::DirectionPBC:
2095             case PullGroupGeometry::Cylinder:
2096             case PullGroupGeometry::Transformation:
2097             case PullGroupGeometry::AngleAxis:
2098                 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2099                 break;
2100             default:
2101                 /* We allow reading of newer tpx files with new pull geometries,
2102                  * but with the same tpx format, with old code. A new geometry
2103                  * only adds a new enum value, which will be out of range for
2104                  * old code. The first place we need to generate an error is
2105                  * here, since the pull code can't handle this.
2106                  * The enum can be used before arriving here only for printing
2107                  * the string corresponding to the geometry, which will result
2108                  * in printing "UNDEFINED".
2109                  */
2110                 gmx_fatal(FARGS,
2111                           "Pull geometry not supported for pull coordinate %d. The geometry enum "
2112                           "%s in the input is larger than that supported by the code (up to %d). "
2113                           "You are probably reading a tpr file generated with a newer version of "
2114                           "GROMACS with an binary from an older version of Gromacs.",
2115                           c + 1,
2116                           enumValueToString(pcrd->params.eGeom),
2117                           static_cast<int>(PullGroupGeometry::Count) - 1);
2118         }
2119
2120         if (pcrd->params.eType == PullingAlgorithm::Constraint)
2121         {
2122             /* Check restrictions of the constraint pull code */
2123             if (pcrd->params.eGeom == PullGroupGeometry::Cylinder
2124                 || pcrd->params.eGeom == PullGroupGeometry::DirectionRelative
2125                 || pcrd->params.eGeom == PullGroupGeometry::Angle
2126                 || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2127                 || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2128             {
2129                 gmx_fatal(FARGS,
2130                           "Pulling of type %s can not be combined with geometry %s. Consider using "
2131                           "pull type %s.",
2132                           enumValueToString(pcrd->params.eType),
2133                           enumValueToString(pcrd->params.eGeom),
2134                           enumValueToString(PullingAlgorithm::Umbrella));
2135             }
2136
2137             GMX_RELEASE_ASSERT(
2138                     !ir->useMts,
2139                     "Constraint pulling can not be combined with multiple time stepping");
2140
2141             pull->bConstraint = TRUE;
2142         }
2143         else
2144         {
2145             pull->bPotential = TRUE;
2146         }
2147
2148         if (pcrd->params.eGeom == PullGroupGeometry::Cylinder)
2149         {
2150             pull->bCylinder = TRUE;
2151         }
2152         else if (pcrd->params.eGeom == PullGroupGeometry::Angle
2153                  || pcrd->params.eGeom == PullGroupGeometry::Dihedral
2154                  || pcrd->params.eGeom == PullGroupGeometry::AngleAxis)
2155         {
2156             pull->bAngle = TRUE;
2157         }
2158
2159         /* We only need to calculate the plain COM of a group
2160          * when it is not only used as a cylinder group.
2161          * Also the absolute reference group 0 needs no COM computation.
2162          */
2163         for (int i = 0; i < pcrd->params.ngroup; i++)
2164         {
2165             int groupIndex = pcrd->params.group[i];
2166             if (groupIndex > 0 && !(pcrd->params.eGeom == PullGroupGeometry::Cylinder && i == 0))
2167             {
2168                 pull->group[groupIndex].needToCalcCom = true;
2169             }
2170         }
2171
2172         /* With non-zero rate the reference value is set at every step */
2173         if (pcrd->params.rate == 0)
2174         {
2175             /* Initialize the constant reference value */
2176             if (pcrd->params.eType != PullingAlgorithm::External)
2177             {
2178                 pcrd->value_ref = sanitizePullCoordReferenceValue(
2179                         pcrd->params,
2180                         pcrd->params.init * pull_conversion_factor_userinput2internal(pcrd->params));
2181             }
2182             else
2183             {
2184                 /* With an external pull potential, the reference value loses
2185                  * it's meaning and should not be used. Setting it to zero
2186                  * makes any terms dependent on it disappear.
2187                  * The only issue this causes is that with cylinder pulling
2188                  * no atoms of the cylinder group within the cylinder radius
2189                  * should be more than half a box length away from the COM of
2190                  * the pull group along the axial direction.
2191                  */
2192                 pcrd->value_ref = 0.0;
2193             }
2194         }
2195
2196         if (pcrd->params.eType == PullingAlgorithm::External)
2197         {
2198             GMX_RELEASE_ASSERT(
2199                     pcrd->params.rate == 0,
2200                     "With an external potential, a pull coordinate should have rate = 0");
2201
2202             /* This external potential needs to be registered later */
2203             pull->numCoordinatesWithExternalPotential++;
2204         }
2205         pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2206     }
2207
2208     pull->numUnregisteredExternalPotentials             = pull->numCoordinatesWithExternalPotential;
2209     pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2210
2211     pull->pbcType = ir->pbcType;
2212     switch (pull->pbcType)
2213     {
2214         case PbcType::No: pull->npbcdim = 0; break;
2215         case PbcType::XY: pull->npbcdim = 2; break;
2216         default: pull->npbcdim = 3; break;
2217     }
2218
2219     if (fplog)
2220     {
2221         gmx_bool bAbs, bCos;
2222
2223         bAbs = FALSE;
2224         for (const pull_coord_work_t& coord : pull->coord)
2225         {
2226             if (pull->group[coord.params.group[0]].params.ind.empty()
2227                 || pull->group[coord.params.group[1]].params.ind.empty())
2228             {
2229                 bAbs = TRUE;
2230             }
2231         }
2232
2233         fprintf(fplog, "\n");
2234         if (pull->bPotential)
2235         {
2236             fprintf(fplog, "Will apply potential COM pulling\n");
2237         }
2238         if (pull->bConstraint)
2239         {
2240             fprintf(fplog, "Will apply constraint COM pulling\n");
2241         }
2242         // Don't include the reference group 0 in output, so we report ngroup-1
2243         int numRealGroups = pull->group.size() - 1;
2244         GMX_RELEASE_ASSERT(numRealGroups > 0,
2245                            "The reference absolute position pull group should always be present");
2246         fprintf(fplog,
2247                 "with %zu pull coordinate%s and %d group%s\n",
2248                 pull->coord.size(),
2249                 pull->coord.size() == 1 ? "" : "s",
2250                 numRealGroups,
2251                 numRealGroups == 1 ? "" : "s");
2252         if (bAbs)
2253         {
2254             fprintf(fplog, "with an absolute reference\n");
2255         }
2256         bCos = FALSE;
2257         // Don't include the reference group 0 in loop
2258         for (size_t g = 1; g < pull->group.size(); g++)
2259         {
2260             if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
2261             {
2262                 /* We are using cosine weighting */
2263                 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2264                 bCos = TRUE;
2265             }
2266         }
2267         if (bCos)
2268         {
2269             please_cite(fplog, "Engin2010");
2270         }
2271     }
2272
2273     pull->bRefAt = FALSE;
2274     pull->cosdim = -1;
2275     for (size_t g = 0; g < pull->group.size(); g++)
2276     {
2277         pull_group_work_t* pgrp;
2278
2279         pgrp = &pull->group[g];
2280         if (!pgrp->params.ind.empty())
2281         {
2282             /* There is an issue when a group is used in multiple coordinates
2283              * and constraints are applied in different dimensions with atoms
2284              * frozen in some, but not all dimensions.
2285              * Since there is only one mass array per group, we can't have
2286              * frozen/non-frozen atoms for different coords at the same time.
2287              * But since this is a very exotic case, we don't check for this.
2288              * A warning is printed in init_pull_group_index.
2289              */
2290
2291             gmx_bool bConstraint;
2292             ivec     pulldim, pulldim_con;
2293
2294             /* Loop over all pull coordinates to see along which dimensions
2295              * this group is pulled and if it is involved in constraints.
2296              */
2297             bConstraint = FALSE;
2298             clear_ivec(pulldim);
2299             clear_ivec(pulldim_con);
2300             for (const pull_coord_work_t& coord : pull->coord)
2301             {
2302                 gmx_bool bGroupUsed = FALSE;
2303                 for (int gi = 0; gi < coord.params.ngroup; gi++)
2304                 {
2305                     if (coord.params.group[gi] == static_cast<int>(g))
2306                     {
2307                         bGroupUsed = TRUE;
2308                     }
2309                 }
2310
2311                 if (bGroupUsed)
2312                 {
2313                     for (int m = 0; m < DIM; m++)
2314                     {
2315                         if (coord.params.dim[m] == 1)
2316                         {
2317                             pulldim[m] = 1;
2318
2319                             if (coord.params.eType == PullingAlgorithm::Constraint)
2320                             {
2321                                 bConstraint    = TRUE;
2322                                 pulldim_con[m] = 1;
2323                             }
2324                         }
2325                     }
2326                 }
2327             }
2328
2329             /* Determine if we need to take PBC into account for calculating
2330              * the COM's of the pull groups.
2331              */
2332             switch (pgrp->epgrppbc)
2333             {
2334                 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
2335                 case epgrppbcCOS:
2336                     if (!pgrp->params.weight.empty())
2337                     {
2338                         gmx_fatal(FARGS,
2339                                   "Pull groups can not have relative weights and cosine weighting "
2340                                   "at same time");
2341                     }
2342                     for (int m = 0; m < DIM; m++)
2343                     {
2344                         if (m < pull->npbcdim && pulldim[m] == 1)
2345                         {
2346                             if (pull->cosdim >= 0 && pull->cosdim != m)
2347                             {
2348                                 gmx_fatal(FARGS,
2349                                           "Can only use cosine weighting with pulling in one "
2350                                           "dimension (use mdp option pull-coord?-dim)");
2351                             }
2352                             pull->cosdim = m;
2353                         }
2354                     }
2355                     break;
2356                 case epgrppbcNONE: break;
2357             }
2358
2359             /* Set the indices */
2360             init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
2361         }
2362         else
2363         {
2364             /* Absolute reference, set the inverse mass to zero.
2365              * This is only relevant (and used) with constraint pulling.
2366              */
2367             pgrp->invtm  = 0;
2368             pgrp->wscale = 1;
2369         }
2370     }
2371
2372     /* If we use cylinder coordinates, do some initialising for them */
2373     if (pull->bCylinder)
2374     {
2375         for (pull_coord_work_t& coord : pull->coord)
2376         {
2377             if (coord.params.eGeom == PullGroupGeometry::Cylinder)
2378             {
2379                 if (pull->group[coord.params.group[0]].params.ind.empty())
2380                 {
2381                     gmx_fatal(FARGS,
2382                               "A cylinder pull group is not supported when using absolute "
2383                               "reference!\n");
2384                 }
2385             }
2386             const auto& group0  = pull->group[coord.params.group[0]];
2387             coord.dynamicGroup0 = std::make_unique<pull_group_work_t>(
2388                     group0.params, group0.atomSet, pull->params.bSetPbcRefToPrevStepCOM, maxNumThreads);
2389         }
2390     }
2391
2392     pull->comSums.resize(maxNumThreads);
2393
2394     comm = &pull->comm;
2395
2396 #if GMX_MPI
2397     /* Use a sub-communicator when we have more than 32 ranks, but not
2398      * when we have an external pull potential, since then the external
2399      * potential provider expects each rank to have the coordinate.
2400      */
2401     comm->bParticipateAll = (cr == nullptr || !haveDDAtomOrdering(*cr) || cr->dd->nnodes <= 32
2402                              || pull->numCoordinatesWithExternalPotential > 0
2403                              || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2404     /* This sub-commicator is not used with comm->bParticipateAll,
2405      * so we can always initialize it to NULL.
2406      */
2407     comm->mpi_comm_com = MPI_COMM_NULL;
2408     comm->nparticipate = 0;
2409     comm->isMasterRank = (cr == nullptr || MASTER(cr));
2410 #else
2411     /* No MPI: 1 rank: all ranks pull */
2412     comm->bParticipateAll = TRUE;
2413     comm->isMasterRank    = true;
2414 #endif
2415     comm->bParticipate = comm->bParticipateAll;
2416     comm->setup_count  = 0;
2417     comm->must_count   = 0;
2418
2419     if (!comm->bParticipateAll && fplog != nullptr)
2420     {
2421         fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2422     }
2423
2424     comm->pbcAtomBuffer.resize(pull->group.size());
2425     comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
2426     if (pull->bCylinder)
2427     {
2428         comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
2429     }
2430
2431     /* We still need to initialize the PBC reference coordinates */
2432     pull->bSetPBCatoms = TRUE;
2433
2434     pull->out_x = nullptr;
2435     pull->out_f = nullptr;
2436
2437     return pull;
2438 }
2439
2440 static void destroy_pull(struct pull_t* pull)
2441 {
2442 #if GMX_MPI
2443     if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2444     {
2445         MPI_Comm_free(&pull->comm.mpi_comm_com);
2446     }
2447 #endif
2448
2449     delete pull;
2450 }
2451
2452 void preparePrevStepPullCom(const t_inputrec*    ir,
2453                             pull_t*              pull_work,
2454                             ArrayRef<const real> masses,
2455                             t_state*             state,
2456                             const t_state*       state_global,
2457                             const t_commrec*     cr,
2458                             bool                 startingFromCheckpoint)
2459 {
2460     if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
2461     {
2462         return;
2463     }
2464     allocStatePrevStepPullCom(state, pull_work);
2465     if (startingFromCheckpoint)
2466     {
2467         if (MASTER(cr))
2468         {
2469             state->pull_com_prev_step = state_global->pull_com_prev_step;
2470         }
2471         if (PAR(cr))
2472         {
2473             /* Only the master rank has the checkpointed COM from the previous step */
2474             gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(),
2475                       &state->pull_com_prev_step[0],
2476                       cr->mpi_comm_mygroup);
2477         }
2478         setPrevStepPullComFromState(pull_work, state);
2479     }
2480     else
2481     {
2482         t_pbc pbc;
2483         set_pbc(&pbc, ir->pbcType, state->box);
2484         initPullComFromPrevStep(
2485                 cr, pull_work, masses, pbc, state->x.arrayRefWithPadding().unpaddedArrayRef());
2486         updatePrevStepPullCom(pull_work, state);
2487     }
2488 }
2489
2490 void finish_pull(struct pull_t* pull)
2491 {
2492     check_external_potential_registration(pull);
2493
2494     if (pull->out_x)
2495     {
2496         gmx_fio_fclose(pull->out_x);
2497     }
2498     if (pull->out_f)
2499     {
2500         gmx_fio_fclose(pull->out_f);
2501     }
2502
2503     destroy_pull(pull);
2504 }
2505
2506 bool pull_have_potential(const pull_t& pull)
2507 {
2508     return pull.bPotential;
2509 }
2510
2511 bool pull_have_constraint(const pull_t& pull)
2512 {
2513     return pull.bConstraint;
2514 }
2515
2516 bool pull_have_constraint(const pull_params_t& pullParameters)
2517 {
2518     for (int c = 0; c < pullParameters.ncoord; c++)
2519     {
2520         if (pullParameters.coord[c].eType == PullingAlgorithm::Constraint)
2521         {
2522             return true;
2523         }
2524     }
2525     return false;
2526 }