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