From: Berk Hess Date: Thu, 7 Mar 2019 09:52:29 +0000 (+0100) Subject: Clean up calc_verlet_buffer_size() X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=415aae89abb57864a1118e287536d654cd412ef5;p=alexxy%2Fgromacs.git Clean up calc_verlet_buffer_size() Use return value instead of passing pointer. Added a separate vsite function for counting the non-linear vsites. Changed pointers to references. Improve variable naming This change is only refactoring. Change-Id: I91b57dee08452d5e44b5b38860ffbd5e01430a3a --- diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index f669c50b44..24fadfbb0c 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -81,6 +81,7 @@ #include "gromacs/mdlib/perf_est.h" #include "gromacs/mdlib/qmmm.h" #include "gromacs/mdlib/sim_util.h" +#include "gromacs/mdlib/vsite.h" #include "gromacs/mdrunutility/mdmodules.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" @@ -1547,8 +1548,6 @@ static void set_verlet_buffer(const gmx_mtop_t *mtop, matrix box, warninp *wi) { - real rlist_1x1; - int n_nonlin_vsite; char warn_buf[STRLEN]; printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, buffer_temp); @@ -1557,17 +1556,18 @@ static void set_verlet_buffer(const gmx_mtop_t *mtop, VerletbufListSetup listSetup1x1; listSetup1x1.cluster_size_i = 1; listSetup1x1.cluster_size_j = 1; - calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1, - buffer_temp, &listSetup1x1, - &n_nonlin_vsite, &rlist_1x1); + const real rlist_1x1 = + calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, + buffer_temp, listSetup1x1); /* Set the pair-list buffer size in ir */ VerletbufListSetup listSetup4x4 = verletbufGetSafeListSetup(ListSetupType::CpuNoSimd); - calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1, - buffer_temp, &listSetup4x4, - &n_nonlin_vsite, &ir->rlist); + ir->rlist = + calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, + buffer_temp, listSetup4x4); + const int n_nonlin_vsite = countNonlinearVsites(*mtop); if (n_nonlin_vsite > 0) { sprintf(warn_buf, "There are %d non-linear virtual site constructions. Their contribution to the energy error is approximated. In most cases this does not affect the error significantly.", n_nonlin_vsite); diff --git a/src/gromacs/mdlib/calc_verletbuf.cpp b/src/gromacs/mdlib/calc_verletbuf.cpp index b494e16652..9183b9c095 100644 --- a/src/gromacs/mdlib/calc_verletbuf.cpp +++ b/src/gromacs/mdlib/calc_verletbuf.cpp @@ -215,12 +215,9 @@ static real getMass(const t_atoms &atoms, static void get_vsite_masses(const gmx_moltype_t &moltype, const gmx_ffparams_t &ffparams, bool setMassesToOne, - gmx::ArrayRef vsite_m, - int *n_nonlin_vsite) + gmx::ArrayRef vsite_m) { - GMX_RELEASE_ASSERT(n_nonlin_vsite, "Expect a valid pointer"); - - *n_nonlin_vsite = 0; + int numNonlinearVsites = 0; /* Check for virtual sites, determine mass from constructing atoms */ for (const auto &ilist : extractILists(moltype.ilist, IF_VSITE)) @@ -283,7 +280,7 @@ static void get_vsite_masses(const gmx_moltype_t &moltype, { vsite_m[a1] = std::min(vsite_m[a1], cam[j]); } - (*n_nonlin_vsite)++; + numNonlinearVsites++; break; } } @@ -322,27 +319,26 @@ static void get_vsite_masses(const gmx_moltype_t &moltype, } } } + + if (debug && numNonlinearVsites > 0) + { + fprintf(debug, "The molecule type has %d non-linear virtual constructions\n", + numNonlinearVsites); + } } static std::vector -get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, - bool setMassesToOne, - int *n_nonlin_vsite) +getVerletBufferAtomtypes(const gmx_mtop_t &mtop, + const bool setMassesToOne) { std::vector att; int ft, i, a1, a2, a3, a; const t_iparams *ip; - int n_nonlin_vsite_mol; - - if (n_nonlin_vsite != nullptr) - { - *n_nonlin_vsite = 0; - } - for (const gmx_molblock_t &molblock : mtop->molblock) + for (const gmx_molblock_t &molblock : mtop.molblock) { int nmol = molblock.nmol; - const gmx_moltype_t &moltype = mtop->moltype[molblock.type]; + const gmx_moltype_t &moltype = mtop.moltype[molblock.type]; const t_atoms *atoms = &moltype.atoms; /* Check for constraints, as they affect the kinetic energy. @@ -359,7 +355,7 @@ get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, for (i = 0; i < il.size(); i += 1+NRAL(ft)) { - ip = &mtop->ffparams.iparams[il.iatoms[i]]; + ip = &mtop.ffparams.iparams[il.iatoms[i]]; a1 = il.iatoms[i+1]; a2 = il.iatoms[i+2]; real mass1 = getMass(*atoms, a1, setMassesToOne); @@ -381,7 +377,7 @@ get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, for (i = 0; i < il.size(); i += 1+NRAL(F_SETTLE)) { - ip = &mtop->ffparams.iparams[il.iatoms[i]]; + ip = &mtop.ffparams.iparams[il.iatoms[i]]; a1 = il.iatoms[i+1]; a2 = il.iatoms[i+2]; a3 = il.iatoms[i+3]; @@ -401,14 +397,9 @@ get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, std::vector vsite_m(atoms->nr); get_vsite_masses(moltype, - mtop->ffparams, + mtop.ffparams, setMassesToOne, - vsite_m, - &n_nonlin_vsite_mol); - if (n_nonlin_vsite != nullptr) - { - *n_nonlin_vsite += nmol*n_nonlin_vsite_mol; - } + vsite_m); for (a = 0; a < atoms->nr; a++) { @@ -854,14 +845,14 @@ static real maxSigma(real kT_fac, return 2*std::sqrt(kT_fac/smallestMass); } -void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, - const t_inputrec *ir, - int nstlist, - int list_lifetime, - real reference_temperature, - const VerletbufListSetup *list_setup, - int *n_nonlin_vsite, - real *rlist) +real +calcVerletBufferSize(const gmx_mtop_t &mtop, + const real boxVolume, + const t_inputrec &ir, + const int nstlist, + const int listLifetime, + real referenceTemperature, + const VerletbufListSetup &listSetup) { double resolution; char *env; @@ -874,24 +865,24 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, real rb, rl; real drift; - if (!EI_DYNAMICS(ir->eI)) + if (!EI_DYNAMICS(ir.eI)) { gmx_incons("Can only determine the Verlet buffer size for integrators that perform dynamics"); } - if (ir->verletbuf_tol <= 0) + if (ir.verletbuf_tol <= 0) { gmx_incons("The Verlet buffer tolerance needs to be larger than zero"); } - if (reference_temperature < 0) + if (referenceTemperature < 0) { /* We use the maximum temperature with multiple T-coupl groups. * We could use a per particle temperature, but since particles * interact, this might underestimate the buffer size. */ - reference_temperature = maxReferenceTemperature(*ir); + referenceTemperature = maxReferenceTemperature(ir); - GMX_RELEASE_ASSERT(reference_temperature >= 0, "Without T-coupling we should not end up here"); + GMX_RELEASE_ASSERT(referenceTemperature >= 0, "Without T-coupling we should not end up here"); } /* Resolution of the buffer size */ @@ -928,14 +919,14 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, */ /* Worst case assumption: HCP packing of particles gives largest distance */ - particle_distance = std::cbrt(boxvol*std::sqrt(2)/mtop->natoms); + particle_distance = std::cbrt(boxVolume*std::sqrt(2)/mtop.natoms); /* TODO: Obtain masses through (future) integrator functionality * to avoid scattering the code with (or forgetting) checks. */ - const bool setMassesToOne = (ir->eI == eiBD && ir->bd_fric > 0); + const bool setMassesToOne = (ir.eI == eiBD && ir.bd_fric > 0); const auto att = - get_verlet_buffer_atomtypes(mtop, setMassesToOne, n_nonlin_vsite); + getVerletBufferAtomtypes(mtop, setMassesToOne); GMX_ASSERT(!att.empty(), "We expect at least one type"); if (debug) @@ -947,45 +938,45 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, pot_derivatives_t ljDisp = { 0, 0, 0 }; pot_derivatives_t ljRep = { 0, 0, 0 }; - real repPow = mtop->ffparams.reppow; + real repPow = mtop.ffparams.reppow; - if (ir->vdwtype == evdwCUT) + if (ir.vdwtype == evdwCUT) { real sw_range, md3_pswf; - switch (ir->vdw_modifier) + switch (ir.vdw_modifier) { case eintmodNONE: case eintmodPOTSHIFT: /* -dV/dr of -r^-6 and r^-reppow */ - ljDisp.md1 = -6*std::pow(ir->rvdw, -7.0); - ljRep.md1 = repPow*std::pow(ir->rvdw, -(repPow + 1)); + ljDisp.md1 = -6*std::pow(ir.rvdw, -7.0); + ljRep.md1 = repPow*std::pow(ir.rvdw, -(repPow + 1)); /* The contribution of the higher derivatives is negligible */ break; case eintmodFORCESWITCH: /* At the cut-off: V=V'=V''=0, so we use only V''' */ - ljDisp.md3 = -md3_force_switch(6.0, ir->rvdw_switch, ir->rvdw); - ljRep.md3 = md3_force_switch(repPow, ir->rvdw_switch, ir->rvdw); + ljDisp.md3 = -md3_force_switch(6.0, ir.rvdw_switch, ir.rvdw); + ljRep.md3 = md3_force_switch(repPow, ir.rvdw_switch, ir.rvdw); break; case eintmodPOTSWITCH: /* At the cut-off: V=V'=V''=0. * V''' is given by the original potential times * the third derivative of the switch function. */ - sw_range = ir->rvdw - ir->rvdw_switch; + sw_range = ir.rvdw - ir.rvdw_switch; md3_pswf = 60.0/gmx::power3(sw_range); - ljDisp.md3 = -std::pow(ir->rvdw, -6.0 )*md3_pswf; - ljRep.md3 = std::pow(ir->rvdw, -repPow)*md3_pswf; + ljDisp.md3 = -std::pow(ir.rvdw, -6.0 )*md3_pswf; + ljRep.md3 = std::pow(ir.rvdw, -repPow)*md3_pswf; break; default: gmx_incons("Unimplemented VdW modifier"); } } - else if (EVDW_PME(ir->vdwtype)) + else if (EVDW_PME(ir.vdwtype)) { - real b = calc_ewaldcoeff_lj(ir->rvdw, ir->ewald_rtol_lj); - real r = ir->rvdw; + real b = calc_ewaldcoeff_lj(ir.rvdw, ir.ewald_rtol_lj); + real r = ir.rvdw; real br = b*r; real br2 = br*br; real br4 = br2*br2; @@ -1001,46 +992,46 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, gmx_fatal(FARGS, "Energy drift calculation is only implemented for plain cut-off Lennard-Jones interactions"); } - elfac = ONE_4PI_EPS0/ir->epsilon_r; + elfac = ONE_4PI_EPS0/ir.epsilon_r; // Determine the 1st and 2nd derivative for the electostatics pot_derivatives_t elec = { 0, 0, 0 }; - if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype)) + if (ir.coulombtype == eelCUT || EEL_RF(ir.coulombtype)) { real eps_rf, k_rf; - if (ir->coulombtype == eelCUT) + if (ir.coulombtype == eelCUT) { eps_rf = 1; k_rf = 0; } else { - eps_rf = ir->epsilon_rf/ir->epsilon_r; + eps_rf = ir.epsilon_rf/ir.epsilon_r; if (eps_rf != 0) { - k_rf = (eps_rf - ir->epsilon_r)/( gmx::power3(ir->rcoulomb) * (2*eps_rf + ir->epsilon_r) ); + k_rf = (eps_rf - ir.epsilon_r)/( gmx::power3(ir.rcoulomb) * (2*eps_rf + ir.epsilon_r) ); } else { /* epsilon_rf = infinity */ - k_rf = 0.5/gmx::power3(ir->rcoulomb); + k_rf = 0.5/gmx::power3(ir.rcoulomb); } } if (eps_rf > 0) { - elec.md1 = elfac*(1.0/gmx::square(ir->rcoulomb) - 2*k_rf*ir->rcoulomb); + elec.md1 = elfac*(1.0/gmx::square(ir.rcoulomb) - 2*k_rf*ir.rcoulomb); } - elec.d2 = elfac*(2.0/gmx::power3(ir->rcoulomb) + 2*k_rf); + elec.d2 = elfac*(2.0/gmx::power3(ir.rcoulomb) + 2*k_rf); } - else if (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD) + else if (EEL_PME(ir.coulombtype) || ir.coulombtype == eelEWALD) { real b, rc, br; - b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol); - rc = ir->rcoulomb; + b = calc_ewaldcoeff_q(ir.rcoulomb, ir.ewald_rtol); + rc = ir.rcoulomb; br = b*rc; elec.md1 = elfac*(b*std::exp(-br*br)*M_2_SQRTPI/rc + std::erfc(br)/(rc*rc)); elec.d2 = elfac/(rc*rc)*(2*b*(1 + br*br)*std::exp(-br*br)*M_2_SQRTPI + 2*std::erfc(br)/rc); @@ -1055,8 +1046,8 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, * For inertial dynamics (not Brownian dynamics) the mass factor * is not included in kT_fac, it is added later. */ - const real kT_fac = displacementVariance(*ir, reference_temperature, - list_lifetime*ir->delta_t); + const real kT_fac = displacementVariance(ir, referenceTemperature, + listLifetime*ir.delta_t); if (debug) { @@ -1075,41 +1066,41 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, { ib = (ib0 + ib1)/2; rb = ib*resolution; - rl = std::max(ir->rvdw, ir->rcoulomb) + rb; + rl = std::max(ir.rvdw, ir.rcoulomb) + rb; /* Calculate the average energy drift at the last step * of the nstlist steps at which the pair-list is used. */ - drift = energyDrift(att, &mtop->ffparams, + drift = energyDrift(att, &mtop.ffparams, kT_fac, &ljDisp, &ljRep, &elec, - ir->rvdw, ir->rcoulomb, - rl, boxvol); + ir.rvdw, ir.rcoulomb, + rl, boxVolume); /* Correct for the fact that we are using a Ni x Nj particle pair list * and not a 1 x 1 particle pair list. This reduces the drift. */ /* We don't have a formula for 8 (yet), use 4 which is conservative */ nb_clust_frac_pairs_not_in_list_at_cutoff = - surface_frac(std::min(list_setup->cluster_size_i, 4), + surface_frac(std::min(listSetup.cluster_size_i, 4), particle_distance, rl)* - surface_frac(std::min(list_setup->cluster_size_j, 4), + surface_frac(std::min(listSetup.cluster_size_j, 4), particle_distance, rl); drift *= nb_clust_frac_pairs_not_in_list_at_cutoff; /* Convert the drift to drift per unit time per atom */ - drift /= nstlist*ir->delta_t*mtop->natoms; + drift /= nstlist*ir.delta_t*mtop.natoms; if (debug) { fprintf(debug, "ib %3d %3d %3d rb %.3f %dx%d fac %.3f drift %.1e\n", ib0, ib, ib1, rb, - list_setup->cluster_size_i, list_setup->cluster_size_j, + listSetup.cluster_size_i, listSetup.cluster_size_j, nb_clust_frac_pairs_not_in_list_at_cutoff, drift); } - if (std::abs(drift) > ir->verletbuf_tol) + if (std::abs(drift) > ir.verletbuf_tol) { ib0 = ib; } @@ -1119,7 +1110,7 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, } } - *rlist = std::max(ir->rvdw, ir->rcoulomb) + ib1*resolution; + return std::max(ir.rvdw, ir.rcoulomb) + ib1*resolution; } /* Returns the pairlist buffer size for use as a minimum buffer size @@ -1363,7 +1354,7 @@ minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop, const bool setMassesToOne = (ir.eI == eiBD && ir.bd_fric > 0); - const auto atomtypes = get_verlet_buffer_atomtypes(&mtop, setMassesToOne, nullptr); + const auto atomtypes = getVerletBufferAtomtypes(mtop, setMassesToOne); const real kT_fac = displacementVariance(ir, temperature, ir.nstlist*ir.delta_t); diff --git a/src/gromacs/mdlib/calc_verletbuf.h b/src/gromacs/mdlib/calc_verletbuf.h index 1a8cdc34ee..e0bc6bf650 100644 --- a/src/gromacs/mdlib/calc_verletbuf.h +++ b/src/gromacs/mdlib/calc_verletbuf.h @@ -90,27 +90,38 @@ enum class ListSetupType */ VerletbufListSetup verletbufGetSafeListSetup(ListSetupType listType); -/* Calculate the non-bonded pair-list buffer size for the Verlet list +/* Returns the non-bonded pair-list radius including computed buffer + * + * Calculate the non-bonded pair-list buffer size for the Verlet list * based on the particle masses, temperature, LJ types, charges * and constraints as well as the non-bonded force behavior at the cut-off. * The pair list update frequency and the list lifetime, which is nstlist-1 * for normal pair-list buffering, are passed separately, as in some cases * we want an estimate for different values than the ones set in the inputrec. - * If reference_temperature < 0, the maximum coupling temperature will be used. + * If referenceTemperature < 0, the maximum coupling temperature will be used. * The target is a maximum average energy jump per atom of - * ir->verletbuf_tol*nstlist*ir->delta_t over the lifetime of the list. - * Returns the number of non-linear virtual sites. For these it's difficult - * to determine their contribution to the drift exaclty, so we approximate. - * Returns the pair-list cut-off. + * inputrec.verletbuf_tol*nstlist*inputrec.delta_t over the lifetime of the list. + * + * \note For non-linear virtual sites it can be problematic to determine their + * contribution to the drift exaclty, so we approximate. + * + * \param[in] mtop The system topology + * \param[in] boxVolume The volume of the unit cell + * \param[in] inputrec The input record + * \param[in] nstlist The pair list update frequency in steps (is not taken from \p inputrec) + * \param[in] listLifetime The lifetime of the pair-list, usually nstlist-1, but could be different for dynamic pruning + * \param[in] referenceTemperature The reference temperature for the ensemble + * \param[in] listSetup The pair-list setup + * \returns The computed pair-list radius including buffer */ -void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, - const t_inputrec *ir, - int nstlist, - int list_lifetime, - real reference_temperature, - const VerletbufListSetup *list_setup, - int *n_nonlin_vsite, - real *rlist); +real +calcVerletBufferSize(const gmx_mtop_t &mtop, + real boxVolume, + const t_inputrec &inputrec, + int nstlist, + int listLifetime, + real referenceTemperature, + const VerletbufListSetup &listSetup); /* Convenience type */ using PartitioningPerMoltype = gmx::ArrayRef; diff --git a/src/gromacs/mdlib/vsite.cpp b/src/gromacs/mdlib/vsite.cpp index cbbae59a48..fda9dd0edd 100644 --- a/src/gromacs/mdlib/vsite.cpp +++ b/src/gromacs/mdlib/vsite.cpp @@ -1844,6 +1844,27 @@ static std::vector atom2cg(const t_block &chargeGroups) return a2cg; } +int countNonlinearVsites(const gmx_mtop_t &mtop) +{ + int numNonlinearVsites = 0; + for (const gmx_molblock_t &molb : mtop.molblock) + { + const gmx_moltype_t &molt = mtop.moltype[molb.type]; + + for (const auto &ilist : extractILists(molt.ilist, IF_VSITE)) + { + if (ilist.functionType != F_VSITE2 && + ilist.functionType != F_VSITE3 && + ilist.functionType != F_VSITEN) + { + numNonlinearVsites += molb.nmol*ilist.iatoms.size()/(1 + NRAL(ilist.functionType)); + } + } + } + + return numNonlinearVsites; +} + int count_intercg_vsites(const gmx_mtop_t *mtop) { int n_intercg_vsite = 0; diff --git a/src/gromacs/mdlib/vsite.h b/src/gromacs/mdlib/vsite.h index b82751206a..b700bf6d4d 100644 --- a/src/gromacs/mdlib/vsite.h +++ b/src/gromacs/mdlib/vsite.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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. @@ -128,6 +128,9 @@ void spread_vsite_f(const gmx_vsite_t *vsite, * as for instance for the PME mesh contribution. */ +/* Return the number of non-linear virtual site constructions in the system */ +int countNonlinearVsites(const gmx_mtop_t &mtop); + int count_intercg_vsites(const gmx_mtop_t *mtop); /* Returns the number of virtual sites that cross charge groups */ diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index d5f6614e63..2d7eceef2c 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -283,8 +283,8 @@ static void prepare_verlet_scheme(FILE *fplog, ListSetupType listType = (makeGpuPairList ? ListSetupType::Gpu : ListSetupType::CpuSimdWhenSupported); VerletbufListSetup listSetup = verletbufGetSafeListSetup(listType); - real rlist_new; - calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1, -1, &listSetup, nullptr, &rlist_new); + const real rlist_new = + calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, -1, listSetup); if (rlist_new != ir->rlist) { diff --git a/src/gromacs/nbnxm/pairlist_tuning.cpp b/src/gromacs/nbnxm/pairlist_tuning.cpp index 26ec903918..97646f6faa 100644 --- a/src/gromacs/nbnxm/pairlist_tuning.cpp +++ b/src/gromacs/nbnxm/pairlist_tuning.cpp @@ -137,7 +137,7 @@ void increaseNstlist(FILE *fp, t_commrec *cr, float listfac_ok, listfac_max; int nstlist_orig, nstlist_prev; - real rlistWithReferenceNstlist, rlist_inc, rlist_ok, rlist_max; + real rlist_inc, rlist_ok, rlist_max; real rlist_new, rlist_prev; size_t nstlist_ind = 0; gmx_bool bBox, bDD, bCont; @@ -247,9 +247,9 @@ void increaseNstlist(FILE *fp, t_commrec *cr, */ nstlist_prev = ir->nstlist; ir->nstlist = nbnxnReferenceNstlist; - calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1, - -1, &listSetup, - nullptr, &rlistWithReferenceNstlist); + const real rlistWithReferenceNstlist = + calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, + -1, listSetup); ir->nstlist = nstlist_prev; /* Determine the pair list size increase due to zero interactions */ @@ -273,7 +273,8 @@ void increaseNstlist(FILE *fp, t_commrec *cr, } /* Set the pair-list buffer size in ir */ - calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1, -1, &listSetup, nullptr, &rlist_new); + rlist_new = + calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, -1, listSetup); /* Does rlist fit in the box? */ bBox = (gmx::square(rlist_new) < max_cutoff2(ir->ePBC, box)); @@ -405,10 +406,10 @@ setDynamicPairlistPruningParameters(const t_inputrec *ir, */ int listLifetime = tunedNstlistPrune - (useGpuList ? 0 : 1); listParams->nstlistPrune = tunedNstlistPrune; - calc_verlet_buffer_size(mtop, det(box), ir, - tunedNstlistPrune, listLifetime, - -1, &listSetup, nullptr, - &listParams->rlistInner); + listParams->rlistInner = + calcVerletBufferSize(*mtop, det(box), *ir, + tunedNstlistPrune, listLifetime, + -1, listSetup); /* On the GPU we apply the dynamic pruning in a rolling fashion * every c_nbnxnGpuRollingListPruningInterval steps, @@ -572,16 +573,17 @@ void setupDynamicPairlistPruning(const gmx::MDLogger &mdlog, } if (supportsDynamicPairlistGenerationInterval(*ir)) { - VerletbufListSetup listSetup1x1 = { 1, 1 }; - real rlistOuter; - real rlistInner; - calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1, - -1, &listSetup1x1, nullptr, &rlistOuter); + const VerletbufListSetup listSetup1x1 = { 1, 1 }; + const real rlistOuter = + calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, + -1, listSetup1x1); + real rlistInner = rlistOuter; if (listParams->useDynamicPruning) { int listLifeTime = listParams->nstlistPrune - (useGpuList ? 0 : 1); - calc_verlet_buffer_size(mtop, det(box), ir, listParams->nstlistPrune, listLifeTime, - -1, &listSetup1x1, nullptr, &rlistInner); + rlistInner = + calcVerletBufferSize(*mtop, det(box), *ir, listParams->nstlistPrune, listLifeTime, + -1, listSetup1x1); } mesg += gmx::formatString("At tolerance %g kJ/mol/ps per atom, equivalent classical 1x1 list would be:\n",