src/gromacs/pulling/pull.cpp

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