From 94209a723061bab561c456c5fe75d077dace8d0a Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Thu, 19 Sep 2019 16:59:26 +0200 Subject: [PATCH] Completely remove charge groups Charge groups remain in the .top and .tpr formats for backward compatibility. Change-Id: I9410d464dae07e0dac5bddfcc0e551c821963547 --- src/gromacs/domdec/domdec.cpp | 5 +- src/gromacs/domdec/domdec_topology.cpp | 178 ++++++-------- src/gromacs/domdec/partition.cpp | 1 - src/gromacs/fileio/tpxio.cpp | 12 +- src/gromacs/gmxana/gmx_saltbr.cpp | 49 +--- src/gromacs/gmxlib/chargegroup.cpp | 311 ------------------------ src/gromacs/gmxlib/chargegroup.h | 69 ------ src/gromacs/gmxpreprocess/grompp.cpp | 93 ------- src/gromacs/gmxpreprocess/grompp_impl.h | 2 - src/gromacs/gmxpreprocess/readir.cpp | 84 ------- src/gromacs/gmxpreprocess/readir.h | 5 - src/gromacs/gmxpreprocess/topio.cpp | 3 +- src/gromacs/gmxpreprocess/toppush.cpp | 23 +- src/gromacs/gmxpreprocess/toppush.h | 2 - src/gromacs/mdlib/constr.cpp | 98 +------- src/gromacs/mdlib/expanded.cpp | 1 - src/gromacs/mdlib/sim_util.cpp | 1 - src/gromacs/mdrun/shellfc.cpp | 23 +- src/gromacs/mdrun/tpi.cpp | 1 - src/gromacs/pbcutil/pbc.cpp | 2 +- src/gromacs/tools/convert_tpr.cpp | 4 +- src/gromacs/topology/mtop_util.cpp | 74 ------ src/gromacs/topology/mtop_util.h | 7 - src/gromacs/topology/topology.cpp | 21 +- src/gromacs/topology/topology.h | 2 - 25 files changed, 110 insertions(+), 961 deletions(-) delete mode 100644 src/gromacs/gmxlib/chargegroup.cpp delete mode 100644 src/gromacs/gmxlib/chargegroup.h diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index 6a3faa45ad..2807a97592 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -2117,8 +2117,9 @@ getSystemInfo(const gmx::MDLogger &mdlog, } else { - systemInfo.haveSplitConstraints = gmx::inter_charge_group_constraints(mtop); - systemInfo.haveSplitSettles = gmx::inter_charge_group_settles(mtop); + systemInfo.haveSplitConstraints = (gmx_mtop_ftype_count(mtop, F_CONSTR) > 0 || + gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0); + systemInfo.haveSplitSettles = (gmx_mtop_ftype_count(mtop, F_SETTLE) > 0); } if (ir.rlist == 0) diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp index 676303f47e..d27b2033c8 100644 --- a/src/gromacs/domdec/domdec_topology.cpp +++ b/src/gromacs/domdec/domdec_topology.cpp @@ -57,7 +57,6 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_network.h" #include "gromacs/domdec/ga2la.h" -#include "gromacs/gmxlib/chargegroup.h" #include "gromacs/gmxlib/network.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/forcerec.h" @@ -128,8 +127,8 @@ struct gmx_reverse_top_t bool bSettle = false; //! \brief All bonded interactions have to be assigned? bool bBCheck = false; - //! \brief Are there bondeds/exclusions between charge-groups? - bool bInterCGInteractions = false; + //! \brief Are there bondeds/exclusions between atoms? + bool bInterAtomicInteractions = false; //! \brief Reverse ilist for all moltypes std::vector ril_mt; //! \brief The size of ril_mt[?].index summed over all entries @@ -629,16 +628,16 @@ static gmx_reverse_top_t make_reverse_top(const gmx_mtop_t *mtop, gmx_bool bFE, rt.bSettle = bSettle; rt.bBCheck = bBCheck; - rt.bInterCGInteractions = mtop->bIntermolecularInteractions; + rt.bInterAtomicInteractions = mtop->bIntermolecularInteractions; rt.ril_mt.resize(mtop->moltype.size()); rt.ril_mt_tot_size = 0; std::vector nint_mt; for (size_t mt = 0; mt < mtop->moltype.size(); mt++) { const gmx_moltype_t &molt = mtop->moltype[mt]; - if (molt.cgs.nr > 1) + if (molt.atoms.nr > 1) { - rt.bInterCGInteractions = true; + rt.bInterAtomicInteractions = true; } /* Make the atom to interaction list for this molecule type */ @@ -1594,7 +1593,7 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, int nbonded_local; gmx_reverse_top_t *rt; - if (dd->reverse_top->bInterCGInteractions) + if (dd->reverse_top->bInterAtomicInteractions) { nzone_bondeds = zones->n; } @@ -1754,7 +1753,7 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones, bRCheckMB = FALSE; bRCheck2B = FALSE; - if (dd->reverse_top->bInterCGInteractions) + if (dd->reverse_top->bInterAtomicInteractions) { /* We need to check to which cell bondeds should be assigned */ rc = dd_cutoff_twobody(dd); @@ -1898,40 +1897,20 @@ static void check_link(t_blocka *link, int cg_gl, int cg_gl_j) } } -/*! \brief Return a vector of the charge group index for all atoms */ -static std::vector make_at2cg(const t_block &cgs) -{ - std::vector at2cg(cgs.index[cgs.nr]); - for (int cg = 0; cg < cgs.nr; cg++) - { - for (int a = cgs.index[cg]; a < cgs.index[cg + 1]; a++) - { - at2cg[a] = cg; - } - } - - return at2cg; -} - t_blocka *makeBondedLinks(const gmx_mtop_t *mtop, cginfo_mb_t *cginfo_mb) { t_blocka *link; cginfo_mb_t *cgi_mb; - /* For each charge group make a list of other charge groups - * in the system that a linked to it via bonded interactions + /* For each atom make a list of other atoms in the system + * that a linked to it via bonded interactions * which are also stored in reverse_top. */ reverse_ilist_t ril_intermol; if (mtop->bIntermolecularInteractions) { - if (ncg_mtop(mtop) < mtop->natoms) - { - gmx_fatal(FARGS, "The combination of intermolecular interactions, charge groups and domain decomposition is not supported. Use cutoff-scheme=Verlet (which removes the charge groups) or run without domain decomposition."); - } - t_atoms atoms; atoms.nr = mtop->natoms; @@ -1945,7 +1924,7 @@ t_blocka *makeBondedLinks(const gmx_mtop_t *mtop, } snew(link, 1); - snew(link->index, ncg_mtop(mtop)+1); + snew(link->index, mtop->natoms + 1); link->nalloc_a = 0; link->a = nullptr; @@ -1960,8 +1939,6 @@ t_blocka *makeBondedLinks(const gmx_mtop_t *mtop, continue; } const gmx_moltype_t &molt = mtop->moltype[molb.type]; - const t_block &cgs = molt.cgs; - std::vector a2c = make_at2cg(cgs); /* Make a reverse ilist in which the interactions are linked * to all atoms, not only the first atom as in gmx_reverse_top. * The constraints are discarded here. @@ -1975,65 +1952,62 @@ t_blocka *makeBondedLinks(const gmx_mtop_t *mtop, int mol; for (mol = 0; mol < (mtop->bIntermolecularInteractions ? molb.nmol : 1); mol++) { - for (int cg = 0; cg < cgs.nr; cg++) + for (int a = 0; a < molt.atoms.nr; a++) { - int cg_gl = cg_offset + cg; + int cg_gl = cg_offset + a; link->index[cg_gl+1] = link->index[cg_gl]; - for (int a = cgs.index[cg]; a < cgs.index[cg + 1]; a++) + int i = ril.index[a]; + while (i < ril.index[a+1]) { - int i = ril.index[a]; - while (i < ril.index[a+1]) + int ftype = ril.il[i++]; + int nral = NRAL(ftype); + /* Skip the ifunc index */ + i++; + for (int j = 0; j < nral; j++) { - int ftype = ril.il[i++]; - int nral = NRAL(ftype); - /* Skip the ifunc index */ - i++; - for (int j = 0; j < nral; j++) + int aj = ril.il[i + j]; + if (aj != a) { - int aj = ril.il[i + j]; - if (a2c[aj] != cg) - { - check_link(link, cg_gl, cg_offset+a2c[aj]); - } + check_link(link, cg_gl, cg_offset + aj); } - i += nral_rt(ftype); } + i += nral_rt(ftype); + } - if (mtop->bIntermolecularInteractions) + if (mtop->bIntermolecularInteractions) + { + int i = ril_intermol.index[a]; + while (i < ril_intermol.index[a+1]) { - int i = ril_intermol.index[a]; - while (i < ril_intermol.index[a+1]) + int ftype = ril_intermol.il[i++]; + int nral = NRAL(ftype); + /* Skip the ifunc index */ + i++; + for (int j = 0; j < nral; j++) { - int ftype = ril_intermol.il[i++]; - int nral = NRAL(ftype); - /* Skip the ifunc index */ - i++; - for (int j = 0; j < nral; j++) - { - /* Here we assume we have no charge groups; - * this has been checked above. - */ - int aj = ril_intermol.il[i + j]; - check_link(link, cg_gl, aj); - } - i += nral_rt(ftype); + /* Here we assume we have no charge groups; + * this has been checked above. + */ + int aj = ril_intermol.il[i + j]; + check_link(link, cg_gl, aj); } + i += nral_rt(ftype); } } if (link->index[cg_gl+1] - link->index[cg_gl] > 0) { - SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]); + SET_CGINFO_BOND_INTER(cgi_mb->cginfo[a]); ncgi++; } } - cg_offset += cgs.nr; + cg_offset += molt.atoms.nr; } - int nlink_mol = link->index[cg_offset] - link->index[cg_offset - cgs.nr]; + int nlink_mol = link->index[cg_offset] - link->index[cg_offset - molt.atoms.nr]; if (debug) { - fprintf(debug, "molecule type '%s' %d cgs has %d cg links through bonded interac.\n", *molt.name, cgs.nr, nlink_mol); + fprintf(debug, "molecule type '%s' %d atoms has %d atom links through bonded interac.\n", *molt.name, molt.atoms.nr, nlink_mol); } if (molb.nmol > mol) @@ -2043,14 +2017,14 @@ t_blocka *makeBondedLinks(const gmx_mtop_t *mtop, srenew(link->a, link->nalloc_a); for (; mol < molb.nmol; mol++) { - for (int cg = 0; cg < cgs.nr; cg++) + for (int a = 0; a < molt.atoms.nr; a++) { - int cg_gl = cg_offset + cg; + int cg_gl = cg_offset + a; link->index[cg_gl + 1] = - link->index[cg_gl + 1 - cgs.nr] + nlink_mol; + link->index[cg_gl + 1 - molt.atoms.nr] + nlink_mol; for (int j = link->index[cg_gl]; j < link->index[cg_gl+1]; j++) { - link->a[j] = link->a[j - nlink_mol] + cgs.nr; + link->a[j] = link->a[j - nlink_mol] + molt.atoms.nr; } if (link->index[cg_gl+1] - link->index[cg_gl] > 0 && cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod) @@ -2059,14 +2033,15 @@ t_blocka *makeBondedLinks(const gmx_mtop_t *mtop, ncgi++; } } - cg_offset += cgs.nr; + cg_offset += molt.atoms.nr; } } } if (debug) { - fprintf(debug, "Of the %d charge groups %d are linked via bonded interactions\n", ncg_mtop(mtop), ncgi); + fprintf(debug, "Of the %d atoms %d are linked via bonded interactions\n", + mtop->natoms, ncgi); } return link; @@ -2094,7 +2069,6 @@ static void update_max_bonded_distance(real r2, int ftype, int a1, int a2, /*! \brief Set the distance, function type and atom indices for the longest distance between charge-groups of molecule type \p molt for two-body and multi-body bonded interactions */ static void bonded_cg_distance_mol(const gmx_moltype_t *molt, - const std::vector &at2cg, gmx_bool bBCheck, gmx_bool bExcl, rvec *cg_cm, bonded_distance_t *bd_2b, bonded_distance_t *bd_mb) @@ -2111,17 +2085,16 @@ static void bonded_cg_distance_mol(const gmx_moltype_t *molt, { for (int ai = 0; ai < nral; ai++) { - int cgi = at2cg[il.iatoms[i+1+ai]]; + int atomI = il.iatoms[i + 1 + ai]; for (int aj = ai + 1; aj < nral; aj++) { - int cgj = at2cg[il.iatoms[i+1+aj]]; - if (cgi != cgj) + int atomJ = il.iatoms[i + 1 + aj]; + if (atomI != atomJ) { - real rij2 = distance2(cg_cm[cgi], cg_cm[cgj]); + real rij2 = distance2(cg_cm[atomI], cg_cm[atomJ]); update_max_bonded_distance(rij2, ftype, - il.iatoms[i+1+ai], - il.iatoms[i+1+aj], + atomI, atomJ, (nral == 2) ? bd_2b : bd_mb); } } @@ -2135,16 +2108,15 @@ static void bonded_cg_distance_mol(const gmx_moltype_t *molt, const t_blocka *excls = &molt->excls; for (int ai = 0; ai < excls->nr; ai++) { - int cgi = at2cg[ai]; for (int j = excls->index[ai]; j < excls->index[ai+1]; j++) { - int cgj = at2cg[excls->a[j]]; - if (cgi != cgj) + int aj = excls->a[j]; + if (ai != aj) { - real rij2 = distance2(cg_cm[cgi], cg_cm[cgj]); + real rij2 = distance2(cg_cm[ai], cg_cm[aj]); /* There is no function type for exclusions, use -1 */ - update_max_bonded_distance(rij2, -1, ai, excls->a[j], bd_2b); + update_max_bonded_distance(rij2, -1, ai, aj, bd_2b); } } } @@ -2215,11 +2187,11 @@ static bool moltypeHasVsite(const gmx_moltype_t &molt) return hasVsite; } -//! Compute charge group centers of mass for molecule \p molt -static void get_cgcm_mol(const gmx_moltype_t *molt, - const gmx_ffparams_t *ffparams, - int ePBC, t_graph *graph, const matrix box, - const rvec *x, rvec *xs, rvec *cg_cm) +//! Returns coordinates not broken over PBC for a molecule +static void getWholeMoleculeCoordinates(const gmx_moltype_t *molt, + const gmx_ffparams_t *ffparams, + int ePBC, t_graph *graph, const matrix box, + const rvec *x, rvec *xs) { int n, i; @@ -2241,7 +2213,7 @@ static void get_cgcm_mol(const gmx_moltype_t *molt, * unchanged, just to be 100% sure that we do not affect * binary reproducibility of simulations. */ - n = molt->cgs.index[molt->cgs.nr]; + n = molt->atoms.nr; for (i = 0; i < n; i++) { copy_rvec(x[i], xs[i]); @@ -2265,8 +2237,6 @@ static void get_cgcm_mol(const gmx_moltype_t *molt, ffparams->iparams.data(), ilist, epbcNONE, TRUE, nullptr, nullptr); } - - calc_cgcm(nullptr, 0, molt->cgs.nr, &molt->cgs, xs, cg_cm); } void dd_bonded_cg_distance(const gmx::MDLogger &mdlog, @@ -2279,7 +2249,7 @@ void dd_bonded_cg_distance(const gmx::MDLogger &mdlog, gmx_bool bExclRequired; int at_offset; t_graph graph; - rvec *xs, *cg_cm; + rvec *xs; bonded_distance_t bd_2b = { 0, -1, -1, -1 }; bonded_distance_t bd_mb = { 0, -1, -1, -1 }; @@ -2291,7 +2261,7 @@ void dd_bonded_cg_distance(const gmx::MDLogger &mdlog, for (const gmx_molblock_t &molb : mtop->molblock) { const gmx_moltype_t &molt = mtop->moltype[molb.type]; - if (molt.cgs.nr == 1 || molb.nmol == 0) + if (molt.atoms.nr == 1 || molb.nmol == 0) { at_offset += molb.nmol*molt.atoms.nr; } @@ -2302,18 +2272,16 @@ void dd_bonded_cg_distance(const gmx::MDLogger &mdlog, mk_graph_moltype(molt, &graph); } - std::vector at2cg = make_at2cg(molt.cgs); snew(xs, molt.atoms.nr); - snew(cg_cm, molt.cgs.nr); for (int mol = 0; mol < molb.nmol; mol++) { - get_cgcm_mol(&molt, &mtop->ffparams, ir->ePBC, &graph, box, - x+at_offset, xs, cg_cm); + getWholeMoleculeCoordinates(&molt, &mtop->ffparams, ir->ePBC, &graph, box, + x+at_offset, xs); bonded_distance_t bd_mol_2b = { 0, -1, -1, -1 }; bonded_distance_t bd_mol_mb = { 0, -1, -1, -1 }; - bonded_cg_distance_mol(&molt, at2cg, bBCheck, bExclRequired, cg_cm, + bonded_cg_distance_mol(&molt, bBCheck, bExclRequired, xs, &bd_mol_2b, &bd_mol_mb); /* Process the mol data adding the atom index offset */ @@ -2328,7 +2296,6 @@ void dd_bonded_cg_distance(const gmx::MDLogger &mdlog, at_offset += molt.atoms.nr; } - sfree(cg_cm); sfree(xs); if (ir->ePBC != epbcNONE) { @@ -2339,11 +2306,6 @@ void dd_bonded_cg_distance(const gmx::MDLogger &mdlog, if (mtop->bIntermolecularInteractions) { - if (ncg_mtop(mtop) < mtop->natoms) - { - gmx_fatal(FARGS, "The combination of intermolecular interactions, charge groups and domain decomposition is not supported. Use cutoff-scheme=Verlet (which removes the charge groups) or run without domain decomposition."); - } - GMX_RELEASE_ASSERT(mtop->intermolecular_ilist, "We should have an ilist when intermolecular interactions are on"); bonded_distance_intermol(*mtop->intermolecular_ilist, diff --git a/src/gromacs/domdec/partition.cpp b/src/gromacs/domdec/partition.cpp index 0508912990..f076b2e2a2 100644 --- a/src/gromacs/domdec/partition.cpp +++ b/src/gromacs/domdec/partition.cpp @@ -61,7 +61,6 @@ #include "gromacs/domdec/localatomsetmanager.h" #include "gromacs/domdec/mdsetup.h" #include "gromacs/ewald/pme.h" -#include "gromacs/gmxlib/chargegroup.h" #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nrnb.h" #include "gromacs/imd/imd.h" diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index 33279bd528..971ebab155 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -2483,7 +2483,17 @@ static void do_moltype(gmx::ISerializer *serializer, do_ilists(serializer, &molt->ilist, file_version); - do_block(serializer, &molt->cgs); + /* TODO: Remove the obsolete charge group index from the file */ + t_block cgs; + cgs.nr = molt->atoms.nr; + cgs.nalloc_index = cgs.nr + 1; + snew(cgs.index, cgs.nalloc_index); + for (int i = 0; i < cgs.nr + 1; i++) + { + cgs.index[i] = i; + } + do_block(serializer, &cgs); + sfree(cgs.index); /* This used to be in the atoms struct */ do_blocka(serializer, &molt->excls); diff --git a/src/gromacs/gmxana/gmx_saltbr.cpp b/src/gromacs/gmxana/gmx_saltbr.cpp index d7f0a8629e..b44d98d539 100644 --- a/src/gromacs/gmxana/gmx_saltbr.cpp +++ b/src/gromacs/gmxana/gmx_saltbr.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -59,28 +59,24 @@ typedef struct { real q; } t_charge; -static t_charge *mk_charge(const t_atoms *atoms, const t_block *cgs, int *nncg) +static t_charge *mk_charge(const t_atoms *atoms, int *nncg) { t_charge *cg = nullptr; char buf[32]; - int i, j, ncg, resnr, anr; + int i, ncg, resnr, anr; real qq; /* Find the charged groups */ ncg = 0; - for (i = 0; (i < cgs->nr); i++) + for (i = 0; i < atoms->nr; i++) { - qq = 0.0; - for (j = cgs->index[i]; (j < cgs->index[i+1]); j++) - { - qq += atoms->atom[j].q; - } + qq = atoms->atom[i].q; if (std::abs(qq) > 1.0e-5) { srenew(cg, ncg+1); cg[ncg].q = qq; cg[ncg].cg = i; - anr = cgs->index[i]; + anr = i; resnr = atoms->atom[anr].resind; sprintf(buf, "%s%d-%d", *(atoms->resinfo[resnr].name), @@ -96,37 +92,13 @@ static t_charge *mk_charge(const t_atoms *atoms, const t_block *cgs, int *nncg) { printf("CG: %10s Q: %6g Atoms:", cg[i].label, cg[i].q); - for (j = cgs->index[cg[i].cg]; (j < cgs->index[cg[i].cg+1]); j++) - { - printf(" %4d", j); - } + printf(" %4d", cg[i].cg); printf("\n"); } return cg; } -static real calc_dist(t_pbc *pbc, rvec x[], const t_block *cgs, int icg, int jcg) -{ - int i, j; - rvec dx; - real d2, mindist2 = 1000; - - for (i = cgs->index[icg]; (i < cgs->index[icg+1]); i++) - { - for (j = cgs->index[jcg]; (j < cgs->index[jcg+1]); j++) - { - pbc_dx(pbc, x[i], x[j], dx); - d2 = norm2(dx); - if (d2 < mindist2) - { - mindist2 = d2; - } - } - } - return std::sqrt(mindist2); -} - int gmx_saltbr(int argc, char *argv[]) { const char *desc[] = { @@ -188,7 +160,7 @@ int gmx_saltbr(int argc, char *argv[]) } top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); - cg = mk_charge(&top->atoms, &(top->cgs), &ncg); + cg = mk_charge(&top->atoms, &ncg); snew(cgdist, ncg); snew(nWithin, ncg); for (i = 0; (i < ncg); i++) @@ -213,8 +185,9 @@ int gmx_saltbr(int argc, char *argv[]) for (j = i+1; (j < ncg); j++) { srenew(cgdist[i][j], teller+1); - cgdist[i][j][teller] = - calc_dist(&pbc, x, &(top->cgs), cg[i].cg, cg[j].cg); + rvec dx; + pbc_dx(&pbc, x[cg[i].cg], x[cg[j].cg], dx); + cgdist[i][j][teller] = norm(dx); if (cgdist[i][j][teller] < truncate) { nWithin[i][j] = 1; diff --git a/src/gromacs/gmxlib/chargegroup.cpp b/src/gromacs/gmxlib/chargegroup.cpp deleted file mode 100644 index dcb2b65646..0000000000 --- a/src/gromacs/gmxlib/chargegroup.cpp +++ /dev/null @@ -1,311 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2018, 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 "chargegroup.h" - -#include - -#include "gromacs/math/vec.h" -#include "gromacs/pbcutil/pbc.h" -#include "gromacs/topology/topology.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/smalloc.h" - - -void calc_chargegroup_radii(const gmx_mtop_t *mtop, rvec *x, - real *rvdw1, real *rvdw2, - real *rcoul1, real *rcoul2) -{ - real r2v1, r2v2, r2c1, r2c2, r2; - int ntype, i, j, m, cg, a_mol, a0, a1, a; - gmx_bool *bLJ; - rvec cen; - - r2v1 = 0; - r2v2 = 0; - r2c1 = 0; - r2c2 = 0; - - ntype = mtop->ffparams.atnr; - snew(bLJ, ntype); - for (i = 0; i < ntype; i++) - { - bLJ[i] = FALSE; - for (j = 0; j < ntype; j++) - { - if (mtop->ffparams.iparams[i*ntype+j].lj.c6 != 0 || - mtop->ffparams.iparams[i*ntype+j].lj.c12 != 0) - { - bLJ[i] = TRUE; - } - } - } - - a_mol = 0; - for (const gmx_molblock_t &molb : mtop->molblock) - { - const gmx_moltype_t *molt = &mtop->moltype[molb.type]; - const t_block *cgs = &molt->cgs; - const t_atom *atom = molt->atoms.atom; - for (m = 0; m < molb.nmol; m++) - { - for (cg = 0; cg < cgs->nr; cg++) - { - a0 = cgs->index[cg]; - a1 = cgs->index[cg+1]; - if (a1 - a0 > 1) - { - clear_rvec(cen); - for (a = a0; a < a1; a++) - { - rvec_inc(cen, x[a_mol+a]); - } - svmul(1.0/(a1-a0), cen, cen); - for (a = a0; a < a1; a++) - { - r2 = distance2(cen, x[a_mol+a]); - if (r2 > r2v2 && (bLJ[atom[a].type ] || - bLJ[atom[a].typeB])) - { - if (r2 > r2v1) - { - r2v2 = r2v1; - r2v1 = r2; - } - else - { - r2v2 = r2; - } - } - if (r2 > r2c2 && - (atom[a].q != 0 || atom[a].qB != 0)) - { - if (r2 > r2c1) - { - r2c2 = r2c1; - r2c1 = r2; - } - else - { - r2c2 = r2; - } - } - } - } - } - a_mol += molt->atoms.nr; - } - } - - sfree(bLJ); - - *rvdw1 = std::sqrt(r2v1); - *rvdw2 = std::sqrt(r2v2); - *rcoul1 = std::sqrt(r2c1); - *rcoul2 = std::sqrt(r2c2); -} - -void calc_cgcm(FILE gmx_unused *fplog, int cg0, int cg1, const t_block *cgs, - rvec pos[], rvec cg_cm[]) -{ - int icg, k, k0, k1, d; - rvec cg; - real nrcg, inv_ncg; - int *cgindex; - -#ifdef DEBUG - fprintf(fplog, "Calculating centre of geometry for charge groups %d to %d\n", - cg0, cg1); -#endif - cgindex = cgs->index; - - /* Compute the center of geometry for all charge groups */ - for (icg = cg0; (icg < cg1); icg++) - { - k0 = cgindex[icg]; - k1 = cgindex[icg+1]; - nrcg = k1-k0; - if (nrcg == 1) - { - copy_rvec(pos[k0], cg_cm[icg]); - } - else - { - inv_ncg = 1.0/nrcg; - - clear_rvec(cg); - for (k = k0; (k < k1); k++) - { - for (d = 0; (d < DIM); d++) - { - cg[d] += pos[k][d]; - } - } - for (d = 0; (d < DIM); d++) - { - cg_cm[icg][d] = inv_ncg*cg[d]; - } - } - } -} - -void put_charge_groups_in_box(FILE gmx_unused *fplog, int cg0, int cg1, - int ePBC, matrix box, const t_block *cgs, - rvec pos[], rvec cg_cm[]) - -{ - int npbcdim, icg, k, k0, k1, d, e; - rvec cg; - real nrcg, inv_ncg; - int *cgindex; - gmx_bool bTric; - - if (ePBC == epbcNONE) - { - gmx_incons("Calling put_charge_groups_in_box for a system without PBC"); - } - -#ifdef DEBUG - fprintf(fplog, "Putting cgs %d to %d in box\n", cg0, cg1); -#endif - cgindex = cgs->index; - - if (ePBC == epbcXY) - { - npbcdim = 2; - } - else - { - npbcdim = 3; - } - - bTric = TRICLINIC(box); - - for (icg = cg0; (icg < cg1); icg++) - { - /* First compute the center of geometry for this charge group */ - k0 = cgindex[icg]; - k1 = cgindex[icg+1]; - nrcg = k1-k0; - - if (nrcg == 1) - { - copy_rvec(pos[k0], cg_cm[icg]); - } - else - { - inv_ncg = 1.0/nrcg; - - clear_rvec(cg); - for (k = k0; (k < k1); k++) - { - for (d = 0; d < DIM; d++) - { - cg[d] += pos[k][d]; - } - } - for (d = 0; d < DIM; d++) - { - cg_cm[icg][d] = inv_ncg*cg[d]; - } - } - /* Now check pbc for this cg */ - if (bTric) - { - for (d = npbcdim-1; d >= 0; d--) - { - while (cg_cm[icg][d] < 0) - { - for (e = d; e >= 0; e--) - { - cg_cm[icg][e] += box[d][e]; - for (k = k0; (k < k1); k++) - { - pos[k][e] += box[d][e]; - } - } - } - while (cg_cm[icg][d] >= box[d][d]) - { - for (e = d; e >= 0; e--) - { - cg_cm[icg][e] -= box[d][e]; - for (k = k0; (k < k1); k++) - { - pos[k][e] -= box[d][e]; - } - } - } - } - } - else - { - for (d = 0; d < npbcdim; d++) - { - while (cg_cm[icg][d] < 0) - { - cg_cm[icg][d] += box[d][d]; - for (k = k0; (k < k1); k++) - { - pos[k][d] += box[d][d]; - } - } - while (cg_cm[icg][d] >= box[d][d]) - { - cg_cm[icg][d] -= box[d][d]; - for (k = k0; (k < k1); k++) - { - pos[k][d] -= box[d][d]; - } - } - } - } -#ifdef DEBUG_PBC - for (d = 0; (d < npbcdim); d++) - { - if ((cg_cm[icg][d] < 0) || (cg_cm[icg][d] >= box[d][d])) - { - gmx_fatal(FARGS, "cg_cm[%d] = %15f %15f %15f\n" - "box = %15f %15f %15f\n", - icg, cg_cm[icg][XX], cg_cm[icg][YY], cg_cm[icg][ZZ], - box[XX][XX], box[YY][YY], box[ZZ][ZZ]); - } - } -#endif - } -} diff --git a/src/gromacs/gmxlib/chargegroup.h b/src/gromacs/gmxlib/chargegroup.h deleted file mode 100644 index 231b190079..0000000000 --- a/src/gromacs/gmxlib/chargegroup.h +++ /dev/null @@ -1,69 +0,0 @@ -/* - * 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) 2010,2014,2015,2018, 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_GMXLIB_CHARGEGROUP_H -#define GMX_GMXLIB_CHARGEGROUP_H - -#include - -#include "gromacs/math/vectypes.h" -#include "gromacs/utility/real.h" - -struct gmx_mtop_t; -struct t_block; - -void calc_chargegroup_radii(const gmx_mtop_t *mtop, rvec *x, - real *rvdw1, real *rvdw2, - real *rcoul1, real *rcoul2); -/* This routine calculates the two largest charge group radii in x, - * separately for VdW and Coulomb interactions. - */ - -void calc_cgcm(FILE *log, int cg0, int cg1, const t_block *cgs, - rvec pos[], rvec cg_cm[]); -/* Routine to compute centers of geometry of charge groups. No periodicity - * is used. - */ - -void put_charge_groups_in_box (FILE *log, int cg0, int cg1, - int ePBC, matrix box, const t_block *cgs, - rvec pos[], - rvec cg_cm[]); -/* This routine puts charge groups in the periodic box, keeping them - * together. - */ - -#endif diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index ba57fa1f0d..11755e631f 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -229,7 +229,6 @@ void InteractionOfType::setForceParameter(int pos, real value) void MoleculeInformation::initMolInfo() { - init_block(&cgs); init_block(&mols); init_blocka(&excls); init_t_atoms(&atoms, 0, FALSE); @@ -243,7 +242,6 @@ void MoleculeInformation::partialCleanUp() void MoleculeInformation::fullCleanUp() { done_atom (&atoms); - done_block(&cgs); done_block(&mols); } @@ -302,70 +300,6 @@ static int check_atom_names(const char *fn1, const char *fn2, return nmismatch; } -static void check_eg_vs_cg(gmx_mtop_t *mtop) -{ - int astart, m, cg, j, firstj; - unsigned char firsteg, eg; - gmx_moltype_t *molt; - - /* Go through all the charge groups and make sure all their - * atoms are in the same energy group. - */ - - astart = 0; - for (const gmx_molblock_t &molb : mtop->molblock) - { - molt = &mtop->moltype[molb.type]; - for (m = 0; m < molb.nmol; m++) - { - for (cg = 0; cg < molt->cgs.nr; cg++) - { - /* Get the energy group of the first atom in this charge group */ - firstj = astart + molt->cgs.index[cg]; - firsteg = getGroupType(mtop->groups, SimulationAtomGroupType::EnergyOutput, firstj); - for (j = molt->cgs.index[cg]+1; j < molt->cgs.index[cg+1]; j++) - { - eg = getGroupType(mtop->groups, SimulationAtomGroupType::EnergyOutput, astart+j); - if (eg != firsteg) - { - gmx_fatal(FARGS, "atoms %d and %d in charge group %d of molecule type '%s' are in different energy groups", - firstj+1, astart+j+1, cg+1, *molt->name); - } - } - } - astart += molt->atoms.nr; - } - } -} - -static void check_cg_sizes(const char *topfn, const t_block *cgs, warninp *wi) -{ - int maxsize, cg; - char warn_buf[STRLEN]; - - maxsize = 0; - for (cg = 0; cg < cgs->nr; cg++) - { - maxsize = std::max(maxsize, cgs->index[cg+1]-cgs->index[cg]); - } - - if (maxsize > MAX_CHARGEGROUP_SIZE) - { - gmx_fatal(FARGS, "The largest charge group contains %d atoms. The maximum is %d.", maxsize, MAX_CHARGEGROUP_SIZE); - } - else if (maxsize > 10) - { - set_warning_line(wi, topfn, -1); - sprintf(warn_buf, - "The largest charge group contains %d atoms.\n" - "Since atoms only see each other when the centers of geometry of the charge groups they belong to are within the cut-off distance, too large charge groups can lead to serious cut-off artifacts.\n" - "For efficiency and accuracy, charge group should consist of a few atoms.\n" - "For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc.", - maxsize); - warning_note(wi, warn_buf); - } -} - static void check_bonds_timestep(const gmx_mtop_t *mtop, double dt, warninp *wi) { /* This check is not intended to ensure accurate integration, @@ -616,7 +550,6 @@ static void molinfo2mtop(gmx::ArrayRef mi, gmx_mtop_t molt.name = mol.name; molt.atoms = mol.atoms; /* ilists are copied later */ - molt.cgs = mol.cgs; molt.excls = mol.excls; pos++; } @@ -1951,15 +1884,6 @@ int gmx_grompp(int argc, char *argv[]) } } - if (ir->cutoff_scheme == ecutsVERLET) - { - fprintf(stderr, "Removing all charge groups because cutoff-scheme=%s\n", - ecutscheme_names[ir->cutoff_scheme]); - - /* Remove all charge groups */ - gmx_mtop_remove_chargegroups(&sys); - } - if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE)) { if (ir->eI == eiCG || ir->eI == eiLBFGS) @@ -2124,11 +2048,6 @@ int gmx_grompp(int argc, char *argv[]) checkForUnboundAtoms(&sys, bVerbose, wi); - for (const gmx_moltype_t &moltype : sys.moltype) - { - check_cg_sizes(ftp2fn(efTOP, NFILE, fnm), &moltype.cgs, wi); - } - if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD) { check_bonds_timestep(&sys, ir->delta_t, wi); @@ -2227,12 +2146,6 @@ int gmx_grompp(int argc, char *argv[]) /* Init the temperature coupling state */ init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */ - if (bVerbose) - { - fprintf(stderr, "Checking consistency between energy and charge groups...\n"); - } - check_eg_vs_cg(&sys); - if (debug) { pr_symtab(debug, 0, "After index", &sys.symtab); @@ -2276,12 +2189,6 @@ int gmx_grompp(int argc, char *argv[]) clear_rvec(state.box[ZZ]); } - if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0) - { - set_warning_line(wi, mdparin, -1); - check_chargegroup_radii(&sys, ir, state.x.rvec_array(), wi); - } - if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype)) { /* Calculate the optimal grid dimensions */ diff --git a/src/gromacs/gmxpreprocess/grompp_impl.h b/src/gromacs/gmxpreprocess/grompp_impl.h index 18411cd7a3..0e38d543e0 100644 --- a/src/gromacs/gmxpreprocess/grompp_impl.h +++ b/src/gromacs/gmxpreprocess/grompp_impl.h @@ -173,8 +173,6 @@ struct MoleculeInformation bool bProcessed = false; //! Atoms in the moelcule. t_atoms atoms; - //! Charge groups in the molecule - t_block cgs; //! Molecules separated in datastructure. t_block mols; //! Exclusions in the molecule. diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index 4acad9cbcc..dbd1c61b27 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -49,7 +49,6 @@ #include "gromacs/awh/read_params.h" #include "gromacs/fileio/readinp.h" #include "gromacs/fileio/warninp.h" -#include "gromacs/gmxlib/chargegroup.h" #include "gromacs/gmxlib/network.h" #include "gromacs/gmxpreprocess/toputil.h" #include "gromacs/math/functions.h" @@ -208,11 +207,6 @@ static void check_nst(const char *desc_nst, int nst, } } -static bool ir_NVE(const t_inputrec *ir) -{ - return (EI_MD(ir->eI) && ir->etc == etcNO); -} - static int lcd(int n1, int n2) { int d, i; @@ -4222,81 +4216,3 @@ void double_check(t_inputrec *ir, matrix box, } } } - -void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir, - rvec *x, - warninp_t wi) -{ - real rvdw1, rvdw2, rcoul1, rcoul2; - char warn_buf[STRLEN]; - - calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2); - - if (rvdw1 > 0) - { - printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n", - rvdw1, rvdw2); - } - if (rcoul1 > 0) - { - printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n", - rcoul1, rcoul2); - } - - if (ir->rlist > 0) - { - if (rvdw1 + rvdw2 > ir->rlist || - rcoul1 + rcoul2 > ir->rlist) - { - sprintf(warn_buf, - "The sum of the two largest charge group radii (%f) " - "is larger than rlist (%f)\n", - std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist); - warning(wi, warn_buf); - } - else - { - /* Here we do not use the zero at cut-off macro, - * since user defined interactions might purposely - * not be zero at the cut-off. - */ - if (ir_vdw_is_zero_at_cutoff(ir) && - rvdw1 + rvdw2 > ir->rlist - ir->rvdw) - { - sprintf(warn_buf, "The sum of the two largest charge group " - "radii (%f) is larger than rlist (%f) - rvdw (%f).\n" - "With exact cut-offs, better performance can be " - "obtained with cutoff-scheme = %s, because it " - "does not use charge groups at all.", - rvdw1+rvdw2, - ir->rlist, ir->rvdw, - ecutscheme_names[ecutsVERLET]); - if (ir_NVE(ir)) - { - warning(wi, warn_buf); - } - else - { - warning_note(wi, warn_buf); - } - } - if (ir_coulomb_is_zero_at_cutoff(ir) && - rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb) - { - sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n" - "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.", - rcoul1+rcoul2, - ir->rlist, ir->rcoulomb, - ecutscheme_names[ecutsVERLET]); - if (ir_NVE(ir)) - { - warning(wi, warn_buf); - } - else - { - warning_note(wi, warn_buf); - } - } - } - } -} diff --git a/src/gromacs/gmxpreprocess/readir.h b/src/gromacs/gmxpreprocess/readir.h index f5375791fb..fff4d452f5 100644 --- a/src/gromacs/gmxpreprocess/readir.h +++ b/src/gromacs/gmxpreprocess/readir.h @@ -111,11 +111,6 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys, warninp_t wi); /* Do even more checks */ -void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir, - rvec *x, - warninp_t wi); -/* Even more checks, charge group radii vs. cut-off's only. */ - void get_ir(const char *mdparin, const char *mdparout, gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts, WriteMdpHeader writeMdpHeader, warninp_t wi); diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp index ac901b2f8c..7d9b44a077 100644 --- a/src/gromacs/gmxpreprocess/topio.cpp +++ b/src/gromacs/gmxpreprocess/topio.cpp @@ -404,7 +404,6 @@ static char **read_topol(const char *infile, const char *outfile, real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */ bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B; double qt = 0, qBt = 0; /* total charge */ - int lastcg = -1; int dcatt = -1, nmol_couple; /* File handling variables */ int status; @@ -726,7 +725,7 @@ static char **read_topol(const char *infile, const char *outfile, break; } case Directive::d_atoms: - push_atom(symtab, &(mi0->cgs), &(mi0->atoms), atypes, pline, &lastcg, wi); + push_atom(symtab, &(mi0->atoms), atypes, pline, wi); break; case Directive::d_pairs: diff --git a/src/gromacs/gmxpreprocess/toppush.cpp b/src/gromacs/gmxpreprocess/toppush.cpp index ff06eb8428..4f2a8d760d 100644 --- a/src/gromacs/gmxpreprocess/toppush.cpp +++ b/src/gromacs/gmxpreprocess/toppush.cpp @@ -1385,32 +1385,17 @@ static void push_atom_now(t_symtab *symtab, t_atoms *at, int atomnr, at->nr++; } -static void push_cg(t_block *block, int *lastindex, int index, int a) -{ - if (((block->nr) && (*lastindex != index)) || (!block->nr)) - { - /* add a new block */ - block->nr++; - srenew(block->index, block->nr+1); - } - block->index[block->nr] = a + 1; - *lastindex = index; -} - -void push_atom(t_symtab *symtab, t_block *cgs, - t_atoms *at, PreprocessingAtomTypes *atypes, char *line, int *lastcg, +void push_atom(t_symtab *symtab, + t_atoms *at, PreprocessingAtomTypes *atypes, char *line, warninp *wi) { - int nr, ptype; + int ptype; int cgnumber, atomnr, type, typeB, nscan; char id[STRLEN], ctype[STRLEN], ctypeB[STRLEN], resnumberic[STRLEN], resname[STRLEN], name[STRLEN], check[STRLEN]; double m, q, mb, qb; real m0, q0, mB, qB; - /* Make a shortcut for writing in this molecule */ - nr = at->nr; - /* Fixed parameters */ if (sscanf(line, "%s%s%s%s%s%d", id, ctype, resnumberic, resname, name, &cgnumber) != 6) @@ -1469,8 +1454,6 @@ void push_atom(t_symtab *symtab, t_block *cgs, } } - push_cg(cgs, lastcg, cgnumber, nr); - push_atom_now(symtab, at, atomnr, atypes->atomNumberFromAtomType(type), type, ctype, ptype, resnumberic, resname, name, m0, q0, typeB, diff --git a/src/gromacs/gmxpreprocess/toppush.h b/src/gromacs/gmxpreprocess/toppush.h index 626fd79f75..8982cfcd52 100644 --- a/src/gromacs/gmxpreprocess/toppush.h +++ b/src/gromacs/gmxpreprocess/toppush.h @@ -86,11 +86,9 @@ void push_nbt(Directive d, t_nbparam **nbt, PreprocessingAtomTypes *atype, warninp *wi); void push_atom(struct t_symtab *symtab, - t_block *cgs, t_atoms *at, PreprocessingAtomTypes *atype, char *line, - int *lastcg, warninp *wi); void push_bond(Directive d, gmx::ArrayRef bondtype, diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index 2e45ed4459..79158dd8b4 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -136,8 +136,6 @@ class Constraints::Impl std::vector at2con_mt; //! A list of atoms to settles for each moleculetype std::vector < std::vector < int>> at2settle_mt; - //! Whether any SETTLES cross charge-group boundaries. - bool bInterCGsettles = false; //! LINCS data. Lincs *lincsd = nullptr; // TODO this should become a unique_ptr //! SHAKE data. @@ -426,7 +424,7 @@ Constraints::Impl::apply(bool bLog, nth = 1; } - /* We do not need full pbc when constraints do not cross charge groups, + /* We do not need full pbc when constraints do not cross update groups * i.e. when dd->constraint_comm==NULL. * Note that PBC for constraints is different from PBC for bondeds. * For constraints there is both forward and backward communication. @@ -1082,7 +1080,7 @@ Constraints::Impl::Impl(const gmx_mtop_t &mtop_p, { if (DOMAINDECOMP(cr) && ddHaveSplitConstraints(*cr->dd)) { - gmx_fatal(FARGS, "SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS"); + gmx_fatal(FARGS, "SHAKE is not supported with domain decomposition and constraint that cross domain boundaries, use LINCS"); } if (nflexcon) { @@ -1102,8 +1100,6 @@ Constraints::Impl::Impl(const gmx_mtop_t &mtop_p, { please_cite(log, "Miyamoto92a"); - bInterCGsettles = inter_charge_group_settles(mtop); - settled = settle_init(mtop); /* Make an atom to settle index for use in domain decomposition */ @@ -1186,96 +1182,6 @@ ArrayRef < const std::vector < int>> Constraints::atom2settle_moltype() const return impl_->at2settle_mt; } - -bool inter_charge_group_constraints(const gmx_mtop_t &mtop) -{ - const gmx_moltype_t *molt; - const t_block *cgs; - int *at2cg, cg, a, ftype, i; - bool bInterCG; - - bInterCG = FALSE; - for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++) - { - molt = &mtop.moltype[mtop.molblock[mb].type]; - - if (molt->ilist[F_CONSTR].size() > 0 || - molt->ilist[F_CONSTRNC].size() > 0 || - molt->ilist[F_SETTLE].size() > 0) - { - cgs = &molt->cgs; - snew(at2cg, molt->atoms.nr); - for (cg = 0; cg < cgs->nr; cg++) - { - for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++) - { - at2cg[a] = cg; - } - } - - for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) - { - const InteractionList &il = molt->ilist[ftype]; - for (i = 0; i < il.size() && !bInterCG; i += 1+NRAL(ftype)) - { - if (at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 2]]) - { - bInterCG = TRUE; - } - } - } - - sfree(at2cg); - } - } - - return bInterCG; -} - -bool inter_charge_group_settles(const gmx_mtop_t &mtop) -{ - const gmx_moltype_t *molt; - const t_block *cgs; - int *at2cg, cg, a, ftype, i; - bool bInterCG; - - bInterCG = FALSE; - for (size_t mb = 0; mb < mtop.molblock.size() && !bInterCG; mb++) - { - molt = &mtop.moltype[mtop.molblock[mb].type]; - - if (molt->ilist[F_SETTLE].size() > 0) - { - cgs = &molt->cgs; - snew(at2cg, molt->atoms.nr); - for (cg = 0; cg < cgs->nr; cg++) - { - for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++) - { - at2cg[a] = cg; - } - } - - for (ftype = F_SETTLE; ftype <= F_SETTLE; ftype++) - { - const InteractionList &il = molt->ilist[ftype]; - for (i = 0; i < il.size() && !bInterCG; i += 1+NRAL(F_SETTLE)) - { - if (at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 2]] || - at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 3]]) - { - bInterCG = TRUE; - } - } - } - - sfree(at2cg); - } - } - - return bInterCG; -} - void do_constrain_first(FILE *fplog, gmx::Constraints *constr, const t_inputrec *ir, const t_mdatoms *md, int natoms, diff --git a/src/gromacs/mdlib/expanded.cpp b/src/gromacs/mdlib/expanded.cpp index 5e0c61b0fd..3f17c1078e 100644 --- a/src/gromacs/mdlib/expanded.cpp +++ b/src/gromacs/mdlib/expanded.cpp @@ -45,7 +45,6 @@ #include "gromacs/fileio/confio.h" #include "gromacs/fileio/gmxfio.h" #include "gromacs/fileio/xtcio.h" -#include "gromacs/gmxlib/chargegroup.h" #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nrnb.h" #include "gromacs/listed_forces/disre.h" diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 069ea8455f..16c891532d 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -53,7 +53,6 @@ #include "gromacs/domdec/partition.h" #include "gromacs/essentialdynamics/edsam.h" #include "gromacs/ewald/pme.h" -#include "gromacs/gmxlib/chargegroup.h" #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nonbonded/nb_free_energy.h" #include "gromacs/gmxlib/nonbonded/nb_kernel.h" diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index 876e4a71f1..42b1accb2b 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -50,7 +50,6 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_struct.h" #include "gromacs/domdec/mdsetup.h" -#include "gromacs/gmxlib/chargegroup.h" #include "gromacs/gmxlib/network.h" #include "gromacs/math/functions.h" #include "gromacs/math/units.h" @@ -309,10 +308,10 @@ gmx_shellfc_t *init_shell_flexcon(FILE *fplog, { gmx_shellfc_t *shfc; t_shell *shell; - int *shell_index = nullptr, *at2cg; + int *shell_index = nullptr; int ns, nshell, nsi; - int i, j, type, a_offset, cg, mol, ftype, nra; + int i, j, type, a_offset, mol, ftype, nra; real qS, alpha; int aS, aN = 0; /* Shell and nucleus */ int bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_ANHARM_POL, F_WATER_POL }; @@ -387,6 +386,7 @@ gmx_shellfc_t *init_shell_flexcon(FILE *fplog, ffparams = &mtop->ffparams; /* Now fill the structures */ + /* TODO: See if we can use update groups that cover shell constructions */ shfc->bInterCG = FALSE; ns = 0; a_offset = 0; @@ -394,18 +394,8 @@ gmx_shellfc_t *init_shell_flexcon(FILE *fplog, { const gmx_molblock_t *molb = &mtop->molblock[mb]; const gmx_moltype_t *molt = &mtop->moltype[molb->type]; - const t_block *cgs = &molt->cgs; - snew(at2cg, molt->atoms.nr); - for (cg = 0; cg < cgs->nr; cg++) - { - for (i = cgs->index[cg]; i < cgs->index[cg+1]; i++) - { - at2cg[i] = cg; - } - } - - const t_atom *atom = molt->atoms.atom; + const t_atom *atom = molt->atoms.atom; for (mol = 0; mol < molb->nmol; mol++) { for (j = 0; (j < NBT); j++) @@ -487,7 +477,7 @@ gmx_shellfc_t *init_shell_flexcon(FILE *fplog, } gmx_fatal(FARGS, "Can not handle more than three bonds per shell\n"); } - if (at2cg[aS] != at2cg[aN]) + if (aS != aN) { /* shell[nsi].bInterCG = TRUE; */ shfc->bInterCG = TRUE; @@ -533,7 +523,6 @@ gmx_shellfc_t *init_shell_flexcon(FILE *fplog, a_offset += molt->atoms.nr; } /* Done with this molecule type */ - sfree(at2cg); } /* Verify whether it's all correct */ @@ -581,7 +570,7 @@ gmx_shellfc_t *init_shell_flexcon(FILE *fplog, { if (fplog) { - fprintf(fplog, "\nNOTE: there are shells that are connected to particles outside their own charge group, will not predict shells positions during the run\n\n"); + fprintf(fplog, "\nNOTE: in the current version shell prediction during the crun is disabled\n\n"); } /* Prediction improves performance, so we should implement either: * 1. communication for the atoms needed for prediction diff --git a/src/gromacs/mdrun/tpi.cpp b/src/gromacs/mdrun/tpi.cpp index b7434f0529..46b185a128 100644 --- a/src/gromacs/mdrun/tpi.cpp +++ b/src/gromacs/mdrun/tpi.cpp @@ -59,7 +59,6 @@ #include "gromacs/fileio/confio.h" #include "gromacs/fileio/trxio.h" #include "gromacs/fileio/xvgr.h" -#include "gromacs/gmxlib/chargegroup.h" #include "gromacs/gmxlib/conformation_utilities.h" #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nrnb.h" diff --git a/src/gromacs/pbcutil/pbc.cpp b/src/gromacs/pbcutil/pbc.cpp index ffa41df880..7cab160d88 100644 --- a/src/gromacs/pbcutil/pbc.cpp +++ b/src/gromacs/pbcutil/pbc.cpp @@ -1585,7 +1585,7 @@ static void low_do_pbc_mtop(FILE *fplog, int ePBC, const matrix box, { const gmx_moltype_t &moltype = mtop->moltype[molb.type]; if (moltype.atoms.nr == 1 || - (!bFirst && moltype.cgs.nr == 1)) + (!bFirst && moltype.atoms.nr == 1)) { /* Just one atom or charge group in the molecule, no PBC required */ as += molb.nmol*moltype.atoms.nr; diff --git a/src/gromacs/tools/convert_tpr.cpp b/src/gromacs/tools/convert_tpr.cpp index d49f581b7e..123f0f52bb 100644 --- a/src/gromacs/tools/convert_tpr.cpp +++ b/src/gromacs/tools/convert_tpr.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -271,7 +271,6 @@ static void reduce_topology_x(int gnx, int index[], bKeep = bKeepIt(gnx, top.atoms.nr, index); invindex = invind(gnx, top.atoms.nr, index); - reduce_block(bKeep, &(top.cgs), "cgs"); reduce_block(bKeep, &(top.mols), "mols"); reduce_blocka(invindex, bKeep, &(top.excls), "excls"); reduce_rvec(gnx, index, x); @@ -301,7 +300,6 @@ static void reduce_topology_x(int gnx, int index[], } } mtop->moltype[0].atoms = top.atoms; - mtop->moltype[0].cgs = top.cgs; mtop->moltype[0].excls = top.excls; mtop->molblock.resize(1); diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp index 2016c631eb..519b1c28b7 100644 --- a/src/gromacs/topology/mtop_util.cpp +++ b/src/gromacs/topology/mtop_util.cpp @@ -169,17 +169,6 @@ void gmx_mtop_count_atomtypes(const gmx_mtop_t *mtop, int state, int typecount[] } } -int ncg_mtop(const gmx_mtop_t *mtop) -{ - int ncg = 0; - for (const gmx_molblock_t &molb : mtop->molblock) - { - ncg += molb.nmol*mtop->moltype[molb.type].cgs.nr; - } - - return ncg; -} - int gmx_mtop_num_molecules(const gmx_mtop_t &mtop) { int numMolecules = 0; @@ -200,23 +189,6 @@ int gmx_mtop_nres(const gmx_mtop_t *mtop) return nres; } -void gmx_mtop_remove_chargegroups(gmx_mtop_t *mtop) -{ - for (gmx_moltype_t &molt : mtop->moltype) - { - t_block &cgs = molt.cgs; - if (cgs.nr < molt.atoms.nr) - { - cgs.nr = molt.atoms.nr; - srenew(cgs.index, cgs.nr + 1); - for (int i = 0; i < cgs.nr + 1; i++) - { - cgs.index[i] = i; - } - } - } -} - AtomIterator::AtomIterator(const gmx_mtop_t &mtop, int globalAtomNumber) : mtop_(&mtop), mblock_(0), atoms_(&mtop.moltype[mtop.molblock[0].type].atoms), @@ -654,32 +626,6 @@ t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop) * The cat routines below are old code from src/kernel/topcat.c */ -static void blockcat(t_block *dest, const t_block *src, int copies) -{ - int i, j, l, nra, size; - - if (src->nr) - { - size = (dest->nr+copies*src->nr+1); - srenew(dest->index, size); - } - - nra = dest->index[dest->nr]; - for (l = dest->nr, j = 0; (j < copies); j++) - { - for (i = 0; (i < src->nr); i++) - { - dest->index[l++] = nra + src->index[i]; - } - if (src->nr > 0) - { - nra += src->index[src->nr]; - } - } - dest->nr += copies*src->nr; - dest->index[dest->nr] = nra; -} - static void blockacat(t_blocka *dest, const t_blocka *src, int copies, int dnum, int snum) { @@ -977,25 +923,6 @@ static void copyAtomtypesFromMtop(const gmx_mtop_t &mtop, } } -/*! \brief Copy cgs from mtop. - * - * Makes a deep copy of cgs(t_block) from gmx_mtop_t. - * Used to initialize legacy topology types. - * - * \param[in] mtop Reference to input mtop. - * \param[in] cgs Pointer to final cgs data structure. - */ -static void copyCgsFromMtop(const gmx_mtop_t &mtop, - t_block *cgs) -{ - init_block(cgs); - for (const gmx_molblock_t &molb : mtop.molblock) - { - const gmx_moltype_t &molt = mtop.moltype[molb.type]; - blockcat(cgs, &molt.cgs, molb.nmol); - } -} - /*! \brief Copy excls from mtop. * * Makes a deep copy of excls(t_blocka) from gmx_mtop_t. @@ -1170,7 +1097,6 @@ static void gen_t_topology(const gmx_mtop_t &mtop, { copyAtomtypesFromMtop(mtop, &top->atomtypes); copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr); - copyCgsFromMtop(mtop, &top->cgs); copyExclsFromMtop(mtop, &top->excls); top->name = mtop.name; diff --git a/src/gromacs/topology/mtop_util.h b/src/gromacs/topology/mtop_util.h index 566e4b1bd2..7f3aa25696 100644 --- a/src/gromacs/topology/mtop_util.h +++ b/src/gromacs/topology/mtop_util.h @@ -69,10 +69,6 @@ gmx_mtop_finalize(gmx_mtop_t *mtop); void gmx_mtop_count_atomtypes(const gmx_mtop_t *mtop, int state, int typecount[]); -/* Returns the total number of charge groups in mtop */ -int -ncg_mtop(const gmx_mtop_t *mtop); - /*!\brief Returns the total number of molecules in mtop * * \param[in] mtop The global topology @@ -82,9 +78,6 @@ int gmx_mtop_num_molecules(const gmx_mtop_t &mtop); /* Returns the total number of residues in mtop. */ int gmx_mtop_nres(const gmx_mtop_t *mtop); -/* Removes the charge groups, i.e. makes single atom charge groups, in mtop */ -void gmx_mtop_remove_chargegroups(gmx_mtop_t *mtop); - class AtomIterator; //! Proxy object returned from AtomIterator diff --git a/src/gromacs/topology/topology.cpp b/src/gromacs/topology/topology.cpp index baa0f5bf26..17c94f6ef6 100644 --- a/src/gromacs/topology/topology.cpp +++ b/src/gromacs/topology/topology.cpp @@ -80,7 +80,6 @@ void init_top(t_topology *top) init_idef(&top->idef); init_atom(&(top->atoms)); init_atomtypes(&(top->atomtypes)); - init_block(&top->cgs); init_block(&top->mols); init_blocka(&top->excls); open_symtab(&top->symtab); @@ -89,7 +88,6 @@ void init_top(t_topology *top) gmx_moltype_t::gmx_moltype_t() : name(nullptr), - cgs(), excls() { init_t_atoms(&atoms, 0, FALSE); @@ -98,7 +96,6 @@ gmx_moltype_t::gmx_moltype_t() : gmx_moltype_t::~gmx_moltype_t() { done_atom(&atoms); - done_block(&cgs); done_blocka(&excls); } @@ -126,7 +123,6 @@ void done_top(t_topology *top) done_atomtypes(&(top->atomtypes)); done_symtab(&(top->symtab)); - done_block(&(top->cgs)); done_block(&(top->mols)); done_blocka(&(top->excls)); } @@ -139,7 +135,6 @@ void done_top_mtop(t_topology *top, gmx_mtop_t *mtop) { done_idef(&top->idef); done_atom(&top->atoms); - done_block(&top->cgs); done_blocka(&top->excls); done_block(&top->mols); done_symtab(&top->symtab); @@ -307,7 +302,6 @@ static void pr_moltype(FILE *fp, int indent, const char *title, pr_indent(fp, indent); fprintf(fp, "name=\"%s\"\n", *(molt->name)); pr_atoms(fp, indent, "atoms", &(molt->atoms), bShowNumbers); - pr_block(fp, indent, "cgs", &molt->cgs, bShowNumbers); pr_blocka(fp, indent, "excls", &molt->excls, bShowNumbers); for (j = 0; (j < F_NRE); j++) { @@ -385,7 +379,6 @@ void pr_top(FILE *fp, int indent, const char *title, const t_topology *top, fprintf(fp, "name=\"%s\"\n", *(top->name)); pr_atoms(fp, indent, "atoms", &(top->atoms), bShowNumbers); pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers); - pr_block(fp, indent, "cgs", &top->cgs, bShowNumbers); pr_block(fp, indent, "mols", &top->mols, bShowNumbers); pr_str(fp, indent, "bIntermolecularInteractions", gmx::boolToString(top->bIntermolecularInteractions)); @@ -476,15 +469,6 @@ static void cmp_cmap(FILE *fp, const gmx_cmap_t *cmap1, const gmx_cmap_t *cmap2, } } -static void cmp_block(FILE *fp, const t_block *b1, const t_block *b2, const char *s) -{ - char buf[32]; - - fprintf(fp, "comparing block %s\n", s); - sprintf(buf, "%s.nr", s); - cmp_int(fp, buf, -1, b1->nr, b2->nr); -} - static void cmp_blocka(FILE *fp, const t_blocka *b1, const t_blocka *b2, const char *s) { char buf[32]; @@ -553,9 +537,7 @@ static void compareMoltypes(FILE *fp, gmx::ArrayRef mt1, gm cmp_str(fp, "Name", i, *mt1[i].name, *mt2[i].name); compareAtoms(fp, &mt1[i].atoms, &mt2[i].atoms, relativeTolerance, absoluteTolerance); compareInteractionLists(fp, &mt1[i].ilist, &mt2[i].ilist); - std::string buf = gmx::formatString("cgs[%d]", i); - cmp_block(fp, &mt1[i].cgs, &mt2[i].cgs, buf.c_str()); - buf = gmx::formatString("excls[%d]", i); + std::string buf = gmx::formatString("excls[%d]", i); cmp_blocka(fp, &mt1[i].excls, &mt2[i].excls, buf.c_str()); } } @@ -691,7 +673,6 @@ void copy_moltype(const gmx_moltype_t *src, gmx_moltype_t *dst) { dst->name = src->name; copy_blocka(&src->excls, &dst->excls); - copy_block(&src->cgs, &dst->cgs); t_atoms *atomsCopy = copy_t_atoms(&src->atoms); dst->atoms = *atomsCopy; sfree(atomsCopy); diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h index 5ed1fc9a24..5b3fc8b103 100644 --- a/src/gromacs/topology/topology.h +++ b/src/gromacs/topology/topology.h @@ -86,7 +86,6 @@ struct gmx_moltype_t char **name; /**< Name of the molecule type */ t_atoms atoms; /**< The atoms in this molecule */ InteractionLists ilist; /**< Interaction list with local indices */ - t_block cgs; /**< The charge groups */ t_blocka excls; /**< The exclusions */ }; @@ -224,7 +223,6 @@ typedef struct t_topology t_idef idef; /* The interaction function definition */ t_atoms atoms; /* The atoms */ t_atomtypes atomtypes; /* Atomtype properties */ - t_block cgs; /* The charge groups */ t_block mols; /* The molecules */ gmx_bool bIntermolecularInteractions; /* Inter.mol. int. ? */ t_blocka excls; /* The exclusions */ -- 2.22.0