2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
52 #include "gromacs/commandline/filenm.h"
53 #include "gromacs/compat/make_unique.h"
54 #include "gromacs/domdec/domdec_struct.h"
55 #include "gromacs/domdec/localatomset.h"
56 #include "gromacs/domdec/localatomsetmanager.h"
57 #include "gromacs/fileio/gmxfio.h"
58 #include "gromacs/fileio/xvgr.h"
59 #include "gromacs/gmxlib/network.h"
60 #include "gromacs/math/functions.h"
61 #include "gromacs/math/units.h"
62 #include "gromacs/math/vec.h"
63 #include "gromacs/math/vectypes.h"
64 #include "gromacs/mdlib/gmx_omp_nthreads.h"
65 #include "gromacs/mdlib/mdrun.h"
66 #include "gromacs/mdtypes/commrec.h"
67 #include "gromacs/mdtypes/forceoutput.h"
68 #include "gromacs/mdtypes/inputrec.h"
69 #include "gromacs/mdtypes/md_enums.h"
70 #include "gromacs/mdtypes/mdatom.h"
71 #include "gromacs/pbcutil/pbc.h"
72 #include "gromacs/topology/mtop_lookup.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/exceptions.h"
76 #include "gromacs/utility/fatalerror.h"
77 #include "gromacs/utility/futil.h"
78 #include "gromacs/utility/gmxassert.h"
79 #include "gromacs/utility/mutex.h"
80 #include "gromacs/utility/pleasecite.h"
81 #include "gromacs/utility/real.h"
82 #include "gromacs/utility/smalloc.h"
83 #include "gromacs/utility/strconvert.h"
85 #include "pull_internal.h"
87 static int groupPbcFromParams(const t_pull_group ¶ms)
94 else if (params.pbcatom >= 0)
104 /* NOTE: The params initialization currently copies pointers.
105 * So the lifetime of the source, currently always inputrec,
106 * should not end before that of this object.
107 * This will be fixed when the pointers are replacted by std::vector.
109 pull_group_work_t::pull_group_work_t(const t_pull_group ¶ms,
110 gmx::LocalAtomSet atomSet) :
112 epgrppbc(groupPbcFromParams(params)),
113 needToCalcCom(false),
123 bool pull_coordinate_is_angletype(const t_pull_coord *pcrd)
125 return (pcrd->eGeom == epullgANGLE ||
126 pcrd->eGeom == epullgDIHEDRAL ||
127 pcrd->eGeom == epullgANGLEAXIS);
130 static bool pull_coordinate_is_directional(const t_pull_coord *pcrd)
132 return (pcrd->eGeom == epullgDIR ||
133 pcrd->eGeom == epullgDIRPBC ||
134 pcrd->eGeom == epullgDIRRELATIVE ||
135 pcrd->eGeom == epullgCYL);
138 const char *pull_coordinate_units(const t_pull_coord *pcrd)
140 return pull_coordinate_is_angletype(pcrd) ? "deg" : "nm";
143 double pull_conversion_factor_userinput2internal(const t_pull_coord *pcrd)
145 if (pull_coordinate_is_angletype(pcrd))
155 double pull_conversion_factor_internal2userinput(const t_pull_coord *pcrd)
157 if (pull_coordinate_is_angletype(pcrd))
167 static std::string append_before_extension(const std::string &pathname,
168 const std::string &to_append)
170 /* Appends to_append before last '.' in pathname */
171 size_t extPos = pathname.find_last_of('.');
172 if (extPos == std::string::npos)
174 return pathname + to_append;
178 return pathname.substr(0, extPos) + to_append +
179 pathname.substr(extPos, std::string::npos);
183 static void pull_print_group_x(FILE *out, const ivec dim,
184 const pull_group_work_t *pgrp)
188 for (m = 0; m < DIM; m++)
192 fprintf(out, "\t%g", pgrp->x[m]);
197 static void pull_print_coord_dr_components(FILE *out, const ivec dim, const dvec dr)
199 for (int m = 0; m < DIM; m++)
203 fprintf(out, "\t%g", dr[m]);
208 static void pull_print_coord_dr(FILE *out, const pull_coord_work_t *pcrd,
209 gmx_bool bPrintRefValue,
210 gmx_bool bPrintComponents)
212 double unit_factor = pull_conversion_factor_internal2userinput(&pcrd->params);
214 fprintf(out, "\t%g", pcrd->spatialData.value*unit_factor);
218 fprintf(out, "\t%g", pcrd->value_ref*unit_factor);
221 if (bPrintComponents)
223 pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr01);
224 if (pcrd->params.ngroup >= 4)
226 pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr23);
228 if (pcrd->params.ngroup >= 6)
230 pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr45);
235 static void pull_print_x(FILE *out, struct pull_t *pull, double t)
237 fprintf(out, "%.4f", t);
239 for (size_t c = 0; c < pull->coord.size(); c++)
241 pull_coord_work_t *pcrd;
243 pcrd = &pull->coord[c];
245 pull_print_coord_dr(out, pcrd,
246 pull->params.bPrintRefValue && pcrd->params.eType != epullEXTERNAL,
247 pull->params.bPrintComp);
249 if (pull->params.bPrintCOM)
251 if (pcrd->params.eGeom == epullgCYL)
253 pull_print_group_x(out, pcrd->params.dim, &pull->dyna[c]);
257 pull_print_group_x(out, pcrd->params.dim,
258 &pull->group[pcrd->params.group[0]]);
260 for (int g = 1; g < pcrd->params.ngroup; g++)
262 pull_print_group_x(out, pcrd->params.dim, &pull->group[pcrd->params.group[g]]);
269 static void pull_print_f(FILE *out, const pull_t *pull, double t)
271 fprintf(out, "%.4f", t);
273 for (const pull_coord_work_t &coord : pull->coord)
275 fprintf(out, "\t%g", coord.scalarForce);
280 void pull_print_output(struct pull_t *pull, int64_t step, double time)
282 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0, "pull_print_output called before all external pull potentials have been applied");
284 if ((pull->params.nstxout != 0) && (step % pull->params.nstxout == 0))
286 pull_print_x(pull->out_x, pull, time);
289 if ((pull->params.nstfout != 0) && (step % pull->params.nstfout == 0))
291 pull_print_f(pull->out_f, pull, time);
295 static void set_legend_for_coord_components(const pull_coord_work_t *pcrd, int coord_index, char **setname, int *nsets_ptr)
297 /* Loop over the distance vectors and print their components. Each vector is made up of two consecutive groups. */
298 for (int g = 0; g < pcrd->params.ngroup; g += 2)
300 /* Loop over the components */
301 for (int m = 0; m < DIM; m++)
303 if (pcrd->params.dim[m])
307 if (g == 0 && pcrd->params.ngroup <= 2)
309 /* For the simplest case we print a simplified legend without group indices, just the cooordinate index
310 and which dimensional component it is. */
311 sprintf(legend, "%d d%c", coord_index + 1, 'X' + m);
315 /* Otherwise, print also the group indices (relative to the coordinate, not the global ones). */
316 sprintf(legend, "%d g %d-%d d%c", coord_index + 1, g + 1, g + 2, 'X' + m);
319 setname[*nsets_ptr] = gmx_strdup(legend);
326 static FILE *open_pull_out(const char *fn, struct pull_t *pull,
327 const gmx_output_env_t *oenv,
329 const ContinuationOptions &continuationOptions)
333 char **setname, buf[50];
335 if (continuationOptions.appendFiles)
337 fp = gmx_fio_fopen(fn, "a+");
341 fp = gmx_fio_fopen(fn, "w+");
344 sprintf(buf, "Position (nm%s)", pull->bAngle ? ", deg" : "");
345 xvgr_header(fp, "Pull COM", "Time (ps)", buf,
350 sprintf(buf, "Force (kJ/mol/nm%s)", pull->bAngle ? ", kJ/mol/rad" : "");
351 xvgr_header(fp, "Pull force", "Time (ps)", buf,
355 /* With default mdp options only the actual coordinate value is printed (1),
356 * but optionally the reference value (+ 1),
357 * the group COMs for all the groups (+ ngroups_max*DIM)
358 * and the components of the distance vectors can be printed (+ (ngroups_max/2)*DIM).
360 snew(setname, pull->coord.size()*(1 + 1 + c_pullCoordNgroupMax*DIM + c_pullCoordNgroupMax/2*DIM));
363 for (size_t c = 0; c < pull->coord.size(); c++)
367 /* The order of this legend should match the order of printing
368 * the data in print_pull_x above.
371 /* The pull coord distance */
372 sprintf(buf, "%zu", c+1);
373 setname[nsets] = gmx_strdup(buf);
375 if (pull->params.bPrintRefValue &&
376 pull->coord[c].params.eType != epullEXTERNAL)
378 sprintf(buf, "%zu ref", c+1);
379 setname[nsets] = gmx_strdup(buf);
382 if (pull->params.bPrintComp)
384 set_legend_for_coord_components(&pull->coord[c], c, setname, &nsets);
387 if (pull->params.bPrintCOM)
389 for (int g = 0; g < pull->coord[c].params.ngroup; g++)
391 /* Legend for reference group position */
392 for (m = 0; m < DIM; m++)
394 if (pull->coord[c].params.dim[m])
396 sprintf(buf, "%zu g %d %c", c+1, g + 1, 'X'+m);
397 setname[nsets] = gmx_strdup(buf);
406 /* For the pull force we always only use one scalar */
407 sprintf(buf, "%zu", c+1);
408 setname[nsets] = gmx_strdup(buf);
414 xvgr_legend(fp, nsets, setname, oenv);
416 for (int c = 0; c < nsets; c++)
426 /* Apply forces in a mass weighted fashion for part of the pull group */
427 static void apply_forces_grp_part(const pull_group_work_t *pgrp,
428 int ind_start, int ind_end,
430 const dvec f_pull, int sign, rvec *f)
432 double inv_wm = pgrp->mwscale;
434 auto localAtomIndices = pgrp->atomSet.localIndex();
435 for (int i = ind_start; i < ind_end; i++)
437 int ii = localAtomIndices[i];
438 double wmass = md->massT[ii];
439 if (!pgrp->localWeights.empty())
441 wmass *= pgrp->localWeights[i];
444 for (int d = 0; d < DIM; d++)
446 f[ii][d] += sign * wmass * f_pull[d] * inv_wm;
451 /* Apply forces in a mass weighted fashion */
452 static void apply_forces_grp(const pull_group_work_t *pgrp,
454 const dvec f_pull, int sign, rvec *f,
457 auto localAtomIndices = pgrp->atomSet.localIndex();
459 if (pgrp->params.nat == 1 && pgrp->atomSet.numAtomsLocal() == 1)
461 /* Only one atom and our rank has this atom: we can skip
462 * the mass weighting, which means that this code also works
463 * for mass=0, e.g. with a virtual site.
465 for (int d = 0; d < DIM; d++)
467 f[localAtomIndices[0]][d] += sign*f_pull[d];
472 if (localAtomIndices.size() <= c_pullMaxNumLocalAtomsSingleThreaded)
474 apply_forces_grp_part(pgrp, 0, localAtomIndices.size(), md, f_pull, sign, f);
478 #pragma omp parallel for num_threads(nthreads) schedule(static)
479 for (int th = 0; th < nthreads; th++)
481 int ind_start = (localAtomIndices.size()*(th + 0))/nthreads;
482 int ind_end = (localAtomIndices.size()*(th + 1))/nthreads;
483 apply_forces_grp_part(pgrp, ind_start, ind_end,
484 md, f_pull, sign, f);
490 /* Apply forces in a mass weighted fashion to a cylinder group */
491 static void apply_forces_cyl_grp(const pull_group_work_t *pgrp,
492 const double dv_corr,
494 const dvec f_pull, double f_scal,
496 int gmx_unused nthreads)
498 double inv_wm = pgrp->mwscale;
500 auto localAtomIndices = pgrp->atomSet.localIndex();
502 /* The cylinder group is always a slab in the system, thus large.
503 * Therefore we always thread-parallelize this group.
505 int numAtomsLocal = localAtomIndices.size();
506 #pragma omp parallel for num_threads(nthreads) schedule(static)
507 for (int i = 0; i < numAtomsLocal; i++)
509 double weight = pgrp->localWeights[i];
514 int ii = localAtomIndices[i];
515 double mass = md->massT[ii];
516 /* The stored axial distance from the cylinder center (dv) needs
517 * to be corrected for an offset (dv_corr), which was unknown when
520 double dv_com = pgrp->dv[i] + dv_corr;
522 /* Here we not only add the pull force working along vec (f_pull),
523 * but also a radial component, due to the dependence of the weights
524 * on the radial distance.
526 for (int m = 0; m < DIM; m++)
528 f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
529 pgrp->mdw[i][m]*dv_com*f_scal);
534 /* Apply torque forces in a mass weighted fashion to the groups that define
535 * the pull vector direction for pull coordinate pcrd.
537 static void apply_forces_vec_torque(const struct pull_t *pull,
538 const pull_coord_work_t *pcrd,
542 const PullCoordSpatialData &spatialData = pcrd->spatialData;
544 /* The component inpr along the pull vector is accounted for in the usual
545 * way. Here we account for the component perpendicular to vec.
548 for (int m = 0; m < DIM; m++)
550 inpr += spatialData.dr01[m]*spatialData.vec[m];
552 /* The torque force works along the component of the distance vector
553 * of between the two "usual" pull groups that is perpendicular to
554 * the pull vector. The magnitude of this force is the "usual" scale force
555 * multiplied with the ratio of the distance between the two "usual" pull
556 * groups and the distance between the two groups that define the vector.
559 for (int m = 0; m < DIM; m++)
561 f_perp[m] = (spatialData.dr01[m] - inpr*spatialData.vec[m])/spatialData.vec_len*pcrd->scalarForce;
564 /* Apply the force to the groups defining the vector using opposite signs */
565 apply_forces_grp(&pull->group[pcrd->params.group[2]], md,
566 f_perp, -1, f, pull->nthreads);
567 apply_forces_grp(&pull->group[pcrd->params.group[3]], md,
568 f_perp, 1, f, pull->nthreads);
571 /* Apply forces in a mass weighted fashion */
572 static void apply_forces_coord(struct pull_t * pull, int coord,
573 const PullCoordVectorForces &forces,
574 const t_mdatoms * md,
577 /* Here it would be more efficient to use one large thread-parallel
578 * region instead of potential parallel regions within apply_forces_grp.
579 * But there could be overlap between pull groups and this would lead
583 const pull_coord_work_t &pcrd = pull->coord[coord];
585 if (pcrd.params.eGeom == epullgCYL)
587 apply_forces_cyl_grp(&pull->dyna[coord], pcrd.spatialData.cyl_dev, md,
588 forces.force01, pcrd.scalarForce, -1, f,
591 /* Sum the force along the vector and the radial force */
593 for (int m = 0; m < DIM; m++)
595 f_tot[m] = forces.force01[m] + pcrd.scalarForce*pcrd.spatialData.ffrad[m];
597 apply_forces_grp(&pull->group[pcrd.params.group[1]], md,
598 f_tot, 1, f, pull->nthreads);
602 if (pcrd.params.eGeom == epullgDIRRELATIVE)
604 /* We need to apply the torque forces to the pull groups
605 * that define the pull vector.
607 apply_forces_vec_torque(pull, &pcrd, md, f);
610 if (pull->group[pcrd.params.group[0]].params.nat > 0)
612 apply_forces_grp(&pull->group[pcrd.params.group[0]], md,
613 forces.force01, -1, f, pull->nthreads);
615 apply_forces_grp(&pull->group[pcrd.params.group[1]], md,
616 forces.force01, 1, f, pull->nthreads);
618 if (pcrd.params.ngroup >= 4)
620 apply_forces_grp(&pull->group[pcrd.params.group[2]], md,
621 forces.force23, -1, f, pull->nthreads);
622 apply_forces_grp(&pull->group[pcrd.params.group[3]], md,
623 forces.force23, 1, f, pull->nthreads);
625 if (pcrd.params.ngroup >= 6)
627 apply_forces_grp(&pull->group[pcrd.params.group[4]], md,
628 forces.force45, -1, f, pull->nthreads);
629 apply_forces_grp(&pull->group[pcrd.params.group[5]], md,
630 forces.force45, 1, f, pull->nthreads);
635 real max_pull_distance2(const pull_coord_work_t *pcrd,
638 /* Note that this maximum distance calculation is more complex than
639 * most other cases in GROMACS, since here we have to take care of
640 * distance calculations that don't involve all three dimensions.
641 * For example, we can use distances that are larger than the
642 * box X and Y dimensions for a box that is elongated along Z.
645 real max_d2 = GMX_REAL_MAX;
647 if (pull_coordinate_is_directional(&pcrd->params))
649 /* Directional pulling along along direction pcrd->vec.
650 * Calculating the exact maximum distance is complex and bug-prone.
651 * So we take a safe approach by not allowing distances that
652 * are larger than half the distance between unit cell faces
653 * along dimensions involved in pcrd->vec.
655 for (int m = 0; m < DIM; m++)
657 if (m < pbc->ndim_ePBC && pcrd->spatialData.vec[m] != 0)
659 real imageDistance2 = gmx::square(pbc->box[m][m]);
660 for (int d = m + 1; d < DIM; d++)
662 imageDistance2 -= gmx::square(pbc->box[d][m]);
664 max_d2 = std::min(max_d2, imageDistance2);
670 /* Distance pulling along dimensions with pcrd->params.dim[d]==1.
671 * We use half the minimum box vector length of the dimensions involved.
672 * This is correct for all cases, except for corner cases with
673 * triclinic boxes where e.g. dim[XX]=1 and dim[YY]=0 with
674 * box[YY][XX]!=0 and box[YY][YY] < box[XX][XX]. But since even
675 * in such corner cases the user could get correct results,
676 * depending on the details of the setup, we avoid further
677 * code complications.
679 for (int m = 0; m < DIM; m++)
681 if (m < pbc->ndim_ePBC && pcrd->params.dim[m] != 0)
683 real imageDistance2 = gmx::square(pbc->box[m][m]);
684 for (int d = 0; d < m; d++)
686 if (pcrd->params.dim[d] != 0)
688 imageDistance2 += gmx::square(pbc->box[m][d]);
691 max_d2 = std::min(max_d2, imageDistance2);
699 /* This function returns the distance based on coordinates xg and xref.
700 * Note that the pull coordinate struct pcrd is not modified.
702 static void low_get_pull_coord_dr(const struct pull_t *pull,
703 const pull_coord_work_t *pcrd,
705 dvec xg, dvec xref, double max_dist2,
708 const pull_group_work_t *pgrp0 = &pull->group[pcrd->params.group[0]];
710 /* Only the first group can be an absolute reference, in that case nat=0 */
711 if (pgrp0->params.nat == 0)
713 for (int m = 0; m < DIM; m++)
715 xref[m] = pcrd->params.origin[m];
720 copy_dvec(xref, xrefr);
722 dvec dref = {0, 0, 0};
723 if (pcrd->params.eGeom == epullgDIRPBC)
725 for (int m = 0; m < DIM; m++)
727 dref[m] = pcrd->value_ref*pcrd->spatialData.vec[m];
729 /* Add the reference position, so we use the correct periodic image */
730 dvec_inc(xrefr, dref);
733 pbc_dx_d(pbc, xg, xrefr, dr);
735 bool directional = pull_coordinate_is_directional(&pcrd->params);
737 for (int m = 0; m < DIM; m++)
739 dr[m] *= pcrd->params.dim[m];
740 if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
745 /* Check if we are close to switching to another periodic image */
746 if (max_dist2 > 0 && dr2 > 0.98*0.98*max_dist2)
748 /* Note that technically there is no issue with switching periodic
749 * image, as pbc_dx_d returns the distance to the closest periodic
750 * image. However in all cases where periodic image switches occur,
751 * the pull results will be useless in practice.
753 gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f).\n%s",
754 pcrd->params.group[0], pcrd->params.group[1],
755 sqrt(dr2), sqrt(0.98*0.98*max_dist2),
756 pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
759 if (pcrd->params.eGeom == epullgDIRPBC)
765 /* This function returns the distance based on the contents of the pull struct.
766 * pull->coord[coord_ind].dr, and potentially vector, are updated.
768 static void get_pull_coord_dr(struct pull_t *pull,
772 pull_coord_work_t *pcrd = &pull->coord[coord_ind];
773 PullCoordSpatialData &spatialData = pcrd->spatialData;
776 if (pcrd->params.eGeom == epullgDIRPBC)
782 md2 = static_cast<double>(max_pull_distance2(pcrd, pbc));
785 if (pcrd->params.eGeom == epullgDIRRELATIVE)
787 /* We need to determine the pull vector */
791 const pull_group_work_t &pgrp2 = pull->group[pcrd->params.group[2]];
792 const pull_group_work_t &pgrp3 = pull->group[pcrd->params.group[3]];
794 pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
796 for (m = 0; m < DIM; m++)
798 vec[m] *= pcrd->params.dim[m];
800 spatialData.vec_len = dnorm(vec);
801 for (m = 0; m < DIM; m++)
803 spatialData.vec[m] = vec[m]/spatialData.vec_len;
807 fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
809 vec[XX], vec[YY], vec[ZZ],
810 spatialData.vec[XX], spatialData.vec[YY], spatialData.vec[ZZ]);
814 pull_group_work_t *pgrp0 = &pull->group[pcrd->params.group[0]];
815 pull_group_work_t *pgrp1 = &pull->group[pcrd->params.group[1]];
817 low_get_pull_coord_dr(pull, pcrd, pbc,
819 pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x,
823 if (pcrd->params.ngroup >= 4)
825 pull_group_work_t *pgrp2, *pgrp3;
826 pgrp2 = &pull->group[pcrd->params.group[2]];
827 pgrp3 = &pull->group[pcrd->params.group[3]];
829 low_get_pull_coord_dr(pull, pcrd, pbc,
835 if (pcrd->params.ngroup >= 6)
837 pull_group_work_t *pgrp4, *pgrp5;
838 pgrp4 = &pull->group[pcrd->params.group[4]];
839 pgrp5 = &pull->group[pcrd->params.group[5]];
841 low_get_pull_coord_dr(pull, pcrd, pbc,
849 /* Modify x so that it is periodic in [-pi, pi)
850 * It is assumed that x is in [-3pi, 3pi) so that x
851 * needs to be shifted by at most one period.
853 static void make_periodic_2pi(double *x)
865 /* This function should always be used to modify pcrd->value_ref */
866 static void low_set_pull_coord_reference_value(pull_coord_work_t *pcrd,
870 GMX_ASSERT(pcrd->params.eType != epullEXTERNAL, "The pull coord reference value should not be used with type external-potential");
872 if (pcrd->params.eGeom == epullgDIST)
876 gmx_fatal(FARGS, "Pull reference distance for coordinate %d (%f) needs to be non-negative", coord_ind + 1, value_ref);
879 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgANGLEAXIS)
881 if (value_ref < 0 || value_ref > M_PI)
883 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));
886 else if (pcrd->params.eGeom == epullgDIHEDRAL)
888 /* Allow pulling to be periodic for dihedral angles by remapping the reference value to the interval [-pi, pi). */
889 make_periodic_2pi(&value_ref);
892 pcrd->value_ref = value_ref;
895 static void update_pull_coord_reference_value(pull_coord_work_t *pcrd, int coord_ind, double t)
897 /* With zero rate the reference value is set initially and doesn't change */
898 if (pcrd->params.rate != 0)
900 double value_ref = (pcrd->params.init + pcrd->params.rate*t)*pull_conversion_factor_userinput2internal(&pcrd->params);
901 low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
905 /* Returns the dihedral angle. Updates the plane normal vectors m, n. */
906 static double get_dihedral_angle_coord(PullCoordSpatialData *spatialData)
909 dvec dr32; /* store instead of dr23? */
911 dsvmul(-1, spatialData->dr23, dr32);
912 dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
913 dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
914 phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
916 /* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
917 * because there r_ij = ri - rj, while here dr01 = r_1 - r_0
918 * Note 2: the angle between the plane normal vectors equals pi only when
919 * both planes coincide. Thus, when phi = pi dr01 will lie in both planes and
920 * we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
922 sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
926 /* Calculates pull->coord[coord_ind].value.
927 * This function also updates pull->coord[coord_ind].dr.
929 static void get_pull_coord_distance(struct pull_t *pull,
933 pull_coord_work_t *pcrd;
936 pcrd = &pull->coord[coord_ind];
938 get_pull_coord_dr(pull, coord_ind, pbc);
940 PullCoordSpatialData &spatialData = pcrd->spatialData;
942 switch (pcrd->params.eGeom)
945 /* Pull along the vector between the com's */
946 spatialData.value = dnorm(spatialData.dr01);
950 case epullgDIRRELATIVE:
953 spatialData.value = 0;
954 for (m = 0; m < DIM; m++)
956 spatialData.value += spatialData.vec[m]*spatialData.dr01[m];
960 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
963 spatialData.value = get_dihedral_angle_coord(&spatialData);
965 case epullgANGLEAXIS:
966 spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
969 gmx_incons("Unsupported pull type in get_pull_coord_distance");
973 /* Returns the deviation from the reference value.
974 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
976 static double get_pull_coord_deviation(struct pull_t *pull,
981 pull_coord_work_t *pcrd;
984 pcrd = &pull->coord[coord_ind];
986 /* Update the reference value before computing the distance,
987 * since it is used in the distance computation with periodic pulling.
989 update_pull_coord_reference_value(pcrd, coord_ind, t);
991 get_pull_coord_distance(pull, coord_ind, pbc);
993 get_pull_coord_distance(pull, coord_ind, pbc);
995 /* Determine the deviation */
996 dev = pcrd->spatialData.value - pcrd->value_ref;
998 /* Check that values are allowed */
999 if (pcrd->params.eGeom == epullgDIST && pcrd->spatialData.value == 0)
1001 /* With no vector we can not determine the direction for the force,
1002 * so we set the force to zero.
1006 else if (pcrd->params.eGeom == epullgDIHEDRAL)
1008 /* The reference value is in [-pi, pi). The coordinate value is in (-pi, pi].
1009 Thus, the unwrapped deviation is here in (-2pi, 2pi].
1010 After making it periodic, the deviation will be in [-pi, pi). */
1011 make_periodic_2pi(&dev);
1017 double get_pull_coord_value(struct pull_t *pull,
1021 get_pull_coord_distance(pull, coord_ind, pbc);
1023 return pull->coord[coord_ind].spatialData.value;
1026 void clear_pull_forces(pull_t *pull)
1028 /* Zeroing the forces is only required for constraint pulling.
1029 * It can happen that multiple constraint steps need to be applied
1030 * and therefore the constraint forces need to be accumulated.
1032 for (pull_coord_work_t &coord : pull->coord)
1034 coord.scalarForce = 0;
1038 /* Apply constraint using SHAKE */
1039 static void do_constraint(struct pull_t *pull, t_pbc *pbc,
1041 gmx_bool bMaster, tensor vir,
1042 double dt, double t)
1045 dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
1046 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
1047 dvec *rnew; /* current 'new' positions of the groups */
1048 double *dr_tot; /* the total update of the coords */
1051 double lambda, rm, invdt = 0;
1052 gmx_bool bConverged_all, bConverged = FALSE;
1053 int niter = 0, ii, j, m, max_iter = 100;
1057 snew(r_ij, pull->coord.size());
1058 snew(dr_tot, pull->coord.size());
1060 snew(rnew, pull->group.size());
1062 /* copy the current unconstrained positions for use in iterations. We
1063 iterate until rinew[i] and rjnew[j] obey the constraints. Then
1064 rinew - pull.x_unc[i] is the correction dr to group i */
1065 for (size_t g = 0; g < pull->group.size(); g++)
1067 copy_dvec(pull->group[g].xp, rnew[g]);
1070 /* Determine the constraint directions from the old positions */
1071 for (size_t c = 0; c < pull->coord.size(); c++)
1073 pull_coord_work_t *pcrd;
1075 pcrd = &pull->coord[c];
1077 if (pcrd->params.eType != epullCONSTRAINT)
1082 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
1083 * We don't modify dr and value anymore, so these values are also used
1086 get_pull_coord_distance(pull, c, pbc);
1088 const PullCoordSpatialData &spatialData = pcrd->spatialData;
1091 fprintf(debug, "Pull coord %zu dr %f %f %f\n",
1092 c, spatialData.dr01[XX], spatialData.dr01[YY], spatialData.dr01[ZZ]);
1095 if (pcrd->params.eGeom == epullgDIR ||
1096 pcrd->params.eGeom == epullgDIRPBC)
1098 /* Select the component along vec */
1100 for (m = 0; m < DIM; m++)
1102 a += spatialData.vec[m]*spatialData.dr01[m];
1104 for (m = 0; m < DIM; m++)
1106 r_ij[c][m] = a*spatialData.vec[m];
1111 copy_dvec(spatialData.dr01, r_ij[c]);
1114 if (dnorm2(r_ij[c]) == 0)
1116 gmx_fatal(FARGS, "Distance for pull coordinate %zu is zero with constraint pulling, which is not allowed.", c + 1);
1120 bConverged_all = FALSE;
1121 while (!bConverged_all && niter < max_iter)
1123 bConverged_all = TRUE;
1125 /* loop over all constraints */
1126 for (size_t c = 0; c < pull->coord.size(); c++)
1128 pull_coord_work_t *pcrd;
1129 pull_group_work_t *pgrp0, *pgrp1;
1132 pcrd = &pull->coord[c];
1134 if (pcrd->params.eType != epullCONSTRAINT)
1139 update_pull_coord_reference_value(pcrd, c, t);
1141 pgrp0 = &pull->group[pcrd->params.group[0]];
1142 pgrp1 = &pull->group[pcrd->params.group[1]];
1144 /* Get the current difference vector */
1145 low_get_pull_coord_dr(pull, pcrd, pbc,
1146 rnew[pcrd->params.group[1]],
1147 rnew[pcrd->params.group[0]],
1152 fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
1155 rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
1157 switch (pcrd->params.eGeom)
1160 if (pcrd->value_ref <= 0)
1162 gmx_fatal(FARGS, "The pull constraint reference distance for group %zu is <= 0 (%f)", c, pcrd->value_ref);
1166 double q, c_a, c_b, c_c;
1168 c_a = diprod(r_ij[c], r_ij[c]);
1169 c_b = diprod(unc_ij, r_ij[c])*2;
1170 c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
1174 q = -0.5*(c_b - std::sqrt(c_b*c_b - 4*c_a*c_c));
1179 q = -0.5*(c_b + std::sqrt(c_b*c_b - 4*c_a*c_c));
1186 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
1187 c_a, c_b, c_c, lambda);
1191 /* The position corrections dr due to the constraints */
1192 dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
1193 dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
1194 dr_tot[c] += -lambda*dnorm(r_ij[c]);
1199 /* A 1-dimensional constraint along a vector */
1201 for (m = 0; m < DIM; m++)
1203 vec[m] = pcrd->spatialData.vec[m];
1204 a += unc_ij[m]*vec[m];
1206 /* Select only the component along the vector */
1207 dsvmul(a, vec, unc_ij);
1208 lambda = a - pcrd->value_ref;
1211 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
1214 /* The position corrections dr due to the constraints */
1215 dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
1216 dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
1217 dr_tot[c] += -lambda;
1220 gmx_incons("Invalid enumeration value for eGeom");
1228 g0 = pcrd->params.group[0];
1229 g1 = pcrd->params.group[1];
1230 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
1231 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
1233 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1234 rnew[g0][0], rnew[g0][1], rnew[g0][2],
1235 rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
1237 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
1238 "", "", "", "", "", "", pcrd->value_ref);
1240 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
1241 dr0[0], dr0[1], dr0[2],
1242 dr1[0], dr1[1], dr1[2],
1246 /* Update the COMs with dr */
1247 dvec_inc(rnew[pcrd->params.group[1]], dr1);
1248 dvec_inc(rnew[pcrd->params.group[0]], dr0);
1251 /* Check if all constraints are fullfilled now */
1252 for (pull_coord_work_t &coord : pull->coord)
1254 if (coord.params.eType != epullCONSTRAINT)
1259 low_get_pull_coord_dr(pull, &coord, pbc,
1260 rnew[coord.params.group[1]],
1261 rnew[coord.params.group[0]],
1264 switch (coord.params.eGeom)
1268 fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
1273 for (m = 0; m < DIM; m++)
1275 vec[m] = coord.spatialData.vec[m];
1277 inpr = diprod(unc_ij, vec);
1278 dsvmul(inpr, vec, unc_ij);
1280 fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
1288 fprintf(debug, "Pull constraint not converged: "
1290 "d_ref = %f, current d = %f\n",
1291 coord.params.group[0], coord.params.group[1],
1292 coord.value_ref, dnorm(unc_ij));
1295 bConverged_all = FALSE;
1300 /* if after all constraints are dealt with and bConverged is still TRUE
1301 we're finished, if not we do another iteration */
1303 if (niter > max_iter)
1305 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
1308 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
1315 /* update atoms in the groups */
1316 for (size_t g = 0; g < pull->group.size(); g++)
1318 const pull_group_work_t *pgrp;
1321 pgrp = &pull->group[g];
1323 /* get the final constraint displacement dr for group g */
1324 dvec_sub(rnew[g], pgrp->xp, dr);
1326 if (dnorm2(dr) == 0)
1328 /* No displacement, no update necessary */
1332 /* update the atom positions */
1333 auto localAtomIndices = pgrp->atomSet.localIndex();
1335 for (gmx::index j = 0; j < localAtomIndices.size(); j++)
1337 ii = localAtomIndices[j];
1338 if (!pgrp->localWeights.empty())
1340 dsvmul(pgrp->wscale*pgrp->localWeights[j], dr, tmp);
1342 for (m = 0; m < DIM; m++)
1348 for (m = 0; m < DIM; m++)
1350 v[ii][m] += invdt*tmp[m];
1356 /* calculate the constraint forces, used for output and virial only */
1357 for (size_t c = 0; c < pull->coord.size(); c++)
1359 pull_coord_work_t *pcrd;
1361 pcrd = &pull->coord[c];
1363 if (pcrd->params.eType != epullCONSTRAINT)
1368 /* Accumulate the forces, in case we have multiple constraint steps */
1369 double force = dr_tot[c]/((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)*dt*dt);
1370 pcrd->scalarForce += force;
1372 if (vir != nullptr && pcrd->params.eGeom != epullgDIRPBC && bMaster)
1376 /* Add the pull contribution to the virial */
1377 /* We have already checked above that r_ij[c] != 0 */
1378 f_invr = pcrd->scalarForce/dnorm(r_ij[c]);
1380 for (j = 0; j < DIM; j++)
1382 for (m = 0; m < DIM; m++)
1384 vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
1390 /* finished! I hope. Give back some memory */
1396 static void add_virial_coord_dr(tensor vir, const dvec dr, const dvec f)
1398 for (int j = 0; j < DIM; j++)
1400 for (int m = 0; m < DIM; m++)
1402 vir[j][m] -= 0.5*f[j]*dr[m];
1407 /* Adds the pull contribution to the virial */
1408 static void add_virial_coord(tensor vir,
1409 const pull_coord_work_t &pcrd,
1410 const PullCoordVectorForces &forces)
1412 if (vir != nullptr && pcrd.params.eGeom != epullgDIRPBC)
1414 /* Add the pull contribution for each distance vector to the virial. */
1415 add_virial_coord_dr(vir, pcrd.spatialData.dr01, forces.force01);
1416 if (pcrd.params.ngroup >= 4)
1418 add_virial_coord_dr(vir, pcrd.spatialData.dr23, forces.force23);
1420 if (pcrd.params.ngroup >= 6)
1422 add_virial_coord_dr(vir, pcrd.spatialData.dr45, forces.force45);
1427 static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t *pcrd,
1428 double dev, real lambda,
1429 real *V, real *dVdl)
1433 k = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
1434 dkdl = pcrd->params.kB - pcrd->params.k;
1436 switch (pcrd->params.eType)
1439 case epullFLATBOTTOM:
1440 case epullFLATBOTTOMHIGH:
1441 /* The only difference between an umbrella and a flat-bottom
1442 * potential is that a flat-bottom is zero above or below
1443 the reference value.
1445 if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0) ||
1446 (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
1451 pcrd->scalarForce = -k*dev;
1452 *V += 0.5* k*gmx::square(dev);
1453 *dVdl += 0.5*dkdl*gmx::square(dev);
1456 pcrd->scalarForce = -k;
1457 *V += k*pcrd->spatialData.value;
1458 *dVdl += dkdl*pcrd->spatialData.value;
1461 gmx_incons("the scalar pull force should not be calculated internally for pull type external");
1463 gmx_incons("Unsupported pull type in do_pull_pot");
1467 static PullCoordVectorForces
1468 calculateVectorForces(const pull_coord_work_t &pcrd)
1470 const t_pull_coord ¶ms = pcrd.params;
1471 const PullCoordSpatialData &spatialData = pcrd.spatialData;
1473 /* The geometry of the coordinate determines how the scalar force relates to the force on each group */
1474 PullCoordVectorForces forces;
1476 if (params.eGeom == epullgDIST)
1478 double invdr01 = spatialData.value > 0 ? 1./spatialData.value : 0.;
1479 for (int m = 0; m < DIM; m++)
1481 forces.force01[m] = pcrd.scalarForce*spatialData.dr01[m]*invdr01;
1484 else if (params.eGeom == epullgANGLE)
1487 double cos_theta, cos_theta2;
1489 cos_theta = cos(spatialData.value);
1490 cos_theta2 = gmx::square(cos_theta);
1492 /* The force at theta = 0, pi is undefined so we should not apply any force.
1493 * cos(theta) = 1 for theta = 0, pi so this is what we check for below.
1494 * Note that dr01 = 0 or dr23 = 0 evaluates to theta = 0 so we do not
1495 * have to check for this before dividing by their norm below.
1499 /* The forces f01 and f23 can be decomposed into one vector parallel to dr01
1500 * and another vector parallel to dr23:
1501 * f01 = -dV/dtheta*(a*dr23/|dr23| - b*dr01*|dr01|)/|dr01|,
1502 * f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
1504 double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1505 double b = a*cos_theta;
1506 double invdr01 = 1./dnorm(spatialData.dr01);
1507 double invdr23 = 1./dnorm(spatialData.dr23);
1508 dvec normalized_dr01, normalized_dr23;
1509 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1510 dsvmul(invdr23, spatialData.dr23, normalized_dr23);
1512 for (int m = 0; m < DIM; m++)
1514 /* Here, f_scal is -dV/dtheta */
1515 forces.force01[m] = pcrd.scalarForce*invdr01*(a*normalized_dr23[m] - b*normalized_dr01[m]);
1516 forces.force23[m] = pcrd.scalarForce*invdr23*(a*normalized_dr01[m] - b*normalized_dr23[m]);
1521 /* No forces to apply for ill-defined cases*/
1522 clear_dvec(forces.force01);
1523 clear_dvec(forces.force23);
1526 else if (params.eGeom == epullgANGLEAXIS)
1528 double cos_theta, cos_theta2;
1530 /* The angle-axis force is exactly the same as the angle force (above) except that in
1531 this case the second vector (dr23) is replaced by the pull vector. */
1532 cos_theta = cos(spatialData.value);
1533 cos_theta2 = gmx::square(cos_theta);
1539 dvec normalized_dr01;
1541 invdr01 = 1./dnorm(spatialData.dr01);
1542 dsvmul(invdr01, spatialData.dr01, normalized_dr01);
1543 a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
1546 for (int m = 0; m < DIM; m++)
1548 forces.force01[m] = pcrd.scalarForce*invdr01*(a*spatialData.vec[m] - b*normalized_dr01[m]);
1553 clear_dvec(forces.force01);
1556 else if (params.eGeom == epullgDIHEDRAL)
1558 double m2, n2, tol, sqrdist_32;
1560 /* Note: there is a small difference here compared to the
1561 dihedral force calculations in the bondeds (ref: Bekker 1994).
1562 There rij = ri - rj, while here dr01 = r1 - r0.
1563 However, all distance vectors occur in form of cross or inner products
1564 so that two signs cancel and we end up with the same expressions.
1565 Also, we treat the more general case of 6 groups (0..5) instead of 4 (i, j, k, l).
1567 m2 = diprod(spatialData.planevec_m, spatialData.planevec_m);
1568 n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
1569 dsvmul(-1, spatialData.dr23, dr32);
1570 sqrdist_32 = diprod(dr32, dr32);
1571 tol = sqrdist_32*GMX_REAL_EPS; /* Avoid tiny angles */
1572 if ((m2 > tol) && (n2 > tol))
1574 double a_01, a_23_01, a_23_45, a_45;
1575 double inv_dist_32, inv_sqrdist_32, dist_32;
1577 inv_dist_32 = gmx::invsqrt(sqrdist_32);
1578 inv_sqrdist_32 = inv_dist_32*inv_dist_32;
1579 dist_32 = sqrdist_32*inv_dist_32;
1581 /* Forces on groups 0, 1 */
1582 a_01 = pcrd.scalarForce*dist_32/m2; /* scalarForce is -dV/dphi */
1583 dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
1585 /* Forces on groups 4, 5 */
1586 a_45 = -pcrd.scalarForce*dist_32/n2;
1587 dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
1589 /* Force on groups 2, 3 (defining the axis) */
1590 a_23_01 = -diprod(spatialData.dr01, dr32)*inv_sqrdist_32;
1591 a_23_45 = -diprod(spatialData.dr45, dr32)*inv_sqrdist_32;
1592 dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
1593 dsvmul(a_23_45, forces.force45, v);
1594 dvec_sub(u, v, forces.force23); /* force on group 3 */
1598 /* No force to apply for ill-defined cases */
1599 clear_dvec(forces.force01);
1600 clear_dvec(forces.force23);
1601 clear_dvec(forces.force45);
1606 for (int m = 0; m < DIM; m++)
1608 forces.force01[m] = pcrd.scalarForce*spatialData.vec[m];
1616 /* We use a global mutex for locking access to the pull data structure
1617 * during registration of external pull potential providers.
1618 * We could use a different, local mutex for each pull object, but the overhead
1619 * is extremely small here and registration is only done during initialization.
1621 static gmx::Mutex registrationMutex;
1623 using Lock = gmx::lock_guard<gmx::Mutex>;
1625 void register_external_pull_potential(struct pull_t *pull,
1627 const char *provider)
1629 GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
1630 GMX_RELEASE_ASSERT(provider != nullptr, "register_external_pull_potential called with NULL as provider name");
1632 if (coord_index < 0 || coord_index >= static_cast<int>(pull->coord.size()))
1634 gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d which is out of the pull coordinate range %d - %zu\n",
1635 provider, coord_index + 1, 1, pull->coord.size());
1638 pull_coord_work_t *pcrd = &pull->coord[coord_index];
1640 if (pcrd->params.eType != epullEXTERNAL)
1642 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'",
1643 provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
1646 GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr, "The external potential provider string for a pull coordinate is NULL");
1648 if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider) != 0)
1650 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'",
1651 provider, coord_index + 1, pcrd->params.externalPotentialProvider);
1654 /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
1655 * pcrd->bExternalPotentialProviderHasBeenRegistered and
1656 * pull->numUnregisteredExternalPotentials.
1658 Lock registrationLock(registrationMutex);
1660 if (pcrd->bExternalPotentialProviderHasBeenRegistered)
1662 gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d more than once",
1663 provider, coord_index + 1);
1666 pcrd->bExternalPotentialProviderHasBeenRegistered = true;
1667 pull->numUnregisteredExternalPotentials--;
1669 GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0, "Negative unregistered potentials, the pull code is inconsistent");
1673 static void check_external_potential_registration(const struct pull_t *pull)
1675 if (pull->numUnregisteredExternalPotentials > 0)
1678 for (c = 0; c < pull->coord.size(); c++)
1680 if (pull->coord[c].params.eType == epullEXTERNAL &&
1681 !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
1687 GMX_RELEASE_ASSERT(c < pull->coord.size(), "Internal inconsistency in the pull potential provider counting");
1689 gmx_fatal(FARGS, "No external provider for external pull potentials have been provided for %d pull coordinates. The first coordinate without provider is number %zu, which expects a module named '%s' to provide the external potential.",
1690 pull->numUnregisteredExternalPotentials,
1691 c + 1, pull->coord[c].params.externalPotentialProvider);
1696 /* Pull takes care of adding the forces of the external potential.
1697 * The external potential module has to make sure that the corresponding
1698 * potential energy is added either to the pull term or to a term
1699 * specific to the external module.
1701 void apply_external_pull_coord_force(struct pull_t *pull,
1704 const t_mdatoms *mdatoms,
1705 gmx::ForceWithVirial *forceWithVirial)
1707 pull_coord_work_t *pcrd;
1709 GMX_ASSERT(coord_index >= 0 && coord_index < static_cast<int>(pull->coord.size()), "apply_external_pull_coord_force called with coord_index out of range");
1711 if (pull->comm.bParticipate)
1713 pcrd = &pull->coord[coord_index];
1715 GMX_RELEASE_ASSERT(pcrd->params.eType == epullEXTERNAL, "The pull force can only be set externally on pull coordinates of external type");
1717 GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered, "apply_external_pull_coord_force called for an unregistered pull coordinate");
1720 pcrd->scalarForce = coord_force;
1722 /* Calculate the forces on the pull groups */
1723 PullCoordVectorForces pullCoordForces = calculateVectorForces(*pcrd);
1725 /* Add the forces for this coordinate to the total virial and force */
1726 if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
1728 matrix virial = { { 0 } };
1729 add_virial_coord(virial, *pcrd, pullCoordForces);
1730 forceWithVirial->addVirialContribution(virial);
1733 apply_forces_coord(pull, coord_index, pullCoordForces, mdatoms,
1734 as_rvec_array(forceWithVirial->force_.data()));
1737 pull->numExternalPotentialsStillToBeAppliedThisStep--;
1740 /* Calculate the pull potential and scalar force for a pull coordinate.
1741 * Returns the vector forces for the pull coordinate.
1743 static PullCoordVectorForces
1744 do_pull_pot_coord(struct pull_t *pull, int coord_ind, t_pbc *pbc,
1745 double t, real lambda,
1746 real *V, tensor vir, real *dVdl)
1748 pull_coord_work_t &pcrd = pull->coord[coord_ind];
1750 assert(pcrd.params.eType != epullCONSTRAINT);
1752 double dev = get_pull_coord_deviation(pull, coord_ind, pbc, t);
1754 calc_pull_coord_scalar_force_and_potential(&pcrd, dev, lambda, V, dVdl);
1756 PullCoordVectorForces pullCoordForces = calculateVectorForces(pcrd);
1758 add_virial_coord(vir, pcrd, pullCoordForces);
1760 return pullCoordForces;
1763 real pull_potential(struct pull_t *pull, const t_mdatoms *md, t_pbc *pbc,
1764 const t_commrec *cr, double t, real lambda,
1765 const rvec *x, gmx::ForceWithVirial *force, real *dvdlambda)
1769 assert(pull != nullptr);
1771 /* Ideally we should check external potential registration only during
1772 * the initialization phase, but that requires another function call
1773 * that should be done exactly in the right order. So we check here.
1775 check_external_potential_registration(pull);
1777 if (pull->comm.bParticipate)
1781 pull_calc_coms(cr, pull, md, pbc, t, x, nullptr);
1783 rvec *f = as_rvec_array(force->force_.data());
1784 matrix virial = { { 0 } };
1785 const bool computeVirial = (force->computeVirial_ && MASTER(cr));
1786 for (size_t c = 0; c < pull->coord.size(); c++)
1788 /* For external potential the force is assumed to be given by an external module by a call to
1789 apply_pull_coord_external_force */
1790 if (pull->coord[c].params.eType == epullCONSTRAINT || pull->coord[c].params.eType == epullEXTERNAL)
1795 PullCoordVectorForces pullCoordForces =
1796 do_pull_pot_coord(pull, c, pbc, t, lambda,
1798 computeVirial ? virial : nullptr,
1801 /* Distribute the force over the atoms in the pulled groups */
1802 apply_forces_coord(pull, c, pullCoordForces, md, f);
1807 force->addVirialContribution(virial);
1812 GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0, "Too few or too many external pull potentials have been applied the previous step");
1813 /* All external pull potentials still need to be applied */
1814 pull->numExternalPotentialsStillToBeAppliedThisStep =
1815 pull->numCoordinatesWithExternalPotential;
1817 return (MASTER(cr) ? V : 0.0);
1820 void pull_constraint(struct pull_t *pull, const t_mdatoms *md, t_pbc *pbc,
1821 const t_commrec *cr, double dt, double t,
1822 rvec *x, rvec *xp, rvec *v, tensor vir)
1824 assert(pull != nullptr);
1826 if (pull->comm.bParticipate)
1828 pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1830 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1834 void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull)
1838 gmx_bool bMustParticipate;
1844 /* We always make the master node participate, such that it can do i/o,
1845 * add the virial and to simplify MC type extensions people might have.
1847 bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
1849 for (pull_group_work_t &group : pull->group)
1851 if (group.epgrppbc == epgrppbcCOS || !group.globalWeights.empty())
1853 group.localWeights.resize(group.atomSet.numAtomsLocal());
1854 for (size_t i = 0; i < group.atomSet.numAtomsLocal(); ++i)
1856 group.localWeights[i] = group.globalWeights[group.atomSet.collectiveIndex()[i]];
1860 GMX_ASSERT(bMustParticipate || dd != nullptr, "Either all ranks (including this rank) participate, or we use DD and need to have access to dd here");
1862 /* We should participate if we have pull or pbc atoms */
1863 if (!bMustParticipate &&
1864 (group.atomSet.numAtomsLocal() > 0 ||
1865 (group.epgrppbc == epgrppbcREFAT &&
1866 group.pbcAtomSet->numAtomsLocal() > 0)))
1868 bMustParticipate = TRUE;
1872 if (!comm->bParticipateAll)
1874 /* Keep currently not required ranks in the communicator
1875 * if they needed to participate up to 20 decompositions ago.
1876 * This avoids frequent rebuilds due to atoms jumping back and forth.
1878 const int64_t history_count = 20;
1879 gmx_bool bWillParticipate;
1882 /* Increase the decomposition counter for the current call */
1883 comm->setup_count++;
1885 if (bMustParticipate)
1887 comm->must_count = comm->setup_count;
1892 (comm->bParticipate &&
1893 comm->must_count >= comm->setup_count - history_count);
1895 if (debug && dd != nullptr)
1897 fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
1898 dd->rank, gmx::boolToString(bMustParticipate), gmx::boolToString(bWillParticipate));
1901 if (bWillParticipate)
1903 /* Count the number of ranks that we want to have participating */
1905 /* Count the number of ranks that need to be added */
1906 count[1] = comm->bParticipate ? 0 : 1;
1914 /* The cost of this global operation will be less that the cost
1915 * of the extra MPI_Comm_split calls that we can avoid.
1917 gmx_sumi(2, count, cr);
1919 /* If we are missing ranks or if we have 20% more ranks than needed
1920 * we make a new sub-communicator.
1922 if (count[1] > 0 || 6*count[0] < 5*comm->nparticipate)
1926 fprintf(debug, "Creating new pull subcommunicator of size %d\n",
1931 if (comm->mpi_comm_com != MPI_COMM_NULL)
1933 MPI_Comm_free(&comm->mpi_comm_com);
1936 /* This might be an extremely expensive operation, so we try
1937 * to avoid this splitting as much as possible.
1939 assert(dd != nullptr);
1940 MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank,
1941 &comm->mpi_comm_com);
1944 comm->bParticipate = bWillParticipate;
1945 comm->nparticipate = count[0];
1949 /* Since the PBC of atoms might have changed, we need to update the PBC */
1950 pull->bSetPBCatoms = TRUE;
1953 static void init_pull_group_index(FILE *fplog, const t_commrec *cr,
1954 int g, pull_group_work_t *pg,
1955 gmx_bool bConstraint, const ivec pulldim_con,
1956 const gmx_mtop_t *mtop,
1957 const t_inputrec *ir, real lambda)
1959 /* With EM and BD there are no masses in the integrator.
1960 * But we still want to have the correct mass-weighted COMs.
1961 * So we store the real masses in the weights.
1963 const bool setWeights = (pg->params.nweight > 0 ||
1964 EI_ENERGY_MINIMIZATION(ir->eI) ||
1967 /* In parallel, store we need to extract localWeights from weights at DD time */
1968 std::vector<real> &weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
1970 const gmx_groups_t *groups = &mtop->groups;
1972 /* Count frozen dimensions and (weighted) mass */
1978 for (int i = 0; i < pg->params.nat; i++)
1980 int ii = pg->params.ind[i];
1981 if (bConstraint && ir->opts.nFreeze)
1983 for (int d = 0; d < DIM; d++)
1985 if (pulldim_con[d] == 1 &&
1986 ir->opts.nFreeze[getGroupType(groups, egcFREEZE, ii)][d])
1992 const t_atom &atom = mtopGetAtomParameters(mtop, ii, &molb);
1994 if (ir->efep == efepNO)
2000 m = (1 - lambda)*atom.m + lambda*atom.mB;
2003 if (pg->params.nweight > 0)
2005 w = pg->params.weight[i];
2011 if (EI_ENERGY_MINIMIZATION(ir->eI))
2013 /* Move the mass to the weight */
2017 else if (ir->eI == eiBD)
2020 if (ir->bd_fric != 0.0f)
2022 mbd = ir->bd_fric*ir->delta_t;
2026 if (groups->grpnr[egcTC] == nullptr)
2028 mbd = ir->delta_t/ir->opts.tau_t[0];
2032 mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
2040 weights.push_back(w);
2049 /* We can have single atom groups with zero mass with potential pulling
2050 * without cosine weighting.
2052 if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
2054 /* With one atom the mass doesn't matter */
2059 gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
2060 pg->params.weight ? " weighted" : "", g);
2066 "Pull group %d: %5d atoms, mass %9.3f",
2067 g, pg->params.nat, tmass);
2068 if (pg->params.weight ||
2069 EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
2071 fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
2073 if (pg->epgrppbc == epgrppbcCOS)
2075 fprintf(fplog, ", cosine weighting will be used");
2077 fprintf(fplog, "\n");
2082 /* A value != 0 signals not frozen, it is updated later */
2088 for (int d = 0; d < DIM; d++)
2090 ndim += pulldim_con[d]*pg->params.nat;
2092 if (fplog && nfrozen > 0 && nfrozen < ndim)
2095 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
2096 " that are subject to pulling are frozen.\n"
2097 " For constraint pulling the whole group will be frozen.\n\n",
2106 init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
2107 const gmx_mtop_t *mtop, const t_commrec *cr, gmx::LocalAtomSetManager *atomSets,
2110 struct pull_t *pull;
2115 /* Copy the pull parameters */
2116 pull->params = *pull_params;
2117 /* Avoid pointer copies */
2118 pull->params.group = nullptr;
2119 pull->params.coord = nullptr;
2121 for (int i = 0; i < pull_params->ngroup; ++i)
2123 pull->group.emplace_back(pull_params->group[i], atomSets->add({pull_params->group[i].ind, pull_params->group[i].ind+pull_params->group[i].nat}));
2126 if (cr != nullptr && DOMAINDECOMP(cr))
2128 /* Set up the global to local atom mapping for PBC atoms */
2129 for (pull_group_work_t &group : pull->group)
2131 if (group.epgrppbc == epgrppbcREFAT)
2133 /* pbcAtomSet consists of a single atom */
2134 group.pbcAtomSet = gmx::compat::make_unique<gmx::LocalAtomSet>(atomSets->add({&group.params.pbcatom, &group.params.pbcatom + 1}));
2139 pull->bPotential = FALSE;
2140 pull->bConstraint = FALSE;
2141 pull->bCylinder = FALSE;
2142 pull->bAngle = FALSE;
2144 GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms");
2146 pull->numCoordinatesWithExternalPotential = 0;
2148 for (int c = 0; c < pull->params.ncoord; c++)
2150 /* Construct a pull coordinate, copying all coordinate parameters */
2151 pull->coord.emplace_back(pull_params->coord[c]);
2153 pull_coord_work_t *pcrd = &pull->coord.back();
2155 switch (pcrd->params.eGeom)
2158 case epullgDIRRELATIVE: /* Direction vector is determined at each step */
2160 case epullgDIHEDRAL:
2165 case epullgANGLEAXIS:
2166 copy_rvec_to_dvec(pull_params->coord[c].vec, pcrd->spatialData.vec);
2169 /* We allow reading of newer tpx files with new pull geometries,
2170 * but with the same tpx format, with old code. A new geometry
2171 * only adds a new enum value, which will be out of range for
2172 * old code. The first place we need to generate an error is
2173 * here, since the pull code can't handle this.
2174 * The enum can be used before arriving here only for printing
2175 * the string corresponding to the geometry, which will result
2176 * in printing "UNDEFINED".
2178 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.",
2179 c + 1, pcrd->params.eGeom, epullgNR - 1);
2182 if (pcrd->params.eType == epullCONSTRAINT)
2184 /* Check restrictions of the constraint pull code */
2185 if (pcrd->params.eGeom == epullgCYL ||
2186 pcrd->params.eGeom == epullgDIRRELATIVE ||
2187 pcrd->params.eGeom == epullgANGLE ||
2188 pcrd->params.eGeom == epullgDIHEDRAL ||
2189 pcrd->params.eGeom == epullgANGLEAXIS)
2191 gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
2192 epull_names[pcrd->params.eType],
2193 epullg_names[pcrd->params.eGeom],
2194 epull_names[epullUMBRELLA]);
2197 pull->bConstraint = TRUE;
2201 pull->bPotential = TRUE;
2204 if (pcrd->params.eGeom == epullgCYL)
2206 pull->bCylinder = TRUE;
2208 else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL || pcrd->params.eGeom == epullgANGLEAXIS)
2210 pull->bAngle = TRUE;
2213 /* We only need to calculate the plain COM of a group
2214 * when it is not only used as a cylinder group.
2215 * Also the absolute reference group 0 needs no COM computation.
2217 for (int i = 0; i < pcrd->params.ngroup; i++)
2219 int groupIndex = pcrd->params.group[i];
2220 if (groupIndex > 0 &&
2221 !(pcrd->params.eGeom == epullgCYL && i == 0))
2223 pull->group[groupIndex].needToCalcCom = true;
2227 /* With non-zero rate the reference value is set at every step */
2228 if (pcrd->params.rate == 0)
2230 /* Initialize the constant reference value */
2231 if (pcrd->params.eType != epullEXTERNAL)
2233 low_set_pull_coord_reference_value(pcrd, c, pcrd->params.init*pull_conversion_factor_userinput2internal(&pcrd->params));
2237 /* With an external pull potential, the reference value loses
2238 * it's meaning and should not be used. Setting it to zero
2239 * makes any terms dependent on it disappear.
2240 * The only issue this causes is that with cylinder pulling
2241 * no atoms of the cylinder group within the cylinder radius
2242 * should be more than half a box length away from the COM of
2243 * the pull group along the axial direction.
2245 pcrd->value_ref = 0.0;
2249 if (pcrd->params.eType == epullEXTERNAL)
2251 GMX_RELEASE_ASSERT(pcrd->params.rate == 0, "With an external potential, a pull coordinate should have rate = 0");
2253 /* This external potential needs to be registered later */
2254 pull->numCoordinatesWithExternalPotential++;
2256 pcrd->bExternalPotentialProviderHasBeenRegistered = false;
2259 pull->numUnregisteredExternalPotentials =
2260 pull->numCoordinatesWithExternalPotential;
2261 pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
2263 pull->ePBC = ir->ePBC;
2266 case epbcNONE: pull->npbcdim = 0; break;
2267 case epbcXY: pull->npbcdim = 2; break;
2268 default: pull->npbcdim = 3; break;
2273 gmx_bool bAbs, bCos;
2276 for (const pull_coord_work_t &coord : pull->coord)
2278 if (pull->group[coord.params.group[0]].params.nat == 0 ||
2279 pull->group[coord.params.group[1]].params.nat == 0)
2285 fprintf(fplog, "\n");
2286 if (pull->bPotential)
2288 fprintf(fplog, "Will apply potential COM pulling\n");
2290 if (pull->bConstraint)
2292 fprintf(fplog, "Will apply constraint COM pulling\n");
2294 // Don't include the reference group 0 in output, so we report ngroup-1
2295 int numRealGroups = pull->group.size() - 1;
2296 GMX_RELEASE_ASSERT(numRealGroups > 0, "The reference absolute position pull group should always be present");
2297 fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n",
2298 pull->coord.size(), pull->coord.size() == 1 ? "" : "s",
2299 numRealGroups, numRealGroups == 1 ? "" : "s");
2302 fprintf(fplog, "with an absolute reference\n");
2305 // Don't include the reference group 0 in loop
2306 for (size_t g = 1; g < pull->group.size(); g++)
2308 if (pull->group[g].params.nat > 1 &&
2309 pull->group[g].params.pbcatom < 0)
2311 /* We are using cosine weighting */
2312 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
2318 please_cite(fplog, "Engin2010");
2322 pull->bRefAt = FALSE;
2324 for (size_t g = 0; g < pull->group.size(); g++)
2326 pull_group_work_t *pgrp;
2328 pgrp = &pull->group[g];
2329 if (pgrp->params.nat > 0)
2331 /* There is an issue when a group is used in multiple coordinates
2332 * and constraints are applied in different dimensions with atoms
2333 * frozen in some, but not all dimensions.
2334 * Since there is only one mass array per group, we can't have
2335 * frozen/non-frozen atoms for different coords at the same time.
2336 * But since this is a very exotic case, we don't check for this.
2337 * A warning is printed in init_pull_group_index.
2340 gmx_bool bConstraint;
2341 ivec pulldim, pulldim_con;
2343 /* Loop over all pull coordinates to see along which dimensions
2344 * this group is pulled and if it is involved in constraints.
2346 bConstraint = FALSE;
2347 clear_ivec(pulldim);
2348 clear_ivec(pulldim_con);
2349 for (const pull_coord_work_t &coord : pull->coord)
2351 gmx_bool bGroupUsed = FALSE;
2352 for (int gi = 0; gi < coord.params.ngroup; gi++)
2354 if (coord.params.group[gi] == static_cast<int>(g))
2362 for (int m = 0; m < DIM; m++)
2364 if (coord.params.dim[m] == 1)
2368 if (coord.params.eType == epullCONSTRAINT)
2378 /* Determine if we need to take PBC into account for calculating
2379 * the COM's of the pull groups.
2381 switch (pgrp->epgrppbc)
2384 pull->bRefAt = TRUE;
2387 if (pgrp->params.weight != nullptr)
2389 gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
2391 for (int m = 0; m < DIM; m++)
2393 if (m < pull->npbcdim && pulldim[m] == 1)
2395 if (pull->cosdim >= 0 && pull->cosdim != m)
2397 gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
2407 /* Set the indices */
2408 init_pull_group_index(fplog, cr, g, pgrp,
2409 bConstraint, pulldim_con,
2414 /* Absolute reference, set the inverse mass to zero.
2415 * This is only relevant (and used) with constraint pulling.
2422 /* If we use cylinder coordinates, do some initialising for them */
2423 if (pull->bCylinder)
2425 for (const pull_coord_work_t &coord : pull->coord)
2427 if (coord.params.eGeom == epullgCYL)
2429 if (pull->group[coord.params.group[0]].params.nat == 0)
2431 gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
2434 const auto &referenceGroup = pull->group[coord.params.group[0]];
2435 pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet);
2439 /* The gmx_omp_nthreads module might not be initialized here, so max(1,) */
2440 pull->nthreads = std::max(1, gmx_omp_nthreads_get(emntDefault));
2441 pull->comSums.resize(pull->nthreads);
2446 /* Use a sub-communicator when we have more than 32 ranks, but not
2447 * when we have an external pull potential, since then the external
2448 * potential provider expects each rank to have the coordinate.
2450 comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) ||
2451 cr->dd->nnodes <= 32 ||
2452 pull->numCoordinatesWithExternalPotential > 0 ||
2453 getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
2454 /* This sub-commicator is not used with comm->bParticipateAll,
2455 * so we can always initialize it to NULL.
2457 comm->mpi_comm_com = MPI_COMM_NULL;
2458 comm->nparticipate = 0;
2459 comm->isMasterRank = (cr == nullptr || MASTER(cr));
2461 /* No MPI: 1 rank: all ranks pull */
2462 comm->bParticipateAll = TRUE;
2463 comm->isMasterRank = true;
2465 comm->bParticipate = comm->bParticipateAll;
2466 comm->setup_count = 0;
2467 comm->must_count = 0;
2469 if (!comm->bParticipateAll && fplog != nullptr)
2471 fprintf(fplog, "Will use a sub-communicator for pull communication\n");
2474 comm->pbcAtomBuffer.resize(pull->group.size());
2475 comm->comBuffer.resize(pull->group.size()*DIM);
2476 if (pull->bCylinder)
2478 comm->cylinderBuffer.resize(pull->coord.size()*c_cylinderBufferStride);
2481 /* We still need to initialize the PBC reference coordinates */
2482 pull->bSetPBCatoms = TRUE;
2484 pull->out_x = nullptr;
2485 pull->out_f = nullptr;
2490 void init_pull_output_files(pull_t *pull,
2492 const t_filenm fnm[],
2493 const gmx_output_env_t *oenv,
2494 const ContinuationOptions &continuationOptions)
2496 /* Check for px and pf filename collision, if we are writing
2498 std::string px_filename, pf_filename;
2499 std::string px_appended, pf_appended;
2502 px_filename = std::string(opt2fn("-px", nfile, fnm));
2503 pf_filename = std::string(opt2fn("-pf", nfile, fnm));
2505 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2507 if ((pull->params.nstxout != 0) &&
2508 (pull->params.nstfout != 0) &&
2509 (px_filename == pf_filename))
2511 if (!opt2bSet("-px", nfile, fnm) && !opt2bSet("-pf", nfile, fnm))
2513 /* We are writing both pull files but neither set directly. */
2516 px_appended = append_before_extension(px_filename, "_pullx");
2517 pf_appended = append_before_extension(pf_filename, "_pullf");
2519 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
2520 pull->out_x = open_pull_out(px_appended.c_str(), pull, oenv,
2521 TRUE, continuationOptions);
2522 pull->out_f = open_pull_out(pf_appended.c_str(), pull, oenv,
2523 FALSE, continuationOptions);
2528 /* If at least one of -px and -pf is set but the filenames are identical: */
2529 gmx_fatal(FARGS, "Identical pull_x and pull_f output filenames %s",
2530 px_filename.c_str());
2533 if (pull->params.nstxout != 0)
2535 pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
2536 TRUE, continuationOptions);
2538 if (pull->params.nstfout != 0)
2540 pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
2541 FALSE, continuationOptions);
2545 static void destroy_pull(struct pull_t *pull)
2548 if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
2550 MPI_Comm_free(&pull->comm.mpi_comm_com);
2557 void finish_pull(struct pull_t *pull)
2559 check_external_potential_registration(pull);
2563 gmx_fio_fclose(pull->out_x);
2567 gmx_fio_fclose(pull->out_f);
2573 gmx_bool pull_have_potential(const struct pull_t *pull)
2575 return pull->bPotential;
2578 gmx_bool pull_have_constraint(const struct pull_t *pull)
2580 return pull->bConstraint;