From ecf42f0e59517fe456d0225c13da523a5213b280 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Fri, 13 Mar 2020 22:49:32 +0100 Subject: [PATCH] Remove use of graph in do_force() Removed all unused code paths using t_graph. TODO: Remove the double coordinate copy in the update code. Change-Id: I9fcba521e85972e7a0d18e348e0330f845c45182 --- src/gromacs/domdec/mdsetup.cpp | 13 - src/gromacs/domdec/mdsetup.h | 5 +- src/gromacs/domdec/partition.cpp | 2 +- src/gromacs/gmxana/gmx_disre.cpp | 19 +- src/gromacs/listed_forces/bonded.cpp | 378 ++++-------------- src/gromacs/listed_forces/bonded.h | 78 ++-- src/gromacs/listed_forces/disre.cpp | 9 - src/gromacs/listed_forces/disre.h | 1 - src/gromacs/listed_forces/listed_forces.cpp | 33 +- src/gromacs/listed_forces/listed_forces.h | 3 - src/gromacs/listed_forces/orires.cpp | 9 - src/gromacs/listed_forces/orires.h | 1 - src/gromacs/listed_forces/pairs.cpp | 44 +- src/gromacs/listed_forces/pairs.h | 2 - src/gromacs/listed_forces/tests/bonded.cpp | 11 +- src/gromacs/mdlib/calcvir.cpp | 96 +---- src/gromacs/mdlib/calcvir.h | 5 +- src/gromacs/mdlib/force.cpp | 74 ++-- src/gromacs/mdlib/force.h | 41 +- src/gromacs/mdlib/sim_util.cpp | 38 +- .../mdlib/tests/leapfrogtestrunners.cpp | 4 +- src/gromacs/mdlib/update.cpp | 21 +- src/gromacs/mdlib/update.h | 2 - src/gromacs/mdlib/vsite.cpp | 271 ++++--------- src/gromacs/mdlib/vsite.h | 2 - src/gromacs/mdrun/md.cpp | 28 +- src/gromacs/mdrun/mimic.cpp | 27 +- src/gromacs/mdrun/minimize.cpp | 61 ++- src/gromacs/mdrun/rerun.cpp | 46 +-- src/gromacs/mdrun/shellfc.cpp | 31 +- src/gromacs/mdrun/shellfc.h | 2 - src/gromacs/mdrun/tpi.cpp | 9 +- src/gromacs/mdtypes/state.cpp | 2 +- src/gromacs/modularsimulator/domdechelper.cpp | 3 +- src/gromacs/modularsimulator/forceelement.cpp | 12 +- .../modularsimulator/shellfcelement.cpp | 10 +- .../modularsimulator/topologyholder.cpp | 6 +- src/gromacs/pbcutil/pbc.cpp | 15 +- src/gromacs/pbcutil/pbc.h | 8 +- 39 files changed, 367 insertions(+), 1055 deletions(-) diff --git a/src/gromacs/domdec/mdsetup.cpp b/src/gromacs/domdec/mdsetup.cpp index b5a0751197..bbe0cc67cf 100644 --- a/src/gromacs/domdec/mdsetup.cpp +++ b/src/gromacs/domdec/mdsetup.cpp @@ -48,7 +48,6 @@ #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/interaction_const.h" #include "gromacs/mdtypes/mdatom.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/topology/topology.h" @@ -65,7 +64,6 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, const gmx_mtop_t& top_global, gmx_localtop_t* top, t_forcerec* fr, - t_graph** graph, gmx::MDAtoms* mdAtoms, gmx::Constraints* constr, gmx_vsite_t* vsite, @@ -115,17 +113,6 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, } } - if (!usingDomDec && ir->pbcType != PbcType::No && !fr->bMolPBC) - { - GMX_ASSERT(graph != nullptr, "We use a graph with PBC (no periodic mols) and without DD"); - - *graph = mk_graph(nullptr, top->idef, 0, top_global.natoms, FALSE, FALSE); - } - else if (graph != nullptr) - { - *graph = nullptr; - } - /* Note that with DD only flexible constraints, not shells, are supported * and these don't require setup in make_local_shells(). * diff --git a/src/gromacs/domdec/mdsetup.h b/src/gromacs/domdec/mdsetup.h index 3752118133..2daa293526 100644 --- a/src/gromacs/domdec/mdsetup.h +++ b/src/gromacs/domdec/mdsetup.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2016,2017,2018,2019, by the GROMACS development team, led by + * Copyright (c) 2016,2017,2018,2019,2020, 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. @@ -50,7 +50,6 @@ struct gmx_shellfc_t; struct gmx_vsite_t; struct t_commrec; struct t_forcerec; -struct t_graph; struct t_inputrec; struct t_mdatoms; @@ -80,7 +79,6 @@ void make_local_shells(const t_commrec* cr, const t_mdatoms* md, gmx_shellfc_t* * \param[in] top_global The global topology * \param[in,out] top The local topology * \param[in,out] fr The force calculation parameter/data record - * \param[out] graph The molecular graph, can be NULL * \param[out] mdAtoms The MD atom data * \param[in,out] constr The constraints handler, can be NULL * \param[in,out] vsite The virtual site data, can be NULL @@ -91,7 +89,6 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, const gmx_mtop_t& top_global, gmx_localtop_t* top, t_forcerec* fr, - t_graph** graph, gmx::MDAtoms* mdAtoms, gmx::Constraints* constr, gmx_vsite_t* vsite, diff --git a/src/gromacs/domdec/partition.cpp b/src/gromacs/domdec/partition.cpp index e83408dfea..7a36e5ae25 100644 --- a/src/gromacs/domdec/partition.cpp +++ b/src/gromacs/domdec/partition.cpp @@ -3155,7 +3155,7 @@ void dd_partition_system(FILE* fplog, comm->atomRanges.end(DDAtomRanges::Type::Constraints), nat_f_novirsum); /* Update atom data for mdatoms and several algorithms */ - mdAlgorithmsSetupAtomData(cr, ir, top_global, top_local, fr, nullptr, mdAtoms, constr, vsite, nullptr); + mdAlgorithmsSetupAtomData(cr, ir, top_global, top_local, fr, mdAtoms, constr, vsite, nullptr); auto mdatoms = mdAtoms->mdatoms(); if (!thisRankHasDuty(cr, DUTY_PME)) diff --git a/src/gromacs/gmxana/gmx_disre.cpp b/src/gromacs/gmxana/gmx_disre.cpp index 397b4c1bfe..baee36ec07 100644 --- a/src/gromacs/gmxana/gmx_disre.cpp +++ b/src/gromacs/gmxana/gmx_disre.cpp @@ -65,7 +65,6 @@ #include "gromacs/mdtypes/md_enums.h" #include "gromacs/mdtypes/mdatom.h" #include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pbcutil/rmpbc.h" #include "gromacs/topology/index.h" @@ -153,7 +152,6 @@ static void check_viol(FILE* log, rvec x[], rvec4 f[], t_pbc* pbc, - t_graph* g, t_dr_result dr[], int clust_id, int isize, @@ -234,7 +232,7 @@ static void check_viol(FILE* log, dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label]; snew(fshift, SHIFTS); - ta_disres(n, &forceatoms[i], forceparams.data(), x, f, fshift, pbc, g, lam, &dvdl, nullptr, + ta_disres(n, &forceatoms[i], forceparams.data(), x, f, fshift, pbc, lam, &dvdl, nullptr, fcd, nullptr); sfree(fshift); viol = fcd->disres.sumviol; @@ -701,7 +699,6 @@ int gmx_disre(int argc, char* argv[]) FILE * out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr; t_fcdata fcd; - t_graph* g; int i, j, kkk; t_trxstatus* status; real t; @@ -774,18 +771,10 @@ int gmx_disre(int argc, char* argv[]) gmx_localtop_t top(topInfo.mtop()->ffparams); gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO); - g = nullptr; pbc_null = nullptr; if (ir->pbcType != PbcType::No) { - if (ir->bPeriodicMols) - { - pbc_null = &pbc; - } - else - { - g = mk_graph(fplog, top.idef, 0, ntopatoms, FALSE, FALSE); - } + pbc_null = &pbc; } if (ftp2bSet(efNDX, NFILE, fnm)) @@ -867,12 +856,12 @@ int gmx_disre(int argc, char* argv[]) } my_clust = clust->inv_clust[j]; range_check(my_clust, 0, clust->clust->nr); - check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, g, dr_clust, + check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, dr_clust, my_clust, isize, index, vvindex, &fcd); } else { - check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, g, &dr, 0, + check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, &dr, 0, isize, index, vvindex, &fcd); } if (bPDB) diff --git a/src/gromacs/listed_forces/bonded.cpp b/src/gromacs/listed_forces/bonded.cpp index d4867af22e..61f6ee32f0 100644 --- a/src/gromacs/listed_forces/bonded.cpp +++ b/src/gromacs/listed_forces/bonded.cpp @@ -65,7 +65,6 @@ #include "gromacs/mdtypes/fcdata.h" #include "gromacs/mdtypes/mdatom.h" #include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pbcutil/pbc_simd.h" #include "gromacs/simd/simd.h" @@ -93,7 +92,6 @@ using BondedFunction = real (*)(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms* md, @@ -205,26 +203,17 @@ namespace /*! \brief Spreads and accumulates the bonded forces to the two atoms and adds the virial contribution when needed * - * When \p g==nullptr, \p shiftIndex is used as the periodic shift. - * When \p g!=nullptr, the graph is used to compute the periodic shift. + * \p shiftIndex is used as the periodic shift. */ template -inline void spreadBondForces(const real bondForce, - const rvec dx, - const int ai, - const int aj, - rvec4* f, - int shiftIndex, - const t_graph* g, - rvec* fshift) +inline void spreadBondForces(const real bondForce, + const rvec dx, + const int ai, + const int aj, + rvec4* f, + int shiftIndex, + rvec* fshift) { - if (computeVirial(flavor) && g) - { - ivec dt; - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt); - shiftIndex = IVEC2IS(dt); - } - for (int m = 0; m < DIM; m++) /* 15 */ { const real fij = bondForce * dx[m]; @@ -257,7 +246,6 @@ real morse_bonds(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -312,8 +300,8 @@ real morse_bonds(int nbonds, *dvdlambda += (cbB - cbA) * omtemp * omtemp - (2 - 2 * omtemp) * omtemp * cb * ((b0B - b0A) * be - (beB - beA) * (dr - b0)); /* 15 */ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 83 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 83 TOTAL */ return vtot; } @@ -326,7 +314,6 @@ real cubic_bonds(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -369,8 +356,8 @@ real cubic_bonds(int nbonds, vtot += vbond; /* 21 */ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 54 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 54 TOTAL */ return vtot; } @@ -382,7 +369,6 @@ real FENE_bonds(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -429,8 +415,8 @@ real FENE_bonds(int nbonds, vtot += vbond; /* 35 */ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 58 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 58 TOTAL */ return vtot; } @@ -468,7 +454,6 @@ real bonds(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -503,8 +488,8 @@ real bonds(int nbonds, vtot += vbond; /* 1*/ fbond *= gmx::invsqrt(dr2); /* 6 */ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 59 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 59 TOTAL */ return vtot; } @@ -516,7 +501,6 @@ real restraint_bonds(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -591,8 +575,8 @@ real restraint_bonds(int nbonds, vtot += vbond; /* 1*/ fbond *= gmx::invsqrt(dr2); /* 6 */ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 59 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 59 TOTAL */ return vtot; } @@ -605,7 +589,6 @@ real polarize(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms* md, @@ -638,8 +621,8 @@ real polarize(int nbonds, vtot += vbond; /* 1*/ fbond *= gmx::invsqrt(dr2); /* 6 */ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 59 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 59 TOTAL */ return vtot; } @@ -651,7 +634,6 @@ real anharm_polarize(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms* md, @@ -693,8 +675,8 @@ real anharm_polarize(int nbonds, fbond *= gmx::invsqrt(dr2); /* 6 */ vtot += vbond; /* 1*/ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 72 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 72 TOTAL */ return vtot; } @@ -706,7 +688,6 @@ real water_pol(int nbonds, rvec4 f[], rvec gmx_unused fshift[], const t_pbc gmx_unused* pbc, - const t_graph gmx_unused* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -718,7 +699,6 @@ real water_pol(int nbonds, * three spatial dimensions in the molecular frame. */ int i, m, aO, aH1, aH2, aD, aS, type, type0, ki; - ivec dt; rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj; real vtot, fij, r_HH, r_OD, r_nW, tx, ty, tz, qS; @@ -790,12 +770,6 @@ real water_pol(int nbonds, kdx[ZZ] = kk[ZZ] * dx[ZZ]; vtot += iprod(dx, kdx); - if (computeVirial(flavor) && g) - { - ivec_sub(SHIFT_IVEC(g, aS), SHIFT_IVEC(g, aD), dt); - ki = IVEC2IS(dt); - } - for (m = 0; (m < DIM); m++) { /* This is a tensor operation but written out for speed */ @@ -858,7 +832,6 @@ real thole_pol(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph gmx_unused* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms* md, @@ -910,7 +883,6 @@ angles(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -920,7 +892,6 @@ angles(int nbonds, int i, ai, aj, ak, t1, t2, type; rvec r_ij, r_kj; real cos_theta, cos_theta2, theta, dVdt, va, vtot; - ivec jt, dt_ij, dt_kj; vtot = 0.0; for (i = 0; i < nbonds;) @@ -970,15 +941,6 @@ angles(int nbonds, } if (computeVirial(flavor)) { - if (g != nullptr) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); @@ -1008,7 +970,6 @@ angles(int nbonds, rvec4 f[], rvec gmx_unused fshift[], const t_pbc* pbc, - const t_graph gmx_unused* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -1156,7 +1117,6 @@ real linear_angles(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -1166,7 +1126,6 @@ real linear_angles(int nbonds, int i, m, ai, aj, ak, t1, t2, type; rvec f_i, f_j, f_k; real L1, kA, kB, aA, aB, dr, dr2, va, vtot, a, b, klin; - ivec jt, dt_ij, dt_kj; rvec r_ij, r_kj, r_ik, dx; L1 = 1 - lambda; @@ -1211,15 +1170,6 @@ real linear_angles(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); @@ -1237,7 +1187,6 @@ urey_bradley(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -1249,7 +1198,6 @@ urey_bradley(int nbonds, real cos_theta, cos_theta2, theta; real dVdt, va, vtot, dr, dr2, vbond, fbond, fik; real kthA, th0A, kUBA, r13A, kthB, th0B, kUBB, r13B; - ivec jt, dt_ij, dt_kj, dt_ik; vtot = 0.0; for (i = 0; (i < nbonds);) @@ -1306,15 +1254,6 @@ urey_bradley(int nbonds, } if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); @@ -1329,11 +1268,6 @@ urey_bradley(int nbonds, vtot += vbond; /* 1*/ fbond *= gmx::invsqrt(dr2); /* 6 */ - if (computeVirial(flavor) && g) - { - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, ak), dt_ik); - ki = IVEC2IS(dt_ik); - } for (m = 0; (m < DIM); m++) /* 15 */ { fik = fbond * r_ik[m]; @@ -1363,7 +1297,6 @@ urey_bradley(int nbonds, rvec4 f[], rvec gmx_unused fshift[], const t_pbc* pbc, - const t_graph gmx_unused* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -1506,7 +1439,6 @@ real quartic_angles(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -1516,7 +1448,6 @@ real quartic_angles(int nbonds, int i, j, ai, aj, ak, t1, t2, type; rvec r_ij, r_kj; real cos_theta, cos_theta2, theta, dt, dVdt, va, dtp, c, vtot; - ivec jt, dt_ij, dt_kj; vtot = 0.0; for (i = 0; (i < nbonds);) @@ -1574,15 +1505,6 @@ real quartic_angles(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); @@ -1714,31 +1636,29 @@ inline void dih_angle_simd(const rvec* x, } // namespace template -void do_dih_fup(int i, - int j, - int k, - int l, - real ddphi, - rvec r_ij, - rvec r_kj, - rvec r_kl, - rvec m, - rvec n, - rvec4 f[], - rvec fshift[], - const t_pbc* pbc, - const t_graph* g, - const rvec x[], - int t1, - int t2, - int t3) +void do_dih_fup(int i, + int j, + int k, + int l, + real ddphi, + rvec r_ij, + rvec r_kj, + rvec r_kl, + rvec m, + rvec n, + rvec4 f[], + rvec fshift[], + const t_pbc* pbc, + const rvec x[], + int t1, + int t2, + int t3) { /* 143 FLOPS */ rvec f_i, f_j, f_k, f_l; rvec uvec, vvec, svec, dx_jl; real iprm, iprn, nrkj, nrkj2, nrkj_1, nrkj_2; real a, b, p, q, toler; - ivec jt, dt_ij, dt_kj, dt_lj; iprm = iprod(m, m); /* 5 */ iprn = iprod(n, n); /* 5 */ @@ -1769,17 +1689,7 @@ void do_dih_fup(int i, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, j), jt); - ivec_sub(SHIFT_IVEC(g, i), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, k), jt, dt_kj); - ivec_sub(SHIFT_IVEC(g, l), jt, dt_lj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - t3 = IVEC2IS(dt_lj); - } - else if (pbc) + if (pbc) { t3 = pbc_rvec_sub(pbc, x[l], x[j], dx_jl); } @@ -1894,7 +1804,6 @@ pdihs(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -1931,9 +1840,9 @@ pdihs(int nbonds, } while (i < nbonds && forceatoms[i + 1] == ai && forceatoms[i + 2] == aj && forceatoms[i + 3] == ak && forceatoms[i + 4] == al); - do_dih_fup(ai, aj, ak, al, ddphi_tot, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, g, x, - t1, t2, t3); /* 112 */ - } /* 223 TOTAL */ + do_dih_fup(ai, aj, ak, al, ddphi_tot, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1, + t2, t3); /* 112 */ + } /* 223 TOTAL */ return vtot; } @@ -1950,7 +1859,6 @@ pdihs(int nbonds, rvec4 f[], rvec gmx_unused fshift[], const t_pbc* pbc, - const t_graph gmx_unused* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -2065,7 +1973,6 @@ rbdihs(int nbonds, rvec4 f[], rvec gmx_unused fshift[], const t_pbc* pbc, - const t_graph gmx_unused* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -2192,7 +2099,6 @@ real idihs(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -2245,7 +2151,7 @@ real idihs(int nbonds, dvdl_term += 0.5 * (kB - kA) * dp2 - kk * dphi0 * dp; - do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, g, x, t1, + do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1, t2, t3); /* 112 */ /* 218 TOTAL */ } @@ -2263,7 +2169,6 @@ real low_angres(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, gmx_bool bZAxis) @@ -2274,7 +2179,6 @@ real low_angres(int nbonds, rvec r_ij, r_kl, f_i, f_k = { 0, 0, 0 }; real st, sth, nrij2, nrkl2, c, cij, ckl; - ivec dt; t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */ vtot = 0.0; @@ -2334,20 +2238,10 @@ real low_angres(int nbonds, if (computeVirial(flavor)) { - if (g) - { - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt); - t1 = IVEC2IS(dt); - } rvec_inc(fshift[t1], f_i); rvec_dec(fshift[CENTRAL], f_i); if (!bZAxis) { - if (g) - { - ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, al), dt); - t2 = IVEC2IS(dt); - } rvec_inc(fshift[t2], f_k); rvec_dec(fshift[CENTRAL], f_k); } @@ -2366,15 +2260,13 @@ real angres(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, t_fcdata gmx_unused* fcd, int gmx_unused* global_atom_index) { - return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, g, lambda, - dvdlambda, FALSE); + return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, lambda, dvdlambda, FALSE); } template @@ -2385,15 +2277,13 @@ real angresz(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, t_fcdata gmx_unused* fcd, int gmx_unused* global_atom_index) { - return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, g, lambda, - dvdlambda, TRUE); + return low_angres(nbonds, forceatoms, forceparams, x, f, fshift, pbc, lambda, dvdlambda, TRUE); } template @@ -2404,7 +2294,6 @@ real dihres(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -2483,8 +2372,8 @@ real dihres(int nbonds, { *dvdlambda += kfac * ddp * ((dphiB - dphiA) - (phi0B - phi0A)); } - do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, g, x, - t1, t2, t3); /* 112 */ + do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1, + t2, t3); /* 112 */ } } return vtot; @@ -2498,7 +2387,6 @@ real unimplemented(int gmx_unused nbonds, rvec4 gmx_unused f[], rvec gmx_unused fshift[], const t_pbc gmx_unused* pbc, - const t_graph gmx_unused* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -2516,7 +2404,6 @@ real restrangles(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -2526,7 +2413,6 @@ real restrangles(int nbonds, int i, d, ai, aj, ak, type, m; int t1, t2; real v, vtot; - ivec jt, dt_ij, dt_kj; rvec f_i, f_j, f_k; double prefactor, ratio_ante, ratio_post; rvec delta_ante, delta_post, vec_temp; @@ -2551,7 +2437,7 @@ real restrangles(int nbonds, real restrangles(int nbonds, const t_iatom forceatoms[],const t_iparams forceparams[], const rvec x[],rvec4 f[],rvec fshift[], - const t_pbc *pbc,const t_graph *g, + const t_pbc *pbc, real gmx_unused lambda,real gmx_unused *dvdlambda, const t_mdatoms gmx_unused *md,t_fcdata gmx_unused *fcd, int gmx_unused *global_atom_index) @@ -2560,7 +2446,6 @@ real restrangles(int nbonds, int t1, t2; rvec r_ij,r_kj; real v, vtot; - ivec jt,dt_ij,dt_kj; rvec f_i, f_j, f_k; real prefactor, ratio_ante, ratio_post; rvec delta_ante, delta_post, vec_temp; @@ -2604,15 +2489,6 @@ real restrangles(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } - rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); @@ -2630,7 +2506,6 @@ real restrdihs(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvlambda, const t_mdatoms gmx_unused* md, @@ -2640,7 +2515,6 @@ real restrdihs(int nbonds, int i, d, type, ai, aj, ak, al; rvec f_i, f_j, f_k, f_l; rvec dx_jl; - ivec jt, dt_ij, dt_kj, dt_lj; int t1, t2, t3; real v, vtot; rvec delta_ante, delta_crnt, delta_post, vec_temp; @@ -2712,17 +2586,7 @@ real restrdihs(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - ivec_sub(SHIFT_IVEC(g, al), jt, dt_lj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - t3 = IVEC2IS(dt_lj); - } - else if (pbc) + if (pbc) { t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl); } @@ -2750,7 +2614,6 @@ real cbtdihs(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -2762,7 +2625,6 @@ real cbtdihs(int nbonds, real v, vtot; rvec vec_temp; rvec f_i, f_j, f_k, f_l; - ivec jt, dt_ij, dt_kj, dt_lj; rvec dx_jl; rvec delta_ante, delta_crnt, delta_post; rvec f_phi_ai, f_phi_aj, f_phi_ak, f_phi_al; @@ -2826,17 +2688,7 @@ real cbtdihs(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - ivec_sub(SHIFT_IVEC(g, al), jt, dt_lj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - t3 = IVEC2IS(dt_lj); - } - else if (pbc) + if (pbc) { t3 = pbc_rvec_sub(pbc, x[al], x[aj], dx_jl); } @@ -2864,7 +2716,6 @@ rbdihs(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -2956,8 +2807,8 @@ rbdihs(int nbonds, ddphi = -ddphi * sin_phi; /* 11 */ - do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, g, x, t1, - t2, t3); /* 112 */ + do_dih_fup(ai, aj, ak, al, ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1, t2, + t3); /* 112 */ vtot += v; } *dvdlambda += dvdl_term; @@ -3008,15 +2859,14 @@ int cmap_setup_grid_index(int ip, int grid_spacing, int* ipm1, int* ipp1, int* i } // namespace -real cmap_dihs(int nbonds, - const t_iatom forceatoms[], - const t_iparams forceparams[], - const gmx_cmap_t* cmap_grid, - const rvec x[], - rvec4 f[], - rvec fshift[], - const struct t_pbc* pbc, - const struct t_graph* g, +real cmap_dihs(int nbonds, + const t_iatom forceatoms[], + const t_iparams forceparams[], + const gmx_cmap_t* cmap_grid, + const rvec x[], + rvec4 f[], + rvec fshift[], + const struct t_pbc* pbc, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -3050,8 +2900,6 @@ real cmap_dihs(int nbonds, rvec a1, b1, a2, b2; rvec f1, g1, h1, f2, g2, h2; rvec dtf1, dtg1, dth1, dtf2, dtg2, dth2; - ivec jt1, dt1_ij, dt1_kj, dt1_lj; - ivec jt2, dt2_ij, dt2_kj, dt2_lj; int loop_index[4][4] = { { 0, 4, 8, 12 }, { 1, 5, 9, 13 }, { 2, 6, 10, 14 }, { 3, 7, 11, 15 } }; @@ -3356,25 +3204,7 @@ real cmap_dihs(int nbonds, /* Shift forces */ if (fshift != nullptr) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, a1j), jt1); - ivec_sub(SHIFT_IVEC(g, a1i), jt1, dt1_ij); - ivec_sub(SHIFT_IVEC(g, a1k), jt1, dt1_kj); - ivec_sub(SHIFT_IVEC(g, a1l), jt1, dt1_lj); - t11 = IVEC2IS(dt1_ij); - t21 = IVEC2IS(dt1_kj); - t31 = IVEC2IS(dt1_lj); - - copy_ivec(SHIFT_IVEC(g, a2j), jt2); - ivec_sub(SHIFT_IVEC(g, a2i), jt2, dt2_ij); - ivec_sub(SHIFT_IVEC(g, a2k), jt2, dt2_kj); - ivec_sub(SHIFT_IVEC(g, a2l), jt2, dt2_lj); - t12 = IVEC2IS(dt2_ij); - t22 = IVEC2IS(dt2_kj); - t32 = IVEC2IS(dt2_lj); - } - else if (pbc) + if (pbc) { t31 = pbc_rvec_sub(pbc, x[a1l], x[a1j], h1); t32 = pbc_rvec_sub(pbc, x[a2l], x[a2j], h2); @@ -3441,7 +3271,6 @@ real g96bonds(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -3468,8 +3297,8 @@ real g96bonds(int nbonds, vtot += 0.5 * vbond; /* 1*/ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 44 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 44 TOTAL */ return vtot; } @@ -3494,7 +3323,6 @@ real g96angles(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -3506,7 +3334,6 @@ real g96angles(int nbonds, real cos_theta, dVdt, va, vtot; real rij_1, rij_2, rkj_1, rkj_2, rijrkj_1; rvec f_i, f_j, f_k; - ivec jt, dt_ij, dt_kj; vtot = 0.0; for (i = 0; (i < nbonds);) @@ -3541,15 +3368,6 @@ real g96angles(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); /* 9 */ @@ -3567,7 +3385,6 @@ real cross_bond_bond(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -3581,7 +3398,6 @@ real cross_bond_bond(int nbonds, rvec r_ij, r_kj; real vtot, vrr, s1, s2, r1, r2, r1e, r2e, krr; rvec f_i, f_j, f_k; - ivec jt, dt_ij, dt_kj; vtot = 0.0; for (i = 0; (i < nbonds);) @@ -3624,15 +3440,6 @@ real cross_bond_bond(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); /* 9 */ @@ -3650,7 +3457,6 @@ real cross_bond_angle(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -3664,7 +3470,6 @@ real cross_bond_angle(int nbonds, rvec r_ij, r_kj, r_ik; real vtot, vrt, s1, s2, s3, r1, r2, r3, r1e, r2e, r3e, krt, k1, k2, k3; rvec f_i, f_j, f_k; - ivec jt, dt_ij, dt_kj; vtot = 0.0; for (i = 0; (i < nbonds);) @@ -3717,15 +3522,6 @@ real cross_bond_angle(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); /* 9 */ @@ -3792,7 +3588,6 @@ real tab_bonds(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -3828,8 +3623,8 @@ real tab_bonds(int nbonds, vtot += vbond; /* 1*/ fbond *= gmx::invsqrt(dr2); /* 6 */ - spreadBondForces(fbond, dx, ai, aj, f, ki, g, fshift); /* 15 */ - } /* 62 TOTAL */ + spreadBondForces(fbond, dx, ai, aj, f, ki, fshift); /* 15 */ + } /* 62 TOTAL */ return vtot; } @@ -3841,7 +3636,6 @@ real tab_angles(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -3851,7 +3645,6 @@ real tab_angles(int nbonds, int i, ai, aj, ak, t1, t2, type, table; rvec r_ij, r_kj; real cos_theta, cos_theta2, theta, dVdt, va, vtot; - ivec jt, dt_ij, dt_kj; vtot = 0.0; for (i = 0; (i < nbonds);) @@ -3899,15 +3692,6 @@ real tab_angles(int nbonds, if (computeVirial(flavor)) { - if (g) - { - copy_ivec(SHIFT_IVEC(g, aj), jt); - - ivec_sub(SHIFT_IVEC(g, ai), jt, dt_ij); - ivec_sub(SHIFT_IVEC(g, ak), jt, dt_kj); - t1 = IVEC2IS(dt_ij); - t2 = IVEC2IS(dt_kj); - } rvec_inc(fshift[t1], f_i); rvec_inc(fshift[CENTRAL], f_j); rvec_inc(fshift[t2], f_k); @@ -3925,7 +3709,6 @@ real tab_dihs(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms gmx_unused* md, @@ -3955,7 +3738,7 @@ real tab_dihs(int nbonds, forceparams[type].tab.kB, phi + M_PI, lambda, &vpd, &ddphi); vtot += vpd; - do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, g, x, t1, + do_dih_fup(ai, aj, ak, al, -ddphi, r_ij, r_kj, r_kl, m, n, f, fshift, pbc, x, t1, t2, t3); /* 112 */ } /* 227 TOTAL */ @@ -4080,25 +3863,24 @@ const gmx::EnumerationArray iatoms = gmx::makeConstArrayRef(ilist.iatoms); - v = calc_one_bond( - thread, ftype, idef, iatoms, idef.numNonperturbedInteractions[ftype], - fr->bondedThreading->workDivision, x, ft, fshift, fr, pbc_null, g, grpp, - nrnb, lambda, dvdlt, md, fcd, stepWork, global_atom_index); + v = calc_one_bond(thread, ftype, idef, iatoms, idef.numNonperturbedInteractions[ftype], + fr->bondedThreading->workDivision, x, ft, fshift, fr, pbc_null, + grpp, nrnb, lambda, dvdlt, md, fcd, stepWork, global_atom_index); epot[ftype] += v; } } @@ -493,7 +490,6 @@ void calc_listed(const t_commrec* cr, const t_forcerec* fr, const struct t_pbc* pbc, const struct t_pbc* pbc_full, - const struct t_graph* g, gmx_enerdata_t* enerd, t_nrnb* nrnb, const real* lambda, @@ -558,9 +554,8 @@ void calc_listed(const t_commrec* cr, /* The dummy array is to have a place to store the dhdl at other values of lambda, which will be thrown away in the end */ real dvdl[efptNR] = { 0 }; - calcBondedForces(idef, x, fr, pbc_null, g, - as_rvec_array(forceWithShiftForces.shiftForces().data()), enerd, nrnb, - lambda, dvdl, md, fcd, stepWork, global_atom_index); + calcBondedForces(idef, x, fr, pbc_null, as_rvec_array(forceWithShiftForces.shiftForces().data()), + enerd, nrnb, lambda, dvdl, md, fcd, stepWork, global_atom_index); wallcycle_sub_stop(wcycle, ewcsLISTED); wallcycle_sub_start(wcycle, ewcsLISTED_BUF_OPS); @@ -593,7 +588,6 @@ void calc_listed_lambda(const InteractionDefinitions& idef, const rvec x[], const t_forcerec* fr, const struct t_pbc* pbc, - const struct t_graph* g, gmx_grppairener_t* grpp, real* epot, gmx::ArrayRef dvdl, @@ -642,7 +636,7 @@ void calc_listed_lambda(const InteractionDefinitions& idef, gmx::StepWorkload tempFlags; tempFlags.computeEnergy = true; v = calc_one_bond(0, ftype, idef, iatomsPerturbed, iatomsPerturbed.ssize(), - workDivision, x, f, fshift, fr, pbc_null, g, grpp, nrnb, lambda, + workDivision, x, f, fshift, fr, pbc_null, grpp, nrnb, lambda, dvdl.data(), md, fcd, tempFlags, global_atom_index); epot[ftype] += v; } @@ -667,7 +661,6 @@ void do_force_listed(struct gmx_wallcycle* wcycle, gmx::ForceOutputs* forceOutputs, const t_forcerec* fr, const struct t_pbc* pbc, - const struct t_graph* graph, gmx_enerdata_t* enerd, t_nrnb* nrnb, const real* lambda, @@ -689,7 +682,7 @@ void do_force_listed(struct gmx_wallcycle* wcycle, set_pbc(&pbc_full, fr->pbcType, box); } calc_listed(cr, ms, wcycle, idef, x, xWholeMolecules, hist, forceOutputs, fr, pbc, &pbc_full, - graph, enerd, nrnb, lambda, md, fcd, global_atom_index, stepWork); + enerd, nrnb, lambda, md, fcd, global_atom_index, stepWork); /* Check if we have to determine energy differences * at foreign lambda's. @@ -715,8 +708,8 @@ void do_force_listed(struct gmx_wallcycle* wcycle, { lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i - 1]); } - calc_listed_lambda(idef, x, fr, pbc, graph, &(enerd->foreign_grpp), - enerd->foreign_term, dvdl, nrnb, lam_i, md, fcd, global_atom_index); + calc_listed_lambda(idef, x, fr, pbc, &(enerd->foreign_grpp), enerd->foreign_term, + dvdl, nrnb, lam_i, md, fcd, global_atom_index); sum_epot(&(enerd->foreign_grpp), enerd->foreign_term); enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT]; for (int j = 0; j < efptNR; j++) diff --git a/src/gromacs/listed_forces/listed_forces.h b/src/gromacs/listed_forces/listed_forces.h index c25e9973c1..67c078430b 100644 --- a/src/gromacs/listed_forces/listed_forces.h +++ b/src/gromacs/listed_forces/listed_forces.h @@ -79,7 +79,6 @@ class InteractionDefinitions; struct t_commrec; struct t_fcdata; struct t_forcerec; -struct t_graph; struct t_lambda; struct t_mdatoms; struct t_nrnb; @@ -101,7 +100,6 @@ using BondedFunction = real (*)(int nbonds, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms* md, @@ -129,7 +127,6 @@ void do_force_listed(struct gmx_wallcycle* wcycle, gmx::ForceOutputs* forceOutputs, const t_forcerec* fr, const struct t_pbc* pbc, - const struct t_graph* graph, gmx_enerdata_t* enerd, t_nrnb* nrnb, const real* lambda, diff --git a/src/gromacs/listed_forces/orires.cpp b/src/gromacs/listed_forces/orires.cpp index fe375477ab..8e257559f4 100644 --- a/src/gromacs/listed_forces/orires.cpp +++ b/src/gromacs/listed_forces/orires.cpp @@ -54,7 +54,6 @@ #include "gromacs/mdtypes/mdatom.h" #include "gromacs/mdtypes/state.h" #include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/mtop_util.h" @@ -653,7 +652,6 @@ real orires(int nfa, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real gmx_unused lambda, real gmx_unused* dvdlambda, const t_mdatoms gmx_unused* md, @@ -661,7 +659,6 @@ real orires(int nfa, int gmx_unused* global_atom_index) { int ex, power, ki = CENTRAL; - ivec dt; real r2, invr, invr2, fc, smooth_fc, dev, devins, pfac; rvec r, Sr, fij; real vtot; @@ -738,12 +735,6 @@ real orires(int nfa, fij[i] = -pfac * dev * (4 * Sr[i] - 2 * (2 + power) * invr2 * iprod(Sr, r) * r[i]); } - if (g) - { - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt); - ki = IVEC2IS(dt); - } - for (int i = 0; i < DIM; i++) { f[ai][i] += fij[i]; diff --git a/src/gromacs/listed_forces/orires.h b/src/gromacs/listed_forces/orires.h index f29557986b..2891ec2dbb 100644 --- a/src/gromacs/listed_forces/orires.h +++ b/src/gromacs/listed_forces/orires.h @@ -116,7 +116,6 @@ real orires(int nfa, rvec4 f[], rvec fshift[], const t_pbc* pbc, - const t_graph* g, real lambda, real* dvdlambda, const t_mdatoms* md, diff --git a/src/gromacs/listed_forces/pairs.cpp b/src/gromacs/listed_forces/pairs.cpp index 489925b962..1f758598b6 100644 --- a/src/gromacs/listed_forces/pairs.cpp +++ b/src/gromacs/listed_forces/pairs.cpp @@ -59,7 +59,6 @@ #include "gromacs/mdtypes/nblist.h" #include "gromacs/mdtypes/simulation_workload.h" #include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pbcutil/pbc_simd.h" #include "gromacs/simd/simd.h" @@ -356,25 +355,23 @@ static real free_energy_evaluate_single(real r2, /*! \brief Calculate pair interactions, supports all types and conditions. */ template -static real do_pairs_general(int ftype, - int nbonds, - const t_iatom iatoms[], - const t_iparams iparams[], - const rvec x[], - rvec4 f[], - rvec fshift[], - const struct t_pbc* pbc, - const struct t_graph* g, - const real* lambda, - real* dvdl, - const t_mdatoms* md, - const t_forcerec* fr, - gmx_grppairener_t* grppener, - int* global_atom_index) +static real do_pairs_general(int ftype, + int nbonds, + const t_iatom iatoms[], + const t_iparams iparams[], + const rvec x[], + rvec4 f[], + rvec fshift[], + const struct t_pbc* pbc, + const real* lambda, + real* dvdl, + const t_mdatoms* md, + const t_forcerec* fr, + gmx_grppairener_t* grppener, + int* global_atom_index) { real qq, c6, c12; rvec dx; - ivec dt; int i, itype, ai, aj, gid; int fshift_index; real r2; @@ -546,12 +543,6 @@ static real do_pairs_general(int ftype, if (computeVirial(flavor)) { - if (g) - { - /* Correct the shift forces using the graph */ - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt); - fshift_index = IVEC2IS(dt); - } if (fshift_index != CENTRAL) { rvec_inc(fshift[fshift_index], dx); @@ -678,7 +669,6 @@ void do_pairs(int ftype, rvec4 f[], rvec fshift[], const struct t_pbc* pbc, - const struct t_graph* g, const real* lambda, real* dvdl, const t_mdatoms* md, @@ -732,13 +722,13 @@ void do_pairs(int ftype, else if (stepWork.computeVirial) { do_pairs_general( - ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, g, lambda, dvdl, md, fr, - grppener, global_atom_index); + ftype, nbonds, iatoms, iparams, x, f, fshift, pbc, lambda, dvdl, md, fr, grppener, + global_atom_index); } else { do_pairs_general(ftype, nbonds, iatoms, iparams, x, f, - fshift, pbc, g, lambda, dvdl, md, fr, + fshift, pbc, lambda, dvdl, md, fr, grppener, global_atom_index); } } diff --git a/src/gromacs/listed_forces/pairs.h b/src/gromacs/listed_forces/pairs.h index e36ac6cd87..af9320e014 100644 --- a/src/gromacs/listed_forces/pairs.h +++ b/src/gromacs/listed_forces/pairs.h @@ -52,7 +52,6 @@ struct gmx_grppairener_t; struct t_forcerec; -struct t_graph; struct t_pbc; namespace gmx @@ -73,7 +72,6 @@ void do_pairs(int ftype, rvec4 f[], rvec fshift[], const struct t_pbc* pbc, - const struct t_graph* g, const real* lambda, real* dvdl, const t_mdatoms* md, diff --git a/src/gromacs/listed_forces/tests/bonded.cpp b/src/gromacs/listed_forces/tests/bonded.cpp index c67f0c6563..f1daf94f52 100644 --- a/src/gromacs/listed_forces/tests/bonded.cpp +++ b/src/gromacs/listed_forces/tests/bonded.cpp @@ -546,12 +546,11 @@ protected: t_mdatoms mdatoms = { 0 }; mdatoms.chargeA = chargeA.data(); OutputQuantities output; - output.energy = calculateSimpleBond( - input_.ftype, iatoms.size(), iatoms.data(), &input_.iparams, - as_rvec_array(x_.data()), output.f, output.fshift, &pbc_, - /* const struct t_graph *g */ nullptr, lambda, &output.dvdlambda, &mdatoms, - /* struct t_fcdata * */ nullptr, ddgatindex.data(), - BondedKernelFlavor::ForcesAndVirialAndEnergy); + output.energy = calculateSimpleBond(input_.ftype, iatoms.size(), iatoms.data(), + &input_.iparams, as_rvec_array(x_.data()), output.f, + output.fshift, &pbc_, lambda, &output.dvdlambda, &mdatoms, + /* struct t_fcdata * */ nullptr, ddgatindex.data(), + BondedKernelFlavor::ForcesAndVirialAndEnergy); // Internal consistency test of both test input // and bonded functions. EXPECT_TRUE((input_.fep || (output.dvdlambda == 0.0))) << "dvdlambda was " << output.dvdlambda; diff --git a/src/gromacs/mdlib/calcvir.cpp b/src/gromacs/mdlib/calcvir.cpp index 0c35b7c918..99398cff84 100644 --- a/src/gromacs/mdlib/calcvir.cpp +++ b/src/gromacs/mdlib/calcvir.cpp @@ -48,7 +48,6 @@ #include "gromacs/math/vectypes.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/utility/gmxassert.h" @@ -143,98 +142,7 @@ void calc_vir(int nxf, const rvec x[], const rvec f[], tensor vir, bool bScrewPB } } - -static void -lo_fcv(int i0, int i1, const real x[], const real f[], tensor vir, const int is[], const real box[], gmx_bool bTriclinic) -{ - int i, i3, tx, ty, tz; - real xx, yy, zz; - real dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0; - - if (bTriclinic) - { - for (i = i0; (i < i1); i++) - { - i3 = DIM * i; - tx = is[i3 + XX]; - ty = is[i3 + YY]; - tz = is[i3 + ZZ]; - - xx = x[i3 + XX] - tx * box[XXXX] - ty * box[YYXX] - tz * box[ZZXX]; - dvxx += xx * f[i3 + XX]; - dvxy += xx * f[i3 + YY]; - dvxz += xx * f[i3 + ZZ]; - - yy = x[i3 + YY] - ty * box[YYYY] - tz * box[ZZYY]; - dvyx += yy * f[i3 + XX]; - dvyy += yy * f[i3 + YY]; - dvyz += yy * f[i3 + ZZ]; - - zz = x[i3 + ZZ] - tz * box[ZZZZ]; - dvzx += zz * f[i3 + XX]; - dvzy += zz * f[i3 + YY]; - dvzz += zz * f[i3 + ZZ]; - } - } - else - { - for (i = i0; (i < i1); i++) - { - i3 = DIM * i; - tx = is[i3 + XX]; - ty = is[i3 + YY]; - tz = is[i3 + ZZ]; - - xx = x[i3 + XX] - tx * box[XXXX]; - dvxx += xx * f[i3 + XX]; - dvxy += xx * f[i3 + YY]; - dvxz += xx * f[i3 + ZZ]; - - yy = x[i3 + YY] - ty * box[YYYY]; - dvyx += yy * f[i3 + XX]; - dvyy += yy * f[i3 + YY]; - dvyz += yy * f[i3 + ZZ]; - - zz = x[i3 + ZZ] - tz * box[ZZZZ]; - dvzx += zz * f[i3 + XX]; - dvzy += zz * f[i3 + YY]; - dvzz += zz * f[i3 + ZZ]; - } - } - - upd_vir(vir[XX], dvxx, dvxy, dvxz); - upd_vir(vir[YY], dvyx, dvyy, dvyz); - upd_vir(vir[ZZ], dvzx, dvzy, dvzz); -} - -void f_calc_vir(int i0, int i1, const rvec x[], const rvec f[], tensor vir, const t_graph* g, const matrix box) +void f_calc_vir(int i0, int i1, const rvec x[], const rvec f[], tensor vir, const matrix box) { - int start, end; - - if (g && g->numNodes() > 0) - { - /* Calculate virial for bonded forces only when they belong to - * this node. - */ - start = std::max(i0, g->at_start); - end = std::min(i1, g->at_end); - lo_fcv(start, end, x[0], f[0], vir, g->ishift[0], box[0], TRICLINIC(box)); - - /* If not all atoms are bonded, calculate their virial contribution - * anyway, without shifting back their coordinates. - * Note the nifty pointer arithmetic... - */ - if (start > i0) - { - calc_vir(start - i0, x + i0, f + i0, vir, FALSE, box); - } - if (end < i1) - { - calc_vir(i1 - end, x + end, f + end, vir, FALSE, box); - } - } - else - { - calc_vir(i1 - i0, x + i0, f + i0, vir, FALSE, box); - } + calc_vir(i1 - i0, x + i0, f + i0, vir, FALSE, box); } diff --git a/src/gromacs/mdlib/calcvir.h b/src/gromacs/mdlib/calcvir.h index c2eedbc289..1569c5bb7b 100644 --- a/src/gromacs/mdlib/calcvir.h +++ b/src/gromacs/mdlib/calcvir.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020, 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. @@ -37,13 +37,12 @@ #include "gromacs/math/vectypes.h" -struct t_graph; struct t_pbc; void calc_vir(int nxf, const rvec x[], const rvec f[], tensor vir, bool bScrewPBC, const matrix box); /* Calculate virial for nxf atoms, and add it to vir */ -void f_calc_vir(int i0, int i1, const rvec x[], const rvec f[], tensor vir, const t_graph* g, const rvec shift_vec[]); +void f_calc_vir(int i0, int i1, const rvec x[], const rvec f[], tensor vir, const rvec shift_vec[]); /* Calculate virial taking periodicity into account */ #endif diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp index fc407a0fe3..591e64f661 100644 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@ -67,13 +67,15 @@ #include "gromacs/mdtypes/mdatom.h" #include "gromacs/mdtypes/simulation_workload.h" #include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/timing/wallcycle.h" #include "gromacs/utility/exceptions.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" +using gmx::ArrayRef; +using gmx::RVec; + static void clearEwaldThreadOutput(ewald_corr_thread_t* ewc_t) { ewc_t->Vcorr_q = 0; @@ -99,29 +101,28 @@ static void reduceEwaldThreadOuput(int nthreads, ewald_corr_thread_t* ewc_t) } } -void do_force_lowlevel(t_forcerec* fr, - const t_inputrec* ir, - const InteractionDefinitions& idef, - const t_commrec* cr, - const gmx_multisim_t* ms, - t_nrnb* nrnb, - gmx_wallcycle_t wcycle, - const t_mdatoms* md, - gmx::ArrayRefWithPadding coordinates, - gmx::ArrayRef xWholeMolecules, - history_t* hist, - gmx::ForceOutputs* forceOutputs, - gmx_enerdata_t* enerd, - t_fcdata* fcd, - const matrix box, - const real* lambda, - const t_graph* graph, - const rvec* mu_tot, - const gmx::StepWorkload& stepWork, - const DDBalanceRegionHandler& ddBalanceRegionHandler) +void do_force_lowlevel(t_forcerec* fr, + const t_inputrec* ir, + const InteractionDefinitions& idef, + const t_commrec* cr, + const gmx_multisim_t* ms, + t_nrnb* nrnb, + gmx_wallcycle_t wcycle, + const t_mdatoms* md, + gmx::ArrayRefWithPadding coordinates, + ArrayRef xWholeMolecules, + history_t* hist, + gmx::ForceOutputs* forceOutputs, + gmx_enerdata_t* enerd, + t_fcdata* fcd, + const matrix box, + const real* lambda, + const rvec* mu_tot, + const gmx::StepWorkload& stepWork, + const DDBalanceRegionHandler& ddBalanceRegionHandler) { // TODO: Replace all uses of x by const coordinates - rvec* x = as_rvec_array(coordinates.paddedArrayRef().data()); + const rvec* x = as_rvec_array(coordinates.paddedArrayRef().data()); auto& forceWithVirial = forceOutputs->forceWithVirial(); @@ -140,33 +141,6 @@ void do_force_lowlevel(t_forcerec* fr, } } - /* Shift the coordinates. Must be done before listed forces and PPPM, - * but is also necessary for SHAKE and update, therefore it can NOT - * go when no listed forces have to be evaluated. - * - * The shifting and PBC code is deliberately not timed, since with - * the Verlet scheme it only takes non-zero time with triclinic - * boxes, and even then the time is around a factor of 100 less - * than the next smallest counter. - */ - - - /* Here sometimes we would not need to shift with NBFonly, - * but we do so anyhow for consistency of the returned coordinates. - */ - if (graph) - { - shift_self(graph, box, x); - if (TRICLINIC(box)) - { - inc_nrnb(nrnb, eNR_SHIFTX, 2 * graph->numNodes()); - } - else - { - inc_nrnb(nrnb, eNR_SHIFTX, graph->numNodes()); - } - } - { t_pbc pbc; @@ -182,7 +156,7 @@ void do_force_lowlevel(t_forcerec* fr, } do_force_listed(wcycle, box, ir->fepvals, cr, ms, idef, x, xWholeMolecules, hist, - forceOutputs, fr, &pbc, graph, enerd, nrnb, lambda, md, fcd, + forceOutputs, fr, &pbc, enerd, nrnb, lambda, md, fcd, DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr, stepWork); } diff --git a/src/gromacs/mdlib/force.h b/src/gromacs/mdlib/force.h index 1c36420304..9d155f369d 100644 --- a/src/gromacs/mdlib/force.h +++ b/src/gromacs/mdlib/force.h @@ -57,7 +57,6 @@ struct pull_t; struct t_commrec; struct t_fcdata; struct t_forcerec; -struct t_graph; struct t_inputrec; struct t_lambda; struct t_mdatoms; @@ -95,7 +94,6 @@ void do_force(FILE* log, gmx_enerdata_t* enerd, t_fcdata* fcd, gmx::ArrayRef lambda, - t_graph* graph, t_forcerec* fr, gmx::MdrunScheduleWorkload* runScheduleWork, const gmx_vsite_t* vsite, @@ -120,26 +118,25 @@ void do_force(FILE* log, * xWholeMolecules only needs to contain whole molecules when orientation * restraints need to be computed and can be empty otherwise. */ -void do_force_lowlevel(t_forcerec* fr, - const t_inputrec* ir, - const InteractionDefinitions& idef, - const t_commrec* cr, - const gmx_multisim_t* ms, - t_nrnb* nrnb, - gmx_wallcycle* wcycle, - const t_mdatoms* md, - gmx::ArrayRefWithPadding coordinates, - gmx::ArrayRef xWholeMolecules, - history_t* hist, - gmx::ForceOutputs* forceOutputs, - gmx_enerdata_t* enerd, - t_fcdata* fcd, - const matrix box, - const real* lambda, - const t_graph* graph, - const rvec* mu_tot, - const gmx::StepWorkload& stepWork, - const DDBalanceRegionHandler& ddBalanceRegionHandler); +void do_force_lowlevel(t_forcerec* fr, + const t_inputrec* ir, + const InteractionDefinitions& idef, + const t_commrec* cr, + const gmx_multisim_t* ms, + t_nrnb* nrnb, + gmx_wallcycle* wcycle, + const t_mdatoms* md, + gmx::ArrayRefWithPadding coordinates, + gmx::ArrayRef xWholeMolecules, + history_t* hist, + gmx::ForceOutputs* forceOutputs, + gmx_enerdata_t* enerd, + t_fcdata* fcd, + const matrix box, + const real* lambda, + const rvec* mu_tot, + const gmx::StepWorkload& stepWork, + const DDBalanceRegionHandler& ddBalanceRegionHandler); /* Call all the force routines */ #endif diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 8cda9b4e9d..c21c809ca2 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -97,7 +97,6 @@ #include "gromacs/nbnxm/nbnxm.h" #include "gromacs/nbnxm/nbnxm_gpu.h" #include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/pull.h" #include "gromacs/pulling/pull_rotation.h" @@ -150,7 +149,6 @@ static void calc_virial(int start, const rvec x[], const gmx::ForceWithShiftForces& forceWithShiftForces, tensor vir_part, - const t_graph* graph, const matrix box, t_nrnb* nrnb, const t_forcerec* fr, @@ -165,7 +163,7 @@ static void calc_virial(int start, * Total virial is computed in global_stat, called from do_md */ const rvec* f = as_rvec_array(forceWithShiftForces.force().data()); - f_calc_vir(start, start + homenr, x, f, vir_part, graph, box); + f_calc_vir(start, start + homenr, x, f, vir_part, box); inc_nrnb(nrnb, eNR_VIRIAL, homenr); if (debug) @@ -288,7 +286,6 @@ static void post_process_forces(const t_commrec* cr, ForceOutputs* forceOutputs, tensor vir_force, const t_mdatoms* mdatoms, - const t_graph* graph, const t_forcerec* fr, const gmx_vsite_t* vsite, const StepWorkload& stepWork) @@ -308,7 +305,7 @@ static void post_process_forces(const t_commrec* cr, */ matrix virial = { { 0 } }; spread_vsite_f(vsite, x, fDirectVir, nullptr, stepWork.computeVirial, virial, nrnb, - top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle); + top->idef, fr->pbcType, fr->bMolPBC, box, cr, wcycle); forceWithVirial.addVirialContribution(virial); } @@ -954,7 +951,6 @@ void do_force(FILE* fplog, gmx_enerdata_t* enerd, t_fcdata* fcd, gmx::ArrayRef lambda, - t_graph* graph, t_forcerec* fr, gmx::MdrunScheduleWorkload* runScheduleWork, const gmx_vsite_t* vsite, @@ -1007,10 +1003,6 @@ void do_force(FILE* fplog, gmx_omp_nthreads_get(emntDefault)); inc_nrnb(nrnb, eNR_SHIFTX, mdatoms->homenr); } - else if (EI_ENERGY_MINIMIZATION(inputrec->eI) && graph) - { - unshift_self(graph, box, as_rvec_array(x.unpaddedArrayRef().data())); - } } nbnxn_atomdata_copy_shiftvec(stepWork.haveDynamicBox, fr->shift_vec, nbv->nbat.get()); @@ -1032,11 +1024,7 @@ void do_force(FILE* fplog, // Otherwise the send will occur after H2D coordinate transfer. if (!thisRankHasDuty(cr, DUTY_PME) && !pmeSendCoordinatesFromGpu) { - /* Send particle coordinates to the pme nodes. - * Since this is only implemented for domain decomposition - * and domain decomposition does not use the graph, - * we do not need to worry about shifting. - */ + /* Send particle coordinates to the pme nodes */ if (!stepWork.doNeighborSearch && simulationWork.useGpuUpdate) { GMX_RELEASE_ASSERT(false, @@ -1116,11 +1104,7 @@ void do_force(FILE* fplog, // Otherwise the send will occur before the H2D coordinate transfer. if (pmeSendCoordinatesFromGpu) { - /* Send particle coordinates to the pme nodes. - * Since this is only implemented for domain decomposition - * and domain decomposition does not use the graph, - * we do not need to worry about shifting. - */ + /* Send particle coordinates to the pme nodes */ gmx_pme_send_coordinates(fr, cr, box, as_rvec_array(x.unpaddedArrayRef().data()), lambda[efptCOUL], lambda[efptVDW], (stepWork.computeVirial || stepWork.computeEnergy), step, simulationWork.useGpuPmePpCommunication, reinitGpuPmePpComms, @@ -1136,11 +1120,6 @@ void do_force(FILE* fplog, /* do gridding for pair search */ if (stepWork.doNeighborSearch) { - if (graph && stepWork.stateChanged) - { - /* Calculate intramolecular shift vectors to make molecules whole */ - mk_mshift(fplog, graph, fr->pbcType, box, as_rvec_array(x.unpaddedArrayRef().data())); - } if (fr->wholeMoleculeTransform && stepWork.stateChanged) { fr->wholeMoleculeTransform->updateForAtomPbcJumps(x.unpaddedArrayRef(), box); @@ -1540,7 +1519,7 @@ void do_force(FILE* fplog, } /* Compute the bonded and non-bonded energies and optionally forces */ do_force_lowlevel(fr, inputrec, top->idef, cr, ms, nrnb, wcycle, mdatoms, x, xWholeMolecules, - hist, &forceOut, enerd, fcd, box, lambda.data(), graph, + hist, &forceOut, enerd, fcd, box, lambda.data(), as_rvec_array(dipoleData.muStateAB), stepWork, ddBalanceRegionHandler); wallcycle_stop(wcycle, ewcFORCE); @@ -1816,15 +1795,14 @@ void do_force(FILE* fplog, { rvec* fshift = as_rvec_array(forceOut.forceWithShiftForces().shiftForces().data()); spread_vsite_f(vsite, as_rvec_array(x.unpaddedArrayRef().data()), f, fshift, FALSE, - nullptr, nrnb, top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle); + nullptr, nrnb, top->idef, fr->pbcType, fr->bMolPBC, box, cr, wcycle); } if (stepWork.computeVirial) { /* Calculation of the virial must be done after vsites! */ calc_virial(0, mdatoms->homenr, as_rvec_array(x.unpaddedArrayRef().data()), - forceOut.forceWithShiftForces(), vir_force, graph, box, nrnb, fr, - inputrec->pbcType); + forceOut.forceWithShiftForces(), vir_force, box, nrnb, fr, inputrec->pbcType); } } @@ -1842,7 +1820,7 @@ void do_force(FILE* fplog, if (stepWork.computeForces) { post_process_forces(cr, step, nrnb, wcycle, top, box, as_rvec_array(x.unpaddedArrayRef().data()), - &forceOut, vir_force, mdatoms, graph, fr, vsite, stepWork); + &forceOut, vir_force, mdatoms, fr, vsite, stepWork); } if (stepWork.computeEnergy) diff --git a/src/gromacs/mdlib/tests/leapfrogtestrunners.cpp b/src/gromacs/mdlib/tests/leapfrogtestrunners.cpp index ec6cff4b78..449dbb1926 100644 --- a/src/gromacs/mdlib/tests/leapfrogtestrunners.cpp +++ b/src/gromacs/mdlib/tests/leapfrogtestrunners.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2019, by the GROMACS development team, led by + * Copyright (c) 2019,2020, 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. @@ -89,7 +89,7 @@ void integrateLeapFrogSimple(LeapFrogTestData* testData, int numSteps) testData->velocityScalingMatrix_, testData->update_.get(), etrtNONE, nullptr, nullptr); finish_update(&testData->inputRecord_, &testData->mdAtoms_, &testData->state_, nullptr, - nullptr, nullptr, testData->update_.get(), nullptr); + testData->update_.get(), nullptr); } auto xp = makeArrayRef(*testData->update_->xp()).subArray(0, testData->numAtoms_); for (int i = 0; i < testData->numAtoms_; i++) diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp index 22fdb3a27d..17089db0a4 100644 --- a/src/gromacs/mdlib/update.cpp +++ b/src/gromacs/mdlib/update.cpp @@ -70,7 +70,6 @@ #include "gromacs/mdtypes/mdatom.h" #include "gromacs/mdtypes/state.h" #include "gromacs/pbcutil/boxutilities.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/pull.h" #include "gromacs/random/tabulatednormaldistribution.h" @@ -1576,8 +1575,6 @@ void update_sd_second_half(int64_t step, void finish_update(const t_inputrec* inputrec, /* input record and box stuff */ const t_mdatoms* md, t_state* state, - const t_graph* graph, - t_nrnb* nrnb, gmx_wallcycle_t wcycle, Update* upd, const gmx::Constraints* constr) @@ -1599,7 +1596,9 @@ void finish_update(const t_inputrec* inputrec, /* input record and box stu * the frozen dimensions. To freeze such degrees of freedom * we copy them back here to later copy them forward. It would * be more elegant and slightly more efficient to copies zero - * times instead of twice, but the graph case below prevents this. + * times instead of twice. + * + * TODO: Now the graph is removed, remove this double copy. */ const ivec* nFreeze = inputrec->opts.nFreeze; bool partialFreezeAndConstraints = false; @@ -1630,19 +1629,7 @@ void finish_update(const t_inputrec* inputrec, /* input record and box stu } } - if (graph && graph->numNodes() > 0) - { - unshift_x(graph, state->box, state->x.rvec_array(), upd->xp()->rvec_array()); - if (TRICLINIC(state->box)) - { - inc_nrnb(nrnb, eNR_SHIFTX, 2 * graph->numNodes()); - } - else - { - inc_nrnb(nrnb, eNR_SHIFTX, graph->numNodes()); - } - } - else + // TODO: Get rid of this copy { auto xp = makeConstArrayRef(*upd->xp()).subArray(0, homenr); auto x = makeArrayRef(state->x).subArray(0, homenr); diff --git a/src/gromacs/mdlib/update.h b/src/gromacs/mdlib/update.h index acd097e3b7..9927a92422 100644 --- a/src/gromacs/mdlib/update.h +++ b/src/gromacs/mdlib/update.h @@ -197,8 +197,6 @@ void update_sd_second_half(int64_t step, void finish_update(const t_inputrec* inputrec, const t_mdatoms* md, t_state* state, - const t_graph* graph, - t_nrnb* nrnb, gmx_wallcycle_t wcycle, gmx::Update* upd, const gmx::Constraints* constr); diff --git a/src/gromacs/mdlib/vsite.cpp b/src/gromacs/mdlib/vsite.cpp index 1973d2088c..bed9b8faa1 100644 --- a/src/gromacs/mdlib/vsite.cpp +++ b/src/gromacs/mdlib/vsite.cpp @@ -55,7 +55,6 @@ #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/mdatom.h" #include "gromacs/pbcutil/ishift.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/timing/wallcycle.h" #include "gromacs/topology/ifunc.h" @@ -669,17 +668,10 @@ void construct_vsites(const gmx_vsite_t* vsite, } } -static void spread_vsite2(const t_iatom ia[], - real a, - const rvec x[], - rvec f[], - rvec fshift[], - const t_pbc* pbc, - const t_graph* g) +static void spread_vsite2(const t_iatom ia[], real a, const rvec x[], rvec f[], rvec fshift[], const t_pbc* pbc) { rvec fi, fj, dx; t_iatom av, ai, aj; - ivec di; int siv, sij; av = ia[1]; @@ -694,14 +686,7 @@ static void spread_vsite2(const t_iatom ia[], rvec_inc(f[aj], fj); /* 6 Flops */ - if (g) - { - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, av), di); - siv = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), di); - sij = IVEC2IS(di); - } - else if (pbc) + if (pbc) { siv = pbc_dx_aiuc(pbc, x[ai], x[av], dx); sij = pbc_dx_aiuc(pbc, x[ai], x[aj], dx); @@ -745,15 +730,14 @@ void constructVsitesGlobal(const gmx_mtop_t& mtop, gmx::ArrayRef x) } } -static void spread_vsite2FD(const t_iatom ia[], - real a, - const rvec x[], - rvec f[], - rvec fshift[], - gmx_bool VirCorr, - matrix dxdf, - const t_pbc* pbc, - const t_graph* g) +static void spread_vsite2FD(const t_iatom ia[], + real a, + const rvec x[], + rvec f[], + rvec fshift[], + gmx_bool VirCorr, + matrix dxdf, + const t_pbc* pbc) { const int av = ia[1]; const int ai = ia[2]; @@ -791,15 +775,7 @@ static void spread_vsite2FD(const t_iatom ia[], if (fshift) { int svi; - if (g) - { - ivec di; - ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di); - svi = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di); - sji = IVEC2IS(di); - } - else if (pbc) + if (pbc) { rvec xvi; svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi); @@ -848,18 +824,10 @@ static void spread_vsite2FD(const t_iatom ia[], /* TOTAL: 38 flops */ } -static void spread_vsite3(const t_iatom ia[], - real a, - real b, - const rvec x[], - rvec f[], - rvec fshift[], - const t_pbc* pbc, - const t_graph* g) +static void spread_vsite3(const t_iatom ia[], real a, real b, const rvec x[], rvec f[], rvec fshift[], const t_pbc* pbc) { rvec fi, fj, fk, dx; int av, ai, aj, ak; - ivec di; int siv, sij, sik; av = ia[1]; @@ -877,16 +845,7 @@ static void spread_vsite3(const t_iatom ia[], rvec_inc(f[ak], fk); /* 9 Flops */ - if (g) - { - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, ia[1]), di); - siv = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), di); - sij = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, ak), di); - sik = IVEC2IS(di); - } - else if (pbc) + if (pbc) { siv = pbc_dx_aiuc(pbc, x[ai], x[av], dx); sij = pbc_dx_aiuc(pbc, x[ai], x[aj], dx); @@ -910,22 +869,20 @@ static void spread_vsite3(const t_iatom ia[], /* TOTAL: 20 flops */ } -static void spread_vsite3FD(const t_iatom ia[], - real a, - real b, - const rvec x[], - rvec f[], - rvec fshift[], - gmx_bool VirCorr, - matrix dxdf, - const t_pbc* pbc, - const t_graph* g) +static void spread_vsite3FD(const t_iatom ia[], + real a, + real b, + const rvec x[], + rvec f[], + rvec fshift[], + gmx_bool VirCorr, + matrix dxdf, + const t_pbc* pbc) { real fproj, a1; rvec xvi, xij, xjk, xix, fv, temp; t_iatom av, ai, aj, ak; int svi, sji, skj; - ivec di; av = ia[1]; ai = ia[2]; @@ -969,16 +926,7 @@ static void spread_vsite3FD(const t_iatom ia[], f[ak][ZZ] += a * temp[ZZ]; /* 19 Flops */ - if (g) - { - ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di); - svi = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di); - sji = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, aj), di); - skj = IVEC2IS(di); - } - else if (pbc) + if (pbc) { svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi); } @@ -1028,22 +976,20 @@ static void spread_vsite3FD(const t_iatom ia[], /* TOTAL: 61 flops */ } -static void spread_vsite3FAD(const t_iatom ia[], - real a, - real b, - const rvec x[], - rvec f[], - rvec fshift[], - gmx_bool VirCorr, - matrix dxdf, - const t_pbc* pbc, - const t_graph* g) +static void spread_vsite3FAD(const t_iatom ia[], + real a, + real b, + const rvec x[], + rvec f[], + rvec fshift[], + gmx_bool VirCorr, + matrix dxdf, + const t_pbc* pbc) { rvec xvi, xij, xjk, xperp, Fpij, Fppp, fv, f1, f2, f3; real a1, b1, c1, c2, invdij, invdij2, invdp, fproj; t_iatom av, ai, aj, ak; int svi, sji, skj, d; - ivec di; av = ia[1]; ai = ia[2]; @@ -1097,16 +1043,7 @@ static void spread_vsite3FAD(const t_iatom ia[], f[ak][ZZ] += f2[ZZ]; /* 30 Flops */ - if (g) - { - ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di); - svi = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di); - sji = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, aj), di); - skj = IVEC2IS(di); - } - else if (pbc) + if (pbc) { svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi); } @@ -1150,22 +1087,20 @@ static void spread_vsite3FAD(const t_iatom ia[], /* TOTAL: 113 flops */ } -static void spread_vsite3OUT(const t_iatom ia[], - real a, - real b, - real c, - const rvec x[], - rvec f[], - rvec fshift[], - gmx_bool VirCorr, - matrix dxdf, - const t_pbc* pbc, - const t_graph* g) +static void spread_vsite3OUT(const t_iatom ia[], + real a, + real b, + real c, + const rvec x[], + rvec f[], + rvec fshift[], + gmx_bool VirCorr, + matrix dxdf, + const t_pbc* pbc) { rvec xvi, xij, xik, fv, fj, fk; real cfx, cfy, cfz; int av, ai, aj, ak; - ivec di; int svi, sji, ski; av = ia[1]; @@ -1200,16 +1135,7 @@ static void spread_vsite3OUT(const t_iatom ia[], rvec_inc(f[ak], fk); /* 15 Flops */ - if (g) - { - ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di); - svi = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di); - sji = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, ai), di); - ski = IVEC2IS(di); - } - else if (pbc) + if (pbc) { svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi); } @@ -1246,22 +1172,20 @@ static void spread_vsite3OUT(const t_iatom ia[], /* TOTAL: 54 flops */ } -static void spread_vsite4FD(const t_iatom ia[], - real a, - real b, - real c, - const rvec x[], - rvec f[], - rvec fshift[], - gmx_bool VirCorr, - matrix dxdf, - const t_pbc* pbc, - const t_graph* g) +static void spread_vsite4FD(const t_iatom ia[], + real a, + real b, + real c, + const rvec x[], + rvec f[], + rvec fshift[], + gmx_bool VirCorr, + matrix dxdf, + const t_pbc* pbc) { real fproj, a1; rvec xvi, xij, xjk, xjl, xix, fv, temp; int av, ai, aj, ak, al; - ivec di; int svi, sji, skj, slj, m; av = ia[1]; @@ -1309,18 +1233,7 @@ static void spread_vsite4FD(const t_iatom ia[], } /* 26 Flops */ - if (g) - { - ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di); - svi = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di); - sji = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, aj), di); - skj = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, al), SHIFT_IVEC(g, aj), di); - slj = IVEC2IS(di); - } - else if (pbc) + if (pbc) { svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi); } @@ -1361,23 +1274,21 @@ static void spread_vsite4FD(const t_iatom ia[], } -static void spread_vsite4FDN(const t_iatom ia[], - real a, - real b, - real c, - const rvec x[], - rvec f[], - rvec fshift[], - gmx_bool VirCorr, - matrix dxdf, - const t_pbc* pbc, - const t_graph* g) +static void spread_vsite4FDN(const t_iatom ia[], + real a, + real b, + real c, + const rvec x[], + rvec f[], + rvec fshift[], + gmx_bool VirCorr, + matrix dxdf, + const t_pbc* pbc) { rvec xvi, xij, xik, xil, ra, rb, rja, rjb, rab, rm, rt; rvec fv, fj, fk, fl; real invrm, denom; real cfx, cfy, cfz; - ivec di; int av, ai, aj, ak, al; int svi, sij, sik, sil; @@ -1478,18 +1389,7 @@ static void spread_vsite4FDN(const t_iatom ia[], rvec_inc(f[al], fl); /* 21 flops */ - if (g) - { - ivec_sub(SHIFT_IVEC(g, av), SHIFT_IVEC(g, ai), di); - svi = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di); - sij = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, ak), SHIFT_IVEC(g, ai), di); - sik = IVEC2IS(di); - ivec_sub(SHIFT_IVEC(g, al), SHIFT_IVEC(g, ai), di); - sil = IVEC2IS(di); - } - else if (pbc) + if (pbc) { svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi); } @@ -1534,13 +1434,11 @@ static int spread_vsiten(const t_iatom ia[], const rvec x[], rvec f[], rvec fshift[], - const t_pbc* pbc, - const t_graph* g) + const t_pbc* pbc) { rvec xv, dx, fi; int n3, av, i, ai; real a; - ivec di; int siv; n3 = 3 * ip[ia[0]].vsiten.n; @@ -1550,12 +1448,7 @@ static int spread_vsiten(const t_iatom ia[], for (i = 0; i < n3; i += 3) { ai = ia[i + 2]; - if (g) - { - ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, av), di); - siv = IVEC2IS(di); - } - else if (pbc) + if (pbc) { siv = pbc_dx_aiuc(pbc, x[ai], xv, dx); } @@ -1597,7 +1490,6 @@ static void spread_vsite_f_thread(const rvec x[], matrix dxdf, ArrayRef ip, ArrayRef ilist, - const t_graph* g, const t_pbc* pbc_null) { const PbcMode pbcMode = getPbcMode(pbc_null); @@ -1636,38 +1528,38 @@ static void spread_vsite_f_thread(const rvec x[], /* Construct the vsite depending on type */ switch (ftype) { - case F_VSITE2: spread_vsite2(ia, a1, x, f, fshift, pbc_null2, g); break; + case F_VSITE2: spread_vsite2(ia, a1, x, f, fshift, pbc_null2); break; case F_VSITE2FD: - spread_vsite2FD(ia, a1, x, f, fshift, VirCorr, dxdf, pbc_null2, g); + spread_vsite2FD(ia, a1, x, f, fshift, VirCorr, dxdf, pbc_null2); break; case F_VSITE3: b1 = ip[tp].vsite.b; - spread_vsite3(ia, a1, b1, x, f, fshift, pbc_null2, g); + spread_vsite3(ia, a1, b1, x, f, fshift, pbc_null2); break; case F_VSITE3FD: b1 = ip[tp].vsite.b; - spread_vsite3FD(ia, a1, b1, x, f, fshift, VirCorr, dxdf, pbc_null2, g); + spread_vsite3FD(ia, a1, b1, x, f, fshift, VirCorr, dxdf, pbc_null2); break; case F_VSITE3FAD: b1 = ip[tp].vsite.b; - spread_vsite3FAD(ia, a1, b1, x, f, fshift, VirCorr, dxdf, pbc_null2, g); + spread_vsite3FAD(ia, a1, b1, x, f, fshift, VirCorr, dxdf, pbc_null2); break; case F_VSITE3OUT: b1 = ip[tp].vsite.b; c1 = ip[tp].vsite.c; - spread_vsite3OUT(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2, g); + spread_vsite3OUT(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2); break; case F_VSITE4FD: b1 = ip[tp].vsite.b; c1 = ip[tp].vsite.c; - spread_vsite4FD(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2, g); + spread_vsite4FD(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2); break; case F_VSITE4FDN: b1 = ip[tp].vsite.b; c1 = ip[tp].vsite.c; - spread_vsite4FDN(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2, g); + spread_vsite4FDN(ia, a1, b1, c1, x, f, fshift, VirCorr, dxdf, pbc_null2); break; - case F_VSITEN: inc = spread_vsiten(ia, ip, x, f, fshift, pbc_null2, g); break; + case F_VSITEN: inc = spread_vsiten(ia, ip, x, f, fshift, pbc_null2); break; default: gmx_fatal(FARGS, "No such vsite type %d in %s, line %d", ftype, __FILE__, __LINE__); } @@ -1707,7 +1599,6 @@ void spread_vsite_f(const gmx_vsite_t* vsite, const InteractionDefinitions& idef, PbcType pbcType, gmx_bool bMolPBC, - const t_graph* g, const matrix box, const t_commrec* cr, gmx_wallcycle* wcycle) @@ -1744,7 +1635,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite, { clear_mat(dxdf); } - spread_vsite_f_thread(x, f, fshift, VirCorr, dxdf, idef.iparams, idef.il, g, pbc_null); + spread_vsite_f_thread(x, f, fshift, VirCorr, dxdf, idef.iparams, idef.il, pbc_null); if (VirCorr) { @@ -1765,7 +1656,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite, clear_mat(vsite->tData[vsite->nthreads]->dxdf); } spread_vsite_f_thread(x, f, fshift, VirCorr, vsite->tData[vsite->nthreads]->dxdf, - idef.iparams, vsite->tData[vsite->nthreads]->ilist, g, pbc_null); + idef.iparams, vsite->tData[vsite->nthreads]->ilist, pbc_null); #pragma omp parallel num_threads(vsite->nthreads) { @@ -1812,7 +1703,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite, copy_rvec(f[idTask->vsite[i]], idTask->force[idTask->vsite[i]]); } spread_vsite_f_thread(x, as_rvec_array(idTask->force.data()), fshift_t, VirCorr, - tData.dxdf, idef.iparams, tData.idTask.ilist, g, pbc_null); + tData.dxdf, idef.iparams, tData.idTask.ilist, pbc_null); /* We need a barrier before reducing forces below * that have been produced by a different thread above. @@ -1851,7 +1742,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite, /* Spread the vsites that spread locally only */ spread_vsite_f_thread(x, f, fshift_t, VirCorr, tData.dxdf, idef.iparams, - tData.ilist, g, pbc_null); + tData.ilist, pbc_null); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } diff --git a/src/gromacs/mdlib/vsite.h b/src/gromacs/mdlib/vsite.h index 291fbc9525..e689405942 100644 --- a/src/gromacs/mdlib/vsite.h +++ b/src/gromacs/mdlib/vsite.h @@ -50,7 +50,6 @@ struct gmx_localtop_t; struct gmx_mtop_t; struct t_commrec; -struct t_graph; struct InteractionList; struct t_mdatoms; struct t_nrnb; @@ -131,7 +130,6 @@ void spread_vsite_f(const gmx_vsite_t* vsite, const InteractionDefinitions& idef, PbcType pbcType, gmx_bool bMolPBC, - const t_graph* g, const matrix box, const t_commrec* cr, gmx_wallcycle* wcycle); diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 91360ccd28..9d341e1eae 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -123,7 +123,6 @@ #include "gromacs/modularsimulator/energyelement.h" #include "gromacs/nbnxm/gpu_data_mgmt.h" #include "gromacs/nbnxm/nbnxm.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/output.h" #include "gromacs/pulling/pull.h" @@ -177,7 +176,6 @@ void gmx::LegacySimulator::do_md() gmx_repl_ex_t repl_ex = nullptr; PaddedHostVector f{}; gmx_global_stat_t gstat; - t_graph* graph = nullptr; gmx_shellfc_t* shellfc; gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition; gmx_bool bTemp, bPres, bTrotter; @@ -342,7 +340,7 @@ void gmx::LegacySimulator::do_md() state = state_global; /* Generate and initialize new topology */ - mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &graph, mdAtoms, constr, vsite, shellfc); + mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, mdAtoms, constr, vsite, shellfc); upd.setNumAtoms(state->natoms); } @@ -388,7 +386,6 @@ void gmx::LegacySimulator::do_md() || (!haveFreeEnergyType(*ir, efptBONDED) && !haveFreeEnergyType(*ir, efptMASS)), "Free energy perturbation of masses and constraints are not supported with the GPU " "update."); - GMX_RELEASE_ASSERT(graph == nullptr, "The graph is not supported with GPU update."); if (constr != nullptr && constr->numConstraintsTotal() > 0) { @@ -839,7 +836,7 @@ void gmx::LegacySimulator::do_md() /* Correct the new box if it is too skewed */ if (inputrecDynamicBox(ir)) { - if (correct_box(fplog, step, state->box, graph)) + if (correct_box(fplog, step, state->box)) { bMasterState = TRUE; // If update is offloaded, it should be informed about the box size change @@ -955,8 +952,8 @@ void gmx::LegacySimulator::do_md() imdSession, pull_work, bNS, force_flags, &top, constr, enerd, fcd, state->natoms, state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(), state->box, state->lambda, &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, graph, - shellfc, fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); + f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, shellfc, + fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); } else { @@ -980,8 +977,8 @@ void gmx::LegacySimulator::do_md() */ do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, graph, - fr, runScheduleWork, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr, + f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, fr, + runScheduleWork, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr, (bNS ? GMX_FORCE_NS : 0) | force_flags, ddBalanceRegionHandler); } @@ -1340,7 +1337,7 @@ void gmx::LegacySimulator::do_md() update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state, cr, nrnb, wcycle, &upd, constr, do_log, do_ene); - finish_update(ir, mdatoms, state, graph, nrnb, wcycle, &upd, constr); + finish_update(ir, mdatoms, state, wcycle, &upd, constr); } if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM) @@ -1371,7 +1368,7 @@ void gmx::LegacySimulator::do_md() * to numerical errors, or are they important * physically? I'm thinking they are just errors, but not completely sure. * For now, will call without actually constraining, constr=NULL*/ - finish_update(ir, mdatoms, state, graph, nrnb, wcycle, &upd, nullptr); + finish_update(ir, mdatoms, state, wcycle, &upd, nullptr); } if (EI_VV(ir->eI)) { @@ -1399,17 +1396,8 @@ void gmx::LegacySimulator::do_md() if (vsite != nullptr) { wallcycle_start(wcycle, ewcVSITECONSTR); - if (graph != nullptr) - { - shift_self(graph, state->box, state->x.rvec_array()); - } construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(), top.idef.iparams, top.idef.il, fr->pbcType, fr->bMolPBC, cr, state->box); - - if (graph != nullptr) - { - unshift_self(graph, state->box, state->x.rvec_array()); - } wallcycle_stop(wcycle, ewcVSITECONSTR); } diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index 6984258b8e..247e25caa6 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -116,7 +116,6 @@ #include "gromacs/mdtypes/state.h" #include "gromacs/mimic/communicator.h" #include "gromacs/mimic/utilities.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/pull.h" #include "gromacs/timing/wallcycle.h" @@ -151,7 +150,6 @@ void gmx::LegacySimulator::do_mimic() rvec mu_tot; PaddedHostVector f{}; gmx_global_stat_t gstat; - t_graph* graph = nullptr; gmx_shellfc_t* shellfc; double cycles; @@ -276,7 +274,7 @@ void gmx::LegacySimulator::do_mimic() /* Copy the pointer to the global state */ state = state_global; - mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &graph, mdAtoms, constr, vsite, shellfc); + mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, mdAtoms, constr, vsite, shellfc); } auto mdatoms = mdAtoms->mdatoms(); @@ -421,8 +419,8 @@ void gmx::LegacySimulator::do_mimic() imdSession, pull_work, bNS, force_flags, &top, constr, enerd, fcd, state->natoms, state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(), state->box, state->lambda, &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, graph, - shellfc, fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); + f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, shellfc, + fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); } else { @@ -435,8 +433,8 @@ void gmx::LegacySimulator::do_mimic() gmx_edsam* ed = nullptr; do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, graph, - fr, runScheduleWork, vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, + f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, fr, + runScheduleWork, vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler); } @@ -455,27 +453,12 @@ void gmx::LegacySimulator::do_mimic() stopHandler->setSignal(); - if (graph) - { - /* Need to unshift here */ - unshift_self(graph, state->box, as_rvec_array(state->x.data())); - } - if (vsite != nullptr) { wallcycle_start(wcycle, ewcVSITECONSTR); - if (graph != nullptr) - { - shift_self(graph, state->box, as_rvec_array(state->x.data())); - } construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()), top.idef.iparams, top.idef.il, fr->pbcType, fr->bMolPBC, cr, state->box); - - if (graph != nullptr) - { - unshift_self(graph, state->box, as_rvec_array(state->x.data())); - } wallcycle_stop(wcycle, ewcVSITECONSTR); } diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index bcda82941b..afe557f276 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -97,7 +97,6 @@ #include "gromacs/mdtypes/mdatom.h" #include "gromacs/mdtypes/mdrunoptions.h" #include "gromacs/mdtypes/state.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/timing/wallcycle.h" #include "gromacs/timing/walltime_accounting.h" @@ -371,7 +370,6 @@ static void init_em(FILE* fplog, gmx_localtop_t* top, t_nrnb* nrnb, t_forcerec* fr, - t_graph** graph, gmx::MDAtoms* mdAtoms, gmx_global_stat_t* gstat, gmx_vsite_t* vsite, @@ -423,8 +421,6 @@ static void init_em(FILE* fplog, imdSession, pull_work, &ems->s, &ems->f, mdAtoms, top, fr, vsite, constr, nrnb, nullptr, FALSE); dd_store_state(cr->dd, &ems->s); - - *graph = nullptr; } else { @@ -434,7 +430,7 @@ static void init_em(FILE* fplog, state_change_natoms(&ems->s, ems->s.natoms); ems->f.resizeWithPadding(ems->s.natoms); - mdAlgorithmsSetupAtomData(cr, ir, *top_global, top, fr, graph, mdAtoms, constr, vsite, + mdAlgorithmsSetupAtomData(cr, ir, *top_global, top, fr, mdAtoms, constr, vsite, shellfc ? *shellfc : nullptr); if (vsite) @@ -794,8 +790,6 @@ public: gmx::Constraints* constr; //! Handles strange things. t_fcdata* fcd; - //! Molecular graph for SHAKE. - t_graph* graph; //! Per-atom data for this domain. gmx::MDAtoms* mdAtoms; //! Handles how to calculate the forces. @@ -851,7 +845,7 @@ void EnergyEvaluator::run(em_state_t* ems, rvec mu_tot, tensor vir, tensor pres, do_force(fplog, cr, ms, inputrec, nullptr, nullptr, imdSession, pull_work, count, nrnb, wcycle, top, ems->s.box, ems->s.x.arrayRefWithPadding(), &ems->s.hist, ems->f.arrayRefWithPadding(), force_vir, mdAtoms->mdatoms(), enerd, fcd, ems->s.lambda, - graph, fr, runScheduleWork, vsite, mu_tot, t, nullptr, + fr, runScheduleWork, vsite, mu_tot, t, nullptr, GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES | GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY | (bNS ? GMX_FORCE_NS : 0), DDBalanceRegionHandler(cr)); @@ -1041,7 +1035,6 @@ void LegacySimulator::do_cg() gmx_localtop_t top(top_global->ffparams); gmx_global_stat_t gstat; - t_graph* graph; double tmp, minstep; real stepsize; real a, b, c, beta = 0.0; @@ -1089,7 +1082,7 @@ void LegacySimulator::do_cg() /* Init em and store the local state in s_min */ init_em(fplog, mdlog, CG, cr, inputrec, imdSession, pull_work, state_global, top_global, s_min, - &top, nrnb, fr, &graph, mdAtoms, &gstat, vsite, constr, nullptr); + &top, nrnb, fr, mdAtoms, &gstat, vsite, constr, nullptr); const bool simulationsShareState = false; gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec, top_global, nullptr, wcycle, @@ -1113,9 +1106,9 @@ void LegacySimulator::do_cg() } EnergyEvaluator energyEvaluator{ - fplog, mdlog, cr, ms, top_global, &top, inputrec, - imdSession, pull_work, nrnb, wcycle, gstat, vsite, constr, - fcd, graph, mdAtoms, fr, runScheduleWork, enerd + fplog, mdlog, cr, ms, top_global, &top, inputrec, imdSession, pull_work, + nrnb, wcycle, gstat, vsite, constr, fcd, mdAtoms, fr, runScheduleWork, + enerd }; /* Call the force routine and some auxiliary (neighboursearching etc.) */ /* do_force always puts the charge groups in the box and shifts again @@ -1640,7 +1633,6 @@ void LegacySimulator::do_lbfgs() em_state_t ems; gmx_localtop_t top(top_global->ffparams); gmx_global_stat_t gstat; - t_graph* graph; int ncorr, nmaxcorr, point, cp, neval, nminstep; double stepsize, step_taken, gpa, gpb, gpc, tmp, minstep; real * rho, *alpha, *p, *s, **dx, **dg; @@ -1703,7 +1695,7 @@ void LegacySimulator::do_lbfgs() /* Init em */ init_em(fplog, mdlog, LBFGS, cr, inputrec, imdSession, pull_work, state_global, top_global, - &ems, &top, nrnb, fr, &graph, mdAtoms, &gstat, vsite, constr, nullptr); + &ems, &top, nrnb, fr, mdAtoms, &gstat, vsite, constr, nullptr); const bool simulationsShareState = false; gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec, top_global, nullptr, wcycle, @@ -1767,9 +1759,9 @@ void LegacySimulator::do_lbfgs() */ neval++; EnergyEvaluator energyEvaluator{ - fplog, mdlog, cr, ms, top_global, &top, inputrec, - imdSession, pull_work, nrnb, wcycle, gstat, vsite, constr, - fcd, graph, mdAtoms, fr, runScheduleWork, enerd + fplog, mdlog, cr, ms, top_global, &top, inputrec, imdSession, pull_work, + nrnb, wcycle, gstat, vsite, constr, fcd, mdAtoms, fr, runScheduleWork, + enerd }; energyEvaluator.run(&ems, mu_tot, vir, pres, -1, TRUE); @@ -2358,7 +2350,6 @@ void LegacySimulator::do_steep() const char* SD = "Steepest Descents"; gmx_localtop_t top(top_global->ffparams); gmx_global_stat_t gstat; - t_graph* graph; real stepsize; real ustep; gmx_bool bDone, bAbort, do_x, do_f; @@ -2384,7 +2375,7 @@ void LegacySimulator::do_steep() /* Init em and store the local state in s_try */ init_em(fplog, mdlog, SD, cr, inputrec, imdSession, pull_work, state_global, top_global, s_try, - &top, nrnb, fr, &graph, mdAtoms, &gstat, vsite, constr, nullptr); + &top, nrnb, fr, mdAtoms, &gstat, vsite, constr, nullptr); const bool simulationsShareState = false; gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec, top_global, nullptr, wcycle, @@ -2414,9 +2405,9 @@ void LegacySimulator::do_steep() sp_header(fplog, SD, inputrec->em_tol, nsteps); } EnergyEvaluator energyEvaluator{ - fplog, mdlog, cr, ms, top_global, &top, inputrec, - imdSession, pull_work, nrnb, wcycle, gstat, vsite, constr, - fcd, graph, mdAtoms, fr, runScheduleWork, enerd + fplog, mdlog, cr, ms, top_global, &top, inputrec, imdSession, pull_work, + nrnb, wcycle, gstat, vsite, constr, fcd, mdAtoms, fr, runScheduleWork, + enerd }; /**** HERE STARTS THE LOOP **** @@ -2593,7 +2584,6 @@ void LegacySimulator::do_nm() int nnodes; gmx_localtop_t top(top_global->ffparams); gmx_global_stat_t gstat; - t_graph* graph; tensor vir, pres; rvec mu_tot = { 0 }; rvec* dfdx; @@ -2630,7 +2620,7 @@ void LegacySimulator::do_nm() /* Init em and store the local state in state_minimum */ init_em(fplog, mdlog, NM, cr, inputrec, imdSession, pull_work, state_global, top_global, - &state_work, &top, nrnb, fr, &graph, mdAtoms, &gstat, vsite, constr, &shellfc); + &state_work, &top, nrnb, fr, mdAtoms, &gstat, vsite, constr, &shellfc); const bool simulationsShareState = false; gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec, top_global, nullptr, wcycle, @@ -2708,9 +2698,9 @@ void LegacySimulator::do_nm() /* Make evaluate_energy do a single node force calculation */ cr->nnodes = 1; EnergyEvaluator energyEvaluator{ - fplog, mdlog, cr, ms, top_global, &top, inputrec, - imdSession, pull_work, nrnb, wcycle, gstat, vsite, constr, - fcd, graph, mdAtoms, fr, runScheduleWork, enerd + fplog, mdlog, cr, ms, top_global, &top, inputrec, imdSession, pull_work, + nrnb, wcycle, gstat, vsite, constr, fcd, mdAtoms, fr, runScheduleWork, + enerd }; energyEvaluator.run(&state_work, mu_tot, vir, pres, -1, TRUE); cr->nnodes = nnodes; @@ -2769,14 +2759,13 @@ void LegacySimulator::do_nm() if (shellfc) { /* Now is the time to relax the shells */ - relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, nullptr, step, inputrec, - imdSession, pull_work, bNS, force_flags, &top, constr, enerd, - fcd, state_work.s.natoms, state_work.s.x.arrayRefWithPadding(), - state_work.s.v.arrayRefWithPadding(), state_work.s.box, - state_work.s.lambda, &state_work.s.hist, - state_work.f.arrayRefWithPadding(), vir, mdatoms, nrnb, - wcycle, graph, shellfc, fr, runScheduleWork, t, mu_tot, - vsite, DDBalanceRegionHandler(nullptr)); + relax_shell_flexcon( + fplog, cr, ms, mdrunOptions.verbose, nullptr, step, inputrec, imdSession, + pull_work, bNS, force_flags, &top, constr, enerd, fcd, state_work.s.natoms, + state_work.s.x.arrayRefWithPadding(), state_work.s.v.arrayRefWithPadding(), + state_work.s.box, state_work.s.lambda, &state_work.s.hist, + state_work.f.arrayRefWithPadding(), vir, mdatoms, nrnb, wcycle, shellfc, + fr, runScheduleWork, t, mu_tot, vsite, DDBalanceRegionHandler(nullptr)); bNS = false; step++; } diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index 41d8c4f12d..3d369223a6 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -113,7 +113,6 @@ #include "gromacs/mdtypes/observableshistory.h" #include "gromacs/mdtypes/state.h" #include "gromacs/mimic/utilities.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pulling/pull.h" #include "gromacs/swap/swapcoords.h" @@ -146,7 +145,6 @@ using gmx::SimulationSignaller; * \param[in] idef Topology parameters, used for constructing vsites * \param[in] timeStep Time step, used for constructing vsites * \param[in] forceRec Force record, used for constructing vsites - * \param[in,out] graph The molecular graph, used for constructing vsites when != nullptr */ static void prepareRerunState(const t_trxframe& rerunFrame, t_state* globalState, @@ -154,8 +152,7 @@ static void prepareRerunState(const t_trxframe& rerunFrame, const gmx_vsite_t* vsite, const InteractionDefinitions& idef, double timeStep, - const t_forcerec& forceRec, - t_graph* graph) + const t_forcerec& forceRec) { auto x = makeArrayRef(globalState->x); auto rerunX = arrayRefFromArray(reinterpret_cast(rerunFrame.x), globalState->natoms); @@ -166,21 +163,9 @@ static void prepareRerunState(const t_trxframe& rerunFrame, { GMX_ASSERT(vsite, "Need valid vsite for constructing vsites"); - if (graph) - { - /* Following is necessary because the graph may get out of sync - * with the coordinates if we only have every N'th coordinate set - */ - mk_mshift(nullptr, graph, forceRec.pbcType, globalState->box, globalState->x.rvec_array()); - shift_self(graph, globalState->box, as_rvec_array(globalState->x.data())); - } construct_vsites(vsite, globalState->x.rvec_array(), timeStep, globalState->v.rvec_array(), idef.iparams, idef.il, forceRec.pbcType, forceRec.bMolPBC, nullptr, globalState->box); - if (graph) - { - unshift_self(graph, globalState->box, globalState->x.rvec_array()); - } } } @@ -205,7 +190,6 @@ void gmx::LegacySimulator::do_rerun() gmx_localtop_t top(top_global->ffparams); PaddedHostVector f{}; gmx_global_stat_t gstat; - t_graph* graph = nullptr; gmx_shellfc_t* shellfc; double cycles; @@ -348,7 +332,7 @@ void gmx::LegacySimulator::do_rerun() /* Copy the pointer to the global state */ state = state_global; - mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &graph, mdAtoms, constr, vsite, shellfc); + mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, mdAtoms, constr, vsite, shellfc); } auto mdatoms = mdAtoms->mdatoms(); @@ -513,8 +497,7 @@ void gmx::LegacySimulator::do_rerun() "decomposition, " "use a single rank"); } - prepareRerunState(rerun_fr, state_global, constructVsites, vsite, top.idef, ir->delta_t, - *fr, graph); + prepareRerunState(rerun_fr, state_global, constructVsites, vsite, top.idef, ir->delta_t, *fr); } isLastStep = isLastStep || stopHandler->stoppingAfterCurrentStep(bNS); @@ -550,8 +533,8 @@ void gmx::LegacySimulator::do_rerun() imdSession, pull_work, bNS, force_flags, &top, constr, enerd, fcd, state->natoms, state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(), state->box, state->lambda, &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, graph, - shellfc, fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); + f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, shellfc, + fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler); } else { @@ -564,8 +547,8 @@ void gmx::LegacySimulator::do_rerun() gmx_edsam* ed = nullptr; do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist, - f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, graph, - fr, runScheduleWork, vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, + f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, fr, + runScheduleWork, vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler); } @@ -584,27 +567,12 @@ void gmx::LegacySimulator::do_rerun() stopHandler->setSignal(); - if (graph) - { - /* Need to unshift here */ - unshift_self(graph, state->box, as_rvec_array(state->x.data())); - } - if (vsite != nullptr) { wallcycle_start(wcycle, ewcVSITECONSTR); - if (graph != nullptr) - { - shift_self(graph, state->box, as_rvec_array(state->x.data())); - } construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()), top.idef.iparams, top.idef.il, fr->pbcType, fr->bMolPBC, cr, state->box); - - if (graph != nullptr) - { - unshift_self(graph, state->box, as_rvec_array(state->x.data())); - } wallcycle_stop(wcycle, ewcVSITECONSTR); } diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index d93e12a8cf..57d2e66ce4 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -70,7 +70,6 @@ #include "gromacs/mdtypes/md_enums.h" #include "gromacs/mdtypes/mdatom.h" #include "gromacs/mdtypes/state.h" -#include "gromacs/pbcutil/mshift.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/mtop_lookup.h" @@ -918,7 +917,6 @@ void relax_shell_flexcon(FILE* fplog, const t_mdatoms* md, t_nrnb* nrnb, gmx_wallcycle_t wcycle, - t_graph* graph, gmx_shellfc_t* shfc, t_forcerec* fr, gmx::MdrunScheduleWorkload* runScheduleWork, @@ -989,17 +987,6 @@ void relax_shell_flexcon(FILE* fplog, */ put_atoms_in_box_omp(fr->pbcType, box, x.subArray(0, md->homenr), gmx_omp_nthreads_get(emntDefault)); - - if (graph) - { - mk_mshift(fplog, graph, fr->pbcType, box, as_rvec_array(x.data())); - } - } - - /* After this all coordinate arrays will contain whole charge groups */ - if (graph) - { - shift_self(graph, box, as_rvec_array(x.data())); } if (nflexcon) @@ -1024,12 +1011,6 @@ void relax_shell_flexcon(FILE* fplog, predict_shells(fplog, x, v, inputrec->delta_t, shells, md->massT, nullptr, bInit); } - /* do_force expected the charge groups to be in the box */ - if (graph) - { - unshift_self(graph, box, as_rvec_array(x.data())); - } - /* Calculate the forces first time around */ if (gmx_debug_at) { @@ -1038,7 +1019,7 @@ void relax_shell_flexcon(FILE* fplog, int shellfc_flags = force_flags | (bVerbose ? GMX_FORCE_ENERGY : 0); do_force(fplog, cr, ms, inputrec, nullptr, enforcedRotation, imdSession, pull_work, mdstep, nrnb, wcycle, top, box, xPadded, hist, forceWithPadding[Min], force_vir, md, enerd, - fcd, lambda, graph, fr, runScheduleWork, vsite, mu_tot, t, nullptr, + fcd, lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, (bDoNS ? GMX_FORCE_NS : 0) | shellfc_flags, ddBalanceRegionHandler); sf_dir = 0; @@ -1118,12 +1099,6 @@ void relax_shell_flexcon(FILE* fplog, /* New positions, Steepest descent */ shell_pos_sd(pos[Min], pos[Try], force[Min], shells, count); - /* do_force expected the charge groups to be in the box */ - if (graph) - { - unshift_self(graph, box, as_rvec_array(pos[Try].data())); - } - if (gmx_debug_at) { pr_rvecs(debug, 0, "RELAX: pos[Min] ", as_rvec_array(pos[Min].data()), homenr); @@ -1132,8 +1107,8 @@ void relax_shell_flexcon(FILE* fplog, /* Try the new positions */ do_force(fplog, cr, ms, inputrec, nullptr, enforcedRotation, imdSession, pull_work, 1, nrnb, wcycle, top, box, posWithPadding[Try], hist, forceWithPadding[Try], force_vir, md, - enerd, fcd, lambda, graph, fr, runScheduleWork, vsite, mu_tot, t, nullptr, - shellfc_flags, ddBalanceRegionHandler); + enerd, fcd, lambda, fr, runScheduleWork, vsite, mu_tot, t, nullptr, shellfc_flags, + ddBalanceRegionHandler); sum_epot(&(enerd->grpp), enerd->term); if (gmx_debug_at) { diff --git a/src/gromacs/mdrun/shellfc.h b/src/gromacs/mdrun/shellfc.h index cf15abbbc9..5a3f31d7e4 100644 --- a/src/gromacs/mdrun/shellfc.h +++ b/src/gromacs/mdrun/shellfc.h @@ -54,7 +54,6 @@ class history_t; struct pull_t; struct t_forcerec; struct t_fcdata; -struct t_graph; struct t_inputrec; class t_state; @@ -105,7 +104,6 @@ void relax_shell_flexcon(FILE* log, const t_mdatoms* md, t_nrnb* nrnb, gmx_wallcycle_t wcycle, - t_graph* graph, gmx_shellfc_t* shfc, t_forcerec* fr, gmx::MdrunScheduleWorkload* runScheduleWork, diff --git a/src/gromacs/mdrun/tpi.cpp b/src/gromacs/mdrun/tpi.cpp index c63a013257..e33081a46a 100644 --- a/src/gromacs/mdrun/tpi.cpp +++ b/src/gromacs/mdrun/tpi.cpp @@ -741,12 +741,7 @@ void LegacySimulator::do_tpi() clear_mat(vir); clear_mat(pres); - /* Calc energy (no forces) on new positions. - * Since we only need the intermolecular energy - * and the RF exclusion terms of the inserted molecule occur - * within a single charge group we can pass NULL for the graph. - * This also avoids shifts that would move charge groups - * out of the box. */ + /* Calc energy (no forces) on new positions. */ /* Make do_force do a single node force calculation */ cr->nnodes = 1; @@ -759,7 +754,7 @@ void LegacySimulator::do_tpi() do_force(fplog, cr, ms, inputrec, nullptr, nullptr, imdSession, pull_work, step, nrnb, wcycle, &top, state_global->box, state_global->x.arrayRefWithPadding(), &state_global->hist, f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, - state_global->lambda, nullptr, fr, runScheduleWork, nullptr, mu_tot, t, nullptr, + state_global->lambda, fr, runScheduleWork, nullptr, mu_tot, t, nullptr, GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY | (bStateChanged ? GMX_FORCE_STATECHANGED : 0), DDBalanceRegionHandler(nullptr)); std::feclearexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); diff --git a/src/gromacs/mdtypes/state.cpp b/src/gromacs/mdtypes/state.cpp index 46711fe941..61c280f997 100644 --- a/src/gromacs/mdtypes/state.cpp +++ b/src/gromacs/mdtypes/state.cpp @@ -260,7 +260,7 @@ t_state::t_state() : void set_box_rel(const t_inputrec* ir, t_state* state) { /* Make sure the box obeys the restrictions before we fix the ratios */ - correct_box(nullptr, 0, state->box, nullptr); + correct_box(nullptr, 0, state->box); clear_mat(state->box_rel); diff --git a/src/gromacs/modularsimulator/domdechelper.cpp b/src/gromacs/modularsimulator/domdechelper.cpp index e731c73bd7..67cc9d52e1 100644 --- a/src/gromacs/modularsimulator/domdechelper.cpp +++ b/src/gromacs/modularsimulator/domdechelper.cpp @@ -138,8 +138,7 @@ void DomDecHelper::run(Step step, Time gmx_unused time) { // TODO: Correcting the box is done here (if using DD) or in ForceElement (non-DD simulations). // Think about unifying this responsibility, could this be done in one place? - t_graph* graph = nullptr; - if (correct_box(fplog_, step, localState->box, graph)) + if (correct_box(fplog_, step, localState->box)) { isMasterState = true; } diff --git a/src/gromacs/modularsimulator/forceelement.cpp b/src/gromacs/modularsimulator/forceelement.cpp index 36f304b3dc..7f29d2f11a 100644 --- a/src/gromacs/modularsimulator/forceelement.cpp +++ b/src/gromacs/modularsimulator/forceelement.cpp @@ -59,7 +59,6 @@ struct gmx_edsam; struct gmx_enfrot; struct gmx_multisim_t; class history_t; -struct t_graph; namespace gmx { @@ -129,17 +128,16 @@ void ForceElement::elementSetup() void ForceElement::run(Step step, Time time, unsigned int flags) { // Disabled functionality - Awh* awh = nullptr; - gmx_edsam* ed = nullptr; - gmx_multisim_t* ms = nullptr; - t_graph* graph = nullptr; + Awh* awh = nullptr; + gmx_edsam* ed = nullptr; + gmx_multisim_t* ms = nullptr; if (!DOMAINDECOMP(cr_) && (flags & GMX_FORCE_NS) && inputrecDynamicBox(inputrec_)) { // TODO: Correcting the box is done in DomDecHelper (if using DD) or here (non-DD simulations). // Think about unifying this responsibility, could this be done in one place? auto box = statePropagatorData_->box(); - correct_box(fplog_, step, box, graph); + correct_box(fplog_, step, box); } /* The coordinates (x) are shifted (to get whole molecules) @@ -159,7 +157,7 @@ void ForceElement::run(Step step, Time time, unsigned int flags) do_force(fplog_, cr_, ms, inputrec_, awh, enforcedRotation_, imdSession_, pull_work_, step, nrnb_, wcycle_, localTopology_, box, x, hist, forces, force_vir, mdAtoms_->mdatoms(), - energyElement_->enerdata(), fcd_, lambda, graph, fr_, runScheduleWork_, vsite_, + energyElement_->enerdata(), fcd_, lambda, fr_, runScheduleWork_, vsite_, energyElement_->muTot(), time, ed, static_cast(flags), ddBalanceRegionHandler_); energyElement_->addToForceVirial(force_vir, step); } diff --git a/src/gromacs/modularsimulator/shellfcelement.cpp b/src/gromacs/modularsimulator/shellfcelement.cpp index 8eca620e4c..5d7b0ffff7 100644 --- a/src/gromacs/modularsimulator/shellfcelement.cpp +++ b/src/gromacs/modularsimulator/shellfcelement.cpp @@ -64,7 +64,6 @@ struct gmx_edsam; struct gmx_enfrot; struct gmx_multisim_t; class history_t; -struct t_graph; namespace gmx { @@ -162,15 +161,14 @@ void ShellFCElement::elementSetup() void ShellFCElement::run(Step step, Time time, bool isNSStep, unsigned int flags) { // Disabled functionality - gmx_multisim_t* ms = nullptr; - t_graph* graph = nullptr; + gmx_multisim_t* ms = nullptr; if (!DOMAINDECOMP(cr_) && isNSStep && inputrecDynamicBox(inputrec_)) { // TODO: Correcting the box is done in DomDecHelper (if using DD) or here (non-DD simulations). // Think about unifying this responsibility, could this be done in one place? auto box = statePropagatorData_->box(); - correct_box(fplog_, step, box, graph); + correct_box(fplog_, step, box); } auto x = statePropagatorData_->positionsView(); @@ -187,8 +185,8 @@ void ShellFCElement::run(Step step, Time time, bool isNSStep, unsigned int flags pull_work_, isNSStep, static_cast(flags), localTopology_, constr_, energyElement_->enerdata(), fcd_, statePropagatorData_->localNumAtoms(), x, v, box, lambda, hist, forces, force_vir, mdAtoms_->mdatoms(), nrnb_, - wcycle_, graph, shellfc_, fr_, runScheduleWork_, time, - energyElement_->muTot(), vsite_, ddBalanceRegionHandler_); + wcycle_, shellfc_, fr_, runScheduleWork_, time, energyElement_->muTot(), + vsite_, ddBalanceRegionHandler_); energyElement_->addToForceVirial(force_vir, step); } diff --git a/src/gromacs/modularsimulator/topologyholder.cpp b/src/gromacs/modularsimulator/topologyholder.cpp index 496b24dd60..d51084a243 100644 --- a/src/gromacs/modularsimulator/topologyholder.cpp +++ b/src/gromacs/modularsimulator/topologyholder.cpp @@ -62,13 +62,11 @@ TopologyHolder::TopologyHolder(const gmx_mtop_t& globalTopology, { if (!DOMAINDECOMP(cr)) { - t_graph* graph = nullptr; // Generate and initialize new topology // Note that most of the data needed for the constructor is used here - // this function should probably be simplified sooner or later. - mdAlgorithmsSetupAtomData(cr, inputrec, globalTopology, localTopology_.get(), fr, &graph, - mdAtoms, constr, vsite, nullptr); - GMX_RELEASE_ASSERT(graph == nullptr, "Graph is not implemented for the modular simulator."); + mdAlgorithmsSetupAtomData(cr, inputrec, globalTopology, localTopology_.get(), fr, mdAtoms, + constr, vsite, nullptr); } } diff --git a/src/gromacs/pbcutil/pbc.cpp b/src/gromacs/pbcutil/pbc.cpp index 61b5cdc066..a97372f88f 100644 --- a/src/gromacs/pbcutil/pbc.cpp +++ b/src/gromacs/pbcutil/pbc.cpp @@ -295,9 +295,9 @@ static int correct_box_elem(FILE* fplog, int step, tensor box, int v, int d) return shift; } -gmx_bool correct_box(FILE* fplog, int step, tensor box, t_graph* graph) +gmx_bool correct_box(FILE* fplog, int step, tensor box) { - int zy, zx, yx, i; + int zy, zx, yx; gmx_bool bCorrected; zy = correct_box_elem(fplog, step, box, ZZ, YY); @@ -306,17 +306,6 @@ gmx_bool correct_box(FILE* fplog, int step, tensor box, t_graph* graph) bCorrected = ((zy != 0) || (zx != 0) || (yx != 0)); - if (bCorrected && graph) - { - /* correct the graph */ - for (i = graph->at_start; i < graph->at_end; i++) - { - graph->ishift[i][YY] -= graph->ishift[i][ZZ] * zy; - graph->ishift[i][XX] -= graph->ishift[i][ZZ] * zx; - graph->ishift[i][XX] -= graph->ishift[i][YY] * yx; - } - } - return bCorrected; } diff --git a/src/gromacs/pbcutil/pbc.h b/src/gromacs/pbcutil/pbc.h index d7c50ff432..1a67c9b98f 100644 --- a/src/gromacs/pbcutil/pbc.h +++ b/src/gromacs/pbcutil/pbc.h @@ -132,8 +132,6 @@ enum ecenterDEF = ecenterTRIC }; -struct t_graph; - /*! \brief Returns the number of dimensions that use pbc * * \param[in] pbcType The periodic boundary condition type @@ -190,16 +188,14 @@ PbcType guessPbcType(const matrix box); /*! \brief Corrects the box if necessary * - * Checks for un-allowed box angles and corrects the box - * and the integer shift vectors in the graph (if \p graph!=NULL) if necessary. + * Checks for un-allowed box angles and corrects the box. * * \param[in] fplog File for debug output * \param[in] step The MD step number * \param[in] box The simulation cell - * \param[in] graph Information about molecular connectivity * \return TRUE when the box was corrected. */ -gmx_bool correct_box(FILE* fplog, int step, tensor box, struct t_graph* graph); +gmx_bool correct_box(FILE* fplog, int step, tensor box); /*! \brief Initiate the periodic boundary condition algorithms. * -- 2.22.0