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/topology/mtop_util.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/gmxassert.h"
65 #include "gromacs/utility/smalloc.h"
67 static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
71 for (m = 0; m < DIM; m++)
75 fprintf(out, "\t%g", pgrp->x[m]);
80 static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd,
81 gmx_bool bPrintRefValue,
82 gmx_bool bPrintComponents)
84 fprintf(out, "\t%g", pcrd->value);
88 fprintf(out, "\t%g", pcrd->value_ref);
95 for (m = 0; m < DIM; m++)
99 fprintf(out, "\t%g", pcrd->dr[m]);
105 static void pull_print_x(FILE *out, t_pull *pull, double t)
110 fprintf(out, "%.4f", t);
112 for (c = 0; c < pull->ncoord; c++)
114 pcrd = &pull->coord[c];
116 if (pull->bPrintCOM1)
118 if (pcrd->eGeom == epullgCYL)
120 pull_print_group_x(out, pcrd->dim, &pull->dyna[c]);
124 pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[0]]);
127 if (pull->bPrintCOM2)
129 pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[1]]);
131 pull_print_coord_dr(out, pcrd,
132 pull->bPrintRefValue, pull->bPrintComp);
137 static void pull_print_f(FILE *out, t_pull *pull, double t)
141 fprintf(out, "%.4f", t);
143 for (c = 0; c < pull->ncoord; c++)
145 fprintf(out, "\t%g", pull->coord[c].f_scal);
150 void pull_print_output(t_pull *pull, gmx_int64_t step, double time)
152 if ((pull->nstxout != 0) && (step % pull->nstxout == 0))
154 pull_print_x(pull->out_x, pull, time);
157 if ((pull->nstfout != 0) && (step % pull->nstfout == 0))
159 pull_print_f(pull->out_f, pull, time);
163 static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv,
164 gmx_bool bCoord, unsigned long Flags)
168 char **setname, buf[10];
170 if (Flags & MD_APPENDFILES)
172 fp = gmx_fio_fopen(fn, "a+");
176 fp = gmx_fio_fopen(fn, "w+");
179 xvgr_header(fp, "Pull COM", "Time (ps)", "Position (nm)",
184 xvgr_header(fp, "Pull force", "Time (ps)", "Force (kJ/mol/nm)",
188 /* With default mdp options only the actual distance is printed,
189 * but optionally 2 COMs, the reference distance and distance
190 * components can also be printed.
192 snew(setname, pull->ncoord*(DIM + DIM + 1 + 1 + DIM));
194 for (c = 0; c < pull->ncoord; c++)
198 /* The order of this legend should match the order of printing
199 * the data in print_pull_x above.
202 if (pull->bPrintCOM1)
204 /* Legend for reference group position */
205 for (m = 0; m < DIM; m++)
207 if (pull->coord[c].dim[m])
209 sprintf(buf, "%d %s%d%c", c+1, "c", 1, 'X'+m);
210 setname[nsets] = gmx_strdup(buf);
215 if (pull->bPrintCOM2)
217 /* Legend for reference group position */
218 for (m = 0; m < DIM; m++)
220 if (pull->coord[c].dim[m])
222 sprintf(buf, "%d %s%d%c", c+1, "c", 2, 'X'+m);
223 setname[nsets] = gmx_strdup(buf);
228 /* The pull coord distance */
229 sprintf(buf, "%d", c+1);
230 setname[nsets] = gmx_strdup(buf);
232 if (pull->bPrintRefValue)
234 sprintf(buf, "%c%d", 'r', c+1);
235 setname[nsets] = gmx_strdup(buf);
238 if (pull->bPrintComp)
240 /* The pull coord distance components */
241 for (m = 0; m < DIM; m++)
243 if (pull->coord[c].dim[m])
245 sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
246 setname[nsets] = gmx_strdup(buf);
254 /* For the pull force we always only use one scalar */
255 sprintf(buf, "%d", c+1);
256 setname[nsets] = gmx_strdup(buf);
262 xvgr_legend(fp, nsets, (const char**)setname, oenv);
264 for (c = 0; c < nsets; c++)
274 /* Apply forces in a mass weighted fashion */
275 static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
276 const dvec f_pull, int sign, rvec *f)
279 double wmass, inv_wm;
281 inv_wm = pgrp->mwscale;
283 if (pgrp->nat == 1 && pgrp->nat_loc == 1)
285 /* Only one atom and our rank has this atom: we can skip
286 * the mass weighting, which means that this code also works
287 * for mass=0, e.g. with a virtual site.
289 for (m = 0; m < DIM; m++)
291 f[pgrp->ind_loc[0]][m] += sign*f_pull[m];
296 for (i = 0; i < pgrp->nat_loc; i++)
298 ii = pgrp->ind_loc[i];
299 wmass = md->massT[ii];
300 if (pgrp->weight_loc)
302 wmass *= pgrp->weight_loc[i];
305 for (m = 0; m < DIM; m++)
307 f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
313 /* Apply forces in a mass weighted fashion to a cylinder group */
314 static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr,
316 const dvec f_pull, double f_scal,
320 double mass, weight, inv_wm, dv_com;
322 inv_wm = pgrp->mwscale;
324 for (i = 0; i < pgrp->nat_loc; i++)
326 ii = pgrp->ind_loc[i];
327 mass = md->massT[ii];
328 weight = pgrp->weight_loc[i];
329 /* The stored axial distance from the cylinder center (dv) needs
330 * to be corrected for an offset (dv_corr), which was unknown when
333 dv_com = pgrp->dv[i] + dv_corr;
335 /* Here we not only add the pull force working along vec (f_pull),
336 * but also a radial component, due to the dependence of the weights
337 * on the radial distance.
339 for (m = 0; m < DIM; m++)
341 f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
342 pgrp->mdw[i][m]*dv_com*f_scal);
347 /* Apply torque forces in a mass weighted fashion to the groups that define
348 * the pull vector direction for pull coordinate pcrd.
350 static void apply_forces_vec_torque(const t_pull *pull,
351 const t_pull_coord *pcrd,
359 /* The component inpr along the pull vector is accounted for in the usual
360 * way. Here we account for the component perpendicular to vec.
363 for (m = 0; m < DIM; m++)
365 inpr += pcrd->dr[m]*pcrd->vec[m];
367 /* The torque force works along the component of the distance vector
368 * of between the two "usual" pull groups that is perpendicular to
369 * the pull vector. The magnitude of this force is the "usual" scale force
370 * multiplied with the ratio of the distance between the two "usual" pull
371 * groups and the distance between the two groups that define the vector.
373 for (m = 0; m < DIM; m++)
375 f_perp[m] = (pcrd->dr[m] - inpr*pcrd->vec[m])/pcrd->vec_len*pcrd->f_scal;
378 /* Apply the force to the groups defining the vector using opposite signs */
379 apply_forces_grp(&pull->group[pcrd->group[2]], md, f_perp, -1, f);
380 apply_forces_grp(&pull->group[pcrd->group[3]], md, f_perp, 1, f);
383 /* Apply forces in a mass weighted fashion */
384 static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
387 const t_pull_coord *pcrd;
389 for (c = 0; c < pull->ncoord; c++)
391 pcrd = &pull->coord[c];
393 if (pcrd->eGeom == epullgCYL)
398 apply_forces_cyl_grp(&pull->dyna[c], pcrd->cyl_dev, md,
399 pcrd->f, pcrd->f_scal, -1, f);
401 /* Sum the force along the vector and the radial force */
402 for (m = 0; m < DIM; m++)
404 f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
406 apply_forces_grp(&pull->group[pcrd->group[1]], md, f_tot, 1, f);
410 if (pcrd->eGeom == epullgDIRRELATIVE)
412 /* We need to apply the torque forces to the pull groups
413 * that define the pull vector.
415 apply_forces_vec_torque(pull, pcrd, md, f);
418 if (pull->group[pcrd->group[0]].nat > 0)
420 apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f);
422 apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
427 static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc)
432 max_d2 = GMX_DOUBLE_MAX;
434 for (m = 0; m < pbc->ndim_ePBC; m++)
436 if (pcrd->dim[m] != 0)
438 max_d2 = std::min(max_d2, static_cast<double>(norm2(pbc->box[m])));
445 /* This function returns the distance based on coordinates xg and xref.
446 * Note that the pull coordinate struct pcrd is not modified.
448 static void low_get_pull_coord_dr(const t_pull *pull,
449 const t_pull_coord *pcrd,
451 dvec xg, dvec xref, double max_dist2,
454 const t_pull_group *pgrp0;
456 dvec xrefr, dref = {0, 0, 0};
459 pgrp0 = &pull->group[pcrd->group[0]];
461 /* Only the first group can be an absolute reference, in that case nat=0 */
464 for (m = 0; m < DIM; m++)
466 xref[m] = pcrd->origin[m];
470 copy_dvec(xref, xrefr);
472 if (pcrd->eGeom == epullgDIRPBC)
474 for (m = 0; m < DIM; m++)
476 dref[m] = pcrd->value_ref*pcrd->vec[m];
478 /* Add the reference position, so we use the correct periodic image */
479 dvec_inc(xrefr, dref);
482 pbc_dx_d(pbc, xg, xrefr, dr);
484 for (m = 0; m < DIM; m++)
486 dr[m] *= pcrd->dim[m];
489 if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
491 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",
492 pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
495 if (pcrd->eGeom == epullgDIRPBC)
501 /* This function returns the distance based on the contents of the pull struct.
502 * pull->coord[coord_ind].dr, and potentially vector, are updated.
504 static void get_pull_coord_dr(t_pull *pull,
511 pcrd = &pull->coord[coord_ind];
513 if (pcrd->eGeom == epullgDIRPBC)
519 md2 = max_pull_distance2(pcrd, pbc);
522 if (pcrd->eGeom == epullgDIRRELATIVE)
524 /* We need to determine the pull vector */
525 const t_pull_group *pgrp2, *pgrp3;
529 pgrp2 = &pull->group[pcrd->group[2]];
530 pgrp3 = &pull->group[pcrd->group[3]];
532 pbc_dx_d(pbc, pgrp3->x, pgrp2->x, vec);
534 for (m = 0; m < DIM; m++)
536 vec[m] *= pcrd->dim[m];
538 pcrd->vec_len = dnorm(vec);
539 for (m = 0; m < DIM; m++)
541 pcrd->vec[m] = vec[m]/pcrd->vec_len;
545 fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
547 vec[XX], vec[YY], vec[ZZ],
548 pcrd->vec[XX], pcrd->vec[YY], pcrd->vec[ZZ]);
552 low_get_pull_coord_dr(pull, pcrd, pbc,
553 pull->group[pcrd->group[1]].x,
554 pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
559 static void update_pull_coord_reference_value(t_pull_coord *pcrd, double t)
561 /* With zero rate the reference value is set initially and doesn't change */
564 pcrd->value_ref = pcrd->init + pcrd->rate*t;
568 /* Calculates pull->coord[coord_ind].value.
569 * This function also updates pull->coord[coord_ind].dr.
571 static void get_pull_coord_distance(t_pull *pull,
578 pcrd = &pull->coord[coord_ind];
580 get_pull_coord_dr(pull, coord_ind, pbc);
585 /* Pull along the vector between the com's */
586 pcrd->value = dnorm(pcrd->dr);
590 case epullgDIRRELATIVE:
594 for (m = 0; m < DIM; m++)
596 pcrd->value += pcrd->vec[m]*pcrd->dr[m];
602 /* Returns the deviation from the reference value.
603 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
605 static double get_pull_coord_deviation(t_pull *pull,
610 static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
611 but is fairly benign */
615 pcrd = &pull->coord[coord_ind];
617 get_pull_coord_distance(pull, coord_ind, pbc);
619 update_pull_coord_reference_value(pcrd, t);
624 /* Pull along the vector between the com's */
625 if (pcrd->value_ref < 0 && !bWarned)
627 fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, pcrd->value_ref);
630 if (pcrd->value == 0)
632 /* With no vector we can not determine the direction for the force,
633 * so we set the force to zero.
639 /* Determine the deviation */
640 dev = pcrd->value - pcrd->value_ref;
645 case epullgDIRRELATIVE:
648 dev = pcrd->value - pcrd->value_ref;
655 void get_pull_coord_value(t_pull *pull,
660 get_pull_coord_distance(pull, coord_ind, pbc);
662 *value = pull->coord[coord_ind].value;
665 void clear_pull_forces(t_pull *pull)
669 /* Zeroing the forces is only required for constraint pulling.
670 * It can happen that multiple constraint steps need to be applied
671 * and therefore the constraint forces need to be accumulated.
673 for (i = 0; i < pull->ncoord; i++)
675 clear_dvec(pull->coord[i].f);
676 pull->coord[i].f_scal = 0;
680 /* Apply constraint using SHAKE */
681 static void do_constraint(t_pull *pull, t_pbc *pbc,
683 gmx_bool bMaster, tensor vir,
687 dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
688 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
689 dvec *rnew; /* current 'new' positions of the groups */
690 double *dr_tot; /* the total update of the coords */
693 double lambda, rm, invdt = 0;
694 gmx_bool bConverged_all, bConverged = FALSE;
695 int niter = 0, g, c, ii, j, m, max_iter = 100;
698 t_pull_group *pgrp0, *pgrp1;
701 snew(r_ij, pull->ncoord);
702 snew(dr_tot, pull->ncoord);
704 snew(rnew, pull->ngroup);
706 /* copy the current unconstrained positions for use in iterations. We
707 iterate until rinew[i] and rjnew[j] obey the constraints. Then
708 rinew - pull.x_unc[i] is the correction dr to group i */
709 for (g = 0; g < pull->ngroup; g++)
711 copy_dvec(pull->group[g].xp, rnew[g]);
714 /* Determine the constraint directions from the old positions */
715 for (c = 0; c < pull->ncoord; c++)
717 pcrd = &pull->coord[c];
719 if (pcrd->eType != epullCONSTRAINT)
724 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
725 * We don't modify dr and value anymore, so these values are also used
728 get_pull_coord_distance(pull, c, pbc);
732 fprintf(debug, "Pull coord %d dr %f %f %f\n",
733 c, pcrd->dr[XX], pcrd->dr[YY], pcrd->dr[ZZ]);
736 if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC)
738 /* Select the component along vec */
740 for (m = 0; m < DIM; m++)
742 a += pcrd->vec[m]*pcrd->dr[m];
744 for (m = 0; m < DIM; m++)
746 r_ij[c][m] = a*pcrd->vec[m];
751 copy_dvec(pcrd->dr, r_ij[c]);
754 if (dnorm2(r_ij[c]) == 0)
756 gmx_fatal(FARGS, "Distance for pull coordinate %d is zero with constraint pulling, which is not allowed.", c + 1);
760 bConverged_all = FALSE;
761 while (!bConverged_all && niter < max_iter)
763 bConverged_all = TRUE;
765 /* loop over all constraints */
766 for (c = 0; c < pull->ncoord; c++)
770 pcrd = &pull->coord[c];
772 if (pcrd->eType != epullCONSTRAINT)
777 update_pull_coord_reference_value(pcrd, t);
779 pgrp0 = &pull->group[pcrd->group[0]];
780 pgrp1 = &pull->group[pcrd->group[1]];
782 /* Get the current difference vector */
783 low_get_pull_coord_dr(pull, pcrd, pbc,
784 rnew[pcrd->group[1]],
785 rnew[pcrd->group[0]],
790 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
793 rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
798 if (pcrd->value_ref <= 0)
800 gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, pcrd->value_ref);
804 double q, c_a, c_b, c_c;
806 c_a = diprod(r_ij[c], r_ij[c]);
807 c_b = diprod(unc_ij, r_ij[c])*2;
808 c_c = diprod(unc_ij, unc_ij) - dsqr(pcrd->value_ref);
812 q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
817 q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
824 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
825 c_a, c_b, c_c, lambda);
829 /* The position corrections dr due to the constraints */
830 dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
831 dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
832 dr_tot[c] += -lambda*dnorm(r_ij[c]);
837 /* A 1-dimensional constraint along a vector */
839 for (m = 0; m < DIM; m++)
841 vec[m] = pcrd->vec[m];
842 a += unc_ij[m]*vec[m];
844 /* Select only the component along the vector */
845 dsvmul(a, vec, unc_ij);
846 lambda = a - pcrd->value_ref;
849 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
852 /* The position corrections dr due to the constraints */
853 dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
854 dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
855 dr_tot[c] += -lambda;
858 gmx_incons("Invalid enumeration value for eGeom");
859 /* Keep static analyzer happy */
871 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
872 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
874 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
875 rnew[g0][0], rnew[g0][1], rnew[g0][2],
876 rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
878 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
879 "", "", "", "", "", "", pcrd->value_ref);
881 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
882 dr0[0], dr0[1], dr0[2],
883 dr1[0], dr1[1], dr1[2],
887 /* Update the COMs with dr */
888 dvec_inc(rnew[pcrd->group[1]], dr1);
889 dvec_inc(rnew[pcrd->group[0]], dr0);
892 /* Check if all constraints are fullfilled now */
893 for (c = 0; c < pull->ncoord; c++)
895 pcrd = &pull->coord[c];
897 if (pcrd->eType != epullCONSTRAINT)
902 low_get_pull_coord_dr(pull, pcrd, pbc,
903 rnew[pcrd->group[1]],
904 rnew[pcrd->group[0]],
911 fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->constr_tol;
916 for (m = 0; m < DIM; m++)
918 vec[m] = pcrd->vec[m];
920 inpr = diprod(unc_ij, vec);
921 dsvmul(inpr, vec, unc_ij);
923 fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->constr_tol;
931 fprintf(debug, "NOT CONVERGED YET: Group %d:"
932 "d_ref = %f, current d = %f\n",
933 g, pcrd->value_ref, dnorm(unc_ij));
936 bConverged_all = FALSE;
941 /* if after all constraints are dealt with and bConverged is still TRUE
942 we're finished, if not we do another iteration */
944 if (niter > max_iter)
946 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
949 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
956 /* update atoms in the groups */
957 for (g = 0; g < pull->ngroup; g++)
959 const t_pull_group *pgrp;
962 pgrp = &pull->group[g];
964 /* get the final constraint displacement dr for group g */
965 dvec_sub(rnew[g], pgrp->xp, dr);
969 /* No displacement, no update necessary */
973 /* update the atom positions */
975 for (j = 0; j < pgrp->nat_loc; j++)
977 ii = pgrp->ind_loc[j];
978 if (pgrp->weight_loc)
980 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
982 for (m = 0; m < DIM; m++)
988 for (m = 0; m < DIM; m++)
990 v[ii][m] += invdt*tmp[m];
996 /* calculate the constraint forces, used for output and virial only */
997 for (c = 0; c < pull->ncoord; c++)
999 pcrd = &pull->coord[c];
1001 if (pcrd->eType != epullCONSTRAINT)
1006 pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
1008 if (vir != NULL && pcrd->eGeom != epullgDIRPBC && bMaster)
1012 /* Add the pull contribution to the virial */
1013 /* We have already checked above that r_ij[c] != 0 */
1014 f_invr = pcrd->f_scal/dnorm(r_ij[c]);
1016 for (j = 0; j < DIM; j++)
1018 for (m = 0; m < DIM; m++)
1020 vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
1026 /* finished! I hope. Give back some memory */
1032 /* Pulling with a harmonic umbrella potential or constant force */
1033 static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
1034 real *V, tensor vir, real *dVdl)
1037 double dev, ndr, invdr = 0;
1041 /* loop over the pull coordinates */
1044 for (c = 0; c < pull->ncoord; c++)
1046 pcrd = &pull->coord[c];
1048 if (pcrd->eType == epullCONSTRAINT)
1053 dev = get_pull_coord_deviation(pull, c, pbc, t);
1055 k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
1056 dkdl = pcrd->kB - pcrd->k;
1058 if (pcrd->eGeom == epullgDIST)
1060 ndr = dnorm(pcrd->dr);
1067 /* With an harmonic umbrella, the force is 0 at r=0,
1068 * so we can set invdr to any value.
1069 * With a constant force, the force at r=0 is not defined,
1070 * so we zero it (this is anyhow a very rare event).
1078 for (m = 0; m < DIM; m++)
1080 ndr += pcrd->vec[m]*pcrd->dr[m];
1084 switch (pcrd->eType)
1087 case epullFLATBOTTOM:
1088 /* The only difference between an umbrella and a flat-bottom
1089 * potential is that a flat-bottom is zero below zero.
1091 if (pcrd->eType == epullFLATBOTTOM && dev < 0)
1096 pcrd->f_scal = -k*dev;
1097 *V += 0.5* k*dsqr(dev);
1098 *dVdl += 0.5*dkdl*dsqr(dev);
1106 gmx_incons("Unsupported pull type in do_pull_pot");
1109 if (pcrd->eGeom == epullgDIST)
1111 for (m = 0; m < DIM; m++)
1113 pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
1118 for (m = 0; m < DIM; m++)
1120 pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
1124 if (vir != NULL && pcrd->eGeom != epullgDIRPBC)
1126 /* Add the pull contribution to the virial */
1127 for (j = 0; j < DIM; j++)
1129 for (m = 0; m < DIM; m++)
1131 vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
1138 real pull_potential(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1139 t_commrec *cr, double t, real lambda,
1140 rvec *x, rvec *f, tensor vir, real *dvdlambda)
1144 pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
1146 do_pull_pot(pull, pbc, t, lambda,
1147 &V, MASTER(cr) ? vir : NULL, &dVdl);
1149 /* Distribute forces over pulled groups */
1150 apply_forces(pull, md, f);
1157 return (MASTER(cr) ? V : 0.0);
1160 void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1161 t_commrec *cr, double dt, double t,
1162 rvec *x, rvec *xp, rvec *v, tensor vir)
1164 pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1166 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1169 static void make_local_pull_group(gmx_ga2la_t ga2la,
1170 t_pull_group *pg, int start, int end)
1175 for (i = 0; i < pg->nat; i++)
1180 if (!ga2la_get_home(ga2la, ii, &ii))
1185 if (ii >= start && ii < end)
1187 /* This is a home atom, add it to the local pull group */
1188 if (pg->nat_loc >= pg->nalloc_loc)
1190 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1191 srenew(pg->ind_loc, pg->nalloc_loc);
1192 if (pg->epgrppbc == epgrppbcCOS || pg->weight)
1194 srenew(pg->weight_loc, pg->nalloc_loc);
1197 pg->ind_loc[pg->nat_loc] = ii;
1200 pg->weight_loc[pg->nat_loc] = pg->weight[i];
1207 void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
1221 for (g = 0; g < pull->ngroup; g++)
1223 make_local_pull_group(ga2la, &pull->group[g],
1227 /* Since the PBC of atoms might have changed, we need to update the PBC */
1228 pull->bSetPBCatoms = TRUE;
1231 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
1232 int g, t_pull_group *pg,
1233 gmx_bool bConstraint, ivec pulldim_con,
1234 gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
1236 int i, ii, d, nfrozen, ndim;
1238 double tmass, wmass, wwmass;
1239 gmx_groups_t *groups;
1240 gmx_mtop_atomlookup_t alook;
1243 if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1245 /* There are no masses in the integrator.
1246 * But we still want to have the correct mass-weighted COMs.
1247 * So we store the real masses in the weights.
1248 * We do not set nweight, so these weights do not end up in the tpx file.
1250 if (pg->nweight == 0)
1252 snew(pg->weight, pg->nat);
1261 pg->weight_loc = NULL;
1265 pg->nat_loc = pg->nat;
1266 pg->ind_loc = pg->ind;
1267 if (pg->epgrppbc == epgrppbcCOS)
1269 snew(pg->weight_loc, pg->nat);
1273 pg->weight_loc = pg->weight;
1277 groups = &mtop->groups;
1279 alook = gmx_mtop_atomlookup_init(mtop);
1285 for (i = 0; i < pg->nat; i++)
1288 gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1289 if (bConstraint && ir->opts.nFreeze)
1291 for (d = 0; d < DIM; d++)
1293 if (pulldim_con[d] == 1 &&
1294 ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1300 if (ir->efep == efepNO)
1306 m = (1 - lambda)*atom->m + lambda*atom->mB;
1308 if (pg->nweight > 0)
1316 if (EI_ENERGY_MINIMIZATION(ir->eI))
1318 /* Move the mass to the weight */
1323 else if (ir->eI == eiBD)
1327 mbd = ir->bd_fric*ir->delta_t;
1331 if (groups->grpnr[egcTC] == NULL)
1333 mbd = ir->delta_t/ir->opts.tau_t[0];
1337 mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1349 gmx_mtop_atomlookup_destroy(alook);
1353 /* We can have single atom groups with zero mass with potential pulling
1354 * without cosine weighting.
1356 if (pg->nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1358 /* With one atom the mass doesn't matter */
1363 gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1364 pg->weight ? " weighted" : "", g);
1370 "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
1371 if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1373 fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1375 if (pg->epgrppbc == epgrppbcCOS)
1377 fprintf(fplog, ", cosine weighting will be used");
1379 fprintf(fplog, "\n");
1384 /* A value != 0 signals not frozen, it is updated later */
1390 for (d = 0; d < DIM; d++)
1392 ndim += pulldim_con[d]*pg->nat;
1394 if (fplog && nfrozen > 0 && nfrozen < ndim)
1397 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1398 " that are subject to pulling are frozen.\n"
1399 " For constraint pulling the whole group will be frozen.\n\n",
1407 void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
1408 gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1409 gmx_bool bOutFile, unsigned long Flags)
1417 pull->bPotential = FALSE;
1418 pull->bConstraint = FALSE;
1419 pull->bCylinder = FALSE;
1420 for (c = 0; c < pull->ncoord; c++)
1423 int calc_com_start, calc_com_end, g;
1425 pcrd = &pull->coord[c];
1427 if (pcrd->eType == epullCONSTRAINT)
1429 /* Check restrictions of the constraint pull code */
1430 if (pcrd->eGeom == epullgCYL ||
1431 pcrd->eGeom == epullgDIRRELATIVE)
1433 gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
1434 epull_names[pcrd->eType],
1435 epullg_names[pcrd->eGeom],
1436 epull_names[epullUMBRELLA]);
1439 pull->bConstraint = TRUE;
1443 pull->bPotential = TRUE;
1446 if (pcrd->eGeom == epullgCYL)
1448 pull->bCylinder = TRUE;
1450 /* We only need to calculate the plain COM of a group
1451 * when it is not only used as a cylinder group.
1453 calc_com_start = (pcrd->eGeom == epullgCYL ? 1 : 0);
1454 calc_com_end = (pcrd->eGeom == epullgDIRRELATIVE ? 4 : 2);
1456 for (g = calc_com_start; g <= calc_com_end; g++)
1458 pull->group[pcrd->group[g]].bCalcCOM = TRUE;
1461 /* With non-zero rate the reference value is set at every step */
1462 if (pcrd->rate == 0)
1464 /* Initialize the constant reference value */
1465 pcrd->value_ref = pcrd->init;
1469 pull->ePBC = ir->ePBC;
1472 case epbcNONE: pull->npbcdim = 0; break;
1473 case epbcXY: pull->npbcdim = 2; break;
1474 default: pull->npbcdim = 3; break;
1479 gmx_bool bAbs, bCos;
1482 for (c = 0; c < pull->ncoord; c++)
1484 if (pull->group[pull->coord[c].group[0]].nat == 0 ||
1485 pull->group[pull->coord[c].group[1]].nat == 0)
1491 fprintf(fplog, "\n");
1492 if (pull->bPotential)
1494 fprintf(fplog, "Will apply potential COM pulling\n");
1496 if (pull->bConstraint)
1498 fprintf(fplog, "Will apply constraint COM pulling\n");
1500 fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
1501 pull->ncoord, pull->ncoord == 1 ? "" : "s",
1502 pull->ngroup, pull->ngroup == 1 ? "" : "s");
1505 fprintf(fplog, "with an absolute reference\n");
1508 for (g = 0; g < pull->ngroup; g++)
1510 if (pull->group[g].nat > 1 &&
1511 pull->group[g].pbcatom < 0)
1513 /* We are using cosine weighting */
1514 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1520 please_cite(fplog, "Engin2010");
1526 pull->dbuf_cyl = NULL;
1527 pull->bRefAt = FALSE;
1529 for (g = 0; g < pull->ngroup; g++)
1531 pgrp = &pull->group[g];
1532 pgrp->epgrppbc = epgrppbcNONE;
1535 /* There is an issue when a group is used in multiple coordinates
1536 * and constraints are applied in different dimensions with atoms
1537 * frozen in some, but not all dimensions.
1538 * Since there is only one mass array per group, we can't have
1539 * frozen/non-frozen atoms for different coords at the same time.
1540 * But since this is a very exotic case, we don't check for this.
1541 * A warning is printed in init_pull_group_index.
1544 gmx_bool bConstraint;
1545 ivec pulldim, pulldim_con;
1547 /* Loop over all pull coordinates to see along which dimensions
1548 * this group is pulled and if it is involved in constraints.
1550 bConstraint = FALSE;
1551 clear_ivec(pulldim);
1552 clear_ivec(pulldim_con);
1553 for (c = 0; c < pull->ncoord; c++)
1555 if (pull->coord[c].group[0] == g ||
1556 pull->coord[c].group[1] == g ||
1557 (pull->coord[c].eGeom == epullgDIRRELATIVE &&
1558 (pull->coord[c].group[2] == g ||
1559 pull->coord[c].group[3] == g)))
1561 for (m = 0; m < DIM; m++)
1563 if (pull->coord[c].dim[m] == 1)
1567 if (pull->coord[c].eType == epullCONSTRAINT)
1577 /* Determine if we need to take PBC into account for calculating
1578 * the COM's of the pull groups.
1580 GMX_RELEASE_ASSERT(pull->npbcdim <= DIM, "npbcdim must be <= the number of dimensions");
1581 for (m = 0; m < pull->npbcdim; m++)
1583 GMX_RELEASE_ASSERT(m <= DIM, "npbcdim must be <= the number of dimensions");
1584 if (pulldim[m] == 1 && pgrp->nat > 1)
1586 if (pgrp->pbcatom >= 0)
1588 pgrp->epgrppbc = epgrppbcREFAT;
1589 pull->bRefAt = TRUE;
1595 gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1597 pgrp->epgrppbc = epgrppbcCOS;
1598 if (pull->cosdim >= 0 && pull->cosdim != m)
1600 gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
1606 /* Set the indices */
1607 init_pull_group_index(fplog, cr, g, pgrp,
1608 bConstraint, pulldim_con,
1613 /* Absolute reference, set the inverse mass to zero.
1614 * This is only relevant (and used) with constraint pulling.
1621 /* If we use cylinder coordinates, do some initialising for them */
1622 if (pull->bCylinder)
1624 snew(pull->dyna, pull->ncoord);
1626 for (c = 0; c < pull->ncoord; c++)
1628 const t_pull_coord *pcrd;
1630 pcrd = &pull->coord[c];
1632 if (pcrd->eGeom == epullgCYL)
1634 if (pull->group[pcrd->group[0]].nat == 0)
1636 gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
1642 /* We still need to initialize the PBC reference coordinates */
1643 pull->bSetPBCatoms = TRUE;
1645 /* Only do I/O when we are doing dynamics and if we are the MASTER */
1650 if (pull->nstxout > 0)
1652 pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
1655 if (pull->nstfout > 0)
1657 pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1663 void finish_pull(t_pull *pull)
1667 gmx_fio_fclose(pull->out_x);
1671 gmx_fio_fclose(pull->out_f);