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