From: Berk Hess Date: Mon, 13 Apr 2015 08:55:04 +0000 (+0200) Subject: Refactored pull data structures X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=66a0ae9d666811dada6df6f39e6f5c8d5f19f11f;p=alexxy%2Fgromacs.git Refactored pull data structures Moved the pull work variables into a new pull_t struct. Only the pull input parameters remain in (renamed) pull_params_t. For the moment struct pull_t pull_work is still present in t_inputrec. This should later be moved to a more appropriate place. Change-Id: I78f6f9736d52514f0327c86e43d4fb2ac88d60cd --- diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index 0b6001e667..c61e2893bc 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -9858,7 +9858,7 @@ void dd_partition_system(FILE *fplog, if (ir->bPull) { /* Update the local pull groups */ - dd_make_local_pull_groups(dd, ir->pull, mdatoms); + dd_make_local_pull_groups(dd, ir->pull_work, mdatoms); } if (ir->bRot) diff --git a/src/gromacs/fileio/tpxio.c b/src/gromacs/fileio/tpxio.c index 84efd2070b..cf77d95e6b 100644 --- a/src/gromacs/fileio/tpxio.c +++ b/src/gromacs/fileio/tpxio.c @@ -632,7 +632,7 @@ static void do_fepvals(t_fileio *fio, t_lambda *fepvals, gmx_bool bRead, int fil } } -static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead, +static void do_pull(t_fileio *fio, pull_params_t *pull, gmx_bool bRead, int file_version, int ePullOld) { int eGeomOld; diff --git a/src/gromacs/gmxlib/txtdump.c b/src/gromacs/gmxlib/txtdump.c index 0732b9dc89..4d04ff8759 100644 --- a/src/gromacs/gmxlib/txtdump.c +++ b/src/gromacs/gmxlib/txtdump.c @@ -760,7 +760,7 @@ static void pr_fepvals(FILE *fp, int indent, t_lambda *fep, gmx_bool bMDPformat) PS("dhdl-derivatives", DHDLDERIVATIVESTYPE(fep->dhdl_derivatives)); }; -static void pr_pull(FILE *fp, int indent, t_pull *pull) +static void pr_pull(FILE *fp, int indent, pull_params_t *pull) { int g; diff --git a/src/gromacs/gmxlib/typedefs.c b/src/gromacs/gmxlib/typedefs.c index ca19c6c39c..12c413d79e 100644 --- a/src/gromacs/gmxlib/typedefs.c +++ b/src/gromacs/gmxlib/typedefs.c @@ -80,21 +80,21 @@ static void done_pull_group(t_pull_group *pgrp) if (pgrp->nat > 0) { sfree(pgrp->ind); - sfree(pgrp->ind_loc); sfree(pgrp->weight); - sfree(pgrp->weight_loc); } } -static void done_pull(t_pull *pull) +static void done_pull_params(pull_params_t *pull) { int i; for (i = 0; i < pull->ngroup+1; i++) { done_pull_group(pull->group); - done_pull_group(pull->dyna); } + + sfree(pull->group); + sfree(pull->coord); } void done_inputrec(t_inputrec *ir) @@ -145,7 +145,7 @@ void done_inputrec(t_inputrec *ir) if (ir->pull) { - done_pull(ir->pull); + done_pull_params(ir->pull); sfree(ir->pull); } } diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c index 9e39ad6ede..3f3e233de5 100644 --- a/src/gromacs/gmxpreprocess/readir.c +++ b/src/gromacs/gmxpreprocess/readir.c @@ -2718,7 +2718,7 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) { t_grpopts *opts; gmx_groups_t *groups; - t_pull *pull; + pull_params_t *pull; int natoms, ai, aj, i, j, d, g, imin, jmin; t_iatom *ia; int *nrdf2, *na_vcm, na_tot; diff --git a/src/gromacs/gmxpreprocess/readir.h b/src/gromacs/gmxpreprocess/readir.h index 49b6de6f73..ea76d63004 100644 --- a/src/gromacs/gmxpreprocess/readir.h +++ b/src/gromacs/gmxpreprocess/readir.h @@ -127,16 +127,16 @@ void do_index(const char* mdparin, /* Routines In readpull.c */ char **read_pullparams(int *ninp_p, t_inpfile **inp, - t_pull *pull, + pull_params_t *pull, warninp_t wi); /* Reads the pull parameters, returns a list of the pull group names */ -void make_pull_groups(t_pull *pull, +void make_pull_groups(pull_params_t *pull, char **pgnames, const t_blocka *grps, char **gnames); /* Process the pull group parameters after reading the index groups */ -void make_pull_coords(t_pull *pull); +void make_pull_coords(pull_params_t *pull); /* Process the pull coordinates after reading the pull groups */ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda, diff --git a/src/gromacs/gmxpreprocess/readpull.c b/src/gromacs/gmxpreprocess/readpull.c index 68752443c2..49793283cc 100644 --- a/src/gromacs/gmxpreprocess/readpull.c +++ b/src/gromacs/gmxpreprocess/readpull.c @@ -186,7 +186,7 @@ static void init_pull_coord(t_pull_coord *pcrd, } char **read_pullparams(int *ninp_p, t_inpfile **inp_p, - t_pull *pull, + pull_params_t *pull, warninp_t wi) { int ninp, i, nchar, nscan, m, idum; @@ -306,7 +306,7 @@ char **read_pullparams(int *ninp_p, t_inpfile **inp_p, return grpbuf; } -void make_pull_groups(t_pull *pull, +void make_pull_groups(pull_params_t *pull, char **pgnames, const t_blocka *grps, char **gnames) { @@ -374,7 +374,7 @@ void make_pull_groups(t_pull *pull, } } -void make_pull_coords(t_pull *pull) +void make_pull_coords(pull_params_t *pull) { int c, d; t_pull_coord *pcrd; @@ -425,26 +425,27 @@ void make_pull_coords(t_pull *pull) void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda, const output_env_t oenv) { - t_mdatoms *md; - t_pull *pull; - t_pbc pbc; - int c; - double t_start; - - init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0); - md = init_mdatoms(NULL, mtop, ir->efep); + pull_params_t *pull; + struct pull_t *pull_work; + t_mdatoms *md; + t_pbc pbc; + int c; + double t_start; + + pull = ir->pull; + pull_work = init_pull(NULL, pull, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0); + md = init_mdatoms(NULL, mtop, ir->efep); atoms2md(mtop, ir, 0, NULL, mtop->natoms, md); if (ir->efep) { update_mdatoms(md, lambda); } - pull = ir->pull; set_pbc(&pbc, ir->ePBC, box); t_start = ir->init_t + ir->init_step*ir->delta_t; - pull_calc_coms(NULL, pull, md, &pbc, t_start, x, NULL); + pull_calc_coms(NULL, pull_work, md, &pbc, t_start, x, NULL); fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n"); for (c = 0; c < pull->ncoord; c++) @@ -469,7 +470,7 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l pcrd->init = 0; } - get_pull_coord_value(pull, c, &pbc, &value); + get_pull_coord_value(pull_work, c, &pbc, &value); fprintf(stderr, " %10.3f nm", value); if (pcrd->bStart) @@ -478,4 +479,6 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l } fprintf(stderr, " %10.3f nm\n", pcrd->init); } + + finish_pull(pull_work); } diff --git a/src/gromacs/legacyheaders/types/inputrec.h b/src/gromacs/legacyheaders/types/inputrec.h index 95cfc402cd..72d85c0119 100644 --- a/src/gromacs/legacyheaders/types/inputrec.h +++ b/src/gromacs/legacyheaders/types/inputrec.h @@ -98,31 +98,12 @@ typedef struct { gmx_bool *bTS; } t_grpopts; -enum { - epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS -}; - typedef struct { int nat; /* Number of atoms in the pull group */ atom_id *ind; /* The global atoms numbers */ - int nat_loc; /* Number of local pull atoms */ - int nalloc_loc; /* Allocation size for ind_loc and weight_loc */ - atom_id *ind_loc; /* Local pull indices */ int nweight; /* The number of weights (0 or nat) */ real *weight; /* Weights (use all 1 when weight==NULL) */ - real *weight_loc; /* Weights for the local indices */ - int epgrppbc; /* The type of pbc for this pull group, see enum above */ atom_id pbcatom; /* The reference atom for pbc (global number) */ - - /* Variables not present in mdp, but used at run time */ - gmx_bool bCalcCOM; /* Calculate COM? Not if only used as cylinder group */ - real mwscale; /* mass*weight scaling factor 1/sum w m */ - real wscale; /* scaling factor for the weights: sum w m/sum w w m */ - real invtm; /* inverse total mass of the group: 1/wscale sum w m */ - dvec *mdw; /* mass*gradient(weight) for atoms */ - double *dv; /* distance to the other group along vec */ - dvec x; /* center of mass before update */ - dvec xp; /* center of mass after update before constraining */ } t_pull_group; typedef struct { @@ -137,16 +118,6 @@ typedef struct { real rate; /* Rate of motion (nm/ps) */ real k; /* force constant */ real kB; /* force constant for state B */ - - /* Variables not present in mdp, but used at run time */ - double value_ref; /* The reference value, usually init+rate*t */ - double value; /* The current value of the coordinate */ - dvec dr; /* The distance from the reference group */ - double vec_len; /* Length of vec for direction-relative */ - dvec ffrad; /* conversion factor from vec to radial force */ - double cyl_dev; /* The deviation from the reference position */ - double f_scal; /* Scalar force for directional pulling */ - dvec f; /* force due to the pulling/constraining */ } t_pull_coord; typedef struct { @@ -229,27 +200,13 @@ typedef struct { gmx_bool bPrintComp; /* Print cartesian components for each coord with geometry=distance */ int nstxout; /* Output frequency for pull x */ int nstfout; /* Output frequency for pull f */ - int ePBC; /* the boundary conditions */ - int npbcdim; /* do pbc in dims 0 <= dim < npbcdim */ - gmx_bool bRefAt; /* do we need reference atoms for a group COM ? */ - int cosdim; /* dimension for cosine weighting, -1 if none */ + t_pull_group *group; /* groups to pull/restrain/etc/ */ t_pull_coord *coord; /* the pull coordinates */ +} pull_params_t; - /* Variables not present in mdp, but used at run time */ - gmx_bool bPotential; /* Are there coordinates with potential? */ - gmx_bool bConstraint; /* Are there constrained coordinates? */ - gmx_bool bCylinder; /* Is group 0 a cylinder group? */ - t_pull_group *dyna; /* dynamic groups for use with local constraints */ - gmx_bool bSetPBCatoms; /* Do we need to set x_pbc for the groups? */ - - rvec *rbuf; /* COM calculation buffer */ - dvec *dbuf; /* COM calculation buffer */ - double *dbuf_cyl; /* cylinder ref. groups COM calculation buffer */ - - FILE *out_x; /* output file for pull data */ - FILE *out_f; /* output file for pull data */ -} t_pull; +/* Abstract type for COM pull caclulations only defined in the pull module */ +struct pull_t; /* Abstract types for enforced rotation only defined in pull_rotation.c */ @@ -460,7 +417,9 @@ typedef struct { real wall_density[2]; /* Number density for walls */ real wall_ewald_zfac; /* Scaling factor for the box for Ewald */ gmx_bool bPull; /* Do we do COM pulling? */ - t_pull *pull; /* The data for center of mass pulling */ + pull_params_t *pull; /* The data for center of mass pulling */ + struct pull_t *pull_work; /* The COM pull force calculation data structure; TODO this pointer should live somewhere else */ + gmx_bool bRot; /* Calculate enforced rotation potential(s)? */ t_rot *rot; /* The data for enforced rotation potentials */ int eSwapCoords; /* Do ion/water position exchanges (CompEL)? */ diff --git a/src/gromacs/mdlib/broadcaststructs.cpp b/src/gromacs/mdlib/broadcaststructs.cpp index c13b7185fe..049fc4cba7 100644 --- a/src/gromacs/mdlib/broadcaststructs.cpp +++ b/src/gromacs/mdlib/broadcaststructs.cpp @@ -493,7 +493,7 @@ static void bc_pull_group(const t_commrec *cr, t_pull_group *pgrp) } } -static void bc_pull(const t_commrec *cr, t_pull *pull) +static void bc_pull(const t_commrec *cr, pull_params_t *pull) { int g; diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index ef655c062d..e4d2f9d3c5 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -614,7 +614,7 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner, if (econq == econqCoord) { - if (ir->bPull && ir->pull->bConstraint) + if (ir->bPull && pull_have_constraint(ir->pull_work)) { if (EI_DYNAMICS(ir->eI)) { @@ -625,7 +625,7 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner, t = ir->init_t; } set_pbc(&pbc, ir->ePBC, box); - pull_constraint(ir->pull, md, &pbc, cr, ir->delta_t, t, x, xprime, v, *vir); + pull_constraint(ir->pull_work, md, &pbc, cr, ir->delta_t, t, x, xprime, v, *vir); } if (constr->ed && delta_step > 0) { @@ -1172,7 +1172,7 @@ gmx_constr_t init_constraints(FILE *fplog, nset = gmx_mtop_ftype_count(mtop, F_SETTLE); if (ncon+nset == 0 && - !(ir->bPull && ir->pull->bConstraint) && + !(ir->bPull && pull_have_constraint(ir->pull_work)) && ed == NULL) { return NULL; diff --git a/src/gromacs/mdlib/mdebin.c b/src/gromacs/mdlib/mdebin.c index d46d50f637..cef361bfd6 100644 --- a/src/gromacs/mdlib/mdebin.c +++ b/src/gromacs/mdlib/mdebin.c @@ -57,6 +57,7 @@ #include "gromacs/math/vec.h" #include "gromacs/mdlib/mdebin_bar.h" #include "gromacs/pbcutil/pbc.h" +#include "gromacs/pulling/pull.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/utility/smalloc.h" @@ -291,7 +292,7 @@ t_mdebin *init_mdebin(ener_file_t fp_ene, } else if (i == F_COM_PULL) { - md->bEner[i] = (ir->bPull && ir->pull->bPotential); + md->bEner[i] = (ir->bPull && pull_have_potential(ir->pull_work)); } else if (i == F_ECONSERVED) { diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index fd4074dde3..08bc5071d3 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -332,7 +332,7 @@ static void pull_potential_wrapper(t_commrec *cr, set_pbc(&pbc, ir->ePBC, box); dvdl = 0; enerd->term[F_COM_PULL] += - pull_potential(ir->pull, mdatoms, &pbc, + pull_potential(ir->pull_work, mdatoms, &pbc, cr, t, lambda[efptRESTRAINT], x, f, vir_force, &dvdl); enerd->dvdl_lin[efptRESTRAINT] += dvdl; wallcycle_stop(wcycle, ewcPULLPOT); @@ -1155,9 +1155,9 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr, clear_rvec(fr->vir_diag_posres); } - if (inputrec->bPull && inputrec->pull->bConstraint) + if (inputrec->bPull && pull_have_constraint(inputrec->pull_work)) { - clear_pull_forces(inputrec->pull); + clear_pull_forces(inputrec->pull_work); } /* We calculate the non-bonded forces, when done on the CPU, here. @@ -1455,7 +1455,7 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr, } } - if (inputrec->bPull && inputrec->pull->bPotential) + if (inputrec->bPull && pull_have_potential(inputrec->pull_work)) { /* Since the COM pulling is always done mass-weighted, no forces are * applied to vsites and this call can be done after vsite spreading. @@ -1798,9 +1798,9 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr, clear_rvec(fr->vir_diag_posres); } - if (inputrec->bPull && inputrec->pull->bConstraint) + if (inputrec->bPull && pull_have_constraint(inputrec->pull_work)) { - clear_pull_forces(inputrec->pull); + clear_pull_forces(inputrec->pull_work); } /* update QMMMrec, if necessary */ @@ -1919,7 +1919,7 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr, } } - if (inputrec->bPull && inputrec->pull->bPotential) + if (inputrec->bPull && pull_have_potential(inputrec->pull_work)) { pull_potential_wrapper(cr, inputrec, box, x, f, vir_force, mdatoms, enerd, lambda, t, diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp index 4635654dc8..3c970da0e4 100644 --- a/src/gromacs/pulling/pull.cpp +++ b/src/gromacs/pulling/pull.cpp @@ -59,12 +59,14 @@ #include "gromacs/legacyheaders/types/commrec.h" #include "gromacs/math/vec.h" #include "gromacs/pbcutil/pbc.h" +#include "gromacs/pulling/pull_internal.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/utility/futil.h" #include "gromacs/utility/gmxassert.h" #include "gromacs/utility/smalloc.h" -static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp) +static void pull_print_group_x(FILE *out, ivec dim, + const pull_group_work_t *pgrp) { int m; @@ -77,7 +79,7 @@ static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp) } } -static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd, +static void pull_print_coord_dr(FILE *out, const pull_coord_work_t *pcrd, gmx_bool bPrintRefValue, gmx_bool bPrintComponents) { @@ -94,7 +96,7 @@ static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd, for (m = 0; m < DIM; m++) { - if (pcrd->dim[m]) + if (pcrd->params.dim[m]) { fprintf(out, "\t%g", pcrd->dr[m]); } @@ -102,39 +104,43 @@ static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd, } } -static void pull_print_x(FILE *out, t_pull *pull, double t) +static void pull_print_x(FILE *out, struct pull_t *pull, double t) { - int c; - t_pull_coord *pcrd; + int c; fprintf(out, "%.4f", t); for (c = 0; c < pull->ncoord; c++) { + pull_coord_work_t *pcrd; + pcrd = &pull->coord[c]; - if (pull->bPrintCOM1) + if (pull->params.bPrintCOM1) { - if (pcrd->eGeom == epullgCYL) + if (pcrd->params.eGeom == epullgCYL) { - pull_print_group_x(out, pcrd->dim, &pull->dyna[c]); + pull_print_group_x(out, pcrd->params.dim, &pull->dyna[c]); } else { - pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[0]]); + pull_print_group_x(out, pcrd->params.dim, + &pull->group[pcrd->params.group[0]]); } } - if (pull->bPrintCOM2) + if (pull->params.bPrintCOM2) { - pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[1]]); + pull_print_group_x(out, pcrd->params.dim, + &pull->group[pcrd->params.group[1]]); } pull_print_coord_dr(out, pcrd, - pull->bPrintRefValue, pull->bPrintComp); + pull->params.bPrintRefValue, + pull->params.bPrintComp); } fprintf(out, "\n"); } -static void pull_print_f(FILE *out, t_pull *pull, double t) +static void pull_print_f(FILE *out, struct pull_t *pull, double t) { int c; @@ -147,20 +153,20 @@ static void pull_print_f(FILE *out, t_pull *pull, double t) fprintf(out, "\n"); } -void pull_print_output(t_pull *pull, gmx_int64_t step, double time) +void pull_print_output(struct pull_t *pull, gmx_int64_t step, double time) { - if ((pull->nstxout != 0) && (step % pull->nstxout == 0)) + if ((pull->params.nstxout != 0) && (step % pull->params.nstxout == 0)) { pull_print_x(pull->out_x, pull, time); } - if ((pull->nstfout != 0) && (step % pull->nstfout == 0)) + if ((pull->params.nstfout != 0) && (step % pull->params.nstfout == 0)) { pull_print_f(pull->out_f, pull, time); } } -static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv, +static FILE *open_pull_out(const char *fn, struct pull_t *pull, const output_env_t oenv, gmx_bool bCoord, unsigned long Flags) { FILE *fp; @@ -199,12 +205,12 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv * the data in print_pull_x above. */ - if (pull->bPrintCOM1) + if (pull->params.bPrintCOM1) { /* Legend for reference group position */ for (m = 0; m < DIM; m++) { - if (pull->coord[c].dim[m]) + if (pull->coord[c].params.dim[m]) { sprintf(buf, "%d %s%d%c", c+1, "c", 1, 'X'+m); setname[nsets] = gmx_strdup(buf); @@ -212,12 +218,12 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv } } } - if (pull->bPrintCOM2) + if (pull->params.bPrintCOM2) { /* Legend for reference group position */ for (m = 0; m < DIM; m++) { - if (pull->coord[c].dim[m]) + if (pull->coord[c].params.dim[m]) { sprintf(buf, "%d %s%d%c", c+1, "c", 2, 'X'+m); setname[nsets] = gmx_strdup(buf); @@ -229,18 +235,18 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv sprintf(buf, "%d", c+1); setname[nsets] = gmx_strdup(buf); nsets++; - if (pull->bPrintRefValue) + if (pull->params.bPrintRefValue) { sprintf(buf, "%c%d", 'r', c+1); setname[nsets] = gmx_strdup(buf); nsets++; } - if (pull->bPrintComp) + if (pull->params.bPrintComp) { /* The pull coord distance components */ for (m = 0; m < DIM; m++) { - if (pull->coord[c].dim[m]) + if (pull->coord[c].params.dim[m]) { sprintf(buf, "%d %s%c", c+1, "d", 'X'+m); setname[nsets] = gmx_strdup(buf); @@ -272,7 +278,8 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv } /* Apply forces in a mass weighted fashion */ -static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md, +static void apply_forces_grp(const pull_group_work_t *pgrp, + const t_mdatoms *md, const dvec f_pull, int sign, rvec *f) { int i, ii, m; @@ -280,7 +287,7 @@ static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md, inv_wm = pgrp->mwscale; - if (pgrp->nat == 1 && pgrp->nat_loc == 1) + if (pgrp->params.nat == 1 && pgrp->nat_loc == 1) { /* Only one atom and our rank has this atom: we can skip * the mass weighting, which means that this code also works @@ -311,7 +318,8 @@ static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md, } /* Apply forces in a mass weighted fashion to a cylinder group */ -static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr, +static void apply_forces_cyl_grp(const pull_group_work_t *pgrp, + double dv_corr, const t_mdatoms *md, const dvec f_pull, double f_scal, int sign, rvec *f) @@ -347,10 +355,10 @@ static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr, /* Apply torque forces in a mass weighted fashion to the groups that define * the pull vector direction for pull coordinate pcrd. */ -static void apply_forces_vec_torque(const t_pull *pull, - const t_pull_coord *pcrd, - const t_mdatoms *md, - rvec *f) +static void apply_forces_vec_torque(const struct pull_t *pull, + const pull_coord_work_t *pcrd, + const t_mdatoms *md, + rvec *f) { double inpr; int m; @@ -376,21 +384,22 @@ static void apply_forces_vec_torque(const t_pull *pull, } /* Apply the force to the groups defining the vector using opposite signs */ - apply_forces_grp(&pull->group[pcrd->group[2]], md, f_perp, -1, f); - apply_forces_grp(&pull->group[pcrd->group[3]], md, f_perp, 1, f); + apply_forces_grp(&pull->group[pcrd->params.group[2]], md, f_perp, -1, f); + apply_forces_grp(&pull->group[pcrd->params.group[3]], md, f_perp, 1, f); } /* Apply forces in a mass weighted fashion */ -static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f) +static void apply_forces(struct pull_t * pull, t_mdatoms * md, rvec *f) { - int c; - const t_pull_coord *pcrd; + int c; for (c = 0; c < pull->ncoord; c++) { + const pull_coord_work_t *pcrd; + pcrd = &pull->coord[c]; - if (pcrd->eGeom == epullgCYL) + if (pcrd->params.eGeom == epullgCYL) { dvec f_tot; int m; @@ -403,11 +412,11 @@ static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f) { f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m]; } - apply_forces_grp(&pull->group[pcrd->group[1]], md, f_tot, 1, f); + apply_forces_grp(&pull->group[pcrd->params.group[1]], md, f_tot, 1, f); } else { - if (pcrd->eGeom == epullgDIRRELATIVE) + if (pcrd->params.eGeom == epullgDIRRELATIVE) { /* We need to apply the torque forces to the pull groups * that define the pull vector. @@ -415,16 +424,17 @@ static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f) apply_forces_vec_torque(pull, pcrd, md, f); } - if (pull->group[pcrd->group[0]].nat > 0) + if (pull->group[pcrd->params.group[0]].params.nat > 0) { - apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f); + apply_forces_grp(&pull->group[pcrd->params.group[0]], md, pcrd->f, -1, f); } - apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f); + apply_forces_grp(&pull->group[pcrd->params.group[1]], md, pcrd->f, 1, f); } } } -static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc) +static double max_pull_distance2(const pull_coord_work_t *pcrd, + const t_pbc *pbc) { double max_d2; int m; @@ -433,7 +443,7 @@ static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc) for (m = 0; m < pbc->ndim_ePBC; m++) { - if (pcrd->dim[m] != 0) + if (pcrd->params.dim[m] != 0) { max_d2 = std::min(max_d2, static_cast(norm2(pbc->box[m]))); } @@ -445,31 +455,31 @@ static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc) /* This function returns the distance based on coordinates xg and xref. * Note that the pull coordinate struct pcrd is not modified. */ -static void low_get_pull_coord_dr(const t_pull *pull, - const t_pull_coord *pcrd, +static void low_get_pull_coord_dr(const struct pull_t *pull, + const pull_coord_work_t *pcrd, const t_pbc *pbc, dvec xg, dvec xref, double max_dist2, dvec dr) { - const t_pull_group *pgrp0; - int m; - dvec xrefr, dref = {0, 0, 0}; - double dr2; + const pull_group_work_t *pgrp0; + int m; + dvec xrefr, dref = {0, 0, 0}; + double dr2; - pgrp0 = &pull->group[pcrd->group[0]]; + pgrp0 = &pull->group[pcrd->params.group[0]]; /* Only the first group can be an absolute reference, in that case nat=0 */ - if (pgrp0->nat == 0) + if (pgrp0->params.nat == 0) { for (m = 0; m < DIM; m++) { - xref[m] = pcrd->origin[m]; + xref[m] = pcrd->params.origin[m]; } } copy_dvec(xref, xrefr); - if (pcrd->eGeom == epullgDIRPBC) + if (pcrd->params.eGeom == epullgDIRPBC) { for (m = 0; m < DIM; m++) { @@ -483,16 +493,17 @@ static void low_get_pull_coord_dr(const t_pull *pull, dr2 = 0; for (m = 0; m < DIM; m++) { - dr[m] *= pcrd->dim[m]; + dr[m] *= pcrd->params.dim[m]; dr2 += dr[m]*dr[m]; } if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2) { 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", - pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2)); + pcrd->params.group[0], pcrd->params.group[1], + sqrt(dr2), sqrt(max_dist2)); } - if (pcrd->eGeom == epullgDIRPBC) + if (pcrd->params.eGeom == epullgDIRPBC) { dvec_inc(dr, dref); } @@ -501,16 +512,16 @@ static void low_get_pull_coord_dr(const t_pull *pull, /* This function returns the distance based on the contents of the pull struct. * pull->coord[coord_ind].dr, and potentially vector, are updated. */ -static void get_pull_coord_dr(t_pull *pull, - int coord_ind, - const t_pbc *pbc) +static void get_pull_coord_dr(struct pull_t *pull, + int coord_ind, + const t_pbc *pbc) { - double md2; - t_pull_coord *pcrd; + double md2; + pull_coord_work_t *pcrd; pcrd = &pull->coord[coord_ind]; - if (pcrd->eGeom == epullgDIRPBC) + if (pcrd->params.eGeom == epullgDIRPBC) { md2 = -1; } @@ -519,21 +530,21 @@ static void get_pull_coord_dr(t_pull *pull, md2 = max_pull_distance2(pcrd, pbc); } - if (pcrd->eGeom == epullgDIRRELATIVE) + if (pcrd->params.eGeom == epullgDIRRELATIVE) { /* We need to determine the pull vector */ - const t_pull_group *pgrp2, *pgrp3; - dvec vec; - int m; + const pull_group_work_t *pgrp2, *pgrp3; + dvec vec; + int m; - pgrp2 = &pull->group[pcrd->group[2]]; - pgrp3 = &pull->group[pcrd->group[3]]; + pgrp2 = &pull->group[pcrd->params.group[2]]; + pgrp3 = &pull->group[pcrd->params.group[3]]; pbc_dx_d(pbc, pgrp3->x, pgrp2->x, vec); for (m = 0; m < DIM; m++) { - vec[m] *= pcrd->dim[m]; + vec[m] *= pcrd->params.dim[m]; } pcrd->vec_len = dnorm(vec); for (m = 0; m < DIM; m++) @@ -550,36 +561,36 @@ static void get_pull_coord_dr(t_pull *pull, } low_get_pull_coord_dr(pull, pcrd, pbc, - pull->group[pcrd->group[1]].x, - pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x, + pull->group[pcrd->params.group[1]].x, + pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->params.group[0]].x, md2, pcrd->dr); } -static void update_pull_coord_reference_value(t_pull_coord *pcrd, double t) +static void update_pull_coord_reference_value(pull_coord_work_t *pcrd, double t) { /* With zero rate the reference value is set initially and doesn't change */ - if (pcrd->rate != 0) + if (pcrd->params.rate != 0) { - pcrd->value_ref = pcrd->init + pcrd->rate*t; + pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t; } } /* Calculates pull->coord[coord_ind].value. * This function also updates pull->coord[coord_ind].dr. */ -static void get_pull_coord_distance(t_pull *pull, - int coord_ind, - const t_pbc *pbc) +static void get_pull_coord_distance(struct pull_t *pull, + int coord_ind, + const t_pbc *pbc) { - t_pull_coord *pcrd; - int m; + pull_coord_work_t *pcrd; + int m; pcrd = &pull->coord[coord_ind]; get_pull_coord_dr(pull, coord_ind, pbc); - switch (pcrd->eGeom) + switch (pcrd->params.eGeom) { case epullgDIST: /* Pull along the vector between the com's */ @@ -602,15 +613,15 @@ static void get_pull_coord_distance(t_pull *pull, /* Returns the deviation from the reference value. * Updates pull->coord[coord_ind].dr, .value and .value_ref. */ -static double get_pull_coord_deviation(t_pull *pull, - int coord_ind, - const t_pbc *pbc, - double t) +static double get_pull_coord_deviation(struct pull_t *pull, + int coord_ind, + const t_pbc *pbc, + double t) { - static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety, - but is fairly benign */ - t_pull_coord *pcrd; - double dev = 0; + static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety, + but is fairly benign */ + pull_coord_work_t *pcrd; + double dev = 0; pcrd = &pull->coord[coord_ind]; @@ -618,7 +629,7 @@ static double get_pull_coord_deviation(t_pull *pull, update_pull_coord_reference_value(pcrd, t); - switch (pcrd->eGeom) + switch (pcrd->params.eGeom) { case epullgDIST: /* Pull along the vector between the com's */ @@ -652,17 +663,17 @@ static double get_pull_coord_deviation(t_pull *pull, return dev; } -void get_pull_coord_value(t_pull *pull, - int coord_ind, - const t_pbc *pbc, - double *value) +void get_pull_coord_value(struct pull_t *pull, + int coord_ind, + const t_pbc *pbc, + double *value) { get_pull_coord_distance(pull, coord_ind, pbc); *value = pull->coord[coord_ind].value; } -void clear_pull_forces(t_pull *pull) +void clear_pull_forces(struct pull_t *pull) { int i; @@ -678,7 +689,7 @@ void clear_pull_forces(t_pull *pull) } /* Apply constraint using SHAKE */ -static void do_constraint(t_pull *pull, t_pbc *pbc, +static void do_constraint(struct pull_t *pull, t_pbc *pbc, rvec *x, rvec *v, gmx_bool bMaster, tensor vir, double dt, double t) @@ -695,8 +706,6 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, int niter = 0, g, c, ii, j, m, max_iter = 100; double a; dvec tmp, tmp3; - t_pull_group *pgrp0, *pgrp1; - t_pull_coord *pcrd; snew(r_ij, pull->ncoord); snew(dr_tot, pull->ncoord); @@ -714,9 +723,11 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, /* Determine the constraint directions from the old positions */ for (c = 0; c < pull->ncoord; c++) { + pull_coord_work_t *pcrd; + pcrd = &pull->coord[c]; - if (pcrd->eType != epullCONSTRAINT) + if (pcrd->params.eType != epullCONSTRAINT) { continue; } @@ -733,7 +744,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, c, pcrd->dr[XX], pcrd->dr[YY], pcrd->dr[ZZ]); } - if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC) + if (pcrd->params.eGeom == epullgDIR || + pcrd->params.eGeom == epullgDIRPBC) { /* Select the component along vec */ a = 0; @@ -765,24 +777,26 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, /* loop over all constraints */ for (c = 0; c < pull->ncoord; c++) { - dvec dr0, dr1; + pull_coord_work_t *pcrd; + pull_group_work_t *pgrp0, *pgrp1; + dvec dr0, dr1; pcrd = &pull->coord[c]; - if (pcrd->eType != epullCONSTRAINT) + if (pcrd->params.eType != epullCONSTRAINT) { continue; } update_pull_coord_reference_value(pcrd, t); - pgrp0 = &pull->group[pcrd->group[0]]; - pgrp1 = &pull->group[pcrd->group[1]]; + pgrp0 = &pull->group[pcrd->params.group[0]]; + pgrp1 = &pull->group[pcrd->params.group[1]]; /* Get the current difference vector */ low_get_pull_coord_dr(pull, pcrd, pbc, - rnew[pcrd->group[1]], - rnew[pcrd->group[0]], + rnew[pcrd->params.group[1]], + rnew[pcrd->params.group[0]], -1, unc_ij); if (debug) @@ -792,7 +806,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, rm = 1.0/(pgrp0->invtm + pgrp1->invtm); - switch (pcrd->eGeom) + switch (pcrd->params.eGeom) { case epullgDIST: if (pcrd->value_ref <= 0) @@ -866,8 +880,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, { int g0, g1; - g0 = pcrd->group[0]; - g1 = pcrd->group[1]; + g0 = pcrd->params.group[0]; + g1 = pcrd->params.group[1]; low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp); low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3); fprintf(debug, @@ -885,30 +899,32 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, } /* END DEBUG */ /* Update the COMs with dr */ - dvec_inc(rnew[pcrd->group[1]], dr1); - dvec_inc(rnew[pcrd->group[0]], dr0); + dvec_inc(rnew[pcrd->params.group[1]], dr1); + dvec_inc(rnew[pcrd->params.group[0]], dr0); } /* Check if all constraints are fullfilled now */ for (c = 0; c < pull->ncoord; c++) { + pull_coord_work_t *pcrd; + pcrd = &pull->coord[c]; - if (pcrd->eType != epullCONSTRAINT) + if (pcrd->params.eType != epullCONSTRAINT) { continue; } low_get_pull_coord_dr(pull, pcrd, pbc, - rnew[pcrd->group[1]], - rnew[pcrd->group[0]], + rnew[pcrd->params.group[1]], + rnew[pcrd->params.group[0]], -1, unc_ij); - switch (pcrd->eGeom) + switch (pcrd->params.eGeom) { case epullgDIST: bConverged = - fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->constr_tol; + fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->params.constr_tol; break; case epullgDIR: case epullgDIRPBC: @@ -920,7 +936,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, inpr = diprod(unc_ij, vec); dsvmul(inpr, vec, unc_ij); bConverged = - fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->constr_tol; + fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->params.constr_tol; break; } @@ -956,8 +972,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, /* update atoms in the groups */ for (g = 0; g < pull->ngroup; g++) { - const t_pull_group *pgrp; - dvec dr; + const pull_group_work_t *pgrp; + dvec dr; pgrp = &pull->group[g]; @@ -996,16 +1012,18 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, /* calculate the constraint forces, used for output and virial only */ for (c = 0; c < pull->ncoord; c++) { + pull_coord_work_t *pcrd; + pcrd = &pull->coord[c]; - if (pcrd->eType != epullCONSTRAINT) + if (pcrd->params.eType != epullCONSTRAINT) { continue; } - pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt); + pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)*dt*dt); - if (vir != NULL && pcrd->eGeom != epullgDIRPBC && bMaster) + if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC && bMaster) { double f_invr; @@ -1030,32 +1048,33 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, } /* Pulling with a harmonic umbrella potential or constant force */ -static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda, +static void do_pull_pot(struct pull_t *pull, t_pbc *pbc, double t, real lambda, real *V, tensor vir, real *dVdl) { - int c, j, m; - double dev, ndr, invdr = 0; - real k, dkdl; - t_pull_coord *pcrd; + int c, j, m; + double dev, ndr, invdr = 0; + real k, dkdl; /* loop over the pull coordinates */ *V = 0; *dVdl = 0; for (c = 0; c < pull->ncoord; c++) { + pull_coord_work_t *pcrd; + pcrd = &pull->coord[c]; - if (pcrd->eType == epullCONSTRAINT) + if (pcrd->params.eType == epullCONSTRAINT) { continue; } dev = get_pull_coord_deviation(pull, c, pbc, t); - k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB; - dkdl = pcrd->kB - pcrd->k; + k = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB; + dkdl = pcrd->params.kB - pcrd->params.k; - if (pcrd->eGeom == epullgDIST) + if (pcrd->params.eGeom == epullgDIST) { ndr = dnorm(pcrd->dr); if (ndr > 0) @@ -1081,14 +1100,14 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda, } } - switch (pcrd->eType) + switch (pcrd->params.eType) { case epullUMBRELLA: case epullFLATBOTTOM: /* The only difference between an umbrella and a flat-bottom * potential is that a flat-bottom is zero below zero. */ - if (pcrd->eType == epullFLATBOTTOM && dev < 0) + if (pcrd->params.eType == epullFLATBOTTOM && dev < 0) { dev = 0; } @@ -1106,7 +1125,7 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda, gmx_incons("Unsupported pull type in do_pull_pot"); } - if (pcrd->eGeom == epullgDIST) + if (pcrd->params.eGeom == epullgDIST) { for (m = 0; m < DIM; m++) { @@ -1121,7 +1140,7 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda, } } - if (vir != NULL && pcrd->eGeom != epullgDIRPBC) + if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC) { /* Add the pull contribution to the virial */ for (j = 0; j < DIM; j++) @@ -1135,12 +1154,14 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda, } } -real pull_potential(t_pull *pull, t_mdatoms *md, t_pbc *pbc, +real pull_potential(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc, t_commrec *cr, double t, real lambda, rvec *x, rvec *f, tensor vir, real *dvdlambda) { real V, dVdl; + assert(pull != NULL); + pull_calc_coms(cr, pull, md, pbc, t, x, NULL); do_pull_pot(pull, pbc, t, lambda, @@ -1157,24 +1178,26 @@ real pull_potential(t_pull *pull, t_mdatoms *md, t_pbc *pbc, return (MASTER(cr) ? V : 0.0); } -void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc, +void pull_constraint(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc, t_commrec *cr, double dt, double t, rvec *x, rvec *xp, rvec *v, tensor vir) { + assert(pull != NULL); + pull_calc_coms(cr, pull, md, pbc, t, x, xp); do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t); } static void make_local_pull_group(gmx_ga2la_t ga2la, - t_pull_group *pg, int start, int end) + pull_group_work_t *pg, int start, int end) { int i, ii; pg->nat_loc = 0; - for (i = 0; i < pg->nat; i++) + for (i = 0; i < pg->params.nat; i++) { - ii = pg->ind[i]; + ii = pg->params.ind[i]; if (ga2la) { if (!ga2la_get_home(ga2la, ii, &ii)) @@ -1189,22 +1212,22 @@ static void make_local_pull_group(gmx_ga2la_t ga2la, { pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1); srenew(pg->ind_loc, pg->nalloc_loc); - if (pg->epgrppbc == epgrppbcCOS || pg->weight) + if (pg->epgrppbc == epgrppbcCOS || pg->params.weight != NULL) { srenew(pg->weight_loc, pg->nalloc_loc); } } pg->ind_loc[pg->nat_loc] = ii; - if (pg->weight) + if (pg->params.weight != NULL) { - pg->weight_loc[pg->nat_loc] = pg->weight[i]; + pg->weight_loc[pg->nat_loc] = pg->params.weight[i]; } pg->nat_loc++; } } } -void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md) +void dd_make_local_pull_groups(gmx_domdec_t *dd, struct pull_t *pull, t_mdatoms *md) { gmx_ga2la_t ga2la; int g; @@ -1229,14 +1252,15 @@ void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md) } static void init_pull_group_index(FILE *fplog, t_commrec *cr, - int g, t_pull_group *pg, + int g, pull_group_work_t *pg, gmx_bool bConstraint, ivec pulldim_con, - gmx_mtop_t *mtop, t_inputrec *ir, real lambda) + const gmx_mtop_t *mtop, + const t_inputrec *ir, real lambda) { int i, ii, d, nfrozen, ndim; real m, w, mbd; double tmass, wmass, wwmass; - gmx_groups_t *groups; + const gmx_groups_t *groups; gmx_mtop_atomlookup_t alook; t_atom *atom; @@ -1247,9 +1271,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, * So we store the real masses in the weights. * We do not set nweight, so these weights do not end up in the tpx file. */ - if (pg->nweight == 0) + if (pg->params.nweight == 0) { - snew(pg->weight, pg->nat); + snew(pg->params.weight, pg->params.nat); } } @@ -1262,15 +1286,15 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, } else { - pg->nat_loc = pg->nat; - pg->ind_loc = pg->ind; + pg->nat_loc = pg->params.nat; + pg->ind_loc = pg->params.ind; if (pg->epgrppbc == epgrppbcCOS) { - snew(pg->weight_loc, pg->nat); + snew(pg->weight_loc, pg->params.nat); } else { - pg->weight_loc = pg->weight; + pg->weight_loc = pg->params.weight; } } @@ -1282,9 +1306,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, tmass = 0; wmass = 0; wwmass = 0; - for (i = 0; i < pg->nat; i++) + for (i = 0; i < pg->params.nat; i++) { - ii = pg->ind[i]; + ii = pg->params.ind[i]; gmx_mtop_atomnr_to_atom(alook, ii, &atom); if (bConstraint && ir->opts.nFreeze) { @@ -1305,9 +1329,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, { m = (1 - lambda)*atom->m + lambda*atom->mB; } - if (pg->nweight > 0) + if (pg->params.nweight > 0) { - w = pg->weight[i]; + w = pg->params.weight[i]; } else { @@ -1316,9 +1340,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, if (EI_ENERGY_MINIMIZATION(ir->eI)) { /* Move the mass to the weight */ - w *= m; - m = 1; - pg->weight[i] = w; + w *= m; + m = 1; + pg->params.weight[i] = w; } else if (ir->eI == eiBD) { @@ -1337,9 +1361,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]]; } } - w *= m/mbd; - m = mbd; - pg->weight[i] = w; + w *= m/mbd; + m = mbd; + pg->params.weight[i] = w; } tmass += m; wmass += m*w; @@ -1353,7 +1377,7 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, /* We can have single atom groups with zero mass with potential pulling * without cosine weighting. */ - if (pg->nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS) + if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS) { /* With one atom the mass doesn't matter */ wwmass = 1; @@ -1361,14 +1385,16 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, else { gmx_fatal(FARGS, "The total%s mass of pull group %d is zero", - pg->weight ? " weighted" : "", g); + pg->params.weight ? " weighted" : "", g); } } if (fplog) { fprintf(fplog, - "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass); - if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD) + "Pull group %d: %5d atoms, mass %9.3f", + g, pg->params.nat, tmass); + if (pg->params.weight || + EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD) { fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass); } @@ -1389,7 +1415,7 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, ndim = 0; for (d = 0; d < DIM; d++) { - ndim += pulldim_con[d]*pg->nat; + ndim += pulldim_con[d]*pg->params.nat; } if (fplog && nfrozen > 0 && nfrozen < ndim) { @@ -1404,35 +1430,92 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, } } -void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], - gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda, - gmx_bool bOutFile, unsigned long Flags) +struct pull_t * +init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir, + int nfile, const t_filenm fnm[], + gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda, + gmx_bool bOutFile, unsigned long Flags) { - t_pull *pull; - t_pull_group *pgrp; - int c, g, m; + struct pull_t *pull; + int g, c, m; + + snew(pull, 1); - pull = ir->pull; + /* Copy the pull parameters */ + pull->params = *pull_params; + /* Avoid pointer copies */ + pull->params.group = NULL; + pull->params.coord = NULL; + + pull->ncoord = pull_params->ncoord; + pull->ngroup = pull_params->ngroup; + snew(pull->coord, pull->ncoord); + snew(pull->group, pull->ngroup); pull->bPotential = FALSE; pull->bConstraint = FALSE; pull->bCylinder = FALSE; + + for (g = 0; g < pull->ngroup; g++) + { + pull_group_work_t *pgrp; + int i; + + pgrp = &pull->group[g]; + + /* Copy the pull group parameters */ + pgrp->params = pull_params->group[g]; + + /* Avoid pointer copies by allocating and copying arrays */ + snew(pgrp->params.ind, pgrp->params.nat); + for (i = 0; i < pgrp->params.nat; i++) + { + pgrp->params.ind[i] = pull_params->group[g].ind[i]; + } + if (pgrp->params.nweight > 0) + { + snew(pgrp->params.ind, pgrp->params.nweight); + for (i = 0; i < pgrp->params.nweight; i++) + { + pgrp->params.weight[i] = pull_params->group[g].weight[i]; + } + } + } + for (c = 0; c < pull->ncoord; c++) { - t_pull_coord *pcrd; - int calc_com_start, calc_com_end, g; + pull_coord_work_t *pcrd; + int calc_com_start, calc_com_end, g; pcrd = &pull->coord[c]; - if (pcrd->eType == epullCONSTRAINT) + /* Copy all pull coordinate parameters */ + pcrd->params = pull_params->coord[c]; + + switch (pcrd->params.eGeom) + { + case epullgDIST: + case epullgDIRRELATIVE: + /* Direction vector is determined at each step */ + break; + case epullgDIR: + case epullgDIRPBC: + case epullgCYL: + copy_rvec(pull_params->coord[c].vec, pcrd->vec); + break; + default: + gmx_incons("Pull geometry not handled"); + } + + if (pcrd->params.eType == epullCONSTRAINT) { /* Check restrictions of the constraint pull code */ - if (pcrd->eGeom == epullgCYL || - pcrd->eGeom == epullgDIRRELATIVE) + if (pcrd->params.eGeom == epullgCYL || + pcrd->params.eGeom == epullgDIRRELATIVE) { gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.", - epull_names[pcrd->eType], - epullg_names[pcrd->eGeom], + epull_names[pcrd->params.eType], + epullg_names[pcrd->params.eGeom], epull_names[epullUMBRELLA]); } @@ -1443,26 +1526,26 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], pull->bPotential = TRUE; } - if (pcrd->eGeom == epullgCYL) + if (pcrd->params.eGeom == epullgCYL) { pull->bCylinder = TRUE; } /* We only need to calculate the plain COM of a group * when it is not only used as a cylinder group. */ - calc_com_start = (pcrd->eGeom == epullgCYL ? 1 : 0); - calc_com_end = (pcrd->eGeom == epullgDIRRELATIVE ? 4 : 2); + calc_com_start = (pcrd->params.eGeom == epullgCYL ? 1 : 0); + calc_com_end = (pcrd->params.eGeom == epullgDIRRELATIVE ? 4 : 2); for (g = calc_com_start; g <= calc_com_end; g++) { - pull->group[pcrd->group[g]].bCalcCOM = TRUE; + pull->group[pcrd->params.group[g]].bCalcCOM = TRUE; } /* With non-zero rate the reference value is set at every step */ - if (pcrd->rate == 0) + if (pcrd->params.rate == 0) { /* Initialize the constant reference value */ - pcrd->value_ref = pcrd->init; + pcrd->value_ref = pcrd->params.init; } } @@ -1481,8 +1564,8 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], bAbs = FALSE; for (c = 0; c < pull->ncoord; c++) { - if (pull->group[pull->coord[c].group[0]].nat == 0 || - pull->group[pull->coord[c].group[1]].nat == 0) + if (pull->group[pull->coord[c].params.group[0]].params.nat == 0 || + pull->group[pull->coord[c].params.group[1]].params.nat == 0) { bAbs = TRUE; } @@ -1507,8 +1590,8 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], bCos = FALSE; for (g = 0; g < pull->ngroup; g++) { - if (pull->group[g].nat > 1 && - pull->group[g].pbcatom < 0) + if (pull->group[g].params.nat > 1 && + pull->group[g].params.pbcatom < 0) { /* We are using cosine weighting */ fprintf(fplog, "Cosine weighting is used for group %d\n", g); @@ -1528,9 +1611,11 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], pull->cosdim = -1; for (g = 0; g < pull->ngroup; g++) { + pull_group_work_t *pgrp; + pgrp = &pull->group[g]; pgrp->epgrppbc = epgrppbcNONE; - if (pgrp->nat > 0) + if (pgrp->params.nat > 0) { /* There is an issue when a group is used in multiple coordinates * and constraints are applied in different dimensions with atoms @@ -1552,19 +1637,19 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], clear_ivec(pulldim_con); for (c = 0; c < pull->ncoord; c++) { - if (pull->coord[c].group[0] == g || - pull->coord[c].group[1] == g || - (pull->coord[c].eGeom == epullgDIRRELATIVE && - (pull->coord[c].group[2] == g || - pull->coord[c].group[3] == g))) + if (pull->coord[c].params.group[0] == g || + pull->coord[c].params.group[1] == g || + (pull->coord[c].params.eGeom == epullgDIRRELATIVE && + (pull->coord[c].params.group[2] == g || + pull->coord[c].params.group[3] == g))) { for (m = 0; m < DIM; m++) { - if (pull->coord[c].dim[m] == 1) + if (pull->coord[c].params.dim[m] == 1) { pulldim[m] = 1; - if (pull->coord[c].eType == epullCONSTRAINT) + if (pull->coord[c].params.eType == epullCONSTRAINT) { bConstraint = TRUE; pulldim_con[m] = 1; @@ -1581,16 +1666,16 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], for (m = 0; m < pull->npbcdim; m++) { GMX_RELEASE_ASSERT(m <= DIM, "npbcdim must be <= the number of dimensions"); - if (pulldim[m] == 1 && pgrp->nat > 1) + if (pulldim[m] == 1 && pgrp->params.nat > 1) { - if (pgrp->pbcatom >= 0) + if (pgrp->params.pbcatom >= 0) { pgrp->epgrppbc = epgrppbcREFAT; pull->bRefAt = TRUE; } else { - if (pgrp->weight) + if (pgrp->params.weight != NULL) { gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time"); } @@ -1625,13 +1710,13 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], for (c = 0; c < pull->ncoord; c++) { - const t_pull_coord *pcrd; + const pull_coord_work_t *pcrd; pcrd = &pull->coord[c]; - if (pcrd->eGeom == epullgCYL) + if (pcrd->params.eGeom == epullgCYL) { - if (pull->group[pcrd->group[0]].nat == 0) + if (pull->group[pcrd->params.group[0]].params.nat == 0) { gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n"); } @@ -1647,20 +1732,65 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], pull->out_f = NULL; if (bOutFile) { - if (pull->nstxout > 0) + if (pull->params.nstxout > 0) { pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE, Flags); } - if (pull->nstfout > 0) + if (pull->params.nstfout > 0) { pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv, FALSE, Flags); } } + + return pull; +} + +static void destroy_pull_group(pull_group_work_t *pgrp) +{ + if (pgrp->ind_loc != pgrp->params.ind) + { + sfree(pgrp->ind_loc); + } + if (pgrp->weight_loc != pgrp->params.weight) + { + sfree(pgrp->weight_loc); + } + sfree(pgrp->mdw); + sfree(pgrp->dv); + + sfree(pgrp->params.ind); + sfree(pgrp->params.weight); +} + +static void destroy_pull(struct pull_t *pull) +{ + int i; + + for (i = 0; i < pull->ngroup; i++) + { + destroy_pull_group(&pull->group[i]); + } + for (i = 0; i < pull->ncoord; i++) + { + if (pull->coord[i].params.eGeom == epullgCYL) + { + destroy_pull_group(&(pull->dyna[i])); + } + } + sfree(pull->group); + sfree(pull->dyna); + sfree(pull->coord); + + sfree(pull->rbuf); + sfree(pull->dbuf); + sfree(pull->dbuf_cyl); + + sfree(pull); } -void finish_pull(t_pull *pull) +void finish_pull(struct pull_t *pull) { if (pull->out_x) { @@ -1670,4 +1800,16 @@ void finish_pull(t_pull *pull) { gmx_fio_fclose(pull->out_f); } + + destroy_pull(pull); +} + +gmx_bool pull_have_potential(const struct pull_t *pull) +{ + return pull->bPotential; +} + +gmx_bool pull_have_constraint(const struct pull_t *pull) +{ + return pull->bConstraint; } diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h index f1f49bd7ac..c73ca640dc 100644 --- a/src/gromacs/pulling/pull.h +++ b/src/gromacs/pulling/pull.h @@ -67,7 +67,7 @@ struct t_pbc; * \param[in] pbc Information structure about periodicity. * \param[out] value The value of the pull coordinate. */ -void get_pull_coord_value(t_pull *pull, +void get_pull_coord_value(struct pull_t *pull, int coord_ind, const struct t_pbc *pbc, double *value); @@ -77,7 +77,7 @@ void get_pull_coord_value(t_pull *pull, * * \param pull The pull group. */ -void clear_pull_forces(t_pull *pull); +void clear_pull_forces(struct pull_t *pull); /*! \brief Determine the COM pull forces and add them to f, return the potential @@ -95,7 +95,7 @@ void clear_pull_forces(t_pull *pull); * * \returns The pull potential energy. */ -real pull_potential(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc, +real pull_potential(struct pull_t *pull, t_mdatoms *md, struct t_pbc *pbc, t_commrec *cr, double t, real lambda, rvec *x, rvec *f, tensor vir, real *dvdlambda); @@ -114,7 +114,7 @@ real pull_potential(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc, * \param[in,out] v Velocities, which may get a pull correction. * \param[in,out] vir The virial, which, if != NULL, gets a pull correction. */ -void pull_constraint(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc, +void pull_constraint(struct pull_t *pull, t_mdatoms *md, struct t_pbc *pbc, t_commrec *cr, double dt, double t, rvec *x, rvec *xp, rvec *v, tensor vir); @@ -127,56 +127,57 @@ void pull_constraint(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc, * \param md All atoms. */ void dd_make_local_pull_groups(gmx_domdec_t *dd, - t_pull *pull, t_mdatoms *md); - - -/*! \brief Get memory and initialize the fields of pull that still need it, and - * do runtype specific initialization. - * - * \param fplog General output file, normally md.log. - * \param ir The inputrec. - * \param nfile Number of files. - * \param fnm Standard filename struct. - * \param mtop The topology of the whole system. - * \param cr Struct for communication info. - * \param oenv Output options. - * \param lambda FEP lambda. - * \param bOutFile Open output files? - * \param Flags Flags passed over from main, used to determine - * whether or not we are appending. + struct pull_t *pull, t_mdatoms *md); + + +/*! \brief Allocate, initialize and return a pull work struct. + * + * \param fplog General output file, normally md.log. + * \param pull_params The pull input parameters containing all pull settings. + * \param ir The inputrec. + * \param nfile Number of files. + * \param fnm Standard filename struct. + * \param mtop The topology of the whole system. + * \param cr Struct for communication info. + * \param oenv Output options. + * \param lambda FEP lambda. + * \param bOutFile Open output files? + * \param Flags Flags passed over from main, used to determine + * whether or not we are appending. */ -void init_pull(FILE *fplog, - t_inputrec *ir, - int nfile, - const t_filenm fnm[], - gmx_mtop_t *mtop, - t_commrec * cr, - const output_env_t oenv, - real lambda, - gmx_bool bOutFile, - unsigned long Flags); +struct pull_t *init_pull(FILE *fplog, + const pull_params_t *pull_params, + const t_inputrec *ir, + int nfile, + const t_filenm fnm[], + gmx_mtop_t *mtop, + t_commrec * cr, + const output_env_t oenv, + real lambda, + gmx_bool bOutFile, + unsigned long Flags); /*! \brief Close the pull output files. * * \param pull The pull group. */ -void finish_pull(t_pull *pull); +void finish_pull(struct pull_t *pull); /*! \brief Print the pull output (x and/or f) * - * \param pull The pull group. + * \param pull The pull data structure. * \param step Time step number. * \param time Time. */ -void pull_print_output(t_pull *pull, gmx_int64_t step, double time); +void pull_print_output(struct pull_t *pull, gmx_int64_t step, double time); /*! \brief Calculates centers of mass all pull groups. * * \param[in] cr Struct for communication info. - * \param[in] pull The pull group. + * \param[in] pull The pull data structure. * \param[in] md All atoms. * \param[in] pbc Information struct about periodicity. * \param[in] t Time, only used for cylinder ref. @@ -185,13 +186,27 @@ void pull_print_output(t_pull *pull, gmx_int64_t step, double time); * */ void pull_calc_coms(t_commrec *cr, - t_pull *pull, + struct pull_t *pull, t_mdatoms *md, struct t_pbc *pbc, double t, rvec x[], rvec *xp); + +/*! \brief Returns if we have pull coordinates with potential pulling. + * + * \param[in] pull The pull data structure. + */ +gmx_bool pull_have_potential(const struct pull_t *pull); + + +/*! \brief Returns if we have pull coordinates with constraint pulling. + * + * \param[in] pull The pull data structure. + */ +gmx_bool pull_have_constraint(const struct pull_t *pull); + #ifdef __cplusplus } #endif diff --git a/src/gromacs/pulling/pull_internal.h b/src/gromacs/pulling/pull_internal.h new file mode 100644 index 0000000000..ad5c10f0e0 --- /dev/null +++ b/src/gromacs/pulling/pull_internal.h @@ -0,0 +1,133 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. + * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, + * and including many others, as listed in the AUTHORS file in the + * top-level source directory and at http://www.gromacs.org. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ + +/*! \libinternal \file + * + * + * \brief + * This file contains datatypes and function declarations for internal + use in the pull code. + * + * \author Berk Hess + * + * \inlibraryapi + */ + +#ifndef GMX_PULLING_PULL_INTERNAL_H +#define GMX_PULLING_PULL_INTERNAL_H + +#include "gromacs/legacyheaders/typedefs.h" + +#ifdef __cplusplus +extern "C" { +#endif + +enum { + epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS +}; + +typedef struct +{ + t_pull_group params; + + gmx_bool bCalcCOM; /* Calculate COM? Not if only used as cylinder group */ + int epgrppbc; /* The type of pbc for this pull group, see enum above */ + + int nat_loc; /* Number of local pull atoms */ + int nalloc_loc; /* Allocation size for ind_loc and weight_loc */ + atom_id *ind_loc; /* Local pull indices */ + real *weight_loc; /* Weights for the local indices */ + + real mwscale; /* mass*weight scaling factor 1/sum w m */ + real wscale; /* scaling factor for the weights: sum w m/sum w w m */ + real invtm; /* inverse total mass of the group: 1/wscale sum w m */ + dvec *mdw; /* mass*gradient(weight) for atoms */ + double *dv; /* distance to the other group along vec */ + dvec x; /* center of mass before update */ + dvec xp; /* center of mass after update before constraining */ +} +pull_group_work_t; + +typedef struct +{ + t_pull_coord params; /* Pull coordinate (constant) parameters */ + + double value_ref; /* The reference value, usually init+rate*t */ + double value; /* The current value of the coordinate */ + dvec dr; /* The distance from the reference group */ + rvec vec; /* The pull direction */ + double vec_len; /* Length of vec for direction-relative */ + dvec ffrad; /* conversion factor from vec to radial force */ + double cyl_dev; /* The deviation from the reference position */ + double f_scal; /* Scalar force for directional pulling */ + dvec f; /* force due to the pulling/constraining */ +} +pull_coord_work_t; + +struct pull_t +{ + pull_params_t params; /* The pull parameters, from inputrec */ + + gmx_bool bPotential; /* Are there coordinates with potential? */ + gmx_bool bConstraint; /* Are there constrained coordinates? */ + + int ePBC; /* the boundary conditions */ + int npbcdim; /* do pbc in dims 0 <= dim < npbcdim */ + gmx_bool bRefAt; /* do we need reference atoms for a group COM ? */ + int cosdim; /* dimension for cosine weighting, -1 if none */ + + int ngroup; /* Number of pull groups */ + int ncoord; /* Number of pull coordinates */ + pull_group_work_t *group; /* The pull group param and work data */ + pull_group_work_t *dyna; /* Dynamic groups for geom=cylinder */ + pull_coord_work_t *coord; /* The pull group param and work data */ + + gmx_bool bCylinder; /* Is group 0 a cylinder group? */ + gmx_bool bSetPBCatoms; /* Do we need to set x_pbc for the groups? */ + + rvec *rbuf; /* COM calculation buffer */ + dvec *dbuf; /* COM calculation buffer */ + double *dbuf_cyl; /* cylinder ref. groups COM calculation buffer */ + FILE *out_x; /* Output file for pull data */ + FILE *out_f; /* Output file for pull data */ +}; + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/gromacs/pulling/pullutil.cpp b/src/gromacs/pulling/pullutil.cpp index 852acc9521..7345aad857 100644 --- a/src/gromacs/pulling/pullutil.cpp +++ b/src/gromacs/pulling/pullutil.cpp @@ -49,11 +49,12 @@ #include "gromacs/math/vec.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/pull.h" +#include "gromacs/pulling/pull_internal.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/futil.h" #include "gromacs/utility/smalloc.h" -static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp, +static void pull_set_pbcatom(t_commrec *cr, pull_group_work_t *pgrp, rvec *x, rvec x_pbc) { @@ -61,7 +62,7 @@ static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp, if (cr != NULL && DOMAINDECOMP(cr)) { - if (ga2la_get_home(cr->dd->ga2la, pgrp->pbcatom, &a)) + if (ga2la_get_home(cr->dd->ga2la, pgrp->params.pbcatom, &a)) { copy_rvec(x[a], x_pbc); } @@ -72,11 +73,11 @@ static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp, } else { - copy_rvec(x[pgrp->pbcatom], x_pbc); + copy_rvec(x[pgrp->params.pbcatom], x_pbc); } } -static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull, +static void pull_set_pbcatoms(t_commrec *cr, struct pull_t *pull, rvec *x, rvec *x_pbc) { @@ -85,7 +86,7 @@ static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull, n = 0; for (g = 0; g < pull->ngroup; g++) { - if (!pull->group[g].bCalcCOM || pull->group[g].pbcatom == -1) + if (!pull->group[g].bCalcCOM || pull->group[g].params.pbcatom == -1) { clear_rvec(x_pbc[g]); } @@ -103,7 +104,7 @@ static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull, } } -static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, +static void make_cyl_refgrps(t_commrec *cr, struct pull_t *pull, t_mdatoms *md, t_pbc *pbc, double t, rvec *x) { /* The size and stride per coord for the reduction buffer */ @@ -111,8 +112,6 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, int c, i, ii, m, start, end; rvec g_x, dx, dir; double inv_cyl_r2; - t_pull_coord *pcrd; - t_pull_group *pref, *pgrp, *pdyna; gmx_ga2la_t ga2la = NULL; if (pull->dbuf_cyl == NULL) @@ -128,13 +127,14 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, start = 0; end = md->homenr; - inv_cyl_r2 = 1/dsqr(pull->cylinder_r); + inv_cyl_r2 = 1/dsqr(pull->params.cylinder_r); /* loop over all groups to make a reference group for each*/ for (c = 0; c < pull->ncoord; c++) { - double sum_a, wmass, wwmass; - dvec radf_fac0, radf_fac1; + pull_coord_work_t *pcrd; + double sum_a, wmass, wwmass; + dvec radf_fac0, radf_fac1; pcrd = &pull->coord[c]; @@ -144,11 +144,13 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, clear_dvec(radf_fac0); clear_dvec(radf_fac1); - if (pcrd->eGeom == epullgCYL) + if (pcrd->params.eGeom == epullgCYL) { + pull_group_work_t *pref, *pgrp, *pdyna; + /* pref will be the same group for all pull coordinates */ - pref = &pull->group[pcrd->group[0]]; - pgrp = &pull->group[pcrd->group[1]]; + pref = &pull->group[pcrd->params.group[0]]; + pgrp = &pull->group[pcrd->params.group[1]]; pdyna = &pull->dyna[c]; copy_rvec(pcrd->vec, dir); pdyna->nat_loc = 0; @@ -158,10 +160,10 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, * we reduced the other group COM over the ranks. This resolves * any PBC issues and we don't need to use a PBC-atom here. */ - if (pcrd->rate != 0) + if (pcrd->params.rate != 0) { /* With rate=0, value_ref is set initially */ - pcrd->value_ref = pcrd->init + pcrd->rate*t; + pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t; } for (m = 0; m < DIM; m++) { @@ -169,12 +171,12 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, } /* loop over all atoms in the main ref group */ - for (i = 0; i < pref->nat; i++) + for (i = 0; i < pref->params.nat; i++) { - ii = pref->ind[i]; + ii = pref->params.ind[i]; if (ga2la) { - if (!ga2la_get_home(ga2la, pref->ind[i], &ii)) + if (!ga2la_get_home(ga2la, pref->params.ind[i], &ii)) { ii = -1; } @@ -260,14 +262,17 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, for (c = 0; c < pull->ncoord; c++) { + pull_coord_work_t *pcrd; + pcrd = &pull->coord[c]; - if (pcrd->eGeom == epullgCYL) + if (pcrd->params.eGeom == epullgCYL) { - double wmass, wwmass, dist; + pull_group_work_t *pdyna, *pgrp; + double wmass, wwmass, dist; pdyna = &pull->dyna[c]; - pgrp = &pull->group[pcrd->group[1]]; + pgrp = &pull->group[pcrd->params.group[1]]; wmass = pull->dbuf_cyl[c*stride+0]; wwmass = pull->dbuf_cyl[c*stride+1]; @@ -327,7 +332,7 @@ static double atan2_0_2pi(double y, double x) /* calculates center of mass of selection index from all coordinates x */ void pull_calc_coms(t_commrec *cr, - t_pull *pull, t_mdatoms *md, t_pbc *pbc, double t, + struct pull_t *pull, t_mdatoms *md, t_pbc *pbc, double t, rvec x[], rvec *xp) { int g; @@ -377,7 +382,7 @@ void pull_calc_coms(t_commrec *cr, for (g = 0; g < pull->ngroup; g++) { - t_pull_group *pgrp; + pull_group_work_t *pgrp; pgrp = &pull->group[g]; @@ -467,7 +472,7 @@ void pull_calc_coms(t_commrec *cr, * Note that with constraint pulling the mass does matter, but * in that case a check group mass != 0 has been done before. */ - if (pgrp->nat == 1 && pgrp->nat_loc == 1 && wmass == 0) + if (pgrp->params.nat == 1 && pgrp->nat_loc == 1 && wmass == 0) { int m; @@ -554,10 +559,10 @@ void pull_calc_coms(t_commrec *cr, for (g = 0; g < pull->ngroup; g++) { - t_pull_group *pgrp; + pull_group_work_t *pgrp; pgrp = &pull->group[g]; - if (pgrp->nat > 0 && pgrp->bCalcCOM) + if (pgrp->params.nat > 0 && pgrp->bCalcCOM) { if (pgrp->epgrppbc != epgrppbcCOS) { diff --git a/src/gromacs/tools/compare.c b/src/gromacs/tools/compare.c index 7be156ab53..fd73cd8c97 100644 --- a/src/gromacs/tools/compare.c +++ b/src/gromacs/tools/compare.c @@ -888,7 +888,7 @@ static void cmp_inputrec(FILE *fp, t_inputrec *ir1, t_inputrec *ir2, real ftol, cmp_cosines(fp, "et", ir1->et, ir2->et, ftol, abstol); } -static void comp_pull_AB(FILE *fp, t_pull *pull, real ftol, real abstol) +static void comp_pull_AB(FILE *fp, pull_params_t *pull, real ftol, real abstol) { int i; diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp index c377038774..a839ee8adb 100644 --- a/src/programs/mdrun/md.cpp +++ b/src/programs/mdrun/md.cpp @@ -1603,7 +1603,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], } if (ir->bPull) { - pull_print_output(ir->pull, step, t); + pull_print_output(ir->pull_work, step, t); } if (do_per_step(step, ir->nstlog)) diff --git a/src/programs/mdrun/runner.cpp b/src/programs/mdrun/runner.cpp index 613611896c..d88c5cd656 100644 --- a/src/programs/mdrun/runner.cpp +++ b/src/programs/mdrun/runner.cpp @@ -1250,8 +1250,10 @@ int mdrunner(gmx_hw_opt_t *hw_opt, if (inputrec->bPull) { /* Initialize pull code */ - init_pull(fplog, inputrec, nfile, fnm, mtop, cr, oenv, inputrec->fepvals->init_lambda, - EI_DYNAMICS(inputrec->eI) && MASTER(cr), Flags); + inputrec->pull_work = + init_pull(fplog, inputrec->pull, inputrec, nfile, fnm, + mtop, cr, oenv, inputrec->fepvals->init_lambda, + EI_DYNAMICS(inputrec->eI) && MASTER(cr), Flags); } if (inputrec->bRot) @@ -1298,7 +1300,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt, if (inputrec->bPull) { - finish_pull(inputrec->pull); + finish_pull(inputrec->pull_work); } if (inputrec->bRot)