src/gromacs/pulling/pull.cpp

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