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