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, 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.
49 #include "gromacs/fileio/filenm.h"
50 #include "gromacs/fileio/gmxfio.h"
51 #include "gromacs/fileio/xvgr.h"
52 #include "gromacs/legacyheaders/copyrite.h"
53 #include "gromacs/legacyheaders/gmx_ga2la.h"
54 #include "gromacs/legacyheaders/macros.h"
55 #include "gromacs/legacyheaders/mdrun.h"
56 #include "gromacs/legacyheaders/names.h"
57 #include "gromacs/legacyheaders/network.h"
58 #include "gromacs/legacyheaders/typedefs.h"
59 #include "gromacs/legacyheaders/types/commrec.h"
60 #include "gromacs/math/vec.h"
61 #include "gromacs/pbcutil/pbc.h"
62 #include "gromacs/pulling/pull_internal.h"
63 #include "gromacs/topology/mtop_util.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/gmxassert.h"
66 #include "gromacs/utility/smalloc.h"
68 static void pull_print_group_x(FILE *out, ivec dim,
69 const pull_group_work_t *pgrp)
73 for (m = 0; m < DIM; m++)
77 fprintf(out, "\t%g", pgrp->x[m]);
82 static void pull_print_coord_dr(FILE *out, const pull_coord_work_t *pcrd,
83 gmx_bool bPrintRefValue,
84 gmx_bool bPrintComponents)
86 fprintf(out, "\t%g", pcrd->value);
90 fprintf(out, "\t%g", pcrd->value_ref);
97 for (m = 0; m < DIM; m++)
99 if (pcrd->params.dim[m])
101 fprintf(out, "\t%g", pcrd->dr[m]);
107 static void pull_print_x(FILE *out, struct pull_t *pull, double t)
111 fprintf(out, "%.4f", t);
113 for (c = 0; c < pull->ncoord; c++)
115 pull_coord_work_t *pcrd;
117 pcrd = &pull->coord[c];
119 if (pull->params.bPrintCOM1)
121 if (pcrd->params.eGeom == epullgCYL)
123 pull_print_group_x(out, pcrd->params.dim, &pull->dyna[c]);
127 pull_print_group_x(out, pcrd->params.dim,
128 &pull->group[pcrd->params.group[0]]);
131 if (pull->params.bPrintCOM2)
133 pull_print_group_x(out, pcrd->params.dim,
134 &pull->group[pcrd->params.group[1]]);
136 pull_print_coord_dr(out, pcrd,
137 pull->params.bPrintRefValue,
138 pull->params.bPrintComp);
143 static void pull_print_f(FILE *out, struct pull_t *pull, double t)
147 fprintf(out, "%.4f", t);
149 for (c = 0; c < pull->ncoord; c++)
151 fprintf(out, "\t%g", pull->coord[c].f_scal);
156 void pull_print_output(struct pull_t *pull, gmx_int64_t step, double time)
158 if ((pull->params.nstxout != 0) && (step % pull->params.nstxout == 0))
160 pull_print_x(pull->out_x, pull, time);
163 if ((pull->params.nstfout != 0) && (step % pull->params.nstfout == 0))
165 pull_print_f(pull->out_f, pull, time);
169 static FILE *open_pull_out(const char *fn, struct pull_t *pull, const output_env_t oenv,
170 gmx_bool bCoord, unsigned long Flags)
174 char **setname, buf[10];
176 if (Flags & MD_APPENDFILES)
178 fp = gmx_fio_fopen(fn, "a+");
182 fp = gmx_fio_fopen(fn, "w+");
185 xvgr_header(fp, "Pull COM", "Time (ps)", "Position (nm)",
190 xvgr_header(fp, "Pull force", "Time (ps)", "Force (kJ/mol/nm)",
194 /* With default mdp options only the actual distance is printed,
195 * but optionally 2 COMs, the reference distance and distance
196 * components can also be printed.
198 snew(setname, pull->ncoord*(DIM + DIM + 1 + 1 + DIM));
200 for (c = 0; c < pull->ncoord; c++)
204 /* The order of this legend should match the order of printing
205 * the data in print_pull_x above.
208 if (pull->params.bPrintCOM1)
210 /* Legend for reference group position */
211 for (m = 0; m < DIM; m++)
213 if (pull->coord[c].params.dim[m])
215 sprintf(buf, "%d %s%d%c", c+1, "c", 1, 'X'+m);
216 setname[nsets] = gmx_strdup(buf);
221 if (pull->params.bPrintCOM2)
223 /* Legend for reference group position */
224 for (m = 0; m < DIM; m++)
226 if (pull->coord[c].params.dim[m])
228 sprintf(buf, "%d %s%d%c", c+1, "c", 2, 'X'+m);
229 setname[nsets] = gmx_strdup(buf);
234 /* The pull coord distance */
235 sprintf(buf, "%d", c+1);
236 setname[nsets] = gmx_strdup(buf);
238 if (pull->params.bPrintRefValue)
240 sprintf(buf, "%c%d", 'r', c+1);
241 setname[nsets] = gmx_strdup(buf);
244 if (pull->params.bPrintComp)
246 /* The pull coord distance components */
247 for (m = 0; m < DIM; m++)
249 if (pull->coord[c].params.dim[m])
251 sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
252 setname[nsets] = gmx_strdup(buf);
260 /* For the pull force we always only use one scalar */
261 sprintf(buf, "%d", c+1);
262 setname[nsets] = gmx_strdup(buf);
268 xvgr_legend(fp, nsets, (const char**)setname, oenv);
270 for (c = 0; c < nsets; c++)
280 /* Apply forces in a mass weighted fashion */
281 static void apply_forces_grp(const pull_group_work_t *pgrp,
283 const dvec f_pull, int sign, rvec *f)
286 double wmass, inv_wm;
288 inv_wm = pgrp->mwscale;
290 if (pgrp->params.nat == 1 && pgrp->nat_loc == 1)
292 /* Only one atom and our rank has this atom: we can skip
293 * the mass weighting, which means that this code also works
294 * for mass=0, e.g. with a virtual site.
296 for (m = 0; m < DIM; m++)
298 f[pgrp->ind_loc[0]][m] += sign*f_pull[m];
303 for (i = 0; i < pgrp->nat_loc; i++)
305 ii = pgrp->ind_loc[i];
306 wmass = md->massT[ii];
307 if (pgrp->weight_loc)
309 wmass *= pgrp->weight_loc[i];
312 for (m = 0; m < DIM; m++)
314 f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
320 /* Apply forces in a mass weighted fashion to a cylinder group */
321 static void apply_forces_cyl_grp(const pull_group_work_t *pgrp,
324 const dvec f_pull, double f_scal,
328 double mass, weight, inv_wm, dv_com;
330 inv_wm = pgrp->mwscale;
332 for (i = 0; i < pgrp->nat_loc; i++)
334 ii = pgrp->ind_loc[i];
335 mass = md->massT[ii];
336 weight = pgrp->weight_loc[i];
337 /* The stored axial distance from the cylinder center (dv) needs
338 * to be corrected for an offset (dv_corr), which was unknown when
341 dv_com = pgrp->dv[i] + dv_corr;
343 /* Here we not only add the pull force working along vec (f_pull),
344 * but also a radial component, due to the dependence of the weights
345 * on the radial distance.
347 for (m = 0; m < DIM; m++)
349 f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
350 pgrp->mdw[i][m]*dv_com*f_scal);
355 /* Apply torque forces in a mass weighted fashion to the groups that define
356 * the pull vector direction for pull coordinate pcrd.
358 static void apply_forces_vec_torque(const struct pull_t *pull,
359 const pull_coord_work_t *pcrd,
367 /* The component inpr along the pull vector is accounted for in the usual
368 * way. Here we account for the component perpendicular to vec.
371 for (m = 0; m < DIM; m++)
373 inpr += pcrd->dr[m]*pcrd->vec[m];
375 /* The torque force works along the component of the distance vector
376 * of between the two "usual" pull groups that is perpendicular to
377 * the pull vector. The magnitude of this force is the "usual" scale force
378 * multiplied with the ratio of the distance between the two "usual" pull
379 * groups and the distance between the two groups that define the vector.
381 for (m = 0; m < DIM; m++)
383 f_perp[m] = (pcrd->dr[m] - inpr*pcrd->vec[m])/pcrd->vec_len*pcrd->f_scal;
386 /* Apply the force to the groups defining the vector using opposite signs */
387 apply_forces_grp(&pull->group[pcrd->params.group[2]], md, f_perp, -1, f);
388 apply_forces_grp(&pull->group[pcrd->params.group[3]], md, f_perp, 1, f);
391 /* Apply forces in a mass weighted fashion */
392 static void apply_forces(struct pull_t * pull, t_mdatoms * md, rvec *f)
396 for (c = 0; c < pull->ncoord; c++)
398 const pull_coord_work_t *pcrd;
400 pcrd = &pull->coord[c];
402 if (pcrd->params.eGeom == epullgCYL)
407 apply_forces_cyl_grp(&pull->dyna[c], pcrd->cyl_dev, md,
408 pcrd->f, pcrd->f_scal, -1, f);
410 /* Sum the force along the vector and the radial force */
411 for (m = 0; m < DIM; m++)
413 f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
415 apply_forces_grp(&pull->group[pcrd->params.group[1]], md, f_tot, 1, f);
419 if (pcrd->params.eGeom == epullgDIRRELATIVE)
421 /* We need to apply the torque forces to the pull groups
422 * that define the pull vector.
424 apply_forces_vec_torque(pull, pcrd, md, f);
427 if (pull->group[pcrd->params.group[0]].params.nat > 0)
429 apply_forces_grp(&pull->group[pcrd->params.group[0]], md, pcrd->f, -1, f);
431 apply_forces_grp(&pull->group[pcrd->params.group[1]], md, pcrd->f, 1, f);
436 static double max_pull_distance2(const pull_coord_work_t *pcrd,
442 max_d2 = GMX_DOUBLE_MAX;
444 for (m = 0; m < pbc->ndim_ePBC; m++)
446 if (pcrd->params.dim[m] != 0)
448 max_d2 = std::min(max_d2, static_cast<double>(norm2(pbc->box[m])));
455 /* This function returns the distance based on coordinates xg and xref.
456 * Note that the pull coordinate struct pcrd is not modified.
458 static void low_get_pull_coord_dr(const struct pull_t *pull,
459 const pull_coord_work_t *pcrd,
461 dvec xg, dvec xref, double max_dist2,
464 const pull_group_work_t *pgrp0;
466 dvec xrefr, dref = {0, 0, 0};
469 pgrp0 = &pull->group[pcrd->params.group[0]];
471 /* Only the first group can be an absolute reference, in that case nat=0 */
472 if (pgrp0->params.nat == 0)
474 for (m = 0; m < DIM; m++)
476 xref[m] = pcrd->params.origin[m];
480 copy_dvec(xref, xrefr);
482 if (pcrd->params.eGeom == epullgDIRPBC)
484 for (m = 0; m < DIM; m++)
486 dref[m] = pcrd->value_ref*pcrd->vec[m];
488 /* Add the reference position, so we use the correct periodic image */
489 dvec_inc(xrefr, dref);
492 pbc_dx_d(pbc, xg, xrefr, dr);
494 for (m = 0; m < DIM; m++)
496 dr[m] *= pcrd->params.dim[m];
499 if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
501 gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f).\nYou might want to consider using \"pull-geometry = direction-periodic\" instead.\n",
502 pcrd->params.group[0], pcrd->params.group[1],
503 sqrt(dr2), sqrt(max_dist2));
506 if (pcrd->params.eGeom == epullgDIRPBC)
512 /* This function returns the distance based on the contents of the pull struct.
513 * pull->coord[coord_ind].dr, and potentially vector, are updated.
515 static void get_pull_coord_dr(struct pull_t *pull,
520 pull_coord_work_t *pcrd;
522 pcrd = &pull->coord[coord_ind];
524 if (pcrd->params.eGeom == epullgDIRPBC)
530 md2 = max_pull_distance2(pcrd, pbc);
533 if (pcrd->params.eGeom == epullgDIRRELATIVE)
535 /* We need to determine the pull vector */
536 const pull_group_work_t *pgrp2, *pgrp3;
540 pgrp2 = &pull->group[pcrd->params.group[2]];
541 pgrp3 = &pull->group[pcrd->params.group[3]];
543 pbc_dx_d(pbc, pgrp3->x, pgrp2->x, vec);
545 for (m = 0; m < DIM; m++)
547 vec[m] *= pcrd->params.dim[m];
549 pcrd->vec_len = dnorm(vec);
550 for (m = 0; m < DIM; m++)
552 pcrd->vec[m] = vec[m]/pcrd->vec_len;
556 fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
558 vec[XX], vec[YY], vec[ZZ],
559 pcrd->vec[XX], pcrd->vec[YY], pcrd->vec[ZZ]);
563 low_get_pull_coord_dr(pull, pcrd, pbc,
564 pull->group[pcrd->params.group[1]].x,
565 pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->params.group[0]].x,
570 static void update_pull_coord_reference_value(pull_coord_work_t *pcrd, double t)
572 /* With zero rate the reference value is set initially and doesn't change */
573 if (pcrd->params.rate != 0)
575 pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t;
579 /* Calculates pull->coord[coord_ind].value.
580 * This function also updates pull->coord[coord_ind].dr.
582 static void get_pull_coord_distance(struct pull_t *pull,
586 pull_coord_work_t *pcrd;
589 pcrd = &pull->coord[coord_ind];
591 get_pull_coord_dr(pull, coord_ind, pbc);
593 switch (pcrd->params.eGeom)
596 /* Pull along the vector between the com's */
597 pcrd->value = dnorm(pcrd->dr);
601 case epullgDIRRELATIVE:
605 for (m = 0; m < DIM; m++)
607 pcrd->value += pcrd->vec[m]*pcrd->dr[m];
613 /* Returns the deviation from the reference value.
614 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
616 static double get_pull_coord_deviation(struct pull_t *pull,
621 static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
622 but is fairly benign */
623 pull_coord_work_t *pcrd;
626 pcrd = &pull->coord[coord_ind];
628 get_pull_coord_distance(pull, coord_ind, pbc);
630 update_pull_coord_reference_value(pcrd, t);
632 switch (pcrd->params.eGeom)
635 /* Pull along the vector between the com's */
636 if (pcrd->value_ref < 0 && !bWarned)
638 fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, pcrd->value_ref);
641 if (pcrd->value == 0)
643 /* With no vector we can not determine the direction for the force,
644 * so we set the force to zero.
650 /* Determine the deviation */
651 dev = pcrd->value - pcrd->value_ref;
656 case epullgDIRRELATIVE:
659 dev = pcrd->value - pcrd->value_ref;
666 void get_pull_coord_value(struct pull_t *pull,
671 get_pull_coord_distance(pull, coord_ind, pbc);
673 *value = pull->coord[coord_ind].value;
676 void clear_pull_forces(struct pull_t *pull)
680 /* Zeroing the forces is only required for constraint pulling.
681 * It can happen that multiple constraint steps need to be applied
682 * and therefore the constraint forces need to be accumulated.
684 for (i = 0; i < pull->ncoord; i++)
686 clear_dvec(pull->coord[i].f);
687 pull->coord[i].f_scal = 0;
691 /* Apply constraint using SHAKE */
692 static void do_constraint(struct pull_t *pull, t_pbc *pbc,
694 gmx_bool bMaster, tensor vir,
698 dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
699 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
700 dvec *rnew; /* current 'new' positions of the groups */
701 double *dr_tot; /* the total update of the coords */
704 double lambda, rm, invdt = 0;
705 gmx_bool bConverged_all, bConverged = FALSE;
706 int niter = 0, g, c, ii, j, m, max_iter = 100;
710 snew(r_ij, pull->ncoord);
711 snew(dr_tot, pull->ncoord);
713 snew(rnew, pull->ngroup);
715 /* copy the current unconstrained positions for use in iterations. We
716 iterate until rinew[i] and rjnew[j] obey the constraints. Then
717 rinew - pull.x_unc[i] is the correction dr to group i */
718 for (g = 0; g < pull->ngroup; g++)
720 copy_dvec(pull->group[g].xp, rnew[g]);
723 /* Determine the constraint directions from the old positions */
724 for (c = 0; c < pull->ncoord; c++)
726 pull_coord_work_t *pcrd;
728 pcrd = &pull->coord[c];
730 if (pcrd->params.eType != epullCONSTRAINT)
735 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
736 * We don't modify dr and value anymore, so these values are also used
739 get_pull_coord_distance(pull, c, pbc);
743 fprintf(debug, "Pull coord %d dr %f %f %f\n",
744 c, pcrd->dr[XX], pcrd->dr[YY], pcrd->dr[ZZ]);
747 if (pcrd->params.eGeom == epullgDIR ||
748 pcrd->params.eGeom == epullgDIRPBC)
750 /* Select the component along vec */
752 for (m = 0; m < DIM; m++)
754 a += pcrd->vec[m]*pcrd->dr[m];
756 for (m = 0; m < DIM; m++)
758 r_ij[c][m] = a*pcrd->vec[m];
763 copy_dvec(pcrd->dr, r_ij[c]);
766 if (dnorm2(r_ij[c]) == 0)
768 gmx_fatal(FARGS, "Distance for pull coordinate %d is zero with constraint pulling, which is not allowed.", c + 1);
772 bConverged_all = FALSE;
773 while (!bConverged_all && niter < max_iter)
775 bConverged_all = TRUE;
777 /* loop over all constraints */
778 for (c = 0; c < pull->ncoord; c++)
780 pull_coord_work_t *pcrd;
781 pull_group_work_t *pgrp0, *pgrp1;
784 pcrd = &pull->coord[c];
786 if (pcrd->params.eType != epullCONSTRAINT)
791 update_pull_coord_reference_value(pcrd, t);
793 pgrp0 = &pull->group[pcrd->params.group[0]];
794 pgrp1 = &pull->group[pcrd->params.group[1]];
796 /* Get the current difference vector */
797 low_get_pull_coord_dr(pull, pcrd, pbc,
798 rnew[pcrd->params.group[1]],
799 rnew[pcrd->params.group[0]],
804 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
807 rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
809 switch (pcrd->params.eGeom)
812 if (pcrd->value_ref <= 0)
814 gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, pcrd->value_ref);
818 double q, c_a, c_b, c_c;
820 c_a = diprod(r_ij[c], r_ij[c]);
821 c_b = diprod(unc_ij, r_ij[c])*2;
822 c_c = diprod(unc_ij, unc_ij) - dsqr(pcrd->value_ref);
826 q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
831 q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
838 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
839 c_a, c_b, c_c, lambda);
843 /* The position corrections dr due to the constraints */
844 dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
845 dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
846 dr_tot[c] += -lambda*dnorm(r_ij[c]);
851 /* A 1-dimensional constraint along a vector */
853 for (m = 0; m < DIM; m++)
855 vec[m] = pcrd->vec[m];
856 a += unc_ij[m]*vec[m];
858 /* Select only the component along the vector */
859 dsvmul(a, vec, unc_ij);
860 lambda = a - pcrd->value_ref;
863 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
866 /* The position corrections dr due to the constraints */
867 dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
868 dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
869 dr_tot[c] += -lambda;
872 gmx_incons("Invalid enumeration value for eGeom");
873 /* Keep static analyzer happy */
883 g0 = pcrd->params.group[0];
884 g1 = pcrd->params.group[1];
885 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
886 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
888 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
889 rnew[g0][0], rnew[g0][1], rnew[g0][2],
890 rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
892 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
893 "", "", "", "", "", "", pcrd->value_ref);
895 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
896 dr0[0], dr0[1], dr0[2],
897 dr1[0], dr1[1], dr1[2],
901 /* Update the COMs with dr */
902 dvec_inc(rnew[pcrd->params.group[1]], dr1);
903 dvec_inc(rnew[pcrd->params.group[0]], dr0);
906 /* Check if all constraints are fullfilled now */
907 for (c = 0; c < pull->ncoord; c++)
909 pull_coord_work_t *pcrd;
911 pcrd = &pull->coord[c];
913 if (pcrd->params.eType != epullCONSTRAINT)
918 low_get_pull_coord_dr(pull, pcrd, pbc,
919 rnew[pcrd->params.group[1]],
920 rnew[pcrd->params.group[0]],
923 switch (pcrd->params.eGeom)
927 fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->params.constr_tol;
932 for (m = 0; m < DIM; m++)
934 vec[m] = pcrd->vec[m];
936 inpr = diprod(unc_ij, vec);
937 dsvmul(inpr, vec, unc_ij);
939 fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->params.constr_tol;
947 fprintf(debug, "NOT CONVERGED YET: Group %d:"
948 "d_ref = %f, current d = %f\n",
949 g, pcrd->value_ref, dnorm(unc_ij));
952 bConverged_all = FALSE;
957 /* if after all constraints are dealt with and bConverged is still TRUE
958 we're finished, if not we do another iteration */
960 if (niter > max_iter)
962 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
965 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
972 /* update atoms in the groups */
973 for (g = 0; g < pull->ngroup; g++)
975 const pull_group_work_t *pgrp;
978 pgrp = &pull->group[g];
980 /* get the final constraint displacement dr for group g */
981 dvec_sub(rnew[g], pgrp->xp, dr);
985 /* No displacement, no update necessary */
989 /* update the atom positions */
991 for (j = 0; j < pgrp->nat_loc; j++)
993 ii = pgrp->ind_loc[j];
994 if (pgrp->weight_loc)
996 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
998 for (m = 0; m < DIM; m++)
1004 for (m = 0; m < DIM; m++)
1006 v[ii][m] += invdt*tmp[m];
1012 /* calculate the constraint forces, used for output and virial only */
1013 for (c = 0; c < pull->ncoord; c++)
1015 pull_coord_work_t *pcrd;
1017 pcrd = &pull->coord[c];
1019 if (pcrd->params.eType != epullCONSTRAINT)
1024 pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)*dt*dt);
1026 if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC && bMaster)
1030 /* Add the pull contribution to the virial */
1031 /* We have already checked above that r_ij[c] != 0 */
1032 f_invr = pcrd->f_scal/dnorm(r_ij[c]);
1034 for (j = 0; j < DIM; j++)
1036 for (m = 0; m < DIM; m++)
1038 vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
1044 /* finished! I hope. Give back some memory */
1050 /* Pulling with a harmonic umbrella potential or constant force */
1051 static void do_pull_pot(struct pull_t *pull, t_pbc *pbc, double t, real lambda,
1052 real *V, tensor vir, real *dVdl)
1055 double dev, ndr, invdr = 0;
1058 /* loop over the pull coordinates */
1061 for (c = 0; c < pull->ncoord; c++)
1063 pull_coord_work_t *pcrd;
1065 pcrd = &pull->coord[c];
1067 if (pcrd->params.eType == epullCONSTRAINT)
1072 dev = get_pull_coord_deviation(pull, c, pbc, t);
1074 k = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
1075 dkdl = pcrd->params.kB - pcrd->params.k;
1077 if (pcrd->params.eGeom == epullgDIST)
1079 ndr = dnorm(pcrd->dr);
1086 /* With an harmonic umbrella, the force is 0 at r=0,
1087 * so we can set invdr to any value.
1088 * With a constant force, the force at r=0 is not defined,
1089 * so we zero it (this is anyhow a very rare event).
1097 for (m = 0; m < DIM; m++)
1099 ndr += pcrd->vec[m]*pcrd->dr[m];
1103 switch (pcrd->params.eType)
1106 case epullFLATBOTTOM:
1107 /* The only difference between an umbrella and a flat-bottom
1108 * potential is that a flat-bottom is zero below zero.
1110 if (pcrd->params.eType == epullFLATBOTTOM && dev < 0)
1115 pcrd->f_scal = -k*dev;
1116 *V += 0.5* k*dsqr(dev);
1117 *dVdl += 0.5*dkdl*dsqr(dev);
1125 gmx_incons("Unsupported pull type in do_pull_pot");
1128 if (pcrd->params.eGeom == epullgDIST)
1130 for (m = 0; m < DIM; m++)
1132 pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
1137 for (m = 0; m < DIM; m++)
1139 pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
1143 if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC)
1145 /* Add the pull contribution to the virial */
1146 for (j = 0; j < DIM; j++)
1148 for (m = 0; m < DIM; m++)
1150 vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
1157 real pull_potential(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
1158 t_commrec *cr, double t, real lambda,
1159 rvec *x, rvec *f, tensor vir, real *dvdlambda)
1163 assert(pull != NULL);
1165 pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
1167 do_pull_pot(pull, pbc, t, lambda,
1168 &V, MASTER(cr) ? vir : NULL, &dVdl);
1170 /* Distribute forces over pulled groups */
1171 apply_forces(pull, md, f);
1178 return (MASTER(cr) ? V : 0.0);
1181 void pull_constraint(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
1182 t_commrec *cr, double dt, double t,
1183 rvec *x, rvec *xp, rvec *v, tensor vir)
1185 assert(pull != NULL);
1187 pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1189 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1192 static void make_local_pull_group(gmx_ga2la_t ga2la,
1193 pull_group_work_t *pg, int start, int end)
1198 for (i = 0; i < pg->params.nat; i++)
1200 ii = pg->params.ind[i];
1203 if (!ga2la_get_home(ga2la, ii, &ii))
1208 if (ii >= start && ii < end)
1210 /* This is a home atom, add it to the local pull group */
1211 if (pg->nat_loc >= pg->nalloc_loc)
1213 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1214 srenew(pg->ind_loc, pg->nalloc_loc);
1215 if (pg->epgrppbc == epgrppbcCOS || pg->params.weight != NULL)
1217 srenew(pg->weight_loc, pg->nalloc_loc);
1220 pg->ind_loc[pg->nat_loc] = ii;
1221 if (pg->params.weight != NULL)
1223 pg->weight_loc[pg->nat_loc] = pg->params.weight[i];
1230 void dd_make_local_pull_groups(gmx_domdec_t *dd, struct pull_t *pull, t_mdatoms *md)
1244 for (g = 0; g < pull->ngroup; g++)
1246 make_local_pull_group(ga2la, &pull->group[g],
1250 /* Since the PBC of atoms might have changed, we need to update the PBC */
1251 pull->bSetPBCatoms = TRUE;
1254 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
1255 int g, pull_group_work_t *pg,
1256 gmx_bool bConstraint, ivec pulldim_con,
1257 const gmx_mtop_t *mtop,
1258 const t_inputrec *ir, real lambda)
1260 int i, ii, d, nfrozen, ndim;
1262 double tmass, wmass, wwmass;
1263 const gmx_groups_t *groups;
1264 gmx_mtop_atomlookup_t alook;
1267 if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1269 /* There are no masses in the integrator.
1270 * But we still want to have the correct mass-weighted COMs.
1271 * So we store the real masses in the weights.
1272 * We do not set nweight, so these weights do not end up in the tpx file.
1274 if (pg->params.nweight == 0)
1276 snew(pg->params.weight, pg->params.nat);
1285 pg->weight_loc = NULL;
1289 pg->nat_loc = pg->params.nat;
1290 pg->ind_loc = pg->params.ind;
1291 if (pg->epgrppbc == epgrppbcCOS)
1293 snew(pg->weight_loc, pg->params.nat);
1297 pg->weight_loc = pg->params.weight;
1301 groups = &mtop->groups;
1303 alook = gmx_mtop_atomlookup_init(mtop);
1309 for (i = 0; i < pg->params.nat; i++)
1311 ii = pg->params.ind[i];
1312 gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1313 if (bConstraint && ir->opts.nFreeze)
1315 for (d = 0; d < DIM; d++)
1317 if (pulldim_con[d] == 1 &&
1318 ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1324 if (ir->efep == efepNO)
1330 m = (1 - lambda)*atom->m + lambda*atom->mB;
1332 if (pg->params.nweight > 0)
1334 w = pg->params.weight[i];
1340 if (EI_ENERGY_MINIMIZATION(ir->eI))
1342 /* Move the mass to the weight */
1345 pg->params.weight[i] = w;
1347 else if (ir->eI == eiBD)
1351 mbd = ir->bd_fric*ir->delta_t;
1355 if (groups->grpnr[egcTC] == NULL)
1357 mbd = ir->delta_t/ir->opts.tau_t[0];
1361 mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1366 pg->params.weight[i] = w;
1373 gmx_mtop_atomlookup_destroy(alook);
1377 /* We can have single atom groups with zero mass with potential pulling
1378 * without cosine weighting.
1380 if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1382 /* With one atom the mass doesn't matter */
1387 gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1388 pg->params.weight ? " weighted" : "", g);
1394 "Pull group %d: %5d atoms, mass %9.3f",
1395 g, pg->params.nat, tmass);
1396 if (pg->params.weight ||
1397 EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1399 fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1401 if (pg->epgrppbc == epgrppbcCOS)
1403 fprintf(fplog, ", cosine weighting will be used");
1405 fprintf(fplog, "\n");
1410 /* A value != 0 signals not frozen, it is updated later */
1416 for (d = 0; d < DIM; d++)
1418 ndim += pulldim_con[d]*pg->params.nat;
1420 if (fplog && nfrozen > 0 && nfrozen < ndim)
1423 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1424 " that are subject to pulling are frozen.\n"
1425 " For constraint pulling the whole group will be frozen.\n\n",
1434 init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
1435 int nfile, const t_filenm fnm[],
1436 gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1437 gmx_bool bOutFile, unsigned long Flags)
1439 struct pull_t *pull;
1444 /* Copy the pull parameters */
1445 pull->params = *pull_params;
1446 /* Avoid pointer copies */
1447 pull->params.group = NULL;
1448 pull->params.coord = NULL;
1450 pull->ncoord = pull_params->ncoord;
1451 pull->ngroup = pull_params->ngroup;
1452 snew(pull->coord, pull->ncoord);
1453 snew(pull->group, pull->ngroup);
1455 pull->bPotential = FALSE;
1456 pull->bConstraint = FALSE;
1457 pull->bCylinder = FALSE;
1459 for (g = 0; g < pull->ngroup; g++)
1461 pull_group_work_t *pgrp;
1464 pgrp = &pull->group[g];
1466 /* Copy the pull group parameters */
1467 pgrp->params = pull_params->group[g];
1469 /* Avoid pointer copies by allocating and copying arrays */
1470 snew(pgrp->params.ind, pgrp->params.nat);
1471 for (i = 0; i < pgrp->params.nat; i++)
1473 pgrp->params.ind[i] = pull_params->group[g].ind[i];
1475 if (pgrp->params.nweight > 0)
1477 snew(pgrp->params.ind, pgrp->params.nweight);
1478 for (i = 0; i < pgrp->params.nweight; i++)
1480 pgrp->params.weight[i] = pull_params->group[g].weight[i];
1485 for (c = 0; c < pull->ncoord; c++)
1487 pull_coord_work_t *pcrd;
1488 int calc_com_start, calc_com_end, g;
1490 pcrd = &pull->coord[c];
1492 /* Copy all pull coordinate parameters */
1493 pcrd->params = pull_params->coord[c];
1495 switch (pcrd->params.eGeom)
1498 case epullgDIRRELATIVE:
1499 /* Direction vector is determined at each step */
1504 copy_rvec(pull_params->coord[c].vec, pcrd->vec);
1507 gmx_incons("Pull geometry not handled");
1510 if (pcrd->params.eType == epullCONSTRAINT)
1512 /* Check restrictions of the constraint pull code */
1513 if (pcrd->params.eGeom == epullgCYL ||
1514 pcrd->params.eGeom == epullgDIRRELATIVE)
1516 gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
1517 epull_names[pcrd->params.eType],
1518 epullg_names[pcrd->params.eGeom],
1519 epull_names[epullUMBRELLA]);
1522 pull->bConstraint = TRUE;
1526 pull->bPotential = TRUE;
1529 if (pcrd->params.eGeom == epullgCYL)
1531 pull->bCylinder = TRUE;
1533 /* We only need to calculate the plain COM of a group
1534 * when it is not only used as a cylinder group.
1536 calc_com_start = (pcrd->params.eGeom == epullgCYL ? 1 : 0);
1537 calc_com_end = (pcrd->params.eGeom == epullgDIRRELATIVE ? 4 : 2);
1539 for (g = calc_com_start; g <= calc_com_end; g++)
1541 pull->group[pcrd->params.group[g]].bCalcCOM = TRUE;
1544 /* With non-zero rate the reference value is set at every step */
1545 if (pcrd->params.rate == 0)
1547 /* Initialize the constant reference value */
1548 pcrd->value_ref = pcrd->params.init;
1552 pull->ePBC = ir->ePBC;
1555 case epbcNONE: pull->npbcdim = 0; break;
1556 case epbcXY: pull->npbcdim = 2; break;
1557 default: pull->npbcdim = 3; break;
1562 gmx_bool bAbs, bCos;
1565 for (c = 0; c < pull->ncoord; c++)
1567 if (pull->group[pull->coord[c].params.group[0]].params.nat == 0 ||
1568 pull->group[pull->coord[c].params.group[1]].params.nat == 0)
1574 fprintf(fplog, "\n");
1575 if (pull->bPotential)
1577 fprintf(fplog, "Will apply potential COM pulling\n");
1579 if (pull->bConstraint)
1581 fprintf(fplog, "Will apply constraint COM pulling\n");
1583 fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
1584 pull->ncoord, pull->ncoord == 1 ? "" : "s",
1585 pull->ngroup, pull->ngroup == 1 ? "" : "s");
1588 fprintf(fplog, "with an absolute reference\n");
1591 for (g = 0; g < pull->ngroup; g++)
1593 if (pull->group[g].params.nat > 1 &&
1594 pull->group[g].params.pbcatom < 0)
1596 /* We are using cosine weighting */
1597 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1603 please_cite(fplog, "Engin2010");
1609 pull->dbuf_cyl = NULL;
1610 pull->bRefAt = FALSE;
1612 for (g = 0; g < pull->ngroup; g++)
1614 pull_group_work_t *pgrp;
1616 pgrp = &pull->group[g];
1617 pgrp->epgrppbc = epgrppbcNONE;
1618 if (pgrp->params.nat > 0)
1620 /* There is an issue when a group is used in multiple coordinates
1621 * and constraints are applied in different dimensions with atoms
1622 * frozen in some, but not all dimensions.
1623 * Since there is only one mass array per group, we can't have
1624 * frozen/non-frozen atoms for different coords at the same time.
1625 * But since this is a very exotic case, we don't check for this.
1626 * A warning is printed in init_pull_group_index.
1629 gmx_bool bConstraint;
1630 ivec pulldim, pulldim_con;
1632 /* Loop over all pull coordinates to see along which dimensions
1633 * this group is pulled and if it is involved in constraints.
1635 bConstraint = FALSE;
1636 clear_ivec(pulldim);
1637 clear_ivec(pulldim_con);
1638 for (c = 0; c < pull->ncoord; c++)
1640 if (pull->coord[c].params.group[0] == g ||
1641 pull->coord[c].params.group[1] == g ||
1642 (pull->coord[c].params.eGeom == epullgDIRRELATIVE &&
1643 (pull->coord[c].params.group[2] == g ||
1644 pull->coord[c].params.group[3] == g)))
1646 for (m = 0; m < DIM; m++)
1648 if (pull->coord[c].params.dim[m] == 1)
1652 if (pull->coord[c].params.eType == epullCONSTRAINT)
1662 /* Determine if we need to take PBC into account for calculating
1663 * the COM's of the pull groups.
1665 GMX_RELEASE_ASSERT(pull->npbcdim <= DIM, "npbcdim must be <= the number of dimensions");
1666 for (m = 0; m < pull->npbcdim; m++)
1668 GMX_RELEASE_ASSERT(m <= DIM, "npbcdim must be <= the number of dimensions");
1669 if (pulldim[m] == 1 && pgrp->params.nat > 1)
1671 if (pgrp->params.pbcatom >= 0)
1673 pgrp->epgrppbc = epgrppbcREFAT;
1674 pull->bRefAt = TRUE;
1678 if (pgrp->params.weight != NULL)
1680 gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1682 pgrp->epgrppbc = epgrppbcCOS;
1683 if (pull->cosdim >= 0 && pull->cosdim != m)
1685 gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
1691 /* Set the indices */
1692 init_pull_group_index(fplog, cr, g, pgrp,
1693 bConstraint, pulldim_con,
1698 /* Absolute reference, set the inverse mass to zero.
1699 * This is only relevant (and used) with constraint pulling.
1706 /* If we use cylinder coordinates, do some initialising for them */
1707 if (pull->bCylinder)
1709 snew(pull->dyna, pull->ncoord);
1711 for (c = 0; c < pull->ncoord; c++)
1713 const pull_coord_work_t *pcrd;
1715 pcrd = &pull->coord[c];
1717 if (pcrd->params.eGeom == epullgCYL)
1719 if (pull->group[pcrd->params.group[0]].params.nat == 0)
1721 gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
1727 /* We still need to initialize the PBC reference coordinates */
1728 pull->bSetPBCatoms = TRUE;
1730 /* Only do I/O when we are doing dynamics and if we are the MASTER */
1735 if (pull->params.nstxout > 0)
1737 pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
1740 if (pull->params.nstfout > 0)
1742 pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1750 static void destroy_pull_group(pull_group_work_t *pgrp)
1752 if (pgrp->ind_loc != pgrp->params.ind)
1754 sfree(pgrp->ind_loc);
1756 if (pgrp->weight_loc != pgrp->params.weight)
1758 sfree(pgrp->weight_loc);
1763 sfree(pgrp->params.ind);
1764 sfree(pgrp->params.weight);
1767 static void destroy_pull(struct pull_t *pull)
1771 for (i = 0; i < pull->ngroup; i++)
1773 destroy_pull_group(&pull->group[i]);
1775 for (i = 0; i < pull->ncoord; i++)
1777 if (pull->coord[i].params.eGeom == epullgCYL)
1779 destroy_pull_group(&(pull->dyna[i]));
1788 sfree(pull->dbuf_cyl);
1793 void finish_pull(struct pull_t *pull)
1797 gmx_fio_fclose(pull->out_x);
1801 gmx_fio_fclose(pull->out_f);
1807 gmx_bool pull_have_potential(const struct pull_t *pull)
1809 return pull->bPotential;
1812 gmx_bool pull_have_constraint(const struct pull_t *pull)
1814 return pull->bConstraint;