src/gromacs/pulling/pull.cpp

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