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