Introduce SimulatorBuilder

[gromacs.git] / src / gromacs / pbcutil / mshift.cpp

/*
 * 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,2016,2017,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 "mshift.h"

#include <cstring>

#include <algorithm>

#include "gromacs/math/vec.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/strconvert.h"
#include "gromacs/utility/stringutil.h"

/************************************************************
 *
 *      S H I F T   U T I L I T I E S
 *
 ************************************************************/


/************************************************************
 *
 *      G R A P H   G E N E R A T I O N    C O D E
 *
 ************************************************************/

static void add_gbond(t_graph *g, int a0, int a1)
{
    int         i;
    int         inda0, inda1;
    gmx_bool    bFound;

    inda0  = a0 - g->at_start;
    inda1  = a1 - g->at_start;
    bFound = FALSE;
    /* Search for a direct edge between a0 and a1.
     * All egdes are bidirectional, so we only need to search one way.
     */
    for (i = 0; (i < g->nedge[inda0] && !bFound); i++)
    {
        bFound = (g->edge[inda0][i] == a1);
    }

    if (!bFound)
    {
        g->edge[inda0][g->nedge[inda0]++] = a1;
        g->edge[inda1][g->nedge[inda1]++] = a0;
    }
}

/* 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.
 */
template <typename T>
static bool mk_igraph(t_graph *g, int ftype, const T &il,
                      int at_start, int at_end,
                      const int *part)
{
    int      i, j, np;
    int      end;
    bool     addedEdge = false;

    end = il.size();

    i = 0;
    while (i < end)
    {
        np = interaction_function[ftype].nratoms;

        if (np > 1 && il.iatoms[i + 1] >= at_start && il.iatoms[i + 1] < at_end)
        {
            if (il.iatoms[i + np] >= at_end)
            {
                gmx_fatal(FARGS,
                          "Molecule in topology has atom numbers below and "
                          "above natoms (%d).\n"
                          "You are probably trying to use a trajectory which does "
                          "not match the first %d atoms of the run input file.\n"
                          "You can make a matching run input file with gmx convert-tpr.",
                          at_end, at_end);
            }
            if (ftype == F_SETTLE)
            {
                /* Bond all the atoms in the settle */
                add_gbond(g, il.iatoms[i + 1], il.iatoms[i + 2]);
                add_gbond(g, il.iatoms[i + 1], il.iatoms[i + 3]);
                addedEdge = true;
            }
            else if (part == nullptr)
            {
                /* Simply add this bond */
                for (j = 1; j < np; j++)
                {
                    add_gbond(g, il.iatoms[i + j], il.iatoms[i + j+1]);
                }
                addedEdge = true;
            }
            else
            {
                /* Add this bond when it connects two unlinked parts of the graph */
                for (j = 1; j < np; j++)
                {
                    if (part[il.iatoms[i + j]] != part[il.iatoms[i + j+1]])
                    {
                        add_gbond(g, il.iatoms[i + j], il.iatoms[i + j+1]);
                        addedEdge = true;
                    }
                }
            }
        }

        i  += np+1;
    }

    return addedEdge;
}

[[noreturn]] static void g_error(int line, const char *file)
{
    gmx_fatal(FARGS, "Trying to print nonexistent graph (file %s, line %d)",
              file, line);
}
#define GCHECK(g) if ((g) == NULL) g_error(__LINE__, __FILE__)

void p_graph(FILE *log, const char *title, t_graph *g)
{
    int         i, j;
    const char *cc[egcolNR] = { "W", "G", "B" };

    GCHECK(g);
    fprintf(log, "graph:  %s\n", title);
    fprintf(log, "nnodes: %d\n", g->nnodes);
    fprintf(log, "nbound: %d\n", g->nbound);
    fprintf(log, "start:  %d\n", g->at_start);
    fprintf(log, "end:    %d\n", g->at_end);
    fprintf(log, " atom shiftx shifty shiftz C nedg    e1    e2 etc.\n");
    for (i = 0; (i < g->nnodes); i++)
    {
        if (g->nedge[i] > 0)
        {
            fprintf(log, "%5d%7d%7d%7d %1s%5d", g->at_start+i+1,
                    g->ishift[g->at_start+i][XX],
                    g->ishift[g->at_start+i][YY],
                    g->ishift[g->at_start+i][ZZ],
                    (g->negc > 0) ? cc[g->egc[i]] : " ",
                    g->nedge[i]);
            for (j = 0; (j < g->nedge[i]); j++)
            {
                fprintf(log, " %5d", g->edge[i][j]+1);
            }
            fprintf(log, "\n");
        }
    }
    fflush(log);
}

template <typename T>
static void calc_1se(t_graph *g, int ftype, const T &il,
                     int nbond[], int at_start, int at_end)
{
    int      k, nratoms, end, j;

    end = il.size();

    for (j = 0; (j < end); j += nratoms + 1)
    {
        nratoms = interaction_function[ftype].nratoms;

        if (ftype == F_SETTLE)
        {
            const int iaa = il.iatoms[j + 1];
            if (iaa >= at_start && iaa < at_end)
            {
                nbond[iaa]              += 2;
                nbond[il.iatoms[j + 2]] += 1;
                nbond[il.iatoms[j + 3]] += 1;
                g->at_start              = std::min(g->at_start, iaa);
                g->at_end                = std::max(g->at_end, iaa+2+1);
            }
        }
        else
        {
            for (k = 1; (k <= nratoms); k++)
            {
                const int iaa = il.iatoms[j + k];
                if (iaa >= at_start && iaa < at_end)
                {
                    g->at_start = std::min(g->at_start, iaa);
                    g->at_end   = std::max(g->at_end,  iaa+1);
                    /* When making the graph we (might) link all atoms in an interaction
                     * sequentially. Therefore the end atoms add 1 to the count,
                     * the middle atoms 2.
                     */
                    if (k == 1 || k == nratoms)
                    {
                        nbond[iaa] += 1;
                    }
                    else
                    {
                        nbond[iaa] += 2;
                    }
                }
            }
        }
    }
}

template <typename T>
static int calc_start_end(FILE *fplog, t_graph *g, const T il[],
                          int at_start, int at_end,
                          int nbond[])
{
    int   i, nnb, nbtot;

    g->at_start = at_end;
    g->at_end   = 0;

    /* First add all the real bonds: they should determine the molecular
     * graph.
     */
    for (i = 0; (i < F_NRE); i++)
    {
        if (interaction_function[i].flags & IF_CHEMBOND)
        {
            calc_1se(g, i, il[i], nbond, at_start, at_end);
        }
    }
    /* Then add all the other interactions in fixed lists, but first
     * check to see what's there already.
     */
    for (i = 0; (i < F_NRE); i++)
    {
        if (!(interaction_function[i].flags & IF_CHEMBOND))
        {
            calc_1se(g, i, il[i], nbond, at_start, at_end);
        }
    }

    nnb   = 0;
    nbtot = 0;
    for (i = g->at_start; (i < g->at_end); i++)
    {
        nbtot += nbond[i];
        nnb    = std::max(nnb, nbond[i]);
    }
    if (fplog)
    {
        fprintf(fplog, "Max number of connections per atom is %d\n", nnb);
        fprintf(fplog, "Total number of connections is %d\n", nbtot);
    }
    return nbtot;
}



static void compact_graph(FILE *fplog, t_graph *g)
{
    int      i, j, n, max_nedge;

    max_nedge = 0;
    n         = 0;
    for (i = 0; i < g->nnodes; i++)
    {
        for (j = 0; j < g->nedge[i]; j++)
        {
            g->edge[0][n++] = g->edge[i][j];
        }
        max_nedge = std::max(max_nedge, g->nedge[i]);
    }
    srenew(g->edge[0], n);
    /* set pointers after srenew because edge[0] might move */
    for (i = 1; i < g->nnodes; i++)
    {
        g->edge[i] = g->edge[i-1] + g->nedge[i-1];
    }

    if (fplog)
    {
        fprintf(fplog, "Max number of graph edges per atom is %d\n",
                max_nedge);
        fprintf(fplog, "Total number of graph edges is %d\n", n);
    }
}

static gmx_bool determine_graph_parts(t_graph *g, int *part)
{
    int      i, e;
    int      nchanged;
    int      at_i, *at_i2;
    gmx_bool bMultiPart;

    /* Initialize the part array with all entries different */
    for (at_i = g->at_start; at_i < g->at_end; at_i++)
    {
        part[at_i] = at_i;
    }

    /* Loop over the graph until the part array is fixed */
    do
    {
        bMultiPart = FALSE;
        nchanged   = 0;
        for (i = 0; (i < g->nnodes); i++)
        {
            at_i  = g->at_start + i;
            at_i2 = g->edge[i];
            for (e = 0; e < g->nedge[i]; e++)
            {
                /* Set part for both nodes to the minimum */
                if (part[at_i2[e]] > part[at_i])
                {
                    part[at_i2[e]] = part[at_i];
                    nchanged++;
                }
                else if (part[at_i2[e]] < part[at_i])
                {
                    part[at_i] = part[at_i2[e]];
                    nchanged++;
                }
            }
            if (part[at_i] != part[g->at_start])
            {
                bMultiPart = TRUE;
            }
        }
        if (debug)
        {
            fprintf(debug, "graph part[] nchanged=%d, bMultiPart=%s\n",
                    nchanged, gmx::boolToString(bMultiPart));
        }
    }
    while (nchanged > 0);

    return bMultiPart;
}

template <typename T>
static void
mk_graph_ilist(FILE *fplog,
               const T *ilist, int at_start, int at_end,
               gmx_bool bShakeOnly, gmx_bool bSettle,
               t_graph *g)
{
    int        *nbond;
    int         i, nbtot;
    gmx_bool    bMultiPart;

    /* The naming is somewhat confusing, but we need g->at0 and g->at1
     * for shifthing coordinates to a new array (not in place) when
     * 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->parts     = t_graph::BondedParts::Single;

    snew(nbond, at_end);
    nbtot = calc_start_end(fplog, g, ilist, at_start, at_end, nbond);

    if (g->at_start >= g->at_end)
    {
        g->at_start = at_start;
        g->at_end   = at_end;
        g->nnodes   = 0;
        g->nbound   = 0;
    }
    else
    {
        g->nnodes = g->at_end - g->at_start;
        snew(g->nedge, g->nnodes);
        snew(g->edge, g->nnodes);
        /* Allocate a single array and set pointers into it */
        snew(g->edge[0], nbtot);
        for (i = 1; (i < g->nnodes); i++)
        {
            g->edge[i] = g->edge[i-1] + nbond[g->at_start+i-1];
        }

        if (!bShakeOnly)
        {
            /* First add all the real bonds: they should determine the molecular
             * graph.
             */
            for (i = 0; (i < F_NRE); i++)
            {
                if (interaction_function[i].flags & IF_CHEMBOND)
                {
                    mk_igraph(g, i, ilist[i], at_start, at_end, nullptr);
                }
            }

            /* Determine of which separated parts the IF_CHEMBOND graph consists.
             * Store the parts in the nbond array.
             */
            bMultiPart = determine_graph_parts(g, nbond);

            if (bMultiPart)
            {
                /* Then add all the other interactions in fixed lists,
                 * 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))
                    {
                        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 */
            compact_graph(fplog, g);
        }
        else
        {
            /* This is a special thing used in splitter.c to generate shake-blocks */
            mk_igraph(g, F_CONSTR, ilist[F_CONSTR], at_start, at_end, nullptr);
            if (bSettle)
            {
                mk_igraph(g, F_SETTLE, ilist[F_SETTLE], at_start, at_end, nullptr);
            }
        }
        g->nbound = 0;
        for (i = 0; (i < g->nnodes); i++)
        {
            if (g->nedge[i] > 0)
            {
                g->nbound++;
            }
        }
    }

    g->negc = 0;
    g->egc  = nullptr;

    sfree(nbond);

    snew(g->ishift, g->at1);

    if (gmx_debug_at)
    {
        p_graph(debug, "graph", g);
    }
}

void mk_graph_moltype(const gmx_moltype_t &moltype,
                      t_graph             *g)
{
    mk_graph_ilist(nullptr, moltype.ilist.data(), 0, moltype.atoms.nr,
                   FALSE, FALSE,
                   g);
}

t_graph *mk_graph(FILE *fplog,
                  const t_idef *idef, int at_start, int at_end,
                  gmx_bool bShakeOnly, gmx_bool bSettle)
{
    t_graph *g;

    snew(g, 1);

    mk_graph_ilist(fplog, idef->il, at_start, at_end, bShakeOnly, bSettle, g);

    return g;
}

void done_graph(t_graph *g)
{
    GCHECK(g);
    if (g->nnodes > 0)
    {
        sfree(g->nedge);
        sfree(g->edge[0]);
        sfree(g->edge);
        sfree(g->egc);
    }
    sfree(g->ishift);
}

/************************************************************
 *
 *      S H I F T   C A L C U L A T I O N   C O D E
 *
 ************************************************************/

static void mk_1shift_tric(int npbcdim, const matrix box, const rvec hbox,
                           const rvec xi, const rvec xj, const int *mi, int *mj)
{
    /* Calculate periodicity for triclinic box... */
    int  m, d;
    rvec dx;

    rvec_sub(xi, xj, dx);

    mj[ZZ] = 0;
    for (m = npbcdim-1; (m >= 0); m--)
    {
        /* If dx < hbox, then xj will be reduced by box, so that
         * xi - xj will be bigger
         */
        if (dx[m] < -hbox[m])
        {
            mj[m] = mi[m]-1;
            for (d = m-1; d >= 0; d--)
            {
                dx[d] += box[m][d];
            }
        }
        else if (dx[m] >= hbox[m])
        {
            mj[m] = mi[m]+1;
            for (d = m-1; d >= 0; d--)
            {
                dx[d] -= box[m][d];
            }
        }
        else
        {
            mj[m] = mi[m];
        }
    }
}

static void mk_1shift(int npbcdim, const rvec hbox, const rvec xi, const rvec xj,
                      const int *mi, int *mj)
{
    /* Calculate periodicity for rectangular box... */
    int  m;
    rvec dx;

    rvec_sub(xi, xj, dx);

    mj[ZZ] = 0;
    for (m = 0; (m < npbcdim); m++)
    {
        /* If dx < hbox, then xj will be reduced by box, so that
         * xi - xj will be bigger
         */
        if (dx[m] < -hbox[m])
        {
            mj[m] = mi[m]-1;
        }
        else if (dx[m] >= hbox[m])
        {
            mj[m] = mi[m]+1;
        }
        else
        {
            mj[m] = mi[m];
        }
    }
}

static void mk_1shift_screw(const matrix box, const rvec hbox,
                            const rvec xi, const rvec xj, const int *mi, int *mj)
{
    /* Calculate periodicity for rectangular box... */
    int  signi, m;
    rvec dx;

    if ((mi[XX] > 0 &&  mi[XX] % 2 == 1) ||
        (mi[XX] < 0 && -mi[XX] % 2 == 1))
    {
        signi = -1;
    }
    else
    {
        signi =  1;
    }

    rvec_sub(xi, xj, dx);

    if (dx[XX] < -hbox[XX])
    {
        mj[XX] = mi[XX] - 1;
    }
    else if (dx[XX] >= hbox[XX])
    {
        mj[XX] = mi[XX] + 1;
    }
    else
    {
        mj[XX] = mi[XX];
    }
    if (mj[XX] != mi[XX])
    {
        /* Rotate */
        dx[YY] = xi[YY] - (box[YY][YY] + box[ZZ][YY] - xj[YY]);
        dx[ZZ] = xi[ZZ] - (box[ZZ][ZZ]               - xj[ZZ]);
    }
    for (m = 1; (m < DIM); m++)
    {
        /* The signs are taken such that we can first shift x and rotate
         * and then shift y and z.
         */
        if (dx[m] < -hbox[m])
        {
            mj[m] = mi[m] - signi;
        }
        else if (dx[m] >= hbox[m])
        {
            mj[m] = mi[m] + signi;
        }
        else
        {
            mj[m] = mi[m];
        }
    }
}

static int mk_grey(egCol egc[], t_graph *g, int *AtomI,
                   int npbcdim, const matrix box, const rvec x[], int *nerror)
{
    int          m, j, ng, ai, aj, g0;
    rvec         dx, hbox;
    gmx_bool     bTriclinic;
    ivec         is_aj;
    t_pbc        pbc;

    for (m = 0; (m < DIM); m++)
    {
        hbox[m] = box[m][m]*0.5;
    }
    bTriclinic = TRICLINIC(box);

    g0 = g->at_start;
    ng = 0;
    ai = g0 + *AtomI;

    /* Loop over all the bonds */
    for (j = 0; (j < g->nedge[ai-g0]); j++)
    {
        aj = g->edge[ai-g0][j];
        /* If there is a white one, make it grey and set pbc */
        if (g->bScrewPBC)
        {
            mk_1shift_screw(box, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
        }
        else if (bTriclinic)
        {
            mk_1shift_tric(npbcdim, box, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
        }
        else
        {
            mk_1shift(npbcdim, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
        }

        if (egc[aj-g0] == egcolWhite)
        {
            if (aj - g0 < *AtomI)
            {
                *AtomI = aj - g0;
            }
            egc[aj-g0] = egcolGrey;

            copy_ivec(is_aj, g->ishift[aj]);

            ng++;
        }
        else if ((is_aj[XX] != g->ishift[aj][XX]) ||
                 (is_aj[YY] != g->ishift[aj][YY]) ||
                 (is_aj[ZZ] != g->ishift[aj][ZZ]))
        {
            if (gmx_debug_at)
            {
                set_pbc(&pbc, -1, box);
                pbc_dx(&pbc, x[ai], x[aj], dx);
                fprintf(debug, "mk_grey: shifts for atom %d due to atom %d\n"
                        "are (%d,%d,%d), should be (%d,%d,%d)\n"
                        "dx = (%g,%g,%g)\n",
                        aj+1, ai+1, is_aj[XX], is_aj[YY], is_aj[ZZ],
                        g->ishift[aj][XX], g->ishift[aj][YY], g->ishift[aj][ZZ],
                        dx[XX], dx[YY], dx[ZZ]);
            }
            (*nerror)++;
        }
    }
    return ng;
}

static int first_colour(int fC, egCol Col, t_graph *g, const egCol egc[])
/* Return the first node with colour Col starting at fC.
 * return -1 if none found.
 */
{
    int i;

    for (i = fC; (i < g->nnodes); i++)
    {
        if ((g->nedge[i] > 0) && (egc[i] == Col))
        {
            return i;
        }
    }

    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[])
{
    static int nerror_tot = 0;
    int        npbcdim;
    int        ng, nnodes, i;
    int        nW, nG, nB; /* Number of Grey, Black, White      */
    int        fW, fG;     /* First of each category    */
    int        nerror = 0;

    g->bScrewPBC = (ePBC == epbcSCREW);

    if (ePBC == epbcXY)
    {
        npbcdim = 2;
    }
    else
    {
        npbcdim = 3;
    }

    GCHECK(g);
    /* This puts everything in the central box, that is does not move it
     * at all. If we return without doing this for a system without bonds
     * (i.e. only settles) all water molecules are moved to the opposite octant
     */
    for (i = g->at0; (i < g->at1); i++)
    {
        g->ishift[i][XX] = g->ishift[i][YY] = g->ishift[i][ZZ] = 0;
    }

    if (!g->nbound)
    {
        return;
    }

    nnodes = g->nnodes;
    if (nnodes > g->negc)
    {
        g->negc = nnodes;
        srenew(g->egc, g->negc);
    }
    memset(g->egc, 0, static_cast<size_t>(nnodes*sizeof(g->egc[0])));

    nW = g->nbound;
    nG = 0;
    nB = 0;

    fW = 0;

    /* We even have a loop invariant:
     * nW+nG+nB == g->nbound
     */
#ifdef DEBUG2
    fprintf(log, "Starting W loop\n");
#endif
    while (nW > 0)
    {
        /* Find the first white, this will allways be a larger
         * number than before, because no nodes are made white
         * in the loop
         */
        if ((fW = first_colour(fW, egcolWhite, g, g->egc)) == -1)
        {
            gmx_fatal(FARGS, "No WHITE nodes found while nW=%d\n", nW);
        }

        /* Make the first white node grey */
        g->egc[fW] = egcolGrey;
        nG++;
        nW--;

        /* Initial value for the first grey */
        fG = fW;
#ifdef DEBUG2
        fprintf(log, "Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
                nW, nG, nB, nW+nG+nB);
#endif
        while (nG > 0)
        {
            if ((fG = first_colour(fG, egcolGrey, g, g->egc)) == -1)
            {
                gmx_fatal(FARGS, "No GREY nodes found while nG=%d\n", nG);
            }

            /* Make the first grey node black */
            g->egc[fG] = egcolBlack;
            nB++;
            nG--;

            /* Make all the neighbours of this black node grey
             * and set their periodicity
             */
            ng = mk_grey(g->egc, g, &fG, npbcdim, box, x, &nerror);
            /* ng is the number of white nodes made grey */
            nG += ng;
            nW -= ng;
        }
    }
    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, "%s", mesg.c_str());
        }

        /* The most likely reason for arriving here is a periodic molecule. */

        nerror_tot++;
        if (nerror_tot <= 100)
        {
            fprintf(stderr, "There were %d inconsistent shifts. Check your topology\n",
                    nerror);
            if (log)
            {
                fprintf(log, "There were %d inconsistent shifts. Check your topology\n",
                        nerror);
            }
        }
        if (nerror_tot == 100)
        {
            fprintf(stderr, "Will stop reporting inconsistent shifts\n");
            if (log)
            {
                fprintf(log, "Will stop reporting inconsistent shifts\n");
            }
        }
    }
}

/************************************************************
 *
 *      A C T U A L   S H I F T   C O D E
 *
 ************************************************************/

void shift_x(const t_graph *g, const matrix box, const rvec x[], rvec x_s[])
{
    ivec    *is;
    int      g0, g1;
    int      j, tx, ty, tz;

    GCHECK(g);
    g0 = g->at_start;
    g1 = g->at_end;
    is = g->ishift;

    for (j = g->at0; j < g0; j++)
    {
        copy_rvec(x[j], x_s[j]);
    }

    if (g->bScrewPBC)
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            if ((tx > 0 && tx % 2 == 1) ||
                (tx < 0 && -tx %2 == 1))
            {
                x_s[j][XX] = x[j][XX] + tx*box[XX][XX];
                x_s[j][YY] = box[YY][YY] + box[ZZ][YY] - x[j][YY];
                x_s[j][ZZ] = box[ZZ][ZZ]               - x[j][ZZ];
            }
            else
            {
                x_s[j][XX] = x[j][XX];
            }
            x_s[j][YY] = x[j][YY] + ty*box[YY][YY] + tz*box[ZZ][YY];
            x_s[j][ZZ] = x[j][ZZ] + tz*box[ZZ][ZZ];
        }
    }
    else if (TRICLINIC(box))
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            x_s[j][XX] = x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
            x_s[j][YY] = x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
            x_s[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
        }
    }
    else
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            x_s[j][XX] = x[j][XX]+tx*box[XX][XX];
            x_s[j][YY] = x[j][YY]+ty*box[YY][YY];
            x_s[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
        }
    }

    for (j = g1; j < g->at1; j++)
    {
        copy_rvec(x[j], x_s[j]);
    }
}

void shift_self(const t_graph *g, const matrix box, rvec x[])
{
    ivec    *is;
    int      g0, g1;
    int      j, tx, ty, tz;

    if (g->bScrewPBC)
    {
        gmx_incons("screw pbc not implemented for shift_self");
    }

    g0 = g->at_start;
    g1 = g->at_end;
    is = g->ishift;

#ifdef DEBUG
    fprintf(stderr, "Shifting atoms %d to %d\n", g0, g0+gn);
#endif
    if (TRICLINIC(box))
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            x[j][XX] = x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
            x[j][YY] = x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
            x[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
        }
    }
    else
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            x[j][XX] = x[j][XX]+tx*box[XX][XX];
            x[j][YY] = x[j][YY]+ty*box[YY][YY];
            x[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
        }
    }
}

void unshift_x(const t_graph *g, const matrix box, rvec x[], const rvec x_s[])
{
    ivec    *is;
    int      g0, g1;
    int      j, tx, ty, tz;

    if (g->bScrewPBC)
    {
        gmx_incons("screw pbc not implemented (yet) for unshift_x");
    }

    g0 = g->at_start;
    g1 = g->at_end;
    is = g->ishift;

    for (j = g->at0; j < g0; j++)
    {
        copy_rvec(x_s[j], x[j]);
    }

    if (TRICLINIC(box))
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            x[j][XX] = x_s[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
            x[j][YY] = x_s[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
            x[j][ZZ] = x_s[j][ZZ]-tz*box[ZZ][ZZ];
        }
    }
    else
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            x[j][XX] = x_s[j][XX]-tx*box[XX][XX];
            x[j][YY] = x_s[j][YY]-ty*box[YY][YY];
            x[j][ZZ] = x_s[j][ZZ]-tz*box[ZZ][ZZ];
        }
    }

    for (j = g1; j < g->at1; j++)
    {
        copy_rvec(x_s[j], x[j]);
    }
}

void unshift_self(const t_graph *g, const matrix box, rvec x[])
{
    ivec *is;
    int   g0, g1;
    int   j, tx, ty, tz;

    if (g->bScrewPBC)
    {
        gmx_incons("screw pbc not implemented for unshift_self");
    }

    g0 = g->at_start;
    g1 = g->at_end;
    is = g->ishift;

    if (TRICLINIC(box))
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            x[j][XX] = x[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
            x[j][YY] = x[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
            x[j][ZZ] = x[j][ZZ]-tz*box[ZZ][ZZ];
        }
    }
    else
    {
        for (j = g0; (j < g1); j++)
        {
            tx = is[j][XX];
            ty = is[j][YY];
            tz = is[j][ZZ];

            x[j][XX] = x[j][XX]-tx*box[XX][XX];
            x[j][YY] = x[j][YY]-ty*box[YY][YY];
            x[j][ZZ] = x[j][ZZ]-tz*box[ZZ][ZZ];
        }
    }
}
#undef GCHECK