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