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