}
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)
#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"
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<reverse_ilist_t> ril_mt;
//! \brief The size of ril_mt[?].index summed over all entries
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<int> 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 */
int nbonded_local;
gmx_reverse_top_t *rt;
- if (dd->reverse_top->bInterCGInteractions)
+ if (dd->reverse_top->bInterAtomicInteractions)
{
nzone_bondeds = zones->n;
}
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);
}
}
-/*! \brief Return a vector of the charge group index for all atoms */
-static std::vector<int> make_at2cg(const t_block &cgs)
-{
- std::vector<int> 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;
}
snew(link, 1);
- snew(link->index, ncg_mtop(mtop)+1);
+ snew(link->index, mtop->natoms + 1);
link->nalloc_a = 0;
link->a = nullptr;
continue;
}
const gmx_moltype_t &molt = mtop->moltype[molb.type];
- const t_block &cgs = molt.cgs;
- std::vector<int> 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.
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)
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)
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;
/*! \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<int> &at2cg,
gmx_bool bBCheck, gmx_bool bExcl, rvec *cg_cm,
bonded_distance_t *bd_2b,
bonded_distance_t *bd_mb)
{
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);
}
}
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);
}
}
}
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;
* 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]);
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,
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 };
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;
}
mk_graph_moltype(molt, &graph);
}
- std::vector<int> 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 */
at_offset += molt.atoms.nr;
}
- sfree(cg_cm);
sfree(xs);
if (ir->ePBC != epbcNONE)
{
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,
#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"
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);
*
* 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.
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),
{
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[] = {
}
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++)
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;
+++ /dev/null
-/*
- * 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 <cmath>
-
-#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
- }
-}
+++ /dev/null
-/*
- * 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 <cstdio>
-
-#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
void MoleculeInformation::initMolInfo()
{
- init_block(&cgs);
init_block(&mols);
init_blocka(&excls);
init_t_atoms(&atoms, 0, FALSE);
void MoleculeInformation::fullCleanUp()
{
done_atom (&atoms);
- done_block(&cgs);
done_block(&mols);
}
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,
molt.name = mol.name;
molt.atoms = mol.atoms;
/* ilists are copied later */
- molt.cgs = mol.cgs;
molt.excls = mol.excls;
pos++;
}
}
}
- 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)
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);
/* 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);
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 */
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.
#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"
}
}
-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;
}
}
}
-
-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);
- }
- }
- }
- }
-}
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);
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;
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:
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)
}
}
- push_cg(cgs, lastcg, cgnumber, nr);
-
push_atom_now(symtab, at, atomnr, atypes->atomNumberFromAtomType(type),
type, ctype, ptype, resnumberic,
resname, name, m0, q0, typeB,
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<InteractionsOfType> bondtype,
std::vector<t_blocka> 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.
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.
{
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)
{
{
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 */
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,
#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"
#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"
#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"
{
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 };
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;
{
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++)
}
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;
a_offset += molt->atoms.nr;
}
/* Done with this molecule type */
- sfree(at2cg);
}
/* Verify whether it's all correct */
{
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
#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"
{
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;
*
* 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.
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);
}
}
mtop->moltype[0].atoms = top.atoms;
- mtop->moltype[0].cgs = top.cgs;
mtop->moltype[0].excls = top.excls;
mtop->molblock.resize(1);
}
}
-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;
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),
* 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)
{
}
}
-/*! \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.
{
copyAtomtypesFromMtop(mtop, &top->atomtypes);
copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr);
- copyCgsFromMtop(mtop, &top->cgs);
copyExclsFromMtop(mtop, &top->excls);
top->name = mtop.name;
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
/* 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
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);
gmx_moltype_t::gmx_moltype_t() :
name(nullptr),
- cgs(),
excls()
{
init_t_atoms(&atoms, 0, FALSE);
gmx_moltype_t::~gmx_moltype_t()
{
done_atom(&atoms);
- done_block(&cgs);
done_blocka(&excls);
}
done_atomtypes(&(top->atomtypes));
done_symtab(&(top->symtab));
- done_block(&(top->cgs));
done_block(&(top->mols));
done_blocka(&(top->excls));
}
{
done_idef(&top->idef);
done_atom(&top->atoms);
- done_block(&top->cgs);
done_blocka(&top->excls);
done_block(&top->mols);
done_symtab(&top->symtab);
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++)
{
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));
}
}
-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];
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());
}
}
{
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);
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 */
};
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 */
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff