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                     /* Keep static analyzer happy */
1172                     clear_dvec(dr0);
1173                     clear_dvec(dr1);
1174             }
1175
1176             /* DEBUG */
1177             if (debug)
1178             {
1179                 int g0, g1;
1180
1181                 g0 = pcrd->params.group[0];
1182                 g1 = pcrd->params.group[1];
1183                 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
1184                 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
1185                 fprintf(debug,
1186                         "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1187                         rnew[g0][0], rnew[g0][1], rnew[g0][2],
1188                         rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
1189                 fprintf(debug,
1190                         "Pull ref %8s %8s %8s   %8s %8s %8s d: %8.5f\n",
1191                         "", "", "", "", "", "", pcrd->value_ref);
1192                 fprintf(debug,
1193                         "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1194                         dr0[0], dr0[1], dr0[2],
1195                         dr1[0], dr1[1], dr1[2],
1196                         dnorm(tmp3));
1197             } /* END DEBUG */
1198
1199             /* Update the COMs with dr */
1200             dvec_inc(rnew[pcrd->params.group[1]], dr1);
1201             dvec_inc(rnew[pcrd->params.group[0]], dr0);
1202         }
1203
1204         /* Check if all constraints are fullfilled now */
1205         for (pull_coord_work_t &coord : pull->coord)
1206         {
1207             if (coord.params.eType != epullCONSTRAINT)
1208             {
1209                 continue;
1210             }
1211
1212             low_get_pull_coord_dr(pull, &coord, pbc,
1213                                   rnew[coord.params.group[1]],
1214                                   rnew[coord.params.group[0]],
1215                                   -1, unc_ij);
1216
1217             switch (coord.params.eGeom)
1218             {
1219                 case epullgDIST:
1220                     bConverged =
1221                         fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1222                     break;
1223                 case epullgDIR:
1224                 case epullgDIRPBC:
1225                 case epullgCYL:
1226                     for (m = 0; m < DIM; m++)
1227                     {
1228                         vec[m] = coord.spatialData.vec[m];
1229                     }
1230                     inpr = diprod(unc_ij, vec);
1231                     dsvmul(inpr, vec, unc_ij);
1232                     bConverged =
1233                         fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1234                     break;
1235             }
1236
1237             if (!bConverged)
1238             {
1239                 if (debug)
1240                 {
1241                     fprintf(debug, "NOT CONVERGED YET: Group %d:"
1242                             "d_ref = %f, current d = %f\n",
1243                             g, coord.value_ref, dnorm(unc_ij));
1244                 }
1245
1246                 bConverged_all = FALSE;
1247             }
1248         }
1249
1250         niter++;
1251         /* if after all constraints are dealt with and bConverged is still TRUE
1252            we're finished, if not we do another iteration */
1253     }
1254     if (niter > max_iter)
1255     {
1256         gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1257     }
1258
1259     /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1260
1261     if (v)
1262     {
1263         invdt = 1/dt;
1264     }
1265
1266     /* update atoms in the groups */
1267     for (g = 0; g < pull->ngroup; g++)
1268     {
1269         const pull_group_work_t *pgrp;
1270         dvec                     dr;
1271
1272         pgrp = &pull->group[g];
1273
1274         /* get the final constraint displacement dr for group g */
1275         dvec_sub(rnew[g], pgrp->xp, dr);
1276
1277         if (dnorm2(dr) == 0)
1278         {
1279             /* No displacement, no update necessary */
1280             continue;
1281         }
1282
1283         /* update the atom positions */
1284         copy_dvec(dr, tmp);
1285         for (j = 0; j < pgrp->nat_loc; j++)
1286         {
1287             ii = pgrp->ind_loc[j];
1288             if (pgrp->weight_loc)
1289             {
1290                 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
1291             }
1292             for (m = 0; m < DIM; m++)
1293             {
1294                 x[ii][m] += tmp[m];
1295             }
1296             if (v)
1297             {
1298                 for (m = 0; m < DIM; m++)
1299                 {
1300                     v[ii][m] += invdt*tmp[m];
1301                 }
1302             }
1303         }
1304     }
1305
1306     /* calculate the constraint forces, used for output and virial only */
1307     for (size_t c = 0; c < pull->coord.size(); c++)
1308     {
1309         pull_coord_work_t *pcrd;
1310
1311         pcrd = &pull->coord[c];
1312
1313         if (pcrd->params.eType != epullCONSTRAINT)
1314         {
1315             continue;
1316         }
1317
1318         /* Accumulate the forces, in case we have multiple constraint steps */
1319         double force = dr_tot[c]/((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)*dt*dt);
1320         pcrd->scalarForce += force;
1321
1322         if (vir != nullptr && pcrd->params.eGeom != epullgDIRPBC && bMaster)
1323         {
1324             double f_invr;
1325
1326             /* Add the pull contribution to the virial */
1327             /* We have already checked above that r_ij[c] != 0 */
1328             f_invr = pcrd->scalarForce/dnorm(r_ij[c]);
1329
1330             for (j = 0; j < DIM; j++)
1331             {
1332                 for (m = 0; m < DIM; m++)
1333                 {
1334                     vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
1335                 }
1336             }
1337         }
1338     }
1339
1340     /* finished! I hope. Give back some memory */
1341     sfree(r_ij);
1342     sfree(dr_tot);
1343     sfree(rnew);
1344 }
1345
1346 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1347 {
1348     for (int j = 0; j < DIM; j++)
1349     {
1350         for (int m = 0; m < DIM; m++)
1351         {
1352             vir[j][m] -= 0.5*f[j]*dr[m];
1353         }
1354     }
1355 }
1356
1357 /* Adds the pull contribution to the virial */
1358 static void add_virial_coord(tensor                       vir,
1359                              const pull_coord_work_t     &pcrd,
1360                              const PullCoordVectorForces &forces)
1361 {
1362     if (vir != nullptr && pcrd.params.eGeom != epullgDIRPBC)
1363     {
1364         /* Add the pull contribution for each distance vector to the virial. */
1365         add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1366         if (pcrd.params.ngroup >= 4)
1367         {
1368             add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1369         }
1370         if (pcrd.params.ngroup >= 6)
1371         {
1372             add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1373         }
1374     }
1375 }
1376
1377 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t *pcrd,
1378                                                        double dev, real lambda,
1379                                                        real *V, real *dVdl)
1380 {
1381     real   k, dkdl;
1382
1383     k    = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
1384     dkdl = pcrd->params.kB - pcrd->params.k;
1385
1386     switch (pcrd->params.eType)
1387     {
1388         case epullUMBRELLA:
1389         case epullFLATBOTTOM:
1390         case epullFLATBOTTOMHIGH:
1391             /* The only difference between an umbrella and a flat-bottom
1392              * potential is that a flat-bottom is zero above or below
1393                the reference value.
1394              */
1395             if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0) ||
1396                 (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
1397             {
1398                 dev = 0;
1399             }
1400
1401             pcrd->scalarForce    = -k*dev;
1402             *V                  += 0.5*   k*gmx::square(dev);
1403             *dVdl               += 0.5*dkdl*gmx::square(dev);
1404             break;
1405         case epullCONST_F:
1406             pcrd->scalarForce    = -k;
1407             *V                  +=    k*pcrd->spatialData.value;
1408             *dVdl               += dkdl*pcrd->spatialData.value;
1409             break;
1410         case epullEXTERNAL:
1411             gmx_incons("the scalar pull force should not be calculated internally for pull type external");
1412             break;
1413         default:
1414             gmx_incons("Unsupported pull type in do_pull_pot");
1415     }
1416 }
1417
1418 static PullCoordVectorForces
1419 calculateVectorForces(const pull_coord_work_t &pcrd)
1420 {
1421     const t_pull_coord         &params      = pcrd.params;
1422     const PullCoordSpatialData &spatialData = pcrd.spatialData;
1423
1424     /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1425     PullCoordVectorForces    forces;
1426
1427     if (params.eGeom == epullgDIST)
1428     {
1429         double invdr01 = spatialData.value > 0 ? 1./spatialData.value : 0.;
1430         for (int m = 0; m < DIM; m++)
1431         {
1432             forces.force01[m] = pcrd.scalarForce*spatialData.dr01[m]*invdr01;
1433         }
1434     }
1435     else if (params.eGeom == epullgANGLE)
1436     {
1437
1438         double cos_theta, cos_theta2;
1439
1440         cos_theta  = cos(spatialData.value);
1441         cos_theta2 = gmx::square(cos_theta);
1442
1443         /* The force at theta = 0, pi is undefined so we should not apply any force.
1444          * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1445          * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1446          * have to check for this before dividing by their norm below.
1447          */
1448         if (cos_theta2 < 1)
1449         {
1450             /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1451              * and another vector parallel to dr23:
1452              * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1453              * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1454              */
1455             double a       = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1456             double b       = a*cos_theta;
1457             double invdr01 = 1./dnorm(spatialData.dr01);
1458             double invdr23 = 1./dnorm(spatialData.dr23);
1459             dvec   normalized_dr01, normalized_dr23;
1460             dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1461             dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1462
1463             for (int m = 0; m < DIM; m++)
1464             {
1465                 /* Here, f_scal is -dV/dtheta */
1466                 forces.force01[m] = pcrd.scalarForce*invdr01*(a*normalized_dr23[m] - b*normalized_dr01[m]);
1467                 forces.force23[m] = pcrd.scalarForce*invdr23*(a*normalized_dr01[m] - b*normalized_dr23[m]);
1468             }
1469         }
1470         else
1471         {
1472             /* No forces to apply for ill-defined cases*/
1473             clear_dvec(forces.force01);
1474             clear_dvec(forces.force23);
1475         }
1476     }
1477     else if (params.eGeom == epullgANGLEAXIS)
1478     {
1479         double cos_theta, cos_theta2;
1480
1481         /* The angle-axis force is exactly the same as the angle force (above) except that in
1482            this case the second vector (dr23) is replaced by the pull vector. */
1483         cos_theta  = cos(spatialData.value);
1484         cos_theta2 = gmx::square(cos_theta);
1485
1486         if (cos_theta2 < 1)
1487         {
1488             double a, b;
1489             double invdr01;
1490             dvec   normalized_dr01;
1491
1492             invdr01 = 1./dnorm(spatialData.dr01);
1493             dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1494             a       = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1495             b       = a*cos_theta;
1496
1497             for (int m = 0; m < DIM; m++)
1498             {
1499                 forces.force01[m] = pcrd.scalarForce*invdr01*(a*spatialData.vec[m] - b*normalized_dr01[m]);
1500             }
1501         }
1502         else
1503         {
1504             clear_dvec(forces.force01);
1505         }
1506     }
1507     else if (params.eGeom == epullgDIHEDRAL)
1508     {
1509         double m2, n2, tol, sqrdist_32;
1510         dvec   dr32;
1511         /* Note: there is a small difference here compared to the
1512            dihedral force calculations in the bondeds (ref: Bekker 1994).
1513            There rij = ri - rj, while here dr01 = r1 - r0.
1514            However, all distance vectors occur in form of cross or inner products
1515            so that two signs cancel and we end up with the same expressions.
1516            Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1517          */
1518         m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1519         n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1520         dsvmul(-1, spatialData.dr23, dr32);
1521         sqrdist_32 = diprod(dr32, dr32);
1522         tol        = sqrdist_32*GMX_REAL_EPS; /* Avoid tiny angles */
1523         if ((m2 > tol) && (n2 > tol))
1524         {
1525             double a_01, a_23_01, a_23_45, a_45;
1526             double inv_dist_32, inv_sqrdist_32, dist_32;
1527             dvec   u, v;
1528             inv_dist_32    = gmx::invsqrt(sqrdist_32);
1529             inv_sqrdist_32 = inv_dist_32*inv_dist_32;
1530             dist_32        = sqrdist_32*inv_dist_32;
1531
1532             /* Forces on groups 0, 1 */
1533             a_01 = pcrd.scalarForce*dist_32/m2;                    /* scalarForce is -dV/dphi */
1534             dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1535
1536             /* Forces on groups 4, 5 */
1537             a_45 = -pcrd.scalarForce*dist_32/n2;
1538             dsvmul(a_45, spatialData.planevec_n, forces.force45);  /* force on group 5 */
1539
1540             /* Force on groups 2, 3 (defining the axis) */
1541             a_23_01 = -diprod(spatialData.dr01, dr32)*inv_sqrdist_32;
1542             a_23_45 = -diprod(spatialData.dr45, dr32)*inv_sqrdist_32;
1543             dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1544             dsvmul(a_23_45, forces.force45, v);
1545             dvec_sub(u, v, forces.force23);      /* force on group 3 */
1546         }
1547         else
1548         {
1549             /* No force to apply for ill-defined cases */
1550             clear_dvec(forces.force01);
1551             clear_dvec(forces.force23);
1552             clear_dvec(forces.force45);
1553         }
1554     }
1555     else
1556     {
1557         for (int m = 0; m < DIM; m++)
1558         {
1559             forces.force01[m] = pcrd.scalarForce*spatialData.vec[m];
1560         }
1561     }
1562
1563     return forces;
1564 }
1565
1566
1567 /* We use a global mutex for locking access to the pull data structure
1568  * during registration of external pull potential providers.
1569  * We could use a different, local mutex for each pull object, but the overhead
1570  * is extremely small here and registration is only done during initialization.
1571  */
1572 static gmx::Mutex registrationMutex;
1573
1574 using Lock = gmx::lock_guard<gmx::Mutex>;
1575
1576 void register_external_pull_potential(struct pull_t *pull,
1577                                       int            coord_index,
1578                                       const char    *provider)
1579 {
1580     GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1581     GMX_RELEASE_ASSERT(provider != nullptr, "register_external_pull_potential called with NULL as provider name");
1582
1583     if (coord_index < 0 || coord_index >= static_cast<int>(pull->coord.size()))
1584     {
1585         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",
1586                   provider, coord_index + 1, 1, pull->coord.size());
1587     }
1588
1589     pull_coord_work_t *pcrd = &pull->coord[coord_index];
1590
1591     if (pcrd->params.eType != epullEXTERNAL)
1592     {
1593         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'",
1594                   provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
1595     }
1596
1597     GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr, "The external potential provider string for a pull coordinate is NULL");
1598
1599     if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider) != 0)
1600     {
1601         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'",
1602                   provider, coord_index + 1, pcrd->params.externalPotentialProvider);
1603     }
1604
1605     /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1606      * pcrd->bExternalPotentialProviderHasBeenRegistered and
1607      * pull->numUnregisteredExternalPotentials.
1608      */
1609     Lock registrationLock(registrationMutex);
1610
1611     if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1612     {
1613         gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d more than once",
1614                   provider, coord_index + 1);
1615     }
1616
1617     pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1618     pull->numUnregisteredExternalPotentials--;
1619
1620     GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0, "Negative unregistered potentials, the pull code is inconsistent");
1621 }
1622
1623
1624 static void check_external_potential_registration(const struct pull_t *pull)
1625 {
1626     if (pull->numUnregisteredExternalPotentials > 0)
1627     {
1628         size_t c;
1629         for (c = 0; c < pull->coord.size(); c++)
1630         {
1631             if (pull->coord[c].params.eType == epullEXTERNAL &&
1632                 !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1633             {
1634                 break;
1635             }
1636         }
1637
1638         GMX_RELEASE_ASSERT(c < pull->coord.size(), "Internal inconsistency in the pull potential provider counting");
1639
1640         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.",
1641                   pull->numUnregisteredExternalPotentials,
1642                   c + 1, pull->coord[c].params.externalPotentialProvider);
1643     }
1644 }
1645
1646
1647 /* Pull takes care of adding the forces of the external potential.
1648  * The external potential module  has to make sure that the corresponding
1649  * potential energy is added either to the pull term or to a term
1650  * specific to the external module.
1651  */
1652 void apply_external_pull_coord_force(struct pull_t        *pull,
1653                                      int                   coord_index,
1654                                      double                coord_force,
1655                                      const t_mdatoms      *mdatoms,
1656                                      gmx::ForceWithVirial *forceWithVirial)
1657 {
1658     pull_coord_work_t *pcrd;
1659
1660     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");
1661
1662     if (pull->comm.bParticipate)
1663     {
1664         pcrd = &pull->coord[coord_index];
1665
1666         GMX_RELEASE_ASSERT(pcrd->params.eType == epullEXTERNAL, "The pull force can only be set externally on pull coordinates of external type");
1667
1668         GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered, "apply_external_pull_coord_force called for an unregistered pull coordinate");
1669
1670         /* Set the force */
1671         pcrd->scalarForce = coord_force;
1672
1673         /* Calculate the forces on the pull groups */
1674         PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1675
1676         /* Add the forces for this coordinate to the total virial and force */
1677         if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1678         {
1679             matrix virial = { { 0 } };
1680             add_virial_coord(virial, *pcrd, pullCoordForces);
1681             forceWithVirial->addVirialContribution(virial);
1682         }
1683
1684         apply_forces_coord(pull, coord_index, pullCoordForces, mdatoms,
1685                            as_rvec_array(forceWithVirial->force_.data()));
1686     }
1687
1688     pull->numExternalPotentialsStillToBeAppliedThisStep--;
1689 }
1690
1691 /* Calculate the pull potential and scalar force for a pull coordinate.
1692  * Returns the vector forces for the pull coordinate.
1693  */
1694 static PullCoordVectorForces
1695 do_pull_pot_coord(struct pull_t *pull, int coord_ind, t_pbc *pbc,
1696                   double t, real lambda,
1697                   real *V, tensor vir, real *dVdl)
1698 {
1699     pull_coord_work_t &pcrd = pull->coord[coord_ind];
1700
1701     assert(pcrd.params.eType != epullCONSTRAINT);
1702
1703     double dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
1704
1705     calc_pull_coord_scalar_force_and_potential(&pcrd, dev, lambda, V, dVdl);
1706
1707     PullCoordVectorForces pullCoordForces = calculateVectorForces(pcrd);
1708
1709     add_virial_coord(vir, pcrd, pullCoordForces);
1710
1711     return pullCoordForces;
1712 }
1713
1714 real pull_potential(struct pull_t *pull, const t_mdatoms *md, t_pbc *pbc,
1715                     const t_commrec *cr, double t, real lambda,
1716                     const rvec *x, gmx::ForceWithVirial *force, real *dvdlambda)
1717 {
1718     real V = 0;
1719
1720     assert(pull != nullptr);
1721
1722     /* Ideally we should check external potential registration only during
1723      * the initialization phase, but that requires another function call
1724      * that should be done exactly in the right order. So we check here.
1725      */
1726     check_external_potential_registration(pull);
1727
1728     if (pull->comm.bParticipate)
1729     {
1730         real dVdl = 0;
1731
1732         pull_calc_coms(cr, pull, md, pbc, t, x, nullptr);
1733
1734         rvec       *f             = as_rvec_array(force->force_.data());
1735         matrix      virial        = { { 0 } };
1736         const bool  computeVirial = (force->computeVirial_ && MASTER(cr));
1737         for (size_t c = 0; c < pull->coord.size(); c++)
1738         {
1739             /* For external potential the force is assumed to be given by an external module by a call to
1740                apply_pull_coord_external_force */
1741             if (pull->coord[c].params.eType == epullCONSTRAINT || pull->coord[c].params.eType == epullEXTERNAL)
1742             {
1743                 continue;
1744             }
1745
1746             PullCoordVectorForces pullCoordForces =
1747                 do_pull_pot_coord(pull, c, pbc, t, lambda,
1748                                   &V,
1749                                   computeVirial ? virial : nullptr,
1750                                   &dVdl);
1751
1752             /* Distribute the force over the atoms in the pulled groups */
1753             apply_forces_coord(pull, c, pullCoordForces, md, f);
1754         }
1755
1756         if (MASTER(cr))
1757         {
1758             force->addVirialContribution(virial);
1759             *dvdlambda += dVdl;
1760         }
1761     }
1762
1763     GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0, "Too few or too many external pull potentials have been applied the previous step");
1764     /* All external pull potentials still need to be applied */
1765     pull->numExternalPotentialsStillToBeAppliedThisStep =
1766         pull->numCoordinatesWithExternalPotential;
1767
1768     return (MASTER(cr) ? V : 0.0);
1769 }
1770
1771 void pull_constraint(struct pull_t *pull, const t_mdatoms *md, t_pbc *pbc,
1772                      const t_commrec *cr, double dt, double t,
1773                      rvec *x, rvec *xp, rvec *v, tensor vir)
1774 {
1775     assert(pull != nullptr);
1776
1777     if (pull->comm.bParticipate)
1778     {
1779         pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1780
1781         do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1782     }
1783 }
1784
1785 static void make_local_pull_group(gmx_ga2la_t *ga2la,
1786                                   pull_group_work_t *pg, int start, int end)
1787 {
1788     int i, ii;
1789
1790     pg->nat_loc = 0;
1791     for (i = 0; i < pg->params.nat; i++)
1792     {
1793         ii = pg->params.ind[i];
1794         if (ga2la)
1795         {
1796             if (!ga2la_get_home(ga2la, ii, &ii))
1797             {
1798                 ii = -1;
1799             }
1800         }
1801         if (ii >= start && ii < end)
1802         {
1803             /* This is a home atom, add it to the local pull group */
1804             if (pg->nat_loc >= pg->nalloc_loc)
1805             {
1806                 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1807                 srenew(pg->ind_loc, pg->nalloc_loc);
1808                 if (pg->epgrppbc == epgrppbcCOS || pg->params.weight != nullptr)
1809                 {
1810                     srenew(pg->weight_loc, pg->nalloc_loc);
1811                 }
1812             }
1813             pg->ind_loc[pg->nat_loc] = ii;
1814             if (pg->params.weight != nullptr)
1815             {
1816                 pg->weight_loc[pg->nat_loc] = pg->params.weight[i];
1817             }
1818             pg->nat_loc++;
1819         }
1820     }
1821 }
1822
1823 void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull, t_mdatoms *md)
1824 {
1825     gmx_domdec_t   *dd;
1826     pull_comm_t    *comm;
1827     gmx_ga2la_t    *ga2la;
1828     gmx_bool        bMustParticipate;
1829     int             g;
1830
1831     dd = cr->dd;
1832
1833     comm = &pull->comm;
1834
1835     if (dd)
1836     {
1837         ga2la = dd->ga2la;
1838     }
1839     else
1840     {
1841         ga2la = nullptr;
1842     }
1843
1844     /* We always make the master node participate, such that it can do i/o,
1845      * add the virial and to simplify MC type extensions people might have.
1846      */
1847     bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1848
1849     for (g = 0; g < pull->ngroup; g++)
1850     {
1851         int a;
1852
1853         make_local_pull_group(ga2la, &pull->group[g],
1854                               0, md->homenr);
1855
1856         GMX_ASSERT(bMustParticipate || dd != nullptr, "Either all ranks (including this rank) participate, or we use DD and need to have access to dd here");
1857
1858         /* We should participate if we have pull or pbc atoms */
1859         if (!bMustParticipate &&
1860             (pull->group[g].nat_loc > 0 ||
1861              (pull->group[g].params.pbcatom >= 0 &&
1862               ga2la_get_home(dd->ga2la, pull->group[g].params.pbcatom, &a))))
1863         {
1864             bMustParticipate = TRUE;
1865         }
1866     }
1867
1868     if (!comm->bParticipateAll)
1869     {
1870         /* Keep currently not required ranks in the communicator
1871          * if they needed to participate up to 20 decompositions ago.
1872          * This avoids frequent rebuilds due to atoms jumping back and forth.
1873          */
1874         const gmx_int64_t history_count = 20;
1875         gmx_bool          bWillParticipate;
1876         int               count[2];
1877
1878         /* Increase the decomposition counter for the current call */
1879         comm->setup_count++;
1880
1881         if (bMustParticipate)
1882         {
1883             comm->must_count = comm->setup_count;
1884         }
1885
1886         bWillParticipate =
1887             bMustParticipate ||
1888             (comm->bParticipate &&
1889              comm->must_count >= comm->setup_count - history_count);
1890
1891         if (debug && dd != nullptr)
1892         {
1893             fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %d, will be part %d\n",
1894                     dd->rank, bMustParticipate, bWillParticipate);
1895         }
1896
1897         if (bWillParticipate)
1898         {
1899             /* Count the number of ranks that we want to have participating */
1900             count[0] = 1;
1901             /* Count the number of ranks that need to be added */
1902             count[1] = comm->bParticipate ? 0 : 1;
1903         }
1904         else
1905         {
1906             count[0] = 0;
1907             count[1] = 0;
1908         }
1909
1910         /* The cost of this global operation will be less that the cost
1911          * of the extra MPI_Comm_split calls that we can avoid.
1912          */
1913         gmx_sumi(2, count, cr);
1914
1915         /* If we are missing ranks or if we have 20% more ranks than needed
1916          * we make a new sub-communicator.
1917          */
1918         if (count[1] > 0 || 6*count[0] < 5*comm->nparticipate)
1919         {
1920             if (debug)
1921             {
1922                 fprintf(debug, "Creating new pull subcommunicator of size %d\n",
1923                         count[0]);
1924             }
1925
1926 #if GMX_MPI
1927             if (comm->mpi_comm_com != MPI_COMM_NULL)
1928             {
1929                 MPI_Comm_free(&comm->mpi_comm_com);
1930             }
1931
1932             /* This might be an extremely expensive operation, so we try
1933              * to avoid this splitting as much as possible.
1934              */
1935             assert(dd != nullptr);
1936             MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank,
1937                            &comm->mpi_comm_com);
1938 #endif
1939
1940             comm->bParticipate = bWillParticipate;
1941             comm->nparticipate = count[0];
1942         }
1943     }
1944
1945     /* Since the PBC of atoms might have changed, we need to update the PBC */
1946     pull->bSetPBCatoms = TRUE;
1947 }
1948
1949 static void init_pull_group_index(FILE *fplog, const t_commrec *cr,
1950                                   int g, pull_group_work_t *pg,
1951                                   gmx_bool bConstraint, ivec pulldim_con,
1952                                   const gmx_mtop_t *mtop,
1953                                   const t_inputrec *ir, real lambda)
1954 {
1955     if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1956     {
1957         /* There are no masses in the integrator.
1958          * But we still want to have the correct mass-weighted COMs.
1959          * So we store the real masses in the weights.
1960          * We do not set nweight, so these weights do not end up in the tpx file.
1961          */
1962         if (pg->params.nweight == 0)
1963         {
1964             snew(pg->params.weight, pg->params.nat);
1965         }
1966     }
1967
1968     if (cr && PAR(cr))
1969     {
1970         pg->nat_loc    = 0;
1971         pg->nalloc_loc = 0;
1972         pg->ind_loc    = nullptr;
1973         pg->weight_loc = nullptr;
1974     }
1975     else
1976     {
1977         pg->nat_loc = pg->params.nat;
1978         pg->ind_loc = pg->params.ind;
1979         if (pg->epgrppbc == epgrppbcCOS)
1980         {
1981             snew(pg->weight_loc, pg->params.nat);
1982         }
1983         else
1984         {
1985             pg->weight_loc = pg->params.weight;
1986         }
1987     }
1988
1989     const gmx_groups_t *groups = &mtop->groups;
1990
1991     /* Count frozen dimensions and (weighted) mass */
1992     int    nfrozen = 0;
1993     double tmass   = 0;
1994     double wmass   = 0;
1995     double wwmass  = 0;
1996     int    molb    = 0;
1997     for (int i = 0; i < pg->params.nat; i++)
1998     {
1999         int ii = pg->params.ind[i];
2000         if (bConstraint && ir->opts.nFreeze)
2001         {
2002             for (int d = 0; d < DIM; d++)
2003             {
2004                 if (pulldim_con[d] == 1 &&
2005                     ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
2006                 {
2007                     nfrozen++;
2008                 }
2009             }
2010         }
2011         const t_atom &atom = mtopGetAtomParameters(mtop, ii, &molb);
2012         real          m;
2013         if (ir->efep == efepNO)
2014         {
2015             m = atom.m;
2016         }
2017         else
2018         {
2019             m = (1 - lambda)*atom.m + lambda*atom.mB;
2020         }
2021         real w;
2022         if (pg->params.nweight > 0)
2023         {
2024             w = pg->params.weight[i];
2025         }
2026         else
2027         {
2028             w = 1;
2029         }
2030         if (EI_ENERGY_MINIMIZATION(ir->eI))
2031         {
2032             /* Move the mass to the weight */
2033             w                   *= m;
2034             m                    = 1;
2035             pg->params.weight[i] = w;
2036         }
2037         else if (ir->eI == eiBD)
2038         {
2039             real mbd;
2040             if (ir->bd_fric)
2041             {
2042                 mbd = ir->bd_fric*ir->delta_t;
2043             }
2044             else
2045             {
2046                 if (groups->grpnr[egcTC] == nullptr)
2047                 {
2048                     mbd = ir->delta_t/ir->opts.tau_t[0];
2049                 }
2050                 else
2051                 {
2052                     mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
2053                 }
2054             }
2055             w                   *= m/mbd;
2056             m                    = mbd;
2057             pg->params.weight[i] = w;
2058         }
2059         tmass  += m;
2060         wmass  += m*w;
2061         wwmass += m*w*w;
2062     }
2063
2064     if (wmass == 0)
2065     {
2066         /* We can have single atom groups with zero mass with potential pulling
2067          * without cosine weighting.
2068          */
2069         if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
2070         {
2071             /* With one atom the mass doesn't matter */
2072             wwmass = 1;
2073         }
2074         else
2075         {
2076             gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
2077                       pg->params.weight ? " weighted" : "", g);
2078         }
2079     }
2080     if (fplog)
2081     {
2082         fprintf(fplog,
2083                 "Pull group %d: %5d atoms, mass %9.3f",
2084                 g, pg->params.nat, tmass);
2085         if (pg->params.weight ||
2086             EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
2087         {
2088             fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
2089         }
2090         if (pg->epgrppbc == epgrppbcCOS)
2091         {
2092             fprintf(fplog, ", cosine weighting will be used");
2093         }
2094         fprintf(fplog, "\n");
2095     }
2096
2097     if (nfrozen == 0)
2098     {
2099         /* A value != 0 signals not frozen, it is updated later */
2100         pg->invtm  = -1.0;
2101     }
2102     else
2103     {
2104         int ndim = 0;
2105         for (int d = 0; d < DIM; d++)
2106         {
2107             ndim += pulldim_con[d]*pg->params.nat;
2108         }
2109         if (fplog && nfrozen > 0 && nfrozen < ndim)
2110         {
2111             fprintf(fplog,
2112                     "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
2113                     "         that are subject to pulling are frozen.\n"
2114                     "         For constraint pulling the whole group will be frozen.\n\n",
2115                     g);
2116         }
2117         pg->invtm  = 0.0;
2118         pg->wscale = 1.0;
2119     }
2120 }
2121
2122 struct pull_t *
2123 init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
2124           const gmx_mtop_t *mtop, const t_commrec *cr,
2125           real lambda)
2126 {
2127     struct pull_t *pull;
2128     pull_comm_t   *comm;
2129
2130     pull = new pull_t;
2131
2132     /* Copy the pull parameters */
2133     pull->params       = *pull_params;
2134     /* Avoid pointer copies */
2135     pull->params.group = nullptr;
2136     pull->params.coord = nullptr;
2137
2138     pull->ngroup       = pull_params->ngroup;
2139     snew(pull->group, pull->ngroup);
2140
2141     pull->bPotential  = FALSE;
2142     pull->bConstraint = FALSE;
2143     pull->bCylinder   = FALSE;
2144     pull->bAngle      = FALSE;
2145
2146     for (int g = 0; g < pull->ngroup; g++)
2147     {
2148         pull_group_work_t *pgrp = &pull->group[g];
2149
2150         /* Copy the pull group parameters */
2151         pgrp->params = pull_params->group[g];
2152
2153         if (g == 0)
2154         {
2155             GMX_RELEASE_ASSERT(pgrp->params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms");
2156             clear_dvec(pgrp->x);
2157             clear_dvec(pgrp->xp);
2158             pgrp->bCalcCOM = FALSE;
2159         }
2160
2161         /* Avoid pointer copies by allocating and copying arrays */
2162         snew(pgrp->params.ind, pgrp->params.nat);
2163         for (int i = 0; i < pgrp->params.nat; i++)
2164         {
2165             pgrp->params.ind[i] = pull_params->group[g].ind[i];
2166         }
2167         if (pgrp->params.nweight > 0)
2168         {
2169             snew(pgrp->params.weight, pgrp->params.nweight);
2170             for (int i = 0; i < pgrp->params.nweight; i++)
2171             {
2172                 pgrp->params.weight[i] = pull_params->group[g].weight[i];
2173             }
2174         }
2175     }
2176
2177     pull->numCoordinatesWithExternalPotential = 0;
2178
2179     for (int c = 0; c < pull->params.ncoord; c++)
2180     {
2181         /* Construct a pull coordinate, copying all coordinate parameters */
2182         pull->coord.emplace_back(pull_params->coord[c]);
2183
2184         pull_coord_work_t *pcrd = &pull->coord.back();
2185
2186         switch (pcrd->params.eGeom)
2187         {
2188             case epullgDIST:
2189             case epullgDIRRELATIVE:  /* Direction vector is determined at each step */
2190             case epullgANGLE:
2191             case epullgDIHEDRAL:
2192                 break;
2193             case epullgDIR:
2194             case epullgDIRPBC:
2195             case epullgCYL:
2196             case epullgANGLEAXIS:
2197                 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2198                 break;
2199             default:
2200                 /* We allow reading of newer tpx files with new pull geometries,
2201                  * but with the same tpx format, with old code. A new geometry
2202                  * only adds a new enum value, which will be out of range for
2203                  * old code. The first place we need to generate an error is
2204                  * here, since the pull code can't handle this.
2205                  * The enum can be used before arriving here only for printing
2206                  * the string corresponding to the geometry, which will result
2207                  * in printing "UNDEFINED".
2208                  */
2209                 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.",
2210                           c + 1, pcrd->params.eGeom, epullgNR - 1);
2211         }
2212
2213         if (pcrd->params.eType == epullCONSTRAINT)
2214         {
2215             /* Check restrictions of the constraint pull code */
2216             if (pcrd->params.eGeom == epullgCYL ||
2217                 pcrd->params.eGeom == epullgDIRRELATIVE ||
2218                 pcrd->params.eGeom == epullgANGLE ||
2219                 pcrd->params.eGeom == epullgDIHEDRAL ||
2220                 pcrd->params.eGeom == epullgANGLEAXIS)
2221             {
2222                 gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
2223                           epull_names[pcrd->params.eType],
2224                           epullg_names[pcrd->params.eGeom],
2225                           epull_names[epullUMBRELLA]);
2226             }
2227
2228             pull->bConstraint = TRUE;
2229         }
2230         else
2231         {
2232             pull->bPotential = TRUE;
2233         }
2234
2235         if (pcrd->params.eGeom == epullgCYL)
2236         {
2237             pull->bCylinder = TRUE;
2238         }
2239         else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL || pcrd->params.eGeom == epullgANGLEAXIS)
2240         {
2241             pull->bAngle = TRUE;
2242         }
2243
2244         /* We only need to calculate the plain COM of a group
2245          * when it is not only used as a cylinder group.
2246          */
2247         int groupRangeStart = (pcrd->params.eGeom == epullgCYL ? 1 : 0);
2248         int groupRangeEnd   = pcrd->params.ngroup + 1;
2249
2250         for (int i = groupRangeStart; i < groupRangeEnd; i++)
2251         {
2252             int groupIndex = pcrd->params.group[i];
2253             if (groupIndex > 0)
2254             {
2255                 pull->group[groupIndex].bCalcCOM = TRUE;
2256             }
2257         }
2258
2259         /* With non-zero rate the reference value is set at every step */
2260         if (pcrd->params.rate == 0)
2261         {
2262             /* Initialize the constant reference value */
2263             if (pcrd->params.eType != epullEXTERNAL)
2264             {
2265                 low_set_pull_coord_reference_value(pcrd, c, pcrd->params.init*pull_conversion_factor_userinput2internal(&pcrd->params));
2266             }
2267             else
2268             {
2269                 /* With an external pull potential, the reference value loses
2270                  * it's meaning and should not be used. Setting it to zero
2271                  * makes any terms dependent on it disappear.
2272                  * The only issue this causes is that with cylinder pulling
2273                  * no atoms of the cylinder group within the cylinder radius
2274                  * should be more than half a box length away from the COM of
2275                  * the pull group along the axial direction.
2276                  */
2277                 pcrd->value_ref = 0.0;
2278             }
2279         }
2280
2281         if (pcrd->params.eType == epullEXTERNAL)
2282         {
2283             GMX_RELEASE_ASSERT(pcrd->params.rate == 0, "With an external potential, a pull coordinate should have rate = 0");
2284
2285             /* This external potential needs to be registered later */
2286             pull->numCoordinatesWithExternalPotential++;
2287         }
2288         pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2289     }
2290
2291     pull->numUnregisteredExternalPotentials             =
2292         pull->numCoordinatesWithExternalPotential;
2293     pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2294
2295     pull->ePBC = ir->ePBC;
2296     switch (pull->ePBC)
2297     {
2298         case epbcNONE: pull->npbcdim = 0; break;
2299         case epbcXY:   pull->npbcdim = 2; break;
2300         default:       pull->npbcdim = 3; break;
2301     }
2302
2303     if (fplog)
2304     {
2305         gmx_bool bAbs, bCos;
2306
2307         bAbs = FALSE;
2308         for (const pull_coord_work_t &coord : pull->coord)
2309         {
2310             if (pull->group[coord.params.group[0]].params.nat == 0 ||
2311                 pull->group[coord.params.group[1]].params.nat == 0)
2312             {
2313                 bAbs = TRUE;
2314             }
2315         }
2316
2317         fprintf(fplog, "\n");
2318         if (pull->bPotential)
2319         {
2320             fprintf(fplog, "Will apply potential COM pulling\n");
2321         }
2322         if (pull->bConstraint)
2323         {
2324             fprintf(fplog, "Will apply constraint COM pulling\n");
2325         }
2326         // Don't include the reference group 0 in output, so we report ngroup-1
2327         GMX_RELEASE_ASSERT(pull->ngroup - 1 > 0, "The reference absolute position pull group should always be present");
2328         fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n",
2329                 pull->coord.size(), pull->coord.size() == 1 ? "" : "s",
2330                 (pull->ngroup - 1), (pull->ngroup - 1) == 1 ? "" : "s");
2331         if (bAbs)
2332         {
2333             fprintf(fplog, "with an absolute reference\n");
2334         }
2335         bCos = FALSE;
2336         // Don't include the reference group 0 in loop
2337         for (int g = 1; g < pull->ngroup; g++)
2338         {
2339             if (pull->group[g].params.nat > 1 &&
2340                 pull->group[g].params.pbcatom < 0)
2341             {
2342                 /* We are using cosine weighting */
2343                 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
2344                 bCos = TRUE;
2345             }
2346         }
2347         if (bCos)
2348         {
2349             please_cite(fplog, "Engin2010");
2350         }
2351     }
2352
2353     pull->bRefAt   = FALSE;
2354     pull->cosdim   = -1;
2355     for (int g = 0; g < pull->ngroup; g++)
2356     {
2357         pull_group_work_t *pgrp;
2358
2359         pgrp           = &pull->group[g];
2360         pgrp->epgrppbc = epgrppbcNONE;
2361         if (pgrp->params.nat > 0)
2362         {
2363             /* There is an issue when a group is used in multiple coordinates
2364              * and constraints are applied in different dimensions with atoms
2365              * frozen in some, but not all dimensions.
2366              * Since there is only one mass array per group, we can't have
2367              * frozen/non-frozen atoms for different coords at the same time.
2368              * But since this is a very exotic case, we don't check for this.
2369              * A warning is printed in init_pull_group_index.
2370              */
2371
2372             gmx_bool bConstraint;
2373             ivec     pulldim, pulldim_con;
2374
2375             /* Loop over all pull coordinates to see along which dimensions
2376              * this group is pulled and if it is involved in constraints.
2377              */
2378             bConstraint = FALSE;
2379             clear_ivec(pulldim);
2380             clear_ivec(pulldim_con);
2381             for (const pull_coord_work_t &coord : pull->coord)
2382             {
2383                 gmx_bool bGroupUsed = FALSE;
2384                 for (int gi = 0; gi < coord.params.ngroup; gi++)
2385                 {
2386                     if (coord.params.group[gi] == g)
2387                     {
2388                         bGroupUsed = TRUE;
2389                     }
2390                 }
2391
2392                 if (bGroupUsed)
2393                 {
2394                     for (int m = 0; m < DIM; m++)
2395                     {
2396                         if (coord.params.dim[m] == 1)
2397                         {
2398                             pulldim[m] = 1;
2399
2400                             if (coord.params.eType == epullCONSTRAINT)
2401                             {
2402                                 bConstraint    = TRUE;
2403                                 pulldim_con[m] = 1;
2404                             }
2405                         }
2406                     }
2407                 }
2408             }
2409
2410             /* Determine if we need to take PBC into account for calculating
2411              * the COM's of the pull groups.
2412              */
2413             GMX_RELEASE_ASSERT(pull->npbcdim <= DIM, "npbcdim must be <= the number of dimensions");
2414             for (int m = 0; m < pull->npbcdim; m++)
2415             {
2416                 GMX_RELEASE_ASSERT(m <= DIM, "npbcdim must be <= the number of dimensions");
2417                 if (pulldim[m] == 1 && pgrp->params.nat > 1)
2418                 {
2419                     if (pgrp->params.pbcatom >= 0)
2420                     {
2421                         pgrp->epgrppbc = epgrppbcREFAT;
2422                         pull->bRefAt   = TRUE;
2423                     }
2424                     else
2425                     {
2426                         if (pgrp->params.weight != nullptr)
2427                         {
2428                             gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
2429                         }
2430                         pgrp->epgrppbc = epgrppbcCOS;
2431                         if (pull->cosdim >= 0 && pull->cosdim != m)
2432                         {
2433                             gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
2434                         }
2435                         pull->cosdim = m;
2436                     }
2437                 }
2438             }
2439             /* Set the indices */
2440             init_pull_group_index(fplog, cr, g, pgrp,
2441                                   bConstraint, pulldim_con,
2442                                   mtop, ir, lambda);
2443         }
2444         else
2445         {
2446             /* Absolute reference, set the inverse mass to zero.
2447              * This is only relevant (and used) with constraint pulling.
2448              */
2449             pgrp->invtm  = 0;
2450             pgrp->wscale = 1;
2451         }
2452     }
2453
2454     /* If we use cylinder coordinates, do some initialising for them */
2455     if (pull->bCylinder)
2456     {
2457         snew(pull->dyna, pull->coord.size());
2458
2459         for (const pull_coord_work_t &coord : pull->coord)
2460         {
2461             if (coord.params.eGeom == epullgCYL)
2462             {
2463                 if (pull->group[coord.params.group[0]].params.nat == 0)
2464                 {
2465                     gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
2466                 }
2467             }
2468         }
2469     }
2470     else
2471     {
2472         pull->dyna = nullptr;
2473     }
2474
2475     /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2476     pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
2477     snew(pull->sum_com, pull->nthreads);
2478
2479     comm = &pull->comm;
2480
2481 #if GMX_MPI
2482     /* Use a sub-communicator when we have more than 32 ranks, but not
2483      * when we have an external pull potential, since then the external
2484      * potential provider expects each rank to have the coordinate.
2485      */
2486     comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) ||
2487                              cr->dd->nnodes <= 32 ||
2488                              pull->numCoordinatesWithExternalPotential > 0 ||
2489                              getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2490     /* This sub-commicator is not used with comm->bParticipateAll,
2491      * so we can always initialize it to NULL.
2492      */
2493     comm->mpi_comm_com    = MPI_COMM_NULL;
2494     comm->nparticipate    = 0;
2495     comm->isMasterRank    = (cr == nullptr || MASTER(cr));
2496 #else
2497     /* No MPI: 1 rank: all ranks pull */
2498     comm->bParticipateAll = TRUE;
2499     comm->isMasterRank    = true;
2500 #endif
2501     comm->bParticipate    = comm->bParticipateAll;
2502     comm->setup_count     = 0;
2503     comm->must_count      = 0;
2504
2505     if (!comm->bParticipateAll && fplog != nullptr)
2506     {
2507         fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2508     }
2509
2510     comm->rbuf     = nullptr;
2511     comm->dbuf     = nullptr;
2512     comm->dbuf_cyl = nullptr;
2513
2514     /* We still need to initialize the PBC reference coordinates */
2515     pull->bSetPBCatoms = TRUE;
2516
2517     pull->out_x = nullptr;
2518     pull->out_f = nullptr;
2519
2520     return pull;
2521 }
2522
2523 void init_pull_output_files(pull_t                    *pull,
2524                             int                        nfile,
2525                             const t_filenm             fnm[],
2526                             const gmx_output_env_t    *oenv,
2527                             const ContinuationOptions &continuationOptions)
2528 {
2529     /* Check for px and pf filename collision, if we are writing
2530        both files */
2531     std::string px_filename, pf_filename;
2532     std::string px_appended, pf_appended;
2533     try
2534     {
2535         px_filename  = std::string(opt2fn("-px", nfile, fnm));
2536         pf_filename  = std::string(opt2fn("-pf", nfile, fnm));
2537     }
2538     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2539
2540     if ((pull->params.nstxout != 0) &&
2541         (pull->params.nstfout != 0) &&
2542         (px_filename == pf_filename))
2543     {
2544         if (!opt2bSet("-px", nfile, fnm) && !opt2bSet("-pf", nfile, fnm))
2545         {
2546             /* We are writing both pull files but neither set directly. */
2547             try
2548             {
2549                 px_appended   = append_before_extension(px_filename, "_pullx");
2550                 pf_appended   = append_before_extension(pf_filename, "_pullf");
2551             }
2552             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2553             pull->out_x = open_pull_out(px_appended.c_str(), pull, oenv,
2554                                         TRUE, continuationOptions);
2555             pull->out_f = open_pull_out(pf_appended.c_str(), pull, oenv,
2556                                         FALSE, continuationOptions);
2557             return;
2558         }
2559         else
2560         {
2561             /* If at least one of -px and -pf is set but the filenames are identical: */
2562             gmx_fatal(FARGS, "Identical pull_x and pull_f output filenames %s",
2563                       px_filename.c_str());
2564         }
2565     }
2566     if (pull->params.nstxout != 0)
2567     {
2568         pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
2569                                     TRUE, continuationOptions);
2570     }
2571     if (pull->params.nstfout != 0)
2572     {
2573         pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
2574                                     FALSE, continuationOptions);
2575     }
2576 }
2577
2578 static void destroy_pull_group(pull_group_work_t *pgrp)
2579 {
2580     if (pgrp->ind_loc != pgrp->params.ind)
2581     {
2582         sfree(pgrp->ind_loc);
2583     }
2584     if (pgrp->weight_loc != pgrp->params.weight)
2585     {
2586         sfree(pgrp->weight_loc);
2587     }
2588     sfree(pgrp->mdw);
2589     sfree(pgrp->dv);
2590
2591     sfree(pgrp->params.ind);
2592     sfree(pgrp->params.weight);
2593 }
2594
2595 static void destroy_pull(struct pull_t *pull)
2596 {
2597     for (int i = 0; i < pull->ngroup; i++)
2598     {
2599         destroy_pull_group(&pull->group[i]);
2600     }
2601     for (size_t c = 0; c < pull->coord.size(); c++)
2602     {
2603         if (pull->coord[c].params.eGeom == epullgCYL)
2604         {
2605             destroy_pull_group(&(pull->dyna[c]));
2606         }
2607     }
2608     sfree(pull->group);
2609     sfree(pull->dyna);
2610
2611 #if GMX_MPI
2612     if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2613     {
2614         MPI_Comm_free(&pull->comm.mpi_comm_com);
2615     }
2616 #endif
2617     sfree(pull->comm.rbuf);
2618     sfree(pull->comm.dbuf);
2619     sfree(pull->comm.dbuf_cyl);
2620
2621     delete pull;
2622 }
2623
2624 void finish_pull(struct pull_t *pull)
2625 {
2626     check_external_potential_registration(pull);
2627
2628     if (pull->out_x)
2629     {
2630         gmx_fio_fclose(pull->out_x);
2631     }
2632     if (pull->out_f)
2633     {
2634         gmx_fio_fclose(pull->out_f);
2635     }
2636
2637     destroy_pull(pull);
2638 }
2639
2640 gmx_bool pull_have_potential(const struct pull_t *pull)
2641 {
2642     return pull->bPotential;
2643 }
2644
2645 gmx_bool pull_have_constraint(const struct pull_t *pull)
2646 {
2647     return pull->bConstraint;
2648 }