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