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