From c73681d683a8e2079ed029a9013f9d690ac65ea5 Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Tue, 20 May 2014 22:32:20 +0200 Subject: [PATCH] Created bonded module Converted files to C++, eliminated unused variables, reduced #include dependencies, used static_cast, made more use of position vectors in x const-correct, removed some gmx_unused arguments that were really unused. Doxygen will follow once there's been some more cleanup and there's a little more internal structure to the module. Left bonded-threading.* files behind, because they don't belong with the bonded-interaction computation machinery. They will move into the bonded module shortly. Noted TODOs for future cleanup. Change-Id: I0fae6f27adcf316cfe673ecf48688cb24dfb3469 --- cmake/gmxGCC44O3BugWorkaround.cmake | 4 +- src/gromacs/CMakeLists.txt | 3 +- src/gromacs/bonded/CMakeLists.txt | 45 +++ .../{gmxlib/bondfree.c => bonded/bonded.cpp} | 300 ++++-------------- .../bondf.h => bonded/bonded.h} | 24 +- .../{gmxlib/restcbt.c => bonded/restcbt.cpp} | 7 +- src/gromacs/{gmxlib => bonded}/restcbt.h | 9 +- src/gromacs/gmxana/anadih.c | 24 +- src/gromacs/gmxana/gmx_dipoles.cpp | 2 +- src/gromacs/gmxana/hxprops.c | 2 +- src/gromacs/gmxana/nrama.c | 6 +- src/gromacs/gmxlib/bonded-threading.cpp | 214 +++++++++++++ src/gromacs/gmxlib/disre.c | 1 - src/gromacs/gmxlib/ifunc.c | 3 +- src/gromacs/gmxlib/nonbonded/nonbonded.c | 2 +- src/gromacs/gmxpreprocess/nm2type.c | 1 - src/gromacs/gmxpreprocess/x2top.c | 2 +- src/gromacs/legacyheaders/bonded-threading.h | 58 ++++ src/gromacs/legacyheaders/force.h | 1 - src/gromacs/mdlib/domdec.c | 3 +- src/gromacs/mdlib/expanded.c | 1 - src/gromacs/mdlib/force.c | 9 +- src/gromacs/mdlib/genborn.c | 3 +- src/gromacs/mdlib/minimize.c | 2 +- src/gromacs/mdlib/qm_gamess.c | 1 - src/gromacs/mdlib/qm_gaussian.c | 1 - src/gromacs/mdlib/qm_mopac.c | 1 - src/gromacs/mdlib/qm_orca.c | 1 - src/gromacs/mdlib/qmmm.c | 1 - src/gromacs/mdlib/sim_util.c | 6 +- src/programs/mdrun/md.cpp | 2 +- 31 files changed, 433 insertions(+), 306 deletions(-) create mode 100644 src/gromacs/bonded/CMakeLists.txt rename src/gromacs/{gmxlib/bondfree.c => bonded/bonded.cpp} (94%) rename src/gromacs/{legacyheaders/bondf.h => bonded/bonded.h} (92%) rename src/gromacs/{gmxlib/restcbt.c => bonded/restcbt.cpp} (99%) rename src/gromacs/{gmxlib => bonded}/restcbt.h (98%) create mode 100644 src/gromacs/gmxlib/bonded-threading.cpp create mode 100644 src/gromacs/legacyheaders/bonded-threading.h diff --git a/cmake/gmxGCC44O3BugWorkaround.cmake b/cmake/gmxGCC44O3BugWorkaround.cmake index 0d6a42da70..6c43ba2a95 100644 --- a/cmake/gmxGCC44O3BugWorkaround.cmake +++ b/cmake/gmxGCC44O3BugWorkaround.cmake @@ -1,7 +1,7 @@ # # This file is part of the GROMACS molecular simulation package. # -# Copyright (c) 2012,2013, by the GROMACS development team, led by +# Copyright (c) 2012,2013,2014, 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. @@ -32,7 +32,7 @@ # To help us fund GROMACS development, we humbly ask that you cite # the research papers on the package. Check out http://www.gromacs.org. -# Due to a bug, gcc 4.4.x crashes when compiling bondfree.c with -O3 and +# Due to a bug, gcc 4.4.x crashes when compiling bonded/bonded.cpp with -O3 and # -fopenmp, but strangely it does not crash with -O2 + all additional options. # -O3 uses. Therefore, for the affected files, when compiling in release mode, # we override -O3 with -O2 and add the additional option. diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt index e9de38158b..e34f0166f8 100644 --- a/src/gromacs/CMakeLists.txt +++ b/src/gromacs/CMakeLists.txt @@ -82,6 +82,7 @@ endif() add_subdirectory(gmxlib) add_subdirectory(mdlib) add_subdirectory(gmxpreprocess) +add_subdirectory(bonded) add_subdirectory(commandline) add_subdirectory(fft) add_subdirectory(linearalgebra) @@ -146,7 +147,7 @@ list(APPEND LIBGROMACS_SOURCES ${GENERATED_VERSION_FILE}) # apply gcc 4.4.x bug workaround if(GMX_USE_GCC44_BUG_WORKAROUND) include(gmxGCC44O3BugWorkaround) - gmx_apply_gcc44_bug_workaround("gmxlib/bondfree.c") + gmx_apply_gcc44_bug_workaround("bonded/bonded.cpp") gmx_apply_gcc44_bug_workaround("mdlib/force.c") gmx_apply_gcc44_bug_workaround("mdlib/constr.c") endif() diff --git a/src/gromacs/bonded/CMakeLists.txt b/src/gromacs/bonded/CMakeLists.txt new file mode 100644 index 0000000000..1542ee9c6e --- /dev/null +++ b/src/gromacs/bonded/CMakeLists.txt @@ -0,0 +1,45 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright (c) 2014, by the GROMACS development team, led by +# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, +# and including many others, as listed in the AUTHORS file in the +# top-level source directory and at http://www.gromacs.org. +# +# GROMACS is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public License +# as published by the Free Software Foundation; either version 2.1 +# of the License, or (at your option) any later version. +# +# GROMACS is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with GROMACS; if not, see +# http://www.gnu.org/licenses, or write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +# +# If you want to redistribute modifications to GROMACS, please +# consider that scientific software is very special. Version +# control is crucial - bugs must be traceable. We will be happy to +# consider code for inclusion in the official distribution, but +# derived work must not be called official GROMACS. Details are found +# in the README & COPYING files - if they are missing, get the +# official version at http://www.gromacs.org. +# +# To help us fund GROMACS development, we humbly ask that you cite +# the research papers on the package. Check out http://www.gromacs.org. + +file(GLOB BONDED_SOURCES *.cpp *.c) +set(LIBGROMACS_SOURCES ${LIBGROMACS_SOURCES} ${BONDED_SOURCES} PARENT_SCOPE) + +set(BONDED_PUBLIC_HEADERS + bonded.h) + +gmx_install_headers(bonded ${BONDED_PUBLIC_HEADERS}) + +if (BUILD_TESTING) +# add_subdirectory(tests) +endif() diff --git a/src/gromacs/gmxlib/bondfree.c b/src/gromacs/bonded/bonded.cpp similarity index 94% rename from src/gromacs/gmxlib/bondfree.c rename to src/gromacs/bonded/bonded.cpp index 55f30ef299..a8757886d4 100644 --- a/src/gromacs/gmxlib/bondfree.c +++ b/src/gromacs/bonded/bonded.cpp @@ -36,16 +36,19 @@ */ #include "gmxpre.h" +#include "bonded.h" + #include "config.h" #include -#include +#include + +#include #include "gromacs/math/units.h" #include "gromacs/math/vec.h" #include "gromacs/math/utilities.h" #include "gromacs/legacyheaders/txtdump.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/ns.h" #include "gromacs/legacyheaders/macros.h" #include "gromacs/legacyheaders/names.h" @@ -53,7 +56,6 @@ #include "gromacs/legacyheaders/orires.h" #include "gromacs/legacyheaders/force.h" #include "gromacs/legacyheaders/nonbonded.h" -#include "restcbt.h" #include "gromacs/pbcutil/ishift.h" #include "gromacs/pbcutil/mshift.h" @@ -64,6 +66,8 @@ #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" +#include "restcbt.h" + /* Find a better place for this? */ const int cmap_coeff_matrix[] = { 1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4, @@ -85,7 +89,9 @@ const int cmap_coeff_matrix[] = { }; - +/* TODO This function should go and live in nonbonded.c where it is + really needed. Here, it only supports giving a fatal error message + with FENE_bonds */ int glatnr(int *global_atom_index, int i) { int atnr; @@ -102,6 +108,7 @@ int glatnr(int *global_atom_index, int i) return atnr; } +/* TODO This kind of code appears in many places. Consolidate it */ static int pbc_rvec_sub(const t_pbc *pbc, const rvec xi, const rvec xj, rvec dx) { if (pbc) @@ -352,7 +359,7 @@ real FENE_bonds(int nbonds, const real half = 0.5; const real one = 1.0; real bm, kb; - real dr, dr2, bm2, omdr2obm2, fbond, vbond, fij, vtot; + real dr2, bm2, omdr2obm2, fbond, vbond, fij, vtot; rvec dx; int i, m, ki, type, ai, aj; ivec dt; @@ -752,7 +759,6 @@ real water_pol(int nbonds, kk[YY] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_y; kk[ZZ] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_z; r_HH = 1.0/forceparams[type0].wpol.rHH; - r_OD = 1.0/forceparams[type0].wpol.rOD; if (debug) { fprintf(debug, "WPOL: qS = %10.5f aS = %5d\n", qS, aS); @@ -876,7 +882,7 @@ static real do_1_thole(const rvec xi, const rvec xj, rvec fi, rvec fj, rvec fshift[], real afac) { rvec r12; - real r12sq, r12_1, r12n, r12bar, v0, v1, fscal, ebar, fff; + real r12sq, r12_1, r12bar, v0, v1, fscal, ebar, fff; int m, t; t = pbc_rvec_sub(pbc, xi, xj, r12); /* 3 */ @@ -915,9 +921,10 @@ real thole_pol(int nbonds, int gmx_unused *global_atom_index) { /* Interaction between two pairs of particles with opposite charge */ - int i, type, a1, da1, a2, da2; - real q1, q2, qq, a, al1, al2, afac; - real V = 0; + int i, type, a1, da1, a2, da2; + real q1, q2, qq, a, al1, al2, afac; + real V = 0; + const real minusOneOnSix = -1.0/6.0; for (i = 0; (i < nbonds); ) { @@ -932,7 +939,7 @@ real thole_pol(int nbonds, al1 = forceparams[type].thole.alpha1; al2 = forceparams[type].thole.alpha2; qq = q1*q2; - afac = a*pow(al1*al2, -1.0/6.0); + afac = a*pow(al1*al2, minusOneOnSix); V += do_1_thole(x[a1], x[a2], f[a1], f[a2], pbc, qq, fshift, afac); V += do_1_thole(x[da1], x[a2], f[da1], f[a2], pbc, -qq, fshift, afac); V += do_1_thole(x[a1], x[da2], f[a1], f[da2], pbc, -qq, fshift, afac); @@ -1717,10 +1724,9 @@ do_dih_fup_noshiftf(int i, int j, int k, int l, real ddphi, rvec m, rvec n, rvec f[]) { rvec f_i, f_j, f_k, f_l; - rvec uvec, vvec, svec, dx_jl; + rvec uvec, vvec, svec; 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 */ @@ -1990,11 +1996,9 @@ pdihs_noener_simd(int nbonds, const int nfa1 = 5; int i, iu, s; int type, ai[GMX_SIMD_REAL_WIDTH], aj[GMX_SIMD_REAL_WIDTH], ak[GMX_SIMD_REAL_WIDTH], al[GMX_SIMD_REAL_WIDTH]; - int t1[GMX_SIMD_REAL_WIDTH], t2[GMX_SIMD_REAL_WIDTH], t3[GMX_SIMD_REAL_WIDTH]; - real ddphi; real dr_array[3*DIM*GMX_SIMD_REAL_WIDTH+GMX_SIMD_REAL_WIDTH], *dr; real buf_array[7*GMX_SIMD_REAL_WIDTH+GMX_SIMD_REAL_WIDTH], *buf; - real *cp, *phi0, *mult, *phi, *p, *q, *sf_i, *msf_l; + real *cp, *phi0, *mult, *p, *q; gmx_simd_real_t phi0_S, phi_S; gmx_simd_real_t mx_S, my_S, mz_S; gmx_simd_real_t nx_S, ny_S, nz_S; @@ -2015,8 +2019,6 @@ pdihs_noener_simd(int nbonds, mult = buf + 2*GMX_SIMD_REAL_WIDTH; p = buf + 3*GMX_SIMD_REAL_WIDTH; q = buf + 4*GMX_SIMD_REAL_WIDTH; - sf_i = buf + 5*GMX_SIMD_REAL_WIDTH; - msf_l = buf + 6*GMX_SIMD_REAL_WIDTH; set_pbc_simd(pbc, &pbc_simd); @@ -2163,7 +2165,7 @@ real idihs(int nbonds, dvdl_term += 0.5*(kB - kA)*dp2 - kk*dphi0*dp; - do_dih_fup(ai, aj, ak, al, (real)(-ddphi), r_ij, r_kj, r_kl, m, n, + 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 */ /* 218 TOTAL */ #ifdef DEBUG @@ -2265,8 +2267,8 @@ real fbposres(int nbonds, int i, ai, m, d, type, npbcdim = 0, fbdim; const t_iparams *pr; real vtot, kk, v; - real ref = 0, dr, dr2, rpot, rfb, rfb2, fact, invdr; - rvec com_sc, rdist, pos, dx, dpdl, fm; + real dr, dr2, rfb, rfb2, fact, invdr; + rvec com_sc, rdist, dx, dpdl, fm; gmx_bool bInvert; npbcdim = ePBC2npbcdim(ePBC); @@ -2381,12 +2383,11 @@ real posres(int nbonds, real lambda, real *dvdlambda, int refcoord_scaling, int ePBC, rvec comA, rvec comB) { - int i, ai, m, d, type, ki, npbcdim = 0; + int i, ai, m, d, type, npbcdim = 0; const t_iparams *pr; real L1; real vtot, kk, fm; - real posA, posB, ref = 0; - rvec comA_sc, comB_sc, rdist, dpdl, pos, dx; + rvec comA_sc, comB_sc, rdist, dpdl, dx; gmx_bool bForceValid = TRUE; if ((f == NULL) || (vir_diag == NULL)) /* should both be null together! */ @@ -2579,7 +2580,7 @@ real dihres(int nbonds, real vtot = 0; int ai, aj, ak, al, i, k, type, t1, t2, t3; real phi0A, phi0B, dphiA, dphiB, kfacA, kfacB, phi0, dphi, kfac; - real phi, ddphi, ddp, ddp2, dp, sign, d2r, fc, L1; + real phi, ddphi, ddp, ddp2, dp, sign, d2r, L1; rvec r_ij, r_kj, r_kl, m, n; L1 = 1.0-lambda; @@ -2686,7 +2687,6 @@ real restrangles(int nbonds, { int i, d, ai, aj, ak, type, m; int t1, t2; - rvec r_ij, r_kj; real v, vtot; ivec jt, dt_ij, dt_kj; rvec f_i, f_j, f_k; @@ -2814,7 +2814,7 @@ real restrdihs(int nbonds, t1 = pbc_rvec_sub(pbc, x[ai], x[aj], vec_temp); pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante); t2 = pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt); - t3 = pbc_rvec_sub(pbc, x[ak], x[al], vec_temp); + pbc_rvec_sub(pbc, x[ak], x[al], vec_temp); pbc_rvec_sub(pbc, x[al], x[ak], delta_post); /* This function computes factors needed for restricted angle potential. @@ -2925,7 +2925,7 @@ real cbtdihs(int nbonds, pbc_rvec_sub(pbc, x[aj], x[ai], delta_ante); t2 = pbc_rvec_sub(pbc, x[ak], x[aj], vec_temp); pbc_rvec_sub(pbc, x[ak], x[aj], delta_crnt); - t3 = pbc_rvec_sub(pbc, x[ak], x[al], vec_temp); + pbc_rvec_sub(pbc, x[ak], x[al], vec_temp); pbc_rvec_sub(pbc, x[al], x[ak], delta_post); /* \brief Compute factors for CBT potential @@ -3155,13 +3155,13 @@ real cmap_dihs(int nbonds, int t11, t21, t31, t12, t22, t32; int iphi1, ip1m1, ip1p1, ip1p2; int iphi2, ip2m1, ip2p1, ip2p2; - int l1, l2, l3, l4; - int pos1, pos2, pos3, pos4, tmp; + int l1, l2, l3; + int pos1, pos2, pos3, pos4; real ty[4], ty1[4], ty2[4], ty12[4], tc[16], tx[16]; - real phi1, psi1, cos_phi1, sin_phi1, sign1, xphi1; - real phi2, psi2, cos_phi2, sin_phi2, sign2, xphi2; - real dx, xx, tt, tu, e, df1, df2, ddf1, ddf2, ddf12, vtot; + real phi1, cos_phi1, sin_phi1, sign1, xphi1; + real phi2, cos_phi2, sin_phi2, sign2, xphi2; + real dx, xx, tt, tu, e, df1, df2, vtot; real ra21, rb21, rg21, rg1, rgr1, ra2r1, rb2r1, rabr1; real ra22, rb22, rg22, rg2, rgr2, ra2r2, rb2r2, rabr2; real fg1, hg1, fga1, hgb1, gaa1, gbb1; @@ -3223,7 +3223,7 @@ real cmap_dihs(int nbonds, b1[1] = r1_kl[2]*r1_kj[0]-r1_kl[0]*r1_kj[2]; b1[2] = r1_kl[0]*r1_kj[1]-r1_kl[1]*r1_kj[0]; /* 9 */ - tmp = pbc_rvec_sub(pbc, x[a1l], x[a1k], h1); + pbc_rvec_sub(pbc, x[a1l], x[a1k], h1); ra21 = iprod(a1, a1); /* 5 */ rb21 = iprod(b1, b1); /* 5 */ @@ -3284,7 +3284,7 @@ real cmap_dihs(int nbonds, b2[1] = r2_kl[2]*r2_kj[0]-r2_kl[0]*r2_kj[2]; b2[2] = r2_kl[0]*r2_kj[1]-r2_kl[1]*r2_kj[0]; /* 9 */ - tmp = pbc_rvec_sub(pbc, x[a2l], x[a2k], h2); + pbc_rvec_sub(pbc, x[a2l], x[a2k], h2); ra22 = iprod(a2, a2); /* 5 */ rb22 = iprod(b2, b2); /* 5 */ @@ -3349,8 +3349,8 @@ real cmap_dihs(int nbonds, dx = 2*M_PI / cmap_grid->grid_spacing; /* Where on the grid are we */ - iphi1 = (int)(xphi1/dx); - iphi2 = (int)(xphi2/dx); + iphi1 = static_cast(xphi1/dx); + iphi2 = static_cast(xphi2/dx); iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1, &ip1p1, &ip1p2); iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1, &ip2p1, &ip2p2); @@ -3415,9 +3415,6 @@ real cmap_dihs(int nbonds, e = 0; df1 = 0; df2 = 0; - ddf1 = 0; - ddf2 = 0; - ddf12 = 0; for (i = 3; i >= 0; i--) { @@ -3428,19 +3425,11 @@ real cmap_dihs(int nbonds, e = tt * e + ((tc[i*4+3]*tu+tc[i*4+2])*tu + tc[i*4+1])*tu+tc[i*4]; df1 = tu * df1 + (3.0*tc[l1]*tt+2.0*tc[l2])*tt+tc[l3]; df2 = tt * df2 + (3.0*tc[i*4+3]*tu+2.0*tc[i*4+2])*tu+tc[i*4+1]; - ddf1 = tu * ddf1 + 2.0*3.0*tc[l1]*tt+2.0*tc[l2]; - ddf2 = tt * ddf2 + 2.0*3.0*tc[4*i+3]*tu+2.0*tc[4*i+2]; } - ddf12 = tc[5] + 2.0*tc[9]*tt + 3.0*tc[13]*tt*tt + 2.0*tu*(tc[6]+2.0*tc[10]*tt+3.0*tc[14]*tt*tt) + - 3.0*tu*tu*(tc[7]+2.0*tc[11]*tt+3.0*tc[15]*tt*tt); - fac = RAD2DEG/dx; df1 = df1 * fac; df2 = df2 * fac; - ddf1 = ddf1 * fac * fac; - ddf2 = ddf2 * fac * fac; - ddf12 = ddf12 * fac * fac; /* CMAP energy */ vtot += e; @@ -3806,7 +3795,7 @@ real cross_bond_angle(int nbonds, /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003) * pp. 842-847 */ - int i, ai, aj, ak, type, m, t1, t2, t3; + int i, ai, aj, ak, type, m, t1, t2; 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; @@ -3827,7 +3816,7 @@ real cross_bond_angle(int nbonds, /* Compute distance vectors ... */ t1 = pbc_rvec_sub(pbc, x[ai], x[aj], r_ij); t2 = pbc_rvec_sub(pbc, x[ak], x[aj], r_kj); - t3 = pbc_rvec_sub(pbc, x[ai], x[ak], r_ik); + pbc_rvec_sub(pbc, x[ai], x[ak], r_ik); /* ... and their lengths */ r1 = norm(r_ij); @@ -3885,7 +3874,7 @@ static real bonded_tab(const char *type, int table_nr, { real k, tabscale, *VFtab, rt, eps, eps2, Yt, Ft, Geps, Heps2, Fp, VV, FF; int n0, nnn; - real v, f, dvdlambda; + real dvdlambda; k = (1.0 - lambda)*kA + lambda*kB; @@ -3893,7 +3882,7 @@ static real bonded_tab(const char *type, int table_nr, VFtab = table->data; rt = r*tabscale; - n0 = rt; + n0 = static_cast(rt); if (n0 >= table->n) { gmx_fatal(FARGS, "A tabulated %s interaction table number %d is out of the table range: r %f, between table indices %d and %d, table length %d", @@ -4020,7 +4009,7 @@ real tab_angles(int nbonds, if (cos_theta2 < 1) { int m; - real snt, st, sth; + real st, sth; real cik, cii, ckk; real nrkj2, nrij2; rvec f_i, f_j, f_k; @@ -4114,11 +4103,9 @@ real tab_dihs(int nbonds, return vtot; } -/* Return if this is a potential calculated in bondfree.c, - * i.e. an interaction that actually calculates a potential and - * works on multiple atoms (not e.g. a connection or a position restraint). - */ -static gmx_inline gmx_bool ftype_is_bonded_potential(int ftype) +/* TODO This function could go away when idef is not a big bucket of + everything. */ +gmx_bool ftype_is_bonded_potential(int ftype) { return (interaction_function[ftype].flags & IF_BOND) && @@ -4126,173 +4113,6 @@ static gmx_inline gmx_bool ftype_is_bonded_potential(int ftype) (ftype < F_GB12 || ftype > F_GB14); } -static void divide_bondeds_over_threads(t_idef *idef, int nthreads) -{ - int ftype; - int nat1; - int t; - int il_nr_thread; - - idef->nthreads = nthreads; - - if (F_NRE*(nthreads+1) > idef->il_thread_division_nalloc) - { - idef->il_thread_division_nalloc = F_NRE*(nthreads+1); - snew(idef->il_thread_division, idef->il_thread_division_nalloc); - } - - for (ftype = 0; ftype < F_NRE; ftype++) - { - if (ftype_is_bonded_potential(ftype)) - { - nat1 = interaction_function[ftype].nratoms + 1; - - for (t = 0; t <= nthreads; t++) - { - /* Divide the interactions equally over the threads. - * When the different types of bonded interactions - * are distributed roughly equally over the threads, - * this should lead to well localized output into - * the force buffer on each thread. - * If this is not the case, a more advanced scheme - * (not implemented yet) will do better. - */ - il_nr_thread = (((idef->il[ftype].nr/nat1)*t)/nthreads)*nat1; - - /* Ensure that distance restraint pairs with the same label - * end up on the same thread. - * This is slighlty tricky code, since the next for iteration - * may have an initial il_nr_thread lower than the final value - * in the previous iteration, but this will anyhow be increased - * to the approriate value again by this while loop. - */ - while (ftype == F_DISRES && - il_nr_thread > 0 && - il_nr_thread < idef->il[ftype].nr && - idef->iparams[idef->il[ftype].iatoms[il_nr_thread]].disres.label == - idef->iparams[idef->il[ftype].iatoms[il_nr_thread-nat1]].disres.label) - { - il_nr_thread += nat1; - } - - idef->il_thread_division[ftype*(nthreads+1)+t] = il_nr_thread; - } - } - } -} - -static unsigned -calc_bonded_reduction_mask(const t_idef *idef, - int shift, - int t, int nt) -{ - unsigned mask; - int ftype, nb, nat1, nb0, nb1, i, a; - - mask = 0; - - for (ftype = 0; ftype < F_NRE; ftype++) - { - if (ftype_is_bonded_potential(ftype)) - { - nb = idef->il[ftype].nr; - if (nb > 0) - { - nat1 = interaction_function[ftype].nratoms + 1; - - /* Divide this interaction equally over the threads. - * This is not stored: should match division in calc_bonds. - */ - nb0 = idef->il_thread_division[ftype*(nt+1)+t]; - nb1 = idef->il_thread_division[ftype*(nt+1)+t+1]; - - for (i = nb0; i < nb1; i += nat1) - { - for (a = 1; a < nat1; a++) - { - mask |= (1U << (idef->il[ftype].iatoms[i+a]>>shift)); - } - } - } - } - } - - return mask; -} - -void setup_bonded_threading(t_forcerec *fr, t_idef *idef) -{ -#define MAX_BLOCK_BITS 32 - int t; - int ctot, c, b; - - assert(fr->nthreads >= 1); - - /* Divide the bonded interaction over the threads */ - divide_bondeds_over_threads(idef, fr->nthreads); - - if (fr->nthreads == 1) - { - fr->red_nblock = 0; - - return; - } - - /* We divide the force array in a maximum of 32 blocks. - * Minimum force block reduction size is 2^6=64. - */ - fr->red_ashift = 6; - while (fr->natoms_force > (int)(MAX_BLOCK_BITS*(1U<red_ashift))) - { - fr->red_ashift++; - } - if (debug) - { - fprintf(debug, "bonded force buffer block atom shift %d bits\n", - fr->red_ashift); - } - - /* Determine to which blocks each thread's bonded force calculation - * contributes. Store this is a mask for each thread. - */ -#pragma omp parallel for num_threads(fr->nthreads) schedule(static) - for (t = 1; t < fr->nthreads; t++) - { - fr->f_t[t].red_mask = - calc_bonded_reduction_mask(idef, fr->red_ashift, t, fr->nthreads); - } - - /* Determine the maximum number of blocks we need to reduce over */ - fr->red_nblock = 0; - ctot = 0; - for (t = 0; t < fr->nthreads; t++) - { - c = 0; - for (b = 0; b < MAX_BLOCK_BITS; b++) - { - if (fr->f_t[t].red_mask & (1U<red_nblock = max(fr->red_nblock, b+1); - c++; - } - } - if (debug) - { - fprintf(debug, "thread %d flags %x count %d\n", - t, fr->f_t[t].red_mask, c); - } - ctot += c; - } - if (debug) - { - fprintf(debug, "Number of blocks to reduce: %d of size %d\n", - fr->red_nblock, 1<red_ashift); - fprintf(debug, "Reduction density %.2f density/#thread %.2f\n", - ctot*(1<red_ashift)/(double)fr->natoms_force, - ctot*(1<red_ashift)/(double)(fr->natoms_force*fr->nthreads)); - } -} - static void zero_thread_forces(f_thread_t *f_t, int n, int nblock, int blocksize) { @@ -4311,7 +4131,7 @@ static void zero_thread_forces(f_thread_t *f_t, int n, if (f_t->red_mask && (1U<f[a]); @@ -4381,7 +4201,7 @@ static void reduce_thread_force_buffer(int n, rvec *f, /* Reduce force buffers for threads that contribute */ a0 = b *block_size; a1 = (b+1)*block_size; - a1 = min(a1, n); + a1 = std::min(a1, n); for (a = a0; a < a1; a++) { for (fb = 0; fb < nfb; fb++) @@ -4452,7 +4272,7 @@ static void reduce_thread_forces(int n, rvec *f, rvec *fshift, static real calc_one_bond(int thread, int ftype, const t_idef *idef, - rvec x[], rvec f[], rvec fshift[], + const rvec x[], rvec f[], rvec fshift[], t_forcerec *fr, const t_pbc *pbc, const t_graph *g, gmx_grppairener_t *grpp, @@ -4489,7 +4309,7 @@ static real calc_one_bond(int thread, { v = cmap_dihs(nbn, iatoms+nb0, idef->iparams, &idef->cmap_grid, - (const rvec*)x, f, fshift, + x, f, fshift, pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd, global_atom_index); } @@ -4500,7 +4320,7 @@ static real calc_one_bond(int thread, /* No energies, shift forces, dvdl */ angles_noener_simd(nbn, idef->il[ftype].iatoms+nb0, idef->iparams, - (const rvec*)x, f, + x, f, pbc, g, lambda[efptFTYPE], md, fcd, global_atom_index); v = 0; @@ -4517,7 +4337,7 @@ static real calc_one_bond(int thread, #endif (nbn, idef->il[ftype].iatoms+nb0, idef->iparams, - (const rvec*)x, f, + x, f, pbc, g, lambda[efptFTYPE], md, fcd, global_atom_index); v = 0; @@ -4526,14 +4346,14 @@ static real calc_one_bond(int thread, { v = interaction_function[ftype].ifunc(nbn, iatoms+nb0, idef->iparams, - (const rvec*)x, f, fshift, + x, f, fshift, pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd, global_atom_index); } } else { - v = do_nonbonded_listed(ftype, nbn, iatoms+nb0, idef->iparams, (const rvec*)x, f, fshift, + v = do_nonbonded_listed(ftype, nbn, iatoms+nb0, idef->iparams, x, f, fshift, pbc, g, lambda, dvdl, md, fr, grpp, global_atom_index); } @@ -4547,22 +4367,20 @@ static real calc_one_bond(int thread, void calc_bonds(const gmx_multisim_t *ms, const t_idef *idef, - rvec x[], history_t *hist, + const rvec x[], history_t *hist, rvec f[], t_forcerec *fr, const t_pbc *pbc, const t_graph *g, gmx_enerdata_t *enerd, t_nrnb *nrnb, real *lambda, const t_mdatoms *md, t_fcdata *fcd, int *global_atom_index, - t_atomtypes gmx_unused *atype, gmx_genborn_t gmx_unused *born, int force_flags) { gmx_bool bCalcEnerVir; int i; - real v, dvdl[efptNR], dvdl_dum[efptNR]; /* The dummy array is to have a place to store the dhdl at other values + real dvdl[efptNR]; /* 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*/ const t_pbc *pbc_null; - char buf[22]; int thread; assert(fr->nthreads == idef->nthreads); @@ -4595,14 +4413,14 @@ void calc_bonds(const gmx_multisim_t *ms, enerd->term[F_ORIRESDEV] = calc_orires_dev(ms, idef->il[F_ORIRES].nr, idef->il[F_ORIRES].iatoms, - idef->iparams, md, (const rvec*)x, + idef->iparams, md, x, pbc_null, fcd, hist); } if (idef->il[F_DISRES].nr) { calc_disres_R_6(idef->il[F_DISRES].nr, idef->il[F_DISRES].iatoms, - idef->iparams, (const rvec*)x, pbc_null, + idef->iparams, x, pbc_null, fcd, hist); #ifdef GMX_MPI if (fcd->disres.nsystems > 1) @@ -4681,7 +4499,7 @@ void calc_bonds(const gmx_multisim_t *ms, } void calc_bonds_lambda(const t_idef *idef, - rvec x[], + const rvec x[], t_forcerec *fr, const t_pbc *pbc, const t_graph *g, gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb, @@ -4690,7 +4508,7 @@ void calc_bonds_lambda(const t_idef *idef, t_fcdata *fcd, int *global_atom_index) { - int i, ftype, nr_nonperturbed, nr; + int ftype, nr_nonperturbed, nr; real v; real dvdl_dum[efptNR]; rvec *f, *fshift; diff --git a/src/gromacs/legacyheaders/bondf.h b/src/gromacs/bonded/bonded.h similarity index 92% rename from src/gromacs/legacyheaders/bondf.h rename to src/gromacs/bonded/bonded.h index 5bd4287505..542f579767 100644 --- a/src/gromacs/legacyheaders/bondf.h +++ b/src/gromacs/bonded/bonded.h @@ -35,8 +35,8 @@ * the research papers on the package. Check out http://www.gromacs.org. */ -#ifndef _bondf_h -#define _bondf_h +#ifndef GMX_BONDED_BONDED_H +#define GMX_BONDED_BONDED_H #include @@ -58,15 +58,20 @@ int glatnr(int *global_atom_index, int i); * When global_atom_index=NULL returns i+1. */ +/*! \brief Return whether this is a potential calculated in + * bonded.cpp, i.e. an interaction that actually calculates a + * potential and works on multiple atoms (not e.g. a connection or a + * position restraint). */ +gmx_bool ftype_is_bonded_potential(int ftype); + void calc_bonds(const gmx_multisim_t *ms, const t_idef *idef, - rvec x[], history_t *hist, + const rvec x[], history_t *hist, rvec f[], t_forcerec *fr, const struct t_pbc *pbc, const struct t_graph *g, gmx_enerdata_t *enerd, t_nrnb *nrnb, real *lambda, const t_mdatoms *md, t_fcdata *fcd, int *ddgatindex, - t_atomtypes *atype, gmx_genborn_t *born, int force_flags); /* * The function calc_bonds() calculates all bonded force interactions. @@ -91,7 +96,7 @@ void calc_bonds(const gmx_multisim_t *ms, */ void calc_bonds_lambda(const t_idef *idef, - rvec x[], + const rvec x[], t_forcerec *fr, const struct t_pbc *pbc, const struct t_graph *g, gmx_grppairener_t *grpp, real *epot, t_nrnb *nrnb, @@ -154,15 +159,8 @@ t_ifunc tab_bonds, tab_angles, tab_dihs; t_ifunc polarize, anharm_polarize, water_pol, thole_pol, angres, angresz, dihres, unimplemented; -/* Divided the bonded interactions over the threads, count=fr->nthreads - * and set up the bonded thread-force buffer reduction. - * This should be called each time the bonded setup changes; - * i.e. at start-up without domain decomposition and at DD. - */ -void setup_bonded_threading(t_forcerec *fr, t_idef *idef); - #ifdef __cplusplus } #endif -#endif /* _bondf_h */ +#endif diff --git a/src/gromacs/gmxlib/restcbt.c b/src/gromacs/bonded/restcbt.cpp similarity index 99% rename from src/gromacs/gmxlib/restcbt.c rename to src/gromacs/bonded/restcbt.cpp index 64c0d00839..45b7e84c01 100644 --- a/src/gromacs/gmxlib/restcbt.c +++ b/src/gromacs/bonded/restcbt.cpp @@ -34,6 +34,8 @@ */ #include "gmxpre.h" +#include "restcbt.h" + #include #include "gromacs/math/units.h" @@ -92,7 +94,7 @@ void compute_factors_restrdihs(int type, const t_iparams forceparams[], real *prefactor_phi, real *v) { - real phi0, sine_phi0, cosine_phi0; + real phi0, cosine_phi0; real k_torsion; real c_self_ante, c_self_crnt, c_self_post; real c_cros_ante, c_cros_acrs, c_cros_post; @@ -100,12 +102,11 @@ void compute_factors_restrdihs(int type, const t_iparams forceparams[], real sine_phi_sq, cosine_phi; real delta_cosine, term_phi_phi0; real ratio_phi_ante, ratio_phi_post; - real cos_phi, norm_phi; + real norm_phi; /* Read parameters phi0 and k_torsion */ phi0 = forceparams[type].pdihs.phiA * DEG2RAD; cosine_phi0 = cos(phi0); - sine_phi0 = sin(phi0); k_torsion = forceparams[type].pdihs.cpA; /* Computation of the cosine of the dihedral angle. The scalar ("dot") product method diff --git a/src/gromacs/gmxlib/restcbt.h b/src/gromacs/bonded/restcbt.h similarity index 98% rename from src/gromacs/gmxlib/restcbt.h rename to src/gromacs/bonded/restcbt.h index 62adbda8e4..9f2490c19b 100644 --- a/src/gromacs/gmxlib/restcbt.h +++ b/src/gromacs/bonded/restcbt.h @@ -48,9 +48,12 @@ * \inlibraryapi */ -#ifndef _restcbt_h -#define _restcbt_h +#ifndef GMX_BONDED_RESTCBT_H +#define GMX_BONDED_RESTCBT_H +#include "gromacs/legacyheaders/types/simple.h" +#include "gromacs/topology/idef.h" +#include "gromacs/math/vec.h" #ifdef __cplusplus extern "C" { @@ -180,4 +183,4 @@ void compute_factors_cbtdihs(int type, const t_iparams forceparams[], } #endif -#endif /* _restcbt_h */ +#endif diff --git a/src/gromacs/gmxana/anadih.c b/src/gromacs/gmxana/anadih.c index 88785cc38a..e3daa2b289 100644 --- a/src/gromacs/gmxana/anadih.c +++ b/src/gromacs/gmxana/anadih.c @@ -45,7 +45,6 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/legacyheaders/macros.h" #include "gromacs/legacyheaders/txtdump.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/fileio/xvgr.h" #include "gromacs/legacyheaders/typedefs.h" #include "gromacs/math/vec.h" @@ -53,6 +52,7 @@ #include "gromacs/fileio/confio.h" #include "gromacs/fileio/trxio.h" +#include "gromacs/bonded/bonded.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/utility/fatalerror.h" @@ -677,7 +677,7 @@ void calc_distribution_props(int nh, int histo[], real start, *S2 = tdc*tdc+tds*tds; } -static void calc_angles(t_pbc *pbc, +static void calc_angles(struct t_pbc *pbc, int n3, atom_id index[], real ang[], rvec x_s[]) { int i, ix, t1, t2; @@ -731,7 +731,7 @@ static real calc_fraction(real angles[], int nangles) } } -static void calc_dihs(t_pbc *pbc, +static void calc_dihs(struct t_pbc *pbc, int n4, atom_id index[], real ang[], rvec x_s[]) { int i, ix, t1, t2, t3; @@ -819,15 +819,15 @@ void read_ang_dih(const char *trj_fn, real *dih[], const output_env_t oenv) { - t_pbc *pbc; - t_trxstatus *status; - int i, angind, natoms, total, teller; - int nangles, n_alloc; - real t, fraction, pifac, aa, angle; - real *angles[2]; - matrix box; - rvec *x; - int cur = 0; + struct t_pbc *pbc; + t_trxstatus *status; + int i, angind, natoms, total, teller; + int nangles, n_alloc; + real t, fraction, pifac, aa, angle; + real *angles[2]; + matrix box; + rvec *x; + int cur = 0; #define prev (1-cur) snew(pbc, 1); diff --git a/src/gromacs/gmxana/gmx_dipoles.cpp b/src/gromacs/gmxana/gmx_dipoles.cpp index 3ca706f979..2382a42d42 100644 --- a/src/gromacs/gmxana/gmx_dipoles.cpp +++ b/src/gromacs/gmxana/gmx_dipoles.cpp @@ -44,7 +44,6 @@ #include "gromacs/legacyheaders/macros.h" #include "gromacs/pbcutil/pbc.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/utility/futil.h" #include "gromacs/legacyheaders/viewit.h" #include "gromacs/legacyheaders/txtdump.h" @@ -60,6 +59,7 @@ #include "gromacs/legacyheaders/copyrite.h" #include "gromacs/fileio/trxio.h" +#include "gromacs/bonded/bonded.h" #include "gromacs/commandline/pargs.h" #include "gromacs/fileio/matio.h" #include "gromacs/fileio/xvgr.h" diff --git a/src/gromacs/gmxana/hxprops.c b/src/gromacs/gmxana/hxprops.c index c789aa964f..c47fca8dbf 100644 --- a/src/gromacs/gmxana/hxprops.c +++ b/src/gromacs/gmxana/hxprops.c @@ -48,8 +48,8 @@ #include "gromacs/topology/index.h" #include "hxprops.h" #include "gromacs/utility/smalloc.h" -#include "gromacs/legacyheaders/bondf.h" +#include "gromacs/bonded/bonded.h" #include "gromacs/utility/fatalerror.h" real ellipticity(int nres, t_bb bb[]) diff --git a/src/gromacs/gmxana/nrama.c b/src/gromacs/gmxana/nrama.c index 46373597c7..edd6d73137 100644 --- a/src/gromacs/gmxana/nrama.c +++ b/src/gromacs/gmxana/nrama.c @@ -42,11 +42,11 @@ #include #include "nrama.h" -#include "gromacs/utility/smalloc.h" -#include "gromacs/legacyheaders/typedefs.h" -#include "gromacs/legacyheaders/bondf.h" + +#include "gromacs/bonded/bonded.h" #include "gromacs/utility/futil.h" #include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/smalloc.h" #include "gromacs/pbcutil/rmpbc.h" static const char *pp_pat[] = { "C", "N", "CA", "C", "N" }; diff --git a/src/gromacs/gmxlib/bonded-threading.cpp b/src/gromacs/gmxlib/bonded-threading.cpp new file mode 100644 index 0000000000..4a69359523 --- /dev/null +++ b/src/gromacs/gmxlib/bonded-threading.cpp @@ -0,0 +1,214 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. + * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, + * and including many others, as listed in the AUTHORS file in the + * top-level source directory and at http://www.gromacs.org. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +#include "gmxpre.h" + +#include "gromacs/legacyheaders/bonded-threading.h" + +#include + +#include + +#include "gromacs/bonded/bonded.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/smalloc.h" + +static void divide_bondeds_over_threads(t_idef *idef, int nthreads) +{ + int ftype; + int nat1; + int t; + int il_nr_thread; + + idef->nthreads = nthreads; + + if (F_NRE*(nthreads+1) > idef->il_thread_division_nalloc) + { + idef->il_thread_division_nalloc = F_NRE*(nthreads+1); + snew(idef->il_thread_division, idef->il_thread_division_nalloc); + } + + for (ftype = 0; ftype < F_NRE; ftype++) + { + if (ftype_is_bonded_potential(ftype)) + { + nat1 = interaction_function[ftype].nratoms + 1; + + for (t = 0; t <= nthreads; t++) + { + /* Divide the interactions equally over the threads. + * When the different types of bonded interactions + * are distributed roughly equally over the threads, + * this should lead to well localized output into + * the force buffer on each thread. + * If this is not the case, a more advanced scheme + * (not implemented yet) will do better. + */ + il_nr_thread = (((idef->il[ftype].nr/nat1)*t)/nthreads)*nat1; + + /* Ensure that distance restraint pairs with the same label + * end up on the same thread. + * This is slighlty tricky code, since the next for iteration + * may have an initial il_nr_thread lower than the final value + * in the previous iteration, but this will anyhow be increased + * to the approriate value again by this while loop. + */ + while (ftype == F_DISRES && + il_nr_thread > 0 && + il_nr_thread < idef->il[ftype].nr && + idef->iparams[idef->il[ftype].iatoms[il_nr_thread]].disres.label == + idef->iparams[idef->il[ftype].iatoms[il_nr_thread-nat1]].disres.label) + { + il_nr_thread += nat1; + } + + idef->il_thread_division[ftype*(nthreads+1)+t] = il_nr_thread; + } + } + } +} + +static unsigned +calc_bonded_reduction_mask(const t_idef *idef, + int shift, + int t, int nt) +{ + unsigned mask; + int ftype, nb, nat1, nb0, nb1, i, a; + + mask = 0; + + for (ftype = 0; ftype < F_NRE; ftype++) + { + if (ftype_is_bonded_potential(ftype)) + { + nb = idef->il[ftype].nr; + if (nb > 0) + { + nat1 = interaction_function[ftype].nratoms + 1; + + /* Divide this interaction equally over the threads. + * This is not stored: should match division in calc_bonds. + */ + nb0 = idef->il_thread_division[ftype*(nt+1)+t]; + nb1 = idef->il_thread_division[ftype*(nt+1)+t+1]; + + for (i = nb0; i < nb1; i += nat1) + { + for (a = 1; a < nat1; a++) + { + mask |= (1U << (idef->il[ftype].iatoms[i+a]>>shift)); + } + } + } + } + } + + return mask; +} + +void setup_bonded_threading(t_forcerec *fr, t_idef *idef) +{ +#define MAX_BLOCK_BITS 32 + int t; + int ctot, c, b; + + assert(fr->nthreads >= 1); + + /* Divide the bonded interaction over the threads */ + divide_bondeds_over_threads(idef, fr->nthreads); + + if (fr->nthreads == 1) + { + fr->red_nblock = 0; + + return; + } + + /* We divide the force array in a maximum of 32 blocks. + * Minimum force block reduction size is 2^6=64. + */ + fr->red_ashift = 6; + while (fr->natoms_force > (int)(MAX_BLOCK_BITS*(1U<red_ashift))) + { + fr->red_ashift++; + } + if (debug) + { + fprintf(debug, "bonded force buffer block atom shift %d bits\n", + fr->red_ashift); + } + + /* Determine to which blocks each thread's bonded force calculation + * contributes. Store this is a mask for each thread. + */ +#pragma omp parallel for num_threads(fr->nthreads) schedule(static) + for (t = 1; t < fr->nthreads; t++) + { + fr->f_t[t].red_mask = + calc_bonded_reduction_mask(idef, fr->red_ashift, t, fr->nthreads); + } + + /* Determine the maximum number of blocks we need to reduce over */ + fr->red_nblock = 0; + ctot = 0; + for (t = 0; t < fr->nthreads; t++) + { + c = 0; + for (b = 0; b < MAX_BLOCK_BITS; b++) + { + if (fr->f_t[t].red_mask & (1U<red_nblock = std::max(fr->red_nblock, b+1); + c++; + } + } + if (debug) + { + fprintf(debug, "thread %d flags %x count %d\n", + t, fr->f_t[t].red_mask, c); + } + ctot += c; + } + if (debug) + { + fprintf(debug, "Number of blocks to reduce: %d of size %d\n", + fr->red_nblock, 1<red_ashift); + fprintf(debug, "Reduction density %.2f density/#thread %.2f\n", + ctot*(1<red_ashift)/(double)fr->natoms_force, + ctot*(1<red_ashift)/(double)(fr->natoms_force*fr->nthreads)); + } +} diff --git a/src/gromacs/gmxlib/disre.c b/src/gromacs/gmxlib/disre.c index 7bcc116fda..764e8e1396 100644 --- a/src/gromacs/gmxlib/disre.c +++ b/src/gromacs/gmxlib/disre.c @@ -46,7 +46,6 @@ #include "gromacs/legacyheaders/types/commrec.h" #include "gromacs/legacyheaders/macros.h" #include "gromacs/utility/futil.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/copyrite.h" #include "gromacs/legacyheaders/disre.h" #include "gromacs/legacyheaders/main.h" diff --git a/src/gromacs/gmxlib/ifunc.c b/src/gromacs/gmxlib/ifunc.c index d85fa96c50..342058cb5a 100644 --- a/src/gromacs/gmxlib/ifunc.c +++ b/src/gromacs/gmxlib/ifunc.c @@ -40,13 +40,12 @@ #include "config.h" +#include "gromacs/bonded/bonded.h" #include "gromacs/legacyheaders/typedefs.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/disre.h" #include "gromacs/legacyheaders/orires.h" #include "gromacs/legacyheaders/genborn.h" - #define def_bonded(str, lstr, nra, nrpa, nrpb, ind, func) \ {str, lstr, (nra), (nrpa), (nrpb), IF_BOND, (ind), (func)} diff --git a/src/gromacs/gmxlib/nonbonded/nonbonded.c b/src/gromacs/gmxlib/nonbonded/nonbonded.c index b6f56da6bd..baaf1184f2 100644 --- a/src/gromacs/gmxlib/nonbonded/nonbonded.c +++ b/src/gromacs/gmxlib/nonbonded/nonbonded.c @@ -43,6 +43,7 @@ #include "thread_mpi/threads.h" +#include "gromacs/bonded/bonded.h" #include "gromacs/legacyheaders/typedefs.h" #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/ns.h" @@ -52,7 +53,6 @@ #include "gromacs/legacyheaders/force.h" #include "gromacs/legacyheaders/names.h" #include "gromacs/legacyheaders/force.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/nrnb.h" #include "gromacs/legacyheaders/nonbonded.h" #include "gromacs/simd/simd.h" diff --git a/src/gromacs/gmxpreprocess/nm2type.c b/src/gromacs/gmxpreprocess/nm2type.c index a29def7c81..e7fe17388d 100644 --- a/src/gromacs/gmxpreprocess/nm2type.c +++ b/src/gromacs/gmxpreprocess/nm2type.c @@ -43,7 +43,6 @@ #include "gromacs/math/utilities.h" #include "gromacs/legacyheaders/macros.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/utility/cstringutil.h" #include "gromacs/utility/smalloc.h" #include "gromacs/fileio/confio.h" diff --git a/src/gromacs/gmxpreprocess/x2top.c b/src/gromacs/gmxpreprocess/x2top.c index 4a720a95cf..4157453f26 100644 --- a/src/gromacs/gmxpreprocess/x2top.c +++ b/src/gromacs/gmxpreprocess/x2top.c @@ -45,7 +45,6 @@ #include "gromacs/legacyheaders/copyrite.h" #include "gromacs/math/utilities.h" #include "gromacs/legacyheaders/macros.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/fileio/gmxfio.h" #include "gromacs/fileio/confio.h" #include "gromacs/math/units.h" @@ -59,6 +58,7 @@ #include "hackblock.h" #include "nm2type.h" +#include "gromacs/bonded/bonded.h" #include "gromacs/commandline/pargs.h" #include "gromacs/math/vec.h" #include "gromacs/pbcutil/pbc.h" diff --git a/src/gromacs/legacyheaders/bonded-threading.h b/src/gromacs/legacyheaders/bonded-threading.h new file mode 100644 index 0000000000..f1ee5c88a5 --- /dev/null +++ b/src/gromacs/legacyheaders/bonded-threading.h @@ -0,0 +1,58 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. + * Copyright (c) 2013,2014, by the GROMACS development team, led by + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, + * and including many others, as listed in the AUTHORS file in the + * top-level source directory and at http://www.gromacs.org. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ + +#ifndef GMX_LEGACYHEADERS_BONDED_THREADING_H +#define GMX_LEGACYHEADERS_BONDED_THREADING_H + +#include "typedefs.h" + +#ifdef __cplusplus +extern "C" { +#endif + +/* Divided the bonded interactions over the threads, count=fr->nthreads + * and set up the bonded thread-force buffer reduction. + * This should be called each time the bonded setup changes; + * i.e. at start-up without domain decomposition and at DD. + */ +void setup_bonded_threading(t_forcerec *fr, t_idef *idef); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/gromacs/legacyheaders/force.h b/src/gromacs/legacyheaders/force.h index c0be55dfd7..2bf5891d27 100644 --- a/src/gromacs/legacyheaders/force.h +++ b/src/gromacs/legacyheaders/force.h @@ -269,7 +269,6 @@ extern void do_force_lowlevel(t_forcerec *fr, t_fcdata *fcd, gmx_localtop_t *top, gmx_genborn_t *born, - t_atomtypes *atype, gmx_bool bBornRadii, matrix box, t_lambda *fepvals, diff --git a/src/gromacs/mdlib/domdec.c b/src/gromacs/mdlib/domdec.c index c871588c01..ff729dc017 100644 --- a/src/gromacs/mdlib/domdec.c +++ b/src/gromacs/mdlib/domdec.c @@ -44,6 +44,7 @@ #include #include +#include "gromacs/bonded/bonded.h" #include "gromacs/legacyheaders/typedefs.h" #include "gromacs/legacyheaders/network.h" #include "gromacs/math/vec.h" @@ -63,7 +64,7 @@ #include "gromacs/legacyheaders/gmx_ga2la.h" #include "gromacs/legacyheaders/macros.h" #include "nbnxn_search.h" -#include "gromacs/legacyheaders/bondf.h" +#include "gromacs/legacyheaders/bonded-threading.h" #include "gromacs/legacyheaders/gmx_omp_nthreads.h" #include "gromacs/legacyheaders/gpu_utils.h" diff --git a/src/gromacs/mdlib/expanded.c b/src/gromacs/mdlib/expanded.c index 78eba4f50a..0da2966100 100644 --- a/src/gromacs/mdlib/expanded.c +++ b/src/gromacs/mdlib/expanded.c @@ -52,7 +52,6 @@ #include "gromacs/math/units.h" #include "gromacs/legacyheaders/mdatoms.h" #include "gromacs/legacyheaders/force.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/pme.h" #include "gromacs/legacyheaders/disre.h" #include "gromacs/legacyheaders/orires.h" diff --git a/src/gromacs/mdlib/force.c b/src/gromacs/mdlib/force.c index b1ab06b301..03f4d9efd8 100644 --- a/src/gromacs/mdlib/force.c +++ b/src/gromacs/mdlib/force.c @@ -50,7 +50,6 @@ #include "gromacs/legacyheaders/network.h" #include "gromacs/legacyheaders/ns.h" #include "gromacs/legacyheaders/nrnb.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/coulomb.h" #include "gromacs/legacyheaders/pme.h" @@ -59,6 +58,7 @@ #include "gromacs/legacyheaders/qmmm.h" #include "gromacs/legacyheaders/gmx_omp_nthreads.h" +#include "gromacs/bonded/bonded.h" #include "gromacs/legacyheaders/types/commrec.h" #include "gromacs/math/vec.h" #include "gromacs/pbcutil/ishift.h" @@ -150,7 +150,6 @@ void do_force_lowlevel(t_forcerec *fr, t_inputrec *ir, t_fcdata *fcd, gmx_localtop_t *top, gmx_genborn_t *born, - t_atomtypes *atype, gmx_bool bBornRadii, matrix box, t_lambda *fepvals, @@ -378,8 +377,8 @@ void do_force_lowlevel(t_forcerec *fr, t_inputrec *ir, { wallcycle_sub_start(wcycle, ewcsBONDED); calc_bonds(cr->ms, - idef, x, hist, f, fr, &pbc, graph, enerd, nrnb, lambda, md, fcd, - DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, atype, born, + idef, (const rvec *) x, hist, f, fr, &pbc, graph, enerd, nrnb, lambda, md, fcd, + DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, flags); /* Check if we have to determine energy differences @@ -399,7 +398,7 @@ void do_force_lowlevel(t_forcerec *fr, t_inputrec *ir, { lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]); } - calc_bonds_lambda(idef, x, fr, &pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md, + calc_bonds_lambda(idef, (const rvec *) x, fr, &pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md, fcd, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL); sum_epot(&(enerd->foreign_grpp), enerd->foreign_term); enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT]; diff --git a/src/gromacs/mdlib/genborn.c b/src/gromacs/mdlib/genborn.c index 53c5128dd3..23caa377f5 100644 --- a/src/gromacs/mdlib/genborn.c +++ b/src/gromacs/mdlib/genborn.c @@ -52,7 +52,6 @@ #include "gromacs/legacyheaders/network.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/legacyheaders/nrnb.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/math/vec.h" #include "gromacs/pbcutil/ishift.h" @@ -89,7 +88,7 @@ typedef struct gbtmpnbls { int list_nalloc; } t_gbtmpnbls; -/* This function is exactly the same as the one in bondfree.c. The reason +/* This function is exactly the same as the one in bonded/bonded.cpp. The reason * it is copied here is that the bonded gb-interactions are evaluated * not in calc_bonds, but rather in calc_gb_forces */ diff --git a/src/gromacs/mdlib/minimize.c b/src/gromacs/mdlib/minimize.c index 83533b671f..5eb03f0c23 100644 --- a/src/gromacs/mdlib/minimize.c +++ b/src/gromacs/mdlib/minimize.c @@ -63,7 +63,7 @@ #include "gromacs/legacyheaders/ns.h" #include "gromacs/topology/mtop_util.h" #include "gromacs/legacyheaders/pme.h" -#include "gromacs/legacyheaders/bondf.h" +#include "gromacs/legacyheaders/bonded-threading.h" #include "gromacs/legacyheaders/gmx_omp_nthreads.h" #include "gromacs/legacyheaders/md_logging.h" diff --git a/src/gromacs/mdlib/qm_gamess.c b/src/gromacs/mdlib/qm_gamess.c index a1cda55476..9fcaed86dc 100644 --- a/src/gromacs/mdlib/qm_gamess.c +++ b/src/gromacs/mdlib/qm_gamess.c @@ -56,7 +56,6 @@ #include "gromacs/legacyheaders/network.h" #include "gromacs/legacyheaders/ns.h" #include "gromacs/legacyheaders/nrnb.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/qmmm.h" #include "gromacs/utility/fatalerror.h" diff --git a/src/gromacs/mdlib/qm_gaussian.c b/src/gromacs/mdlib/qm_gaussian.c index a24a67f325..1092a65046 100644 --- a/src/gromacs/mdlib/qm_gaussian.c +++ b/src/gromacs/mdlib/qm_gaussian.c @@ -56,7 +56,6 @@ #include "gromacs/legacyheaders/network.h" #include "gromacs/legacyheaders/ns.h" #include "gromacs/legacyheaders/nrnb.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/qmmm.h" #include "gromacs/utility/fatalerror.h" diff --git a/src/gromacs/mdlib/qm_mopac.c b/src/gromacs/mdlib/qm_mopac.c index 8e9f6a4e58..7f1bea8a83 100644 --- a/src/gromacs/mdlib/qm_mopac.c +++ b/src/gromacs/mdlib/qm_mopac.c @@ -56,7 +56,6 @@ #include "gromacs/legacyheaders/network.h" #include "gromacs/legacyheaders/ns.h" #include "gromacs/legacyheaders/nrnb.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/qmmm.h" #include "gromacs/utility/fatalerror.h" diff --git a/src/gromacs/mdlib/qm_orca.c b/src/gromacs/mdlib/qm_orca.c index c48b419b01..3b967db012 100644 --- a/src/gromacs/mdlib/qm_orca.c +++ b/src/gromacs/mdlib/qm_orca.c @@ -54,7 +54,6 @@ #include "gromacs/legacyheaders/network.h" #include "gromacs/legacyheaders/ns.h" #include "gromacs/legacyheaders/nrnb.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/qmmm.h" #include "gromacs/utility/fatalerror.h" diff --git a/src/gromacs/mdlib/qmmm.c b/src/gromacs/mdlib/qmmm.c index f6679a552a..ccb566ae9c 100644 --- a/src/gromacs/mdlib/qmmm.c +++ b/src/gromacs/mdlib/qmmm.c @@ -55,7 +55,6 @@ #include "gromacs/legacyheaders/network.h" #include "gromacs/legacyheaders/ns.h" #include "gromacs/legacyheaders/nrnb.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/legacyheaders/qmmm.h" #include "gromacs/legacyheaders/typedefs.h" diff --git a/src/gromacs/mdlib/sim_util.c b/src/gromacs/mdlib/sim_util.c index 07f89c5114..c7dc5d22b1 100644 --- a/src/gromacs/mdlib/sim_util.c +++ b/src/gromacs/mdlib/sim_util.c @@ -60,7 +60,6 @@ #include "gromacs/math/units.h" #include "gromacs/legacyheaders/mdatoms.h" #include "gromacs/legacyheaders/force.h" -#include "gromacs/legacyheaders/bondf.h" #include "gromacs/legacyheaders/pme.h" #include "gromacs/legacyheaders/disre.h" #include "gromacs/legacyheaders/orires.h" @@ -80,6 +79,7 @@ #include "../gmxlib/nonbonded/nb_kernel.h" #include "../gmxlib/nonbonded/nb_free_energy.h" +#include "gromacs/bonded/bonded.h" #include "gromacs/legacyheaders/types/commrec.h" #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_data_mgmt.h" #include "gromacs/pbcutil/ishift.h" @@ -1335,7 +1335,7 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr, do_force_lowlevel(fr, inputrec, &(top->idef), cr, nrnb, wcycle, mdatoms, x, hist, f, bSepLRF ? fr->f_twin : f, enerd, fcd, top, fr->born, - &(top->atomtypes), bBornRadii, box, + bBornRadii, box, inputrec->fepvals, lambda, graph, &(top->excls), fr->mu_tot, flags, &cycles_pme); @@ -1884,7 +1884,7 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr, do_force_lowlevel(fr, inputrec, &(top->idef), cr, nrnb, wcycle, mdatoms, x, hist, f, bSepLRF ? fr->f_twin : f, enerd, fcd, top, fr->born, - &(top->atomtypes), bBornRadii, box, + bBornRadii, box, inputrec->fepvals, lambda, graph, &(top->excls), fr->mu_tot, flags, diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp index 259484cc18..8c994cb7c2 100644 --- a/src/programs/mdrun/md.cpp +++ b/src/programs/mdrun/md.cpp @@ -74,7 +74,7 @@ #include "gromacs/legacyheaders/txtdump.h" #include "gromacs/utility/cstringutil.h" #include "pme_loadbal.h" -#include "gromacs/legacyheaders/bondf.h" +#include "gromacs/legacyheaders/bonded-threading.h" #include "membed.h" #include "gromacs/legacyheaders/types/nlistheuristics.h" #include "gromacs/legacyheaders/types/iteratedconstraints.h" -- 2.22.0