#include "gromacs/topology/ifunc.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringutil.h"
/************************************************************
*
}
}
-static void mk_igraph(t_graph *g, int ftype, const t_ilist *il,
+/* Make the actual graph for an ilist, returns whether an edge was added.
+ *
+ * When a non-null part array is supplied with part indices for each atom,
+ * edges are only added when atoms have a different part index.
+ */
+static bool mk_igraph(t_graph *g, int ftype, const t_ilist *il,
int at_start, int at_end,
const int *part)
{
t_iatom *ia;
int i, j, np;
int end;
+ bool addedEdge = false;
end = il->nr;
ia = il->iatoms;
{
np = interaction_function[ftype].nratoms;
- if (ia[1] >= at_start && ia[1] < at_end)
+ if (np > 1 && ia[1] >= at_start && ia[1] < at_end)
{
if (ia[np] >= at_end)
{
/* Bond all the atoms in the settle */
add_gbond(g, ia[1], ia[2]);
add_gbond(g, ia[1], ia[3]);
+ addedEdge = true;
}
else if (part == nullptr)
{
{
add_gbond(g, ia[j], ia[j+1]);
}
+ addedEdge = true;
}
else
{
if (part[ia[j]] != part[ia[j+1]])
{
add_gbond(g, ia[j], ia[j+1]);
+ addedEdge = true;
}
}
}
ia += np+1;
i += np+1;
}
+
+ return addedEdge;
}
gmx_noreturn static void g_error(int line, const char *file)
* some atoms are not connected by the graph, which runs from
* g->at_start (>= g->at0) to g->at_end (<= g->at1).
*/
- g->at0 = at_start;
- g->at1 = at_end;
+ g->at0 = at_start;
+ g->at1 = at_end;
+ g->parts = t_graph::BondedParts::Single;
snew(nbond, at_end);
nbtot = calc_start_end(fplog, g, ilist, at_start, at_end, nbond);
* but only when they connect parts of the graph
* that are not connected through IF_CHEMBOND interactions.
*/
+ bool addedEdge = false;
for (i = 0; (i < F_NRE); i++)
{
if (!(interaction_function[i].flags & IF_CHEMBOND))
{
- mk_igraph(g, i, &(ilist[i]), at_start, at_end, nbond);
+ bool addedEdgeForType =
+ mk_igraph(g, i, &(ilist[i]), at_start, at_end, nbond);
+ addedEdge = (addedEdge || addedEdgeForType);
}
}
+
+ if (addedEdge)
+ {
+ g->parts = t_graph::BondedParts::MultipleConnected;
+ }
+ else
+ {
+ g->parts = t_graph::BondedParts::MultipleDisconnected;
+ }
}
/* Removed all the unused space from the edge array */
return -1;
}
+/* Returns the maximum length of the graph edges for coordinates x */
+static real maxEdgeLength(const t_graph g,
+ int ePBC,
+ const matrix box,
+ const rvec x[])
+{
+ t_pbc pbc;
+
+ set_pbc(&pbc, ePBC, box);
+
+ real maxEdgeLength2 = 0;
+
+ for (int node = 0; node < g.nnodes; node++)
+ {
+ for (int edge = 0; edge < g.nedge[node]; edge++)
+ {
+ int nodeJ = g.edge[node][edge];
+ rvec dx;
+ pbc_dx(&pbc, x[g.at0 + node], x[g.at0 + nodeJ], dx);
+ maxEdgeLength2 = std::max(maxEdgeLength2, norm2(dx));
+ }
+ }
+
+ return std::sqrt(maxEdgeLength2);
+}
+
void mk_mshift(FILE *log, t_graph *g, int ePBC,
const matrix box, const rvec x[])
{
}
if (nerror > 0)
{
+ /* We use a threshold of 0.25*boxSize for generating a fatal error
+ * to allow removing PBC for systems with periodic molecules.
+ *
+ * TODO: Consider a better solution for systems with periodic
+ * molecules. Ideally analysis tools should only ask to make
+ * non-periodic molecules whole in a system with periodic mols.
+ */
+ constexpr real c_relativeDistanceThreshold = 0.25;
+
+ int numPbcDimensions = ePBC2npbcdim(ePBC);
+ GMX_RELEASE_ASSERT(numPbcDimensions > 0, "Expect PBC with graph");
+ real minBoxSize = norm(box[XX]);
+ for (int d = 1; d < numPbcDimensions; d++)
+ {
+ minBoxSize = std::min(minBoxSize, norm(box[d]));
+ }
+ real maxDistance = maxEdgeLength(*g, ePBC, box, x);
+ if (maxDistance >= c_relativeDistanceThreshold*minBoxSize)
+ {
+ std::string mesg = gmx::formatString("There are inconsistent shifts over periodic boundaries in a molecule type consisting of %d atoms. The longest distance involved in such interactions is %.3f nm which is %s half the box length.",
+ g->at1 - g->at0, maxDistance,
+ maxDistance >= 0.5*minBoxSize ? "above" : "close to");
+
+ switch (g->parts)
+ {
+ case t_graph::BondedParts::MultipleConnected:
+ /* Ideally we should check if the long distances are
+ * actually between the parts, but that would require
+ * a lot of extra code.
+ */
+ mesg += " This molecule type consists of muliple parts, e.g. monomers, that are connected by interactions that are not chemical bonds, e.g. restraints. Such systems can not be treated. The only solution is increasing the box size.";
+ break;
+ default:
+ mesg += " Either you have excessively large distances between atoms in bonded interactions or your system is exploding.";
+ }
+ gmx_fatal(FARGS, mesg.c_str());
+ }
+
+ /* The most likely reason for arriving here is a periodic molecule. */
+
nerror_tot++;
if (nerror_tot <= 100)
{
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff