src/gromacs/pulling/pull.cpp

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