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.
48 #include "gromacs/fileio/filenm.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/legacyheaders/copyrite.h"
52 #include "gromacs/legacyheaders/gmx_ga2la.h"
53 #include "gromacs/legacyheaders/macros.h"
54 #include "gromacs/legacyheaders/mdrun.h"
55 #include "gromacs/legacyheaders/names.h"
56 #include "gromacs/legacyheaders/network.h"
57 #include "gromacs/legacyheaders/typedefs.h"
58 #include "gromacs/legacyheaders/types/commrec.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/topology/mtop_util.h"
62 #include "gromacs/utility/futil.h"
63 #include "gromacs/utility/smalloc.h"
65 static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
69 for (m = 0; m < DIM; m++)
73 fprintf(out, "\t%g", pgrp->x[m]);
78 static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd,
79 gmx_bool bPrintRefValue,
80 gmx_bool bPrintComponents)
82 fprintf(out, "\t%g", pcrd->value);
86 fprintf(out, "\t%g", pcrd->value_ref);
93 for (m = 0; m < DIM; m++)
97 fprintf(out, "\t%g", pcrd->dr[m]);
103 static void pull_print_x(FILE *out, t_pull *pull, double t)
108 fprintf(out, "%.4f", t);
110 for (c = 0; c < pull->ncoord; c++)
112 pcrd = &pull->coord[c];
114 if (pull->bPrintCOM1)
116 if (pcrd->eGeom == epullgCYL)
118 pull_print_group_x(out, pcrd->dim, &pull->dyna[c]);
122 pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[0]]);
125 if (pull->bPrintCOM2)
127 pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[1]]);
129 pull_print_coord_dr(out, pcrd,
130 pull->bPrintRefValue, pull->bPrintComp);
135 static void pull_print_f(FILE *out, t_pull *pull, double t)
139 fprintf(out, "%.4f", t);
141 for (c = 0; c < pull->ncoord; c++)
143 fprintf(out, "\t%g", pull->coord[c].f_scal);
148 void pull_print_output(t_pull *pull, gmx_int64_t step, double time)
150 if ((pull->nstxout != 0) && (step % pull->nstxout == 0))
152 pull_print_x(pull->out_x, pull, time);
155 if ((pull->nstfout != 0) && (step % pull->nstfout == 0))
157 pull_print_f(pull->out_f, pull, time);
161 static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv,
162 gmx_bool bCoord, unsigned long Flags)
166 char **setname, buf[10];
168 if (Flags & MD_APPENDFILES)
170 fp = gmx_fio_fopen(fn, "a+");
174 fp = gmx_fio_fopen(fn, "w+");
177 xvgr_header(fp, "Pull COM", "Time (ps)", "Position (nm)",
182 xvgr_header(fp, "Pull force", "Time (ps)", "Force (kJ/mol/nm)",
186 /* With default mdp options only the actual distance is printed,
187 * but optionally 2 COMs, the reference distance and distance
188 * components can also be printed.
190 snew(setname, pull->ncoord*(DIM + DIM + 1 + 1 + DIM));
192 for (c = 0; c < pull->ncoord; c++)
196 /* The order of this legend should match the order of printing
197 * the data in print_pull_x above.
200 if (pull->bPrintCOM1)
202 /* Legend for reference group position */
203 for (m = 0; m < DIM; m++)
205 if (pull->coord[c].dim[m])
207 sprintf(buf, "%d %s%d%c", c+1, "c", 1, 'X'+m);
208 setname[nsets] = gmx_strdup(buf);
213 if (pull->bPrintCOM2)
215 /* Legend for reference group position */
216 for (m = 0; m < DIM; m++)
218 if (pull->coord[c].dim[m])
220 sprintf(buf, "%d %s%d%c", c+1, "c", 2, 'X'+m);
221 setname[nsets] = gmx_strdup(buf);
226 /* The pull coord distance */
227 sprintf(buf, "%d", c+1);
228 setname[nsets] = gmx_strdup(buf);
230 if (pull->bPrintRefValue)
232 sprintf(buf, "%c%d", 'r', c+1);
233 setname[nsets] = gmx_strdup(buf);
236 if (pull->bPrintComp)
238 /* The pull coord distance components */
239 for (m = 0; m < DIM; m++)
241 if (pull->coord[c].dim[m])
243 sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
244 setname[nsets] = gmx_strdup(buf);
252 /* For the pull force we always only use one scalar */
253 sprintf(buf, "%d", c+1);
254 setname[nsets] = gmx_strdup(buf);
260 xvgr_legend(fp, nsets, (const char**)setname, oenv);
262 for (c = 0; c < nsets; c++)
272 /* Apply forces in a mass weighted fashion */
273 static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
274 const dvec f_pull, int sign, rvec *f)
277 double wmass, inv_wm;
279 inv_wm = pgrp->mwscale;
281 if (pgrp->nat == 1 && pgrp->nat_loc == 1)
283 /* Only one atom and our rank has this atom: we can skip
284 * the mass weighting, which means that this code also works
285 * for mass=0, e.g. with a virtual site.
287 for (m = 0; m < DIM; m++)
289 f[pgrp->ind_loc[0]][m] += sign*f_pull[m];
294 for (i = 0; i < pgrp->nat_loc; i++)
296 ii = pgrp->ind_loc[i];
297 wmass = md->massT[ii];
298 if (pgrp->weight_loc)
300 wmass *= pgrp->weight_loc[i];
303 for (m = 0; m < DIM; m++)
305 f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
311 /* Apply forces in a mass weighted fashion to a cylinder group */
312 static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr,
314 const dvec f_pull, double f_scal,
318 double mass, weight, inv_wm, dv_com;
320 inv_wm = pgrp->mwscale;
322 for (i = 0; i < pgrp->nat_loc; i++)
324 ii = pgrp->ind_loc[i];
325 mass = md->massT[ii];
326 weight = pgrp->weight_loc[i];
327 /* The stored axial distance from the cylinder center (dv) needs
328 * to be corrected for an offset (dv_corr), which was unknown when
331 dv_com = pgrp->dv[i] + dv_corr;
333 /* Here we not only add the pull force working along vec (f_pull),
334 * but also a radial component, due to the dependence of the weights
335 * on the radial distance.
337 for (m = 0; m < DIM; m++)
339 f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
340 pgrp->mdw[i][m]*dv_com*f_scal);
345 /* Apply torque forces in a mass weighted fashion to the groups that define
346 * the pull vector direction for pull coordinate pcrd.
348 static void apply_forces_vec_torque(const t_pull *pull,
349 const t_pull_coord *pcrd,
357 /* The component inpr along the pull vector is accounted for in the usual
358 * way. Here we account for the component perpendicular to vec.
361 for (m = 0; m < DIM; m++)
363 inpr += pcrd->dr[m]*pcrd->vec[m];
365 /* The torque force works along the component of the distance vector
366 * of between the two "usual" pull groups that is perpendicular to
367 * the pull vector. The magnitude of this force is the "usual" scale force
368 * multiplied with the ratio of the distance between the two "usual" pull
369 * groups and the distance between the two groups that define the vector.
371 for (m = 0; m < DIM; m++)
373 f_perp[m] = (pcrd->dr[m] - inpr*pcrd->vec[m])/pcrd->vec_len*pcrd->f_scal;
376 /* Apply the force to the groups defining the vector using opposite signs */
377 apply_forces_grp(&pull->group[pcrd->group[2]], md, f_perp, -1, f);
378 apply_forces_grp(&pull->group[pcrd->group[3]], md, f_perp, 1, f);
381 /* Apply forces in a mass weighted fashion */
382 static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
385 const t_pull_coord *pcrd;
387 for (c = 0; c < pull->ncoord; c++)
389 pcrd = &pull->coord[c];
391 if (pcrd->eGeom == epullgCYL)
396 apply_forces_cyl_grp(&pull->dyna[c], pcrd->cyl_dev, md,
397 pcrd->f, pcrd->f_scal, -1, f);
399 /* Sum the force along the vector and the radial force */
400 for (m = 0; m < DIM; m++)
402 f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
404 apply_forces_grp(&pull->group[pcrd->group[1]], md, f_tot, 1, f);
408 if (pcrd->eGeom == epullgDIRRELATIVE)
410 /* We need to apply the torque forces to the pull groups
411 * that define the pull vector.
413 apply_forces_vec_torque(pull, pcrd, md, f);
416 if (pull->group[pcrd->group[0]].nat > 0)
418 apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f);
420 apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
425 static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc)
430 max_d2 = GMX_DOUBLE_MAX;
432 for (m = 0; m < pbc->ndim_ePBC; m++)
434 if (pcrd->dim[m] != 0)
436 max_d2 = std::min(max_d2, static_cast<double>(norm2(pbc->box[m])));
443 /* This function returns the distance based on coordinates xg and xref.
444 * Note that the pull coordinate struct pcrd is not modified.
446 static void low_get_pull_coord_dr(const t_pull *pull,
447 const t_pull_coord *pcrd,
449 dvec xg, dvec xref, double max_dist2,
452 const t_pull_group *pgrp0;
454 dvec xrefr, dref = {0, 0, 0};
457 pgrp0 = &pull->group[pcrd->group[0]];
459 /* Only the first group can be an absolute reference, in that case nat=0 */
462 for (m = 0; m < DIM; m++)
464 xref[m] = pcrd->origin[m];
468 copy_dvec(xref, xrefr);
470 if (pcrd->eGeom == epullgDIRPBC)
472 for (m = 0; m < DIM; m++)
474 dref[m] = pcrd->value_ref*pcrd->vec[m];
476 /* Add the reference position, so we use the correct periodic image */
477 dvec_inc(xrefr, dref);
480 pbc_dx_d(pbc, xg, xrefr, dr);
482 for (m = 0; m < DIM; m++)
484 dr[m] *= pcrd->dim[m];
487 if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
489 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",
490 pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
493 if (pcrd->eGeom == epullgDIRPBC)
499 /* This function returns the distance based on the contents of the pull struct.
500 * pull->coord[coord_ind].dr, and potentially vector, are updated.
502 static void get_pull_coord_dr(t_pull *pull,
509 pcrd = &pull->coord[coord_ind];
511 if (pcrd->eGeom == epullgDIRPBC)
517 md2 = max_pull_distance2(pcrd, pbc);
520 if (pcrd->eGeom == epullgDIRRELATIVE)
522 /* We need to determine the pull vector */
523 const t_pull_group *pgrp2, *pgrp3;
527 pgrp2 = &pull->group[pcrd->group[2]];
528 pgrp3 = &pull->group[pcrd->group[3]];
530 pbc_dx_d(pbc, pgrp3->x, pgrp2->x, vec);
532 for (m = 0; m < DIM; m++)
534 vec[m] *= pcrd->dim[m];
536 pcrd->vec_len = dnorm(vec);
537 for (m = 0; m < DIM; m++)
539 pcrd->vec[m] = vec[m]/pcrd->vec_len;
543 fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
545 vec[XX], vec[YY], vec[ZZ],
546 pcrd->vec[XX], pcrd->vec[YY], pcrd->vec[ZZ]);
550 low_get_pull_coord_dr(pull, pcrd, pbc,
551 pull->group[pcrd->group[1]].x,
552 pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
557 static void update_pull_coord_reference_value(t_pull_coord *pcrd, double t)
559 /* With zero rate the reference value is set initially and doesn't change */
562 pcrd->value_ref = pcrd->init + pcrd->rate*t;
566 /* Calculates pull->coord[coord_ind].value.
567 * This function also updates pull->coord[coord_ind].dr.
569 static void get_pull_coord_distance(t_pull *pull,
576 pcrd = &pull->coord[coord_ind];
578 get_pull_coord_dr(pull, coord_ind, pbc);
583 /* Pull along the vector between the com's */
584 pcrd->value = dnorm(pcrd->dr);
588 case epullgDIRRELATIVE:
592 for (m = 0; m < DIM; m++)
594 pcrd->value += pcrd->vec[m]*pcrd->dr[m];
600 /* Returns the deviation from the reference value.
601 * Updates pull->coord[coord_ind].dr, .value and .value_ref.
603 static double get_pull_coord_deviation(t_pull *pull,
608 static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
609 but is fairly benign */
613 pcrd = &pull->coord[coord_ind];
615 get_pull_coord_distance(pull, coord_ind, pbc);
617 update_pull_coord_reference_value(pcrd, t);
622 /* Pull along the vector between the com's */
623 if (pcrd->value_ref < 0 && !bWarned)
625 fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, pcrd->value_ref);
628 if (pcrd->value == 0)
630 /* With no vector we can not determine the direction for the force,
631 * so we set the force to zero.
637 /* Determine the deviation */
638 dev = pcrd->value - pcrd->value_ref;
643 case epullgDIRRELATIVE:
646 dev = pcrd->value - pcrd->value_ref;
653 void get_pull_coord_value(t_pull *pull,
658 get_pull_coord_distance(pull, coord_ind, pbc);
660 *value = pull->coord[coord_ind].value;
663 void clear_pull_forces(t_pull *pull)
667 /* Zeroing the forces is only required for constraint pulling.
668 * It can happen that multiple constraint steps need to be applied
669 * and therefore the constraint forces need to be accumulated.
671 for (i = 0; i < pull->ncoord; i++)
673 clear_dvec(pull->coord[i].f);
674 pull->coord[i].f_scal = 0;
678 /* Apply constraint using SHAKE */
679 static void do_constraint(t_pull *pull, t_pbc *pbc,
681 gmx_bool bMaster, tensor vir,
685 dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
686 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
687 dvec *rnew; /* current 'new' positions of the groups */
688 double *dr_tot; /* the total update of the coords */
691 double lambda, rm, invdt = 0;
692 gmx_bool bConverged_all, bConverged = FALSE;
693 int niter = 0, g, c, ii, j, m, max_iter = 100;
696 t_pull_group *pgrp0, *pgrp1;
699 snew(r_ij, pull->ncoord);
700 snew(dr_tot, pull->ncoord);
702 snew(rnew, pull->ngroup);
704 /* copy the current unconstrained positions for use in iterations. We
705 iterate until rinew[i] and rjnew[j] obey the constraints. Then
706 rinew - pull.x_unc[i] is the correction dr to group i */
707 for (g = 0; g < pull->ngroup; g++)
709 copy_dvec(pull->group[g].xp, rnew[g]);
712 /* Determine the constraint directions from the old positions */
713 for (c = 0; c < pull->ncoord; c++)
715 pcrd = &pull->coord[c];
717 if (pcrd->eType != epullCONSTRAINT)
722 /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
723 * We don't modify dr and value anymore, so these values are also used
726 get_pull_coord_distance(pull, c, pbc);
730 fprintf(debug, "Pull coord %d dr %f %f %f\n",
731 c, pcrd->dr[XX], pcrd->dr[YY], pcrd->dr[ZZ]);
734 if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC)
736 /* Select the component along vec */
738 for (m = 0; m < DIM; m++)
740 a += pcrd->vec[m]*pcrd->dr[m];
742 for (m = 0; m < DIM; m++)
744 r_ij[c][m] = a*pcrd->vec[m];
749 copy_dvec(pcrd->dr, r_ij[c]);
752 if (dnorm2(r_ij[c]) == 0)
754 gmx_fatal(FARGS, "Distance for pull coordinate %d is zero with constraint pulling, which is not allowed.", c + 1);
758 bConverged_all = FALSE;
759 while (!bConverged_all && niter < max_iter)
761 bConverged_all = TRUE;
763 /* loop over all constraints */
764 for (c = 0; c < pull->ncoord; c++)
768 pcrd = &pull->coord[c];
770 if (pcrd->eType != epullCONSTRAINT)
775 update_pull_coord_reference_value(pcrd, t);
777 pgrp0 = &pull->group[pcrd->group[0]];
778 pgrp1 = &pull->group[pcrd->group[1]];
780 /* Get the current difference vector */
781 low_get_pull_coord_dr(pull, pcrd, pbc,
782 rnew[pcrd->group[1]],
783 rnew[pcrd->group[0]],
788 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
791 rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
796 if (pcrd->value_ref <= 0)
798 gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, pcrd->value_ref);
802 double q, c_a, c_b, c_c;
804 c_a = diprod(r_ij[c], r_ij[c]);
805 c_b = diprod(unc_ij, r_ij[c])*2;
806 c_c = diprod(unc_ij, unc_ij) - dsqr(pcrd->value_ref);
810 q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
815 q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
822 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
823 c_a, c_b, c_c, lambda);
827 /* The position corrections dr due to the constraints */
828 dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
829 dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
830 dr_tot[c] += -lambda*dnorm(r_ij[c]);
835 /* A 1-dimensional constraint along a vector */
837 for (m = 0; m < DIM; m++)
839 vec[m] = pcrd->vec[m];
840 a += unc_ij[m]*vec[m];
842 /* Select only the component along the vector */
843 dsvmul(a, vec, unc_ij);
844 lambda = a - pcrd->value_ref;
847 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
850 /* The position corrections dr due to the constraints */
851 dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
852 dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
853 dr_tot[c] += -lambda;
856 gmx_incons("Invalid enumeration value for eGeom");
857 /* Keep static analyzer happy */
869 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
870 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
872 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
873 rnew[g0][0], rnew[g0][1], rnew[g0][2],
874 rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
876 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
877 "", "", "", "", "", "", pcrd->value_ref);
879 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
880 dr0[0], dr0[1], dr0[2],
881 dr1[0], dr1[1], dr1[2],
885 /* Update the COMs with dr */
886 dvec_inc(rnew[pcrd->group[1]], dr1);
887 dvec_inc(rnew[pcrd->group[0]], dr0);
890 /* Check if all constraints are fullfilled now */
891 for (c = 0; c < pull->ncoord; c++)
893 pcrd = &pull->coord[c];
895 if (pcrd->eType != epullCONSTRAINT)
900 low_get_pull_coord_dr(pull, pcrd, pbc,
901 rnew[pcrd->group[1]],
902 rnew[pcrd->group[0]],
909 fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->constr_tol;
914 for (m = 0; m < DIM; m++)
916 vec[m] = pcrd->vec[m];
918 inpr = diprod(unc_ij, vec);
919 dsvmul(inpr, vec, unc_ij);
921 fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->constr_tol;
929 fprintf(debug, "NOT CONVERGED YET: Group %d:"
930 "d_ref = %f, current d = %f\n",
931 g, pcrd->value_ref, dnorm(unc_ij));
934 bConverged_all = FALSE;
939 /* if after all constraints are dealt with and bConverged is still TRUE
940 we're finished, if not we do another iteration */
942 if (niter > max_iter)
944 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
947 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
954 /* update atoms in the groups */
955 for (g = 0; g < pull->ngroup; g++)
957 const t_pull_group *pgrp;
960 pgrp = &pull->group[g];
962 /* get the final constraint displacement dr for group g */
963 dvec_sub(rnew[g], pgrp->xp, dr);
967 /* No displacement, no update necessary */
971 /* update the atom positions */
973 for (j = 0; j < pgrp->nat_loc; j++)
975 ii = pgrp->ind_loc[j];
976 if (pgrp->weight_loc)
978 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
980 for (m = 0; m < DIM; m++)
986 for (m = 0; m < DIM; m++)
988 v[ii][m] += invdt*tmp[m];
994 /* calculate the constraint forces, used for output and virial only */
995 for (c = 0; c < pull->ncoord; c++)
997 pcrd = &pull->coord[c];
999 if (pcrd->eType != epullCONSTRAINT)
1004 pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
1006 if (vir != NULL && pcrd->eGeom != epullgDIRPBC && bMaster)
1010 /* Add the pull contribution to the virial */
1011 /* We have already checked above that r_ij[c] != 0 */
1012 f_invr = pcrd->f_scal/dnorm(r_ij[c]);
1014 for (j = 0; j < DIM; j++)
1016 for (m = 0; m < DIM; m++)
1018 vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
1024 /* finished! I hope. Give back some memory */
1030 /* Pulling with a harmonic umbrella potential or constant force */
1031 static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
1032 real *V, tensor vir, real *dVdl)
1035 double dev, ndr, invdr = 0;
1039 /* loop over the pull coordinates */
1042 for (c = 0; c < pull->ncoord; c++)
1044 pcrd = &pull->coord[c];
1046 if (pcrd->eType == epullCONSTRAINT)
1051 dev = get_pull_coord_deviation(pull, c, pbc, t);
1053 k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
1054 dkdl = pcrd->kB - pcrd->k;
1056 if (pcrd->eGeom == epullgDIST)
1058 ndr = dnorm(pcrd->dr);
1065 /* With an harmonic umbrella, the force is 0 at r=0,
1066 * so we can set invdr to any value.
1067 * With a constant force, the force at r=0 is not defined,
1068 * so we zero it (this is anyhow a very rare event).
1076 for (m = 0; m < DIM; m++)
1078 ndr += pcrd->vec[m]*pcrd->dr[m];
1082 switch (pcrd->eType)
1085 case epullFLATBOTTOM:
1086 /* The only difference between an umbrella and a flat-bottom
1087 * potential is that a flat-bottom is zero below zero.
1089 if (pcrd->eType == epullFLATBOTTOM && dev < 0)
1094 pcrd->f_scal = -k*dev;
1095 *V += 0.5* k*dsqr(dev);
1096 *dVdl += 0.5*dkdl*dsqr(dev);
1104 gmx_incons("Unsupported pull type in do_pull_pot");
1107 if (pcrd->eGeom == epullgDIST)
1109 for (m = 0; m < DIM; m++)
1111 pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
1116 for (m = 0; m < DIM; m++)
1118 pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
1122 if (vir != NULL && pcrd->eGeom != epullgDIRPBC)
1124 /* Add the pull contribution to the virial */
1125 for (j = 0; j < DIM; j++)
1127 for (m = 0; m < DIM; m++)
1129 vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
1136 real pull_potential(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1137 t_commrec *cr, double t, real lambda,
1138 rvec *x, rvec *f, tensor vir, real *dvdlambda)
1142 pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
1144 do_pull_pot(pull, pbc, t, lambda,
1145 &V, MASTER(cr) ? vir : NULL, &dVdl);
1147 /* Distribute forces over pulled groups */
1148 apply_forces(pull, md, f);
1155 return (MASTER(cr) ? V : 0.0);
1158 void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
1159 t_commrec *cr, double dt, double t,
1160 rvec *x, rvec *xp, rvec *v, tensor vir)
1162 pull_calc_coms(cr, pull, md, pbc, t, x, xp);
1164 do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
1167 static void make_local_pull_group(gmx_ga2la_t ga2la,
1168 t_pull_group *pg, int start, int end)
1173 for (i = 0; i < pg->nat; i++)
1178 if (!ga2la_get_home(ga2la, ii, &ii))
1183 if (ii >= start && ii < end)
1185 /* This is a home atom, add it to the local pull group */
1186 if (pg->nat_loc >= pg->nalloc_loc)
1188 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
1189 srenew(pg->ind_loc, pg->nalloc_loc);
1190 if (pg->epgrppbc == epgrppbcCOS || pg->weight)
1192 srenew(pg->weight_loc, pg->nalloc_loc);
1195 pg->ind_loc[pg->nat_loc] = ii;
1198 pg->weight_loc[pg->nat_loc] = pg->weight[i];
1205 void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
1219 for (g = 0; g < pull->ngroup; g++)
1221 make_local_pull_group(ga2la, &pull->group[g],
1225 /* Since the PBC of atoms might have changed, we need to update the PBC */
1226 pull->bSetPBCatoms = TRUE;
1229 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
1230 int g, t_pull_group *pg,
1231 gmx_bool bConstraint, ivec pulldim_con,
1232 gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
1234 int i, ii, d, nfrozen, ndim;
1236 double tmass, wmass, wwmass;
1237 gmx_groups_t *groups;
1238 gmx_mtop_atomlookup_t alook;
1241 if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1243 /* There are no masses in the integrator.
1244 * But we still want to have the correct mass-weighted COMs.
1245 * So we store the real masses in the weights.
1246 * We do not set nweight, so these weights do not end up in the tpx file.
1248 if (pg->nweight == 0)
1250 snew(pg->weight, pg->nat);
1259 pg->weight_loc = NULL;
1263 pg->nat_loc = pg->nat;
1264 pg->ind_loc = pg->ind;
1265 if (pg->epgrppbc == epgrppbcCOS)
1267 snew(pg->weight_loc, pg->nat);
1271 pg->weight_loc = pg->weight;
1275 groups = &mtop->groups;
1277 alook = gmx_mtop_atomlookup_init(mtop);
1283 for (i = 0; i < pg->nat; i++)
1286 gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1287 if (bConstraint && ir->opts.nFreeze)
1289 for (d = 0; d < DIM; d++)
1291 if (pulldim_con[d] == 1 &&
1292 ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1298 if (ir->efep == efepNO)
1304 m = (1 - lambda)*atom->m + lambda*atom->mB;
1306 if (pg->nweight > 0)
1314 if (EI_ENERGY_MINIMIZATION(ir->eI))
1316 /* Move the mass to the weight */
1321 else if (ir->eI == eiBD)
1325 mbd = ir->bd_fric*ir->delta_t;
1329 if (groups->grpnr[egcTC] == NULL)
1331 mbd = ir->delta_t/ir->opts.tau_t[0];
1335 mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1347 gmx_mtop_atomlookup_destroy(alook);
1351 /* We can have single atom groups with zero mass with potential pulling
1352 * without cosine weighting.
1354 if (pg->nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
1356 /* With one atom the mass doesn't matter */
1361 gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1362 pg->weight ? " weighted" : "", g);
1368 "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
1369 if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1371 fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1373 if (pg->epgrppbc == epgrppbcCOS)
1375 fprintf(fplog, ", cosine weighting will be used");
1377 fprintf(fplog, "\n");
1382 /* A value != 0 signals not frozen, it is updated later */
1388 for (d = 0; d < DIM; d++)
1390 ndim += pulldim_con[d]*pg->nat;
1392 if (fplog && nfrozen > 0 && nfrozen < ndim)
1395 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1396 " that are subject to pulling are frozen.\n"
1397 " For constraint pulling the whole group will be frozen.\n\n",
1405 void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
1406 gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1407 gmx_bool bOutFile, unsigned long Flags)
1415 pull->bPotential = FALSE;
1416 pull->bConstraint = FALSE;
1417 pull->bCylinder = FALSE;
1418 for (c = 0; c < pull->ncoord; c++)
1421 int calc_com_start, calc_com_end, g;
1423 pcrd = &pull->coord[c];
1425 if (pcrd->eType == epullCONSTRAINT)
1427 /* Check restrictions of the constraint pull code */
1428 if (pcrd->eGeom == epullgCYL ||
1429 pcrd->eGeom == epullgDIRRELATIVE)
1431 gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
1432 epull_names[pcrd->eType],
1433 epullg_names[pcrd->eGeom],
1434 epull_names[epullUMBRELLA]);
1437 pull->bConstraint = TRUE;
1441 pull->bPotential = TRUE;
1444 if (pcrd->eGeom == epullgCYL)
1446 pull->bCylinder = TRUE;
1448 /* We only need to calculate the plain COM of a group
1449 * when it is not only used as a cylinder group.
1451 calc_com_start = (pcrd->eGeom == epullgCYL ? 1 : 0);
1452 calc_com_end = (pcrd->eGeom == epullgDIRRELATIVE ? 4 : 2);
1454 for (g = calc_com_start; g <= calc_com_end; g++)
1456 pull->group[pcrd->group[g]].bCalcCOM = TRUE;
1459 /* With non-zero rate the reference value is set at every step */
1460 if (pcrd->rate == 0)
1462 /* Initialize the constant reference value */
1463 pcrd->value_ref = pcrd->init;
1467 pull->ePBC = ir->ePBC;
1470 case epbcNONE: pull->npbcdim = 0; break;
1471 case epbcXY: pull->npbcdim = 2; break;
1472 default: pull->npbcdim = 3; break;
1477 gmx_bool bAbs, bCos;
1480 for (c = 0; c < pull->ncoord; c++)
1482 if (pull->group[pull->coord[c].group[0]].nat == 0 ||
1483 pull->group[pull->coord[c].group[1]].nat == 0)
1489 fprintf(fplog, "\n");
1490 if (pull->bPotential)
1492 fprintf(fplog, "Will apply potential COM pulling\n");
1494 if (pull->bConstraint)
1496 fprintf(fplog, "Will apply constraint COM pulling\n");
1498 fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
1499 pull->ncoord, pull->ncoord == 1 ? "" : "s",
1500 pull->ngroup, pull->ngroup == 1 ? "" : "s");
1503 fprintf(fplog, "with an absolute reference\n");
1506 for (g = 0; g < pull->ngroup; g++)
1508 if (pull->group[g].nat > 1 &&
1509 pull->group[g].pbcatom < 0)
1511 /* We are using cosine weighting */
1512 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1518 please_cite(fplog, "Engin2010");
1524 pull->dbuf_cyl = NULL;
1525 pull->bRefAt = FALSE;
1527 for (g = 0; g < pull->ngroup; g++)
1529 pgrp = &pull->group[g];
1530 pgrp->epgrppbc = epgrppbcNONE;
1533 /* There is an issue when a group is used in multiple coordinates
1534 * and constraints are applied in different dimensions with atoms
1535 * frozen in some, but not all dimensions.
1536 * Since there is only one mass array per group, we can't have
1537 * frozen/non-frozen atoms for different coords at the same time.
1538 * But since this is a very exotic case, we don't check for this.
1539 * A warning is printed in init_pull_group_index.
1542 gmx_bool bConstraint;
1543 ivec pulldim, pulldim_con;
1545 /* Loop over all pull coordinates to see along which dimensions
1546 * this group is pulled and if it is involved in constraints.
1548 bConstraint = FALSE;
1549 clear_ivec(pulldim);
1550 clear_ivec(pulldim_con);
1551 for (c = 0; c < pull->ncoord; c++)
1553 if (pull->coord[c].group[0] == g ||
1554 pull->coord[c].group[1] == g ||
1555 (pull->coord[c].eGeom == epullgDIRRELATIVE &&
1556 (pull->coord[c].group[2] == g ||
1557 pull->coord[c].group[3] == g)))
1559 for (m = 0; m < DIM; m++)
1561 if (pull->coord[c].dim[m] == 1)
1565 if (pull->coord[c].eType == epullCONSTRAINT)
1575 /* Determine if we need to take PBC into account for calculating
1576 * the COM's of the pull groups.
1578 for (m = 0; m < pull->npbcdim; m++)
1580 if (pulldim[m] == 1 && pgrp->nat > 1)
1582 if (pgrp->pbcatom >= 0)
1584 pgrp->epgrppbc = epgrppbcREFAT;
1585 pull->bRefAt = TRUE;
1591 gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1593 pgrp->epgrppbc = epgrppbcCOS;
1594 if (pull->cosdim >= 0 && pull->cosdim != m)
1596 gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
1602 /* Set the indices */
1603 init_pull_group_index(fplog, cr, g, pgrp,
1604 bConstraint, pulldim_con,
1609 /* Absolute reference, set the inverse mass to zero.
1610 * This is only relevant (and used) with constraint pulling.
1617 /* If we use cylinder coordinates, do some initialising for them */
1618 if (pull->bCylinder)
1620 snew(pull->dyna, pull->ncoord);
1622 for (c = 0; c < pull->ncoord; c++)
1624 const t_pull_coord *pcrd;
1626 pcrd = &pull->coord[c];
1628 if (pcrd->eGeom == epullgCYL)
1630 if (pull->group[pcrd->group[0]].nat == 0)
1632 gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
1638 /* We still need to initialize the PBC reference coordinates */
1639 pull->bSetPBCatoms = TRUE;
1641 /* Only do I/O when we are doing dynamics and if we are the MASTER */
1646 if (pull->nstxout > 0)
1648 pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
1651 if (pull->nstfout > 0)
1653 pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1659 void finish_pull(t_pull *pull)
1663 gmx_fio_fclose(pull->out_x);
1667 gmx_fio_fclose(pull->out_f);