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