src/gromacs/pulling/pull.cpp

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