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