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, 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.
46 #include "gromacs/utility/futil.h"
47 #include "gromacs/legacyheaders/typedefs.h"
48 #include "gromacs/legacyheaders/types/commrec.h"
49 #include "gromacs/legacyheaders/network.h"
51 #include "gromacs/legacyheaders/names.h"
52 #include "gromacs/pbcutil/pbc.h"
53 #include "gromacs/topology/mtop_util.h"
54 #include "gromacs/legacyheaders/mdrun.h"
55 #include "gromacs/legacyheaders/gmx_ga2la.h"
56 #include "gromacs/legacyheaders/copyrite.h"
57 #include "gromacs/legacyheaders/macros.h"
59 #include "gromacs/fileio/filenm.h"
60 #include "gromacs/fileio/gmxfio.h"
61 #include "gromacs/fileio/xvgr.h"
62 #include "gromacs/math/vec.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, ivec dim, const t_pull_coord *pcrd)
82 for (m = 0; m < DIM; m++)
86 fprintf(out, "\t%g", pcrd->dr[m]);
91 static void pull_print_x(FILE *out, t_pull *pull, double t)
94 const t_pull_coord *pcrd;
96 fprintf(out, "%.4f", t);
98 for (c = 0; c < pull->ncoord; c++)
100 pcrd = &pull->coord[c];
106 pull_print_group_x(out, pull->dim, &pull->dyna[c]);
110 pull_print_group_x(out, pull->dim, &pull->group[pcrd->group[0]]);
113 pull_print_coord_dr(out, pull->dim, pcrd);
118 static void pull_print_f(FILE *out, t_pull *pull, double t)
122 fprintf(out, "%.4f", t);
124 for (c = 0; c < pull->ncoord; c++)
126 fprintf(out, "\t%g", pull->coord[c].f_scal);
131 void pull_print_output(t_pull *pull, gmx_int64_t step, double time)
133 if ((pull->nstxout != 0) && (step % pull->nstxout == 0))
135 pull_print_x(pull->out_x, pull, time);
138 if ((pull->nstfout != 0) && (step % pull->nstfout == 0))
140 pull_print_f(pull->out_f, pull, time);
144 static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv,
145 gmx_bool bCoord, unsigned long Flags)
149 char **setname, buf[10];
151 if (Flags & MD_APPENDFILES)
153 fp = gmx_fio_fopen(fn, "a+");
157 fp = gmx_fio_fopen(fn, "w+");
160 xvgr_header(fp, "Pull COM", "Time (ps)", "Position (nm)",
165 xvgr_header(fp, "Pull force", "Time (ps)", "Force (kJ/mol/nm)",
169 snew(setname, 2*pull->ncoord*DIM);
171 for (c = 0; c < pull->ncoord; c++)
177 for (m = 0; m < DIM; m++)
181 sprintf(buf, "%d %s%c", c+1, "c", 'X'+m);
182 setname[nsets] = gmx_strdup(buf);
187 for (m = 0; m < DIM; m++)
191 sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
192 setname[nsets] = gmx_strdup(buf);
199 sprintf(buf, "%d", c+1);
200 setname[nsets] = gmx_strdup(buf);
206 xvgr_legend(fp, nsets, (const char**)setname, oenv);
208 for (c = 0; c < nsets; c++)
218 /* Apply forces in a mass weighted fashion */
219 static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
220 const dvec f_pull, int sign, rvec *f)
223 double wmass, inv_wm;
225 inv_wm = pgrp->wscale*pgrp->invtm;
227 for (i = 0; i < pgrp->nat_loc; i++)
229 ii = pgrp->ind_loc[i];
230 wmass = md->massT[ii];
231 if (pgrp->weight_loc)
233 wmass *= pgrp->weight_loc[i];
236 for (m = 0; m < DIM; m++)
238 f[ii][m] += sign * wmass * f_pull[m] * inv_wm;
243 /* Apply forces in a mass weighted fashion */
244 static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
247 const t_pull_coord *pcrd;
249 for (c = 0; c < pull->ncoord; c++)
251 pcrd = &pull->coord[c];
255 apply_forces_grp(&pull->dyna[c], md, pcrd->f, -1, f);
259 if (pull->group[pcrd->group[0]].nat > 0)
261 apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f);
264 apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
268 static double max_pull_distance2(const t_pull *pull, const t_pbc *pbc)
273 max_d2 = GMX_DOUBLE_MAX;
275 if (pull->eGeom != epullgDIRPBC)
277 for (m = 0; m < pbc->ndim_ePBC; m++)
279 if (pull->dim[m] != 0)
281 max_d2 = min(max_d2, norm2(pbc->box[m]));
289 static void low_get_pull_coord_dr(const t_pull *pull,
290 const t_pull_coord *pcrd,
291 const t_pbc *pbc, double t,
292 dvec xg, dvec xref, double max_dist2,
295 const t_pull_group *pgrp0, *pgrp1;
297 dvec xrefr, dref = {0, 0, 0};
300 pgrp0 = &pull->group[pcrd->group[0]];
301 pgrp1 = &pull->group[pcrd->group[1]];
303 /* Only the first group can be an absolute reference, in that case nat=0 */
306 for (m = 0; m < DIM; m++)
308 xref[m] = pcrd->origin[m];
312 copy_dvec(xref, xrefr);
314 if (pull->eGeom == epullgDIRPBC)
316 for (m = 0; m < DIM; m++)
318 dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m];
320 /* Add the reference position, so we use the correct periodic image */
321 dvec_inc(xrefr, dref);
324 pbc_dx_d(pbc, xg, xrefr, dr);
326 for (m = 0; m < DIM; m++)
328 dr[m] *= pull->dim[m];
331 if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
333 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",
334 pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
337 if (pull->eGeom == epullgDIRPBC)
343 static void get_pull_coord_dr(const t_pull *pull,
345 const t_pbc *pbc, double t,
349 const t_pull_coord *pcrd;
351 if (pull->eGeom == epullgDIRPBC)
357 md2 = max_pull_distance2(pull, pbc);
360 pcrd = &pull->coord[coord_ind];
362 low_get_pull_coord_dr(pull, pcrd, pbc, t,
363 pull->group[pcrd->group[1]].x,
364 PULL_CYL(pull) ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
369 void get_pull_coord_distance(const t_pull *pull,
371 const t_pbc *pbc, double t,
372 dvec dr, double *dev)
374 static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
375 but is fairly benign */
376 const t_pull_coord *pcrd;
378 double ref, drs, inpr;
380 pcrd = &pull->coord[coord_ind];
382 get_pull_coord_dr(pull, coord_ind, pbc, t, dr);
384 ref = pcrd->init + pcrd->rate*t;
389 /* Pull along the vector between the com's */
390 if (ref < 0 && !bWarned)
392 fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref);
398 /* With no vector we can not determine the direction for the force,
399 * so we set the force to zero.
405 /* Determine the deviation */
414 for (m = 0; m < DIM; m++)
416 inpr += pcrd->vec[m]*dr[m];
423 void clear_pull_forces(t_pull *pull)
427 /* Zeroing the forces is only required for constraint pulling.
428 * It can happen that multiple constraint steps need to be applied
429 * and therefore the constraint forces need to be accumulated.
431 for (i = 0; i < pull->ncoord; i++)
433 clear_dvec(pull->coord[i].f);
434 pull->coord[i].f_scal = 0;
438 /* Apply constraint using SHAKE */
439 static void do_constraint(t_pull *pull, t_pbc *pbc,
441 gmx_bool bMaster, tensor vir,
445 dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
446 dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
447 dvec *rnew; /* current 'new' positions of the groups */
448 double *dr_tot; /* the total update of the coords */
452 double lambda, rm, mass, invdt = 0;
453 gmx_bool bConverged_all, bConverged = FALSE;
454 int niter = 0, g, c, ii, j, m, max_iter = 100;
456 dvec f; /* the pull force */
458 t_pull_group *pdyna, *pgrp0, *pgrp1;
461 snew(r_ij, pull->ncoord);
462 snew(dr_tot, pull->ncoord);
464 snew(rnew, pull->ngroup);
466 /* copy the current unconstrained positions for use in iterations. We
467 iterate until rinew[i] and rjnew[j] obey the constraints. Then
468 rinew - pull.x_unc[i] is the correction dr to group i */
469 for (g = 0; g < pull->ngroup; g++)
471 copy_dvec(pull->group[g].xp, rnew[g]);
475 /* There is only one pull coordinate and reference group */
476 copy_dvec(pull->dyna[0].xp, rnew[pull->coord[0].group[0]]);
479 /* Determine the constraint directions from the old positions */
480 for (c = 0; c < pull->ncoord; c++)
482 get_pull_coord_dr(pull, c, pbc, t, r_ij[c]);
483 /* Store the difference vector at time t for printing */
484 copy_dvec(r_ij[c], pull->coord[c].dr);
487 fprintf(debug, "Pull coord %d dr %f %f %f\n",
488 c, r_ij[c][XX], r_ij[c][YY], r_ij[c][ZZ]);
491 if (pull->eGeom == epullgDIR || pull->eGeom == epullgDIRPBC)
493 /* Select the component along vec */
495 for (m = 0; m < DIM; m++)
497 a += pull->coord[c].vec[m]*r_ij[c][m];
499 for (m = 0; m < DIM; m++)
501 r_ij[c][m] = a*pull->coord[c].vec[m];
505 if (dnorm2(r_ij[c]) == 0)
507 gmx_fatal(FARGS, "Distance for pull coordinate %d is zero with constraint pulling, which is not allowed.", c + 1);
511 bConverged_all = FALSE;
512 while (!bConverged_all && niter < max_iter)
514 bConverged_all = TRUE;
516 /* loop over all constraints */
517 for (c = 0; c < pull->ncoord; c++)
521 pcrd = &pull->coord[c];
522 pgrp0 = &pull->group[pcrd->group[0]];
523 pgrp1 = &pull->group[pcrd->group[1]];
525 /* Get the current difference vector */
526 low_get_pull_coord_dr(pull, pcrd, pbc, t,
527 rnew[pcrd->group[1]],
528 rnew[pcrd->group[0]],
531 ref = pcrd->init + pcrd->rate*t;
535 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
538 rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
545 gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref);
549 double q, c_a, c_b, c_c;
551 c_a = diprod(r_ij[c], r_ij[c]);
552 c_b = diprod(unc_ij, r_ij[c])*2;
553 c_c = diprod(unc_ij, unc_ij) - dsqr(ref);
557 q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
562 q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
569 "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
570 c_a, c_b, c_c, lambda);
574 /* The position corrections dr due to the constraints */
575 dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
576 dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
577 dr_tot[c] += -lambda*dnorm(r_ij[c]);
582 /* A 1-dimensional constraint along a vector */
584 for (m = 0; m < DIM; m++)
586 vec[m] = pcrd->vec[m];
587 a += unc_ij[m]*vec[m];
589 /* Select only the component along the vector */
590 dsvmul(a, vec, unc_ij);
594 fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
597 /* The position corrections dr due to the constraints */
598 dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
599 dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
600 dr_tot[c] += -lambda;
611 low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp);
612 low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3);
614 "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
615 rnew[g0][0], rnew[g0][1], rnew[g0][2],
616 rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
618 "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
619 "", "", "", "", "", "", ref);
621 "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
622 dr0[0], dr0[1], dr0[2],
623 dr1[0], dr1[1], dr1[2],
627 /* Update the COMs with dr */
628 dvec_inc(rnew[pcrd->group[1]], dr1);
629 dvec_inc(rnew[pcrd->group[0]], dr0);
632 /* Check if all constraints are fullfilled now */
633 for (c = 0; c < pull->ncoord; c++)
635 pcrd = &pull->coord[c];
636 ref = pcrd->init + pcrd->rate*t;
638 low_get_pull_coord_dr(pull, pcrd, pbc, t,
639 rnew[pcrd->group[1]],
640 rnew[pcrd->group[0]],
646 bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol;
651 for (m = 0; m < DIM; m++)
653 vec[m] = pcrd->vec[m];
655 inpr = diprod(unc_ij, vec);
656 dsvmul(inpr, vec, unc_ij);
658 fabs(diprod(unc_ij, vec) - ref) < pull->constr_tol;
666 fprintf(debug, "NOT CONVERGED YET: Group %d:"
667 "d_ref = %f, current d = %f\n",
668 g, ref, dnorm(unc_ij));
671 bConverged_all = FALSE;
676 /* if after all constraints are dealt with and bConverged is still TRUE
677 we're finished, if not we do another iteration */
679 if (niter > max_iter)
681 gmx_fatal(FARGS, "Too many iterations for constraint run: %d", niter);
684 /* DONE ITERATING, NOW UPDATE COORDINATES AND CALC. CONSTRAINT FORCES */
691 /* update atoms in the groups */
692 for (g = 0; g < pull->ngroup; g++)
694 const t_pull_group *pgrp;
697 if (PULL_CYL(pull) && g == pull->coord[0].group[0])
699 pgrp = &pull->dyna[0];
703 pgrp = &pull->group[g];
706 /* get the final constraint displacement dr for group g */
707 dvec_sub(rnew[g], pgrp->xp, dr);
708 /* select components of dr */
709 for (m = 0; m < DIM; m++)
711 dr[m] *= pull->dim[m];
714 /* update the atom positions */
716 for (j = 0; j < pgrp->nat_loc; j++)
718 ii = pgrp->ind_loc[j];
719 if (pgrp->weight_loc)
721 dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
723 for (m = 0; m < DIM; m++)
729 for (m = 0; m < DIM; m++)
731 v[ii][m] += invdt*tmp[m];
737 /* calculate the constraint forces, used for output and virial only */
738 for (c = 0; c < pull->ncoord; c++)
740 pcrd = &pull->coord[c];
741 pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
747 /* Add the pull contribution to the virial */
748 /* We have already checked above that r_ij[c] != 0 */
749 f_invr = pcrd->f_scal/dnorm(r_ij[c]);
751 for (j = 0; j < DIM; j++)
753 for (m = 0; m < DIM; m++)
755 vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
761 /* finished! I hope. Give back some memory */
767 /* Pulling with a harmonic umbrella potential or constant force */
768 static void do_pull_pot(int ePull,
769 t_pull *pull, t_pbc *pbc, double t, real lambda,
770 real *V, tensor vir, real *dVdl)
773 double dev, ndr, invdr;
777 /* loop over the pull coordinates */
780 for (c = 0; c < pull->ncoord; c++)
782 pcrd = &pull->coord[c];
784 get_pull_coord_distance(pull, c, pbc, t, pcrd->dr, &dev);
786 k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
787 dkdl = pcrd->kB - pcrd->k;
792 ndr = dnorm(pcrd->dr);
799 /* With an harmonic umbrella, the force is 0 at r=0,
800 * so we can set invdr to any value.
801 * With a constant force, the force at r=0 is not defined,
802 * so we zero it (this is anyhow a very rare event).
806 if (ePull == epullUMBRELLA)
808 pcrd->f_scal = -k*dev;
809 *V += 0.5* k*dsqr(dev);
810 *dVdl += 0.5*dkdl*dsqr(dev);
818 for (m = 0; m < DIM; m++)
820 pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
826 if (ePull == epullUMBRELLA)
828 pcrd->f_scal = -k*dev;
829 *V += 0.5* k*dsqr(dev);
830 *dVdl += 0.5*dkdl*dsqr(dev);
835 for (m = 0; m < DIM; m++)
837 ndr += pcrd->vec[m]*pcrd->dr[m];
843 for (m = 0; m < DIM; m++)
845 pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
852 /* Add the pull contribution to the virial */
853 for (j = 0; j < DIM; j++)
855 for (m = 0; m < DIM; m++)
857 vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
864 real pull_potential(int ePull, t_pull *pull, t_mdatoms *md, t_pbc *pbc,
865 t_commrec *cr, double t, real lambda,
866 rvec *x, rvec *f, tensor vir, real *dvdlambda)
870 pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
872 do_pull_pot(ePull, pull, pbc, t, lambda,
873 &V, pull->bVirial && MASTER(cr) ? vir : NULL, &dVdl);
875 /* Distribute forces over pulled groups */
876 apply_forces(pull, md, f);
883 return (MASTER(cr) ? V : 0.0);
886 void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
887 t_commrec *cr, double dt, double t,
888 rvec *x, rvec *xp, rvec *v, tensor vir)
890 pull_calc_coms(cr, pull, md, pbc, t, x, xp);
892 do_constraint(pull, pbc, xp, v, pull->bVirial && MASTER(cr), vir, dt, t);
895 static void make_local_pull_group(gmx_ga2la_t ga2la,
896 t_pull_group *pg, int start, int end)
901 for (i = 0; i < pg->nat; i++)
906 if (!ga2la_get_home(ga2la, ii, &ii))
911 if (ii >= start && ii < end)
913 /* This is a home atom, add it to the local pull group */
914 if (pg->nat_loc >= pg->nalloc_loc)
916 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
917 srenew(pg->ind_loc, pg->nalloc_loc);
918 if (pg->epgrppbc == epgrppbcCOS || pg->weight)
920 srenew(pg->weight_loc, pg->nalloc_loc);
923 pg->ind_loc[pg->nat_loc] = ii;
926 pg->weight_loc[pg->nat_loc] = pg->weight[i];
933 void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
947 for (g = 0; g < pull->ngroup; g++)
949 make_local_pull_group(ga2la, &pull->group[g],
953 /* Since the PBC of atoms might have changed, we need to update the PBC */
954 pull->bSetPBCatoms = TRUE;
957 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
958 int g, t_pull_group *pg, ivec pulldims,
959 gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
961 int i, ii, d, nfrozen, ndim;
963 double tmass, wmass, wwmass;
964 gmx_groups_t *groups;
965 gmx_mtop_atomlookup_t alook;
968 if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
970 /* There are no masses in the integrator.
971 * But we still want to have the correct mass-weighted COMs.
972 * So we store the real masses in the weights.
973 * We do not set nweight, so these weights do not end up in the tpx file.
975 if (pg->nweight == 0)
977 snew(pg->weight, pg->nat);
986 pg->weight_loc = NULL;
990 pg->nat_loc = pg->nat;
991 pg->ind_loc = pg->ind;
992 if (pg->epgrppbc == epgrppbcCOS)
994 snew(pg->weight_loc, pg->nat);
998 pg->weight_loc = pg->weight;
1002 groups = &mtop->groups;
1004 alook = gmx_mtop_atomlookup_init(mtop);
1010 for (i = 0; i < pg->nat; i++)
1013 gmx_mtop_atomnr_to_atom(alook, ii, &atom);
1014 if (ir->opts.nFreeze)
1016 for (d = 0; d < DIM; d++)
1018 if (pulldims[d] && ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
1024 if (ir->efep == efepNO)
1030 m = (1 - lambda)*atom->m + lambda*atom->mB;
1032 if (pg->nweight > 0)
1040 if (EI_ENERGY_MINIMIZATION(ir->eI))
1042 /* Move the mass to the weight */
1047 else if (ir->eI == eiBD)
1051 mbd = ir->bd_fric*ir->delta_t;
1055 if (groups->grpnr[egcTC] == NULL)
1057 mbd = ir->delta_t/ir->opts.tau_t[0];
1061 mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
1073 gmx_mtop_atomlookup_destroy(alook);
1077 gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
1078 pg->weight ? " weighted" : "", g);
1083 "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
1084 if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
1086 fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
1088 if (pg->epgrppbc == epgrppbcCOS)
1090 fprintf(fplog, ", cosine weighting will be used");
1092 fprintf(fplog, "\n");
1097 /* A value > 0 signals not frozen, it is updated later */
1103 for (d = 0; d < DIM; d++)
1105 ndim += pulldims[d]*pg->nat;
1107 if (fplog && nfrozen > 0 && nfrozen < ndim)
1110 "\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
1111 " that are subject to pulling are frozen.\n"
1112 " For pulling the whole group will be frozen.\n\n",
1120 void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
1121 gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
1122 gmx_bool bOutFile, unsigned long Flags)
1126 int c, g, start = 0, end = 0, m;
1130 pull->ePBC = ir->ePBC;
1133 case epbcNONE: pull->npbcdim = 0; break;
1134 case epbcXY: pull->npbcdim = 2; break;
1135 default: pull->npbcdim = 3; break;
1140 gmx_bool bAbs, bCos;
1143 for (c = 0; c < pull->ncoord; c++)
1145 if (pull->group[pull->coord[c].group[0]].nat == 0 ||
1146 pull->group[pull->coord[c].group[1]].nat == 0)
1152 fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
1153 EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom));
1154 fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
1155 pull->ncoord, pull->ncoord == 1 ? "" : "s",
1156 pull->ngroup, pull->ngroup == 1 ? "" : "s");
1159 fprintf(fplog, "with an absolute reference\n");
1162 for (g = 0; g < pull->ngroup; g++)
1164 if (pull->group[g].nat > 1 &&
1165 pull->group[g].pbcatom < 0)
1167 /* We are using cosine weighting */
1168 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
1174 please_cite(fplog, "Engin2010");
1178 /* We always add the virial contribution,
1179 * except for geometry = direction_periodic where this is impossible.
1181 pull->bVirial = (pull->eGeom != epullgDIRPBC);
1182 if (getenv("GMX_NO_PULLVIR") != NULL)
1186 fprintf(fplog, "Found env. var., will not add the virial contribution of the COM pull forces\n");
1188 pull->bVirial = FALSE;
1193 pull->dbuf_cyl = NULL;
1194 pull->bRefAt = FALSE;
1196 for (g = 0; g < pull->ngroup; g++)
1198 pgrp = &pull->group[g];
1199 pgrp->epgrppbc = epgrppbcNONE;
1202 /* Determine if we need to take PBC into account for calculating
1203 * the COM's of the pull groups.
1205 for (m = 0; m < pull->npbcdim; m++)
1207 if (pull->dim[m] && pgrp->nat > 1)
1209 if (pgrp->pbcatom >= 0)
1211 pgrp->epgrppbc = epgrppbcREFAT;
1212 pull->bRefAt = TRUE;
1218 gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
1220 pgrp->epgrppbc = epgrppbcCOS;
1221 if (pull->cosdim >= 0 && pull->cosdim != m)
1223 gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull_dim)");
1229 /* Set the indices */
1230 init_pull_group_index(fplog, cr, g, pgrp, pull->dim, mtop, ir, lambda);
1231 if (PULL_CYL(pull) && pgrp->invtm == 0)
1233 gmx_fatal(FARGS, "Can not have frozen atoms in a cylinder pull group");
1238 /* Absolute reference, set the inverse mass to zero */
1244 /* if we use dynamic reference groups, do some initialising for them */
1247 if (ir->ePull == epullCONSTRAINT && pull->ncoord > 1)
1249 /* We can't easily update the single reference group with multiple
1250 * constraints. This would require recalculating COMs.
1252 gmx_fatal(FARGS, "Constraint COM pulling supports only one coordinate with geometry=cylinder, you can use umbrella pulling with multiple coordinates");
1255 for (c = 0; c < pull->ncoord; c++)
1257 if (pull->group[pull->coord[c].group[0]].nat == 0)
1259 gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n");
1263 snew(pull->dyna, pull->ncoord);
1266 /* We still need to initialize the PBC reference coordinates */
1267 pull->bSetPBCatoms = TRUE;
1269 /* Only do I/O when we are doing dynamics and if we are the MASTER */
1274 if (pull->nstxout > 0)
1276 pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE, Flags);
1278 if (pull->nstfout > 0)
1280 pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
1286 void finish_pull(t_pull *pull)
1290 gmx_fio_fclose(pull->out_x);
1294 gmx_fio_fclose(pull->out_f);