Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / gmxlib / splitter.c
index 80f12625643954a0e071e1a087e2628a3bd60c8b..cb2eb90f6b976fb595d82462a4baaf291760d9c1 100644 (file)
-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
  *
- *
- *                This source code is part of
- *
- *                 G   R   O   M   A   C   S
- *
- *          GROningen MAchine for Chemical Simulations
- *
- *                        VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
  * of the License, or (at your option) any later version.
  *
- * If you want to redistribute modifications, 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 www.gromacs.org.
+ * 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.
  *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * 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.
  *
- * For more info, check our website at http://www.gromacs.org
+ * 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.
  *
- * And Hey:
- * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "gromacs/legacyheaders/splitter.h"
 
 #include <assert.h>
-#include <stdio.h>
+#include <stdlib.h>
 #include <string.h>
 
-#include "sysstuff.h"
-#include "macros.h"
-#include "smalloc.h"
-#include "typedefs.h"
-#include "mshift.h"
-#include "invblock.h"
-#include "txtdump.h"
-#include <math.h>
-#include "gmx_fatal.h"
-#include "splitter.h"
-
-typedef struct {
-    int      nr;
-    t_iatom *ia;
-} t_sf;
-
-static t_sf *init_sf(int nr)
-{
-    t_sf *sf;
-    int   i;
-
-    snew(sf, nr);
-    for (i = 0; (i < nr); i++)
-    {
-        sf[i].nr = 0;
-        sf[i].ia = NULL;
-    }
-
-    return sf;
-}
-
-static void done_sf(int nr, t_sf *sf)
-{
-    int i;
-
-    for (i = 0; (i < nr); i++)
-    {
-        sf[i].nr = 0;
-        sfree(sf[i].ia);
-        sf[i].ia = NULL;
-    }
-    sfree(sf);
-}
-
-static void push_sf(t_sf *sf, int nr, t_iatom ia[])
-{
-    int i;
-
-    srenew(sf->ia, sf->nr+nr);
-    for (i = 0; (i < nr); i++)
-    {
-        sf->ia[sf->nr+i] = ia[i];
-    }
-    sf->nr += nr;
-}
-
-static int min_nodeid(int nr, atom_id list[], int hid[])
-{
-    int i, nodeid, minnodeid;
-
-    if (nr <= 0)
-    {
-        gmx_incons("Invalid node number");
-    }
-    minnodeid = hid[list[0]];
-    for (i = 1; (i < nr); i++)
-    {
-        if ((nodeid = hid[list[i]]) < minnodeid)
-        {
-            minnodeid = nodeid;
-        }
-    }
-
-    return minnodeid;
-}
-
-
-static void split_force2(t_inputrec *ir, int nnodes, int hid[], int ftype, t_ilist *ilist,
-                         int *multinr,
-                         int *constr_min_nodeid, int * constr_max_nodeid,
-                         int *left_range, int *right_range)
-{
-    int      i, j, k, type, nodeid, nratoms, tnr;
-    int      nvsite_constr;
-    t_iatom  ai, aj;
-    int      node_low_ai, node_low_aj, node_high_ai, node_high_aj;
-    int      node_low, node_high;
-    int      nodei, nodej;
-    t_sf    *sf;
-    int      nextra;
-
-    sf = init_sf(nnodes);
-
-    node_high = node_low = 0;
-    nextra    = 0;
-
-    /* Walk along all the bonded forces, find the appropriate node
-     * to calc it on, and add it to that nodes list.
-     */
-    for (i = 0; i < ilist->nr; i += (1+nratoms))
-    {
-        type    = ilist->iatoms[i];
-        nratoms = interaction_function[ftype].nratoms;
-
-        if (ftype == F_CONSTR)
-        {
-            ai = ilist->iatoms[i+1];
-            aj = ilist->iatoms[i+2];
-
-            nodei  = hid[ai];
-            nodej  = hid[aj];
-            nodeid = nodei;
-
-            if (ir->eConstrAlg == econtLINCS)
-            {
-                node_low_ai   = constr_min_nodeid[ai];
-                node_low_aj   = constr_min_nodeid[aj];
-                node_high_ai  = constr_max_nodeid[ai];
-                node_high_aj  = constr_max_nodeid[aj];
-
-                node_low      = min(node_low_ai, node_low_aj);
-                node_high     = max(node_high_ai, node_high_aj);
-
-                if (node_high-nodei > 1 || nodei-node_low > 1 ||
-                    node_high-nodej > 1 || nodej-node_low > 1)
-                {
-                    gmx_fatal(FARGS, "Constraint dependencies further away than next-neighbor\n"
-                              "in particle decomposition. Constraint between atoms %d--%d evaluated\n"
-                              "on node %d and %d, but atom %d has connections within %d bonds (lincs_order)\n"
-                              "of node %d, and atom %d has connections within %d bonds of node %d.\n"
-                              "Reduce the # nodes, lincs_order, or\n"
-                              "try domain decomposition.", ai, aj, nodei, nodej, ai, ir->nProjOrder, node_low, aj, ir->nProjOrder, node_high);
-                }
-
-                if (node_low < nodei || node_low < nodej)
-                {
-                    right_range[node_low] = max(right_range[node_low], aj);
-                }
-                if (node_high > nodei || node_high > nodej)
-                {
-                    left_range[node_high] = min(left_range[node_high], ai);
-                }
-            }
-            else
-            {
-                /* Shake */
-                if (hid[ilist->iatoms[i+2]] != nodei)
-                {
-                    gmx_fatal(FARGS, "Shake block crossing node boundaries\n"
-                              "constraint between atoms (%d,%d) (try LINCS instead!)",
-                              ilist->iatoms[i+1]+1, ilist->iatoms[i+2]+1);
-                }
-            }
-        }
-        else if (ftype == F_SETTLE)
-        {
-            /* Only the first particle is stored for settles ... */
-            ai     = ilist->iatoms[i+1];
-            nodeid = hid[ai];
-            if (nodeid != hid[ilist->iatoms[i+2]] ||
-                nodeid != hid[ilist->iatoms[i+3]])
-            {
-                gmx_fatal(FARGS, "Settle block crossing node boundaries\n"
-                          "constraint between atoms %d, %d, %d)",
-                          ai, ilist->iatoms[i+2], ilist->iatoms[i+3]);
-            }
-        }
-        else if (interaction_function[ftype].flags & IF_VSITE)
-        {
-            /* Virtual sites are special, since we need to pre-communicate
-             * their coordinates to construct vsites before then main
-             * coordinate communication.
-             * Vsites can have constructing atoms both larger and smaller than themselves.
-             * To minimize communication and book-keeping, each vsite is constructed on
-             * the home node of the atomnr of the vsite.
-             * Since the vsite coordinates too have to be communicated to the next node,
-             * we need to
-             *
-             * 1. Pre-communicate coordinates of constructing atoms
-             * 2. Construct the vsite
-             * 3. Perform main coordinate communication
-             *
-             * Note that this has change from gromacs 4.0 and earlier, where the vsite
-             * was constructed on the home node of the lowest index of any of the constructing
-             * atoms and the vsite itself.
-             */
-
-            if (ftype == F_VSITE2)
-            {
-                nvsite_constr = 2;
-            }
-            else if (ftype == F_VSITE4FD || ftype == F_VSITE4FDN)
-            {
-                nvsite_constr = 4;
-            }
-            else
-            {
-                nvsite_constr = 3;
-            }
-
-            /* Vsites are constructed on the home node of the actual site to save communication
-             * and simplify the book-keeping.
-             */
-            nodeid = hid[ilist->iatoms[i+1]];
-
-            for (k = 2; k < nvsite_constr+2; k++)
-            {
-                if (hid[ilist->iatoms[i+k]] < (nodeid-1) ||
-                    hid[ilist->iatoms[i+k]] > (nodeid+1))
-                {
-                    gmx_fatal(FARGS, "Virtual site %d and its constructing"
-                              " atoms are not on the same or adjacent\n"
-                              " nodes. This is necessary to avoid a lot\n"
-                              " of extra communication. The easiest way"
-                              " to ensure this is to place virtual sites\n"
-                              " close to the constructing atoms.\n"
-                              " Sorry, but you will have to rework your topology!\n",
-                              ilist->iatoms[i+1]);
-                }
-            }
-        }
-        else
-        {
-            nodeid = min_nodeid(nratoms, &ilist->iatoms[i+1], hid);
-        }
-
-        if (ftype == F_CONSTR && ir->eConstrAlg == econtLINCS)
-        {
-            push_sf(&(sf[nodeid]), nratoms+1, &(ilist->iatoms[i]));
-
-            if (node_low < nodeid)
-            {
-                push_sf(&(sf[node_low]), nratoms+1, &(ilist->iatoms[i]));
-                nextra += nratoms+1;
-            }
-            if (node_high > nodeid)
-            {
-                push_sf(&(sf[node_high]), nratoms+1, &(ilist->iatoms[i]));
-                nextra += nratoms+1;
-            }
-        }
-        else
-        {
-            push_sf(&(sf[nodeid]), nratoms+1, &(ilist->iatoms[i]));
-        }
-    }
-
-    if (nextra > 0)
-    {
-        ilist->nr += nextra;
-        srenew(ilist->iatoms, ilist->nr);
-    }
-
-    tnr = 0;
-    for (nodeid = 0; (nodeid < nnodes); nodeid++)
-    {
-        for (i = 0; (i < sf[nodeid].nr); i++)
-        {
-            ilist->iatoms[tnr++] = sf[nodeid].ia[i];
-        }
-
-        multinr[nodeid]  = (nodeid == 0) ? 0 : multinr[nodeid-1];
-        multinr[nodeid] += sf[nodeid].nr;
-    }
-
-    if (tnr != ilist->nr)
-    {
-        gmx_incons("Splitting forces over processors");
-    }
-
-    done_sf(nnodes, sf);
-}
-
-static int *home_index(int nnodes, t_block *cgs, int *multinr)
-{
-    /* This routine determines the node id for each particle */
-    int *hid;
-    int  nodeid, j0, j1, j, k;
-
-    snew(hid, cgs->index[cgs->nr]);
-    /* Initiate to -1 to make it possible to check afterwards,
-     * all hid's should be set in the loop below
-     */
-    for (k = 0; (k < cgs->index[cgs->nr]); k++)
-    {
-        hid[k] = -1;
-    }
-
-    /* loop over nodes */
-    for (nodeid = 0; (nodeid < nnodes); nodeid++)
-    {
-        j0 = (nodeid == 0) ? 0 : multinr[nodeid-1];
-        j1 = multinr[nodeid];
-
-        /* j0 and j1 are the boundariesin the index array */
-        for (j = j0; (j < j1); j++)
-        {
-            for (k = cgs->index[j]; (k < cgs->index[j+1]); k++)
-            {
-                hid[k] = nodeid;
-            }
-        }
-    }
-    /* Now verify that all hid's are not -1 */
-    for (k = 0; (k < cgs->index[cgs->nr]); k++)
-    {
-        if (hid[k] == -1)
-        {
-            gmx_fatal(FARGS, "hid[%d] = -1, cgs->nr = %d, natoms = %d",
-                      k, cgs->nr, cgs->index[cgs->nr]);
-        }
-    }
-
-    return hid;
-}
-
-typedef struct {
-    int atom, ic, is;
-} t_border;
-
-void set_bor(t_border *b, int atom, int ic, int is)
-{
-    if (debug)
-    {
-        fprintf(debug, "border @ atom %5d [ ic = %5d,  is = %5d ]\n", atom, ic, is);
-    }
-    b->atom = atom;
-    b->ic   = ic;
-    b->is   = is;
-}
-
-static gmx_bool is_bor(atom_id ai[], int i)
-{
-    return ((ai[i] != ai[i-1]) || ((ai[i] == NO_ATID) && (ai[i-1] == NO_ATID)));
-}
-
-static t_border *mk_border(FILE *fp, int natom, atom_id *invcgs,
-                           atom_id *invshk, int *nb)
-{
-    t_border *bor;
-    atom_id  *sbor, *cbor;
-    int       i, j, is, ic, ns, nc, nbor;
-
-    if (debug)
-    {
-        for (i = 0; (i < natom); i++)
-        {
-            fprintf(debug, "atom: %6d  cgindex: %6d  shkindex: %6d\n",
-                    i, invcgs[i], invshk[i]);
-        }
-    }
-
-    snew(sbor, natom+1);
-    snew(cbor, natom+1);
-    ns = nc = 1;
-    for (i = 1; (i < natom); i++)
-    {
-        if (is_bor(invcgs, i))
-        {
-            cbor[nc++] = i;
-        }
-        if (is_bor(invshk, i))
-        {
-            sbor[ns++] = i;
-        }
-    }
-    sbor[ns] = 0;
-    cbor[nc] = 0;
-    if (fp)
-    {
-        fprintf(fp, "There are %d charge group borders", nc);
-        if (invshk != NULL)
-        {
-            fprintf(fp, " and %d shake borders", ns);
-        }
-        fprintf(fp, ".\n");
-    }
-    snew(bor, max(nc, ns));
-    ic = is = nbor = 0;
-    while ((ic < nc) || (is < ns))
-    {
-        if (sbor[is] == cbor[ic])
-        {
-            set_bor(&(bor[nbor]), cbor[ic], ic, is);
-            nbor++;
-            if (ic < nc)
-            {
-                ic++;
-            }
-            if (is < ns)
-            {
-                is++;
-            }
-        }
-        else if (cbor[ic] > sbor[is])
-        {
-            if (is == ns)
-            {
-                set_bor(&(bor[nbor]), cbor[ic], ic, is);
-                nbor++;
-                if (ic < nc)
-                {
-                    ic++;
-                }
-            }
-            else if (is < ns)
-            {
-                is++;
-            }
-        }
-        else if (ic < nc)
-        {
-            ic++;
-        }
-        else
-        {
-            is++; /*gmx_fatal(FARGS,"Can't happen is=%d, ic=%d (%s, %d)",
-                     is,ic,__FILE__,__LINE__);*/
-        }
-    }
-    if (fp)
-    {
-        fprintf(fp, "There are %d total borders\n", nbor);
-    }
-
-    if (debug)
-    {
-        fprintf(debug, "There are %d actual bor entries\n", nbor);
-        for (i = 0; (i < nbor); i++)
-        {
-            fprintf(debug, "bor[%5d] = atom: %d  ic: %d  is: %d\n", i,
-                    bor[i].atom, bor[i].ic, bor[i].is);
-        }
-    }
-
-    *nb = nbor;
-
-    return bor;
-}
-
-static void split_blocks(FILE *fp, t_inputrec *ir, int nnodes,
-                         t_block *cgs, t_blocka *sblock, real capacity[],
-                         int *multinr_cgs)
-{
-    int         natoms, *maxatom;
-    int         i, ii, ai, b0, b1;
-    int         nodeid, last_shk, nbor;
-    t_border   *border;
-    double      tload, tcap;
-
-    gmx_bool    bSHK;
-    atom_id    *shknum, *cgsnum;
-
-    natoms = cgs->index[cgs->nr];
-
-    if (NULL != debug)
-    {
-        pr_block(debug, 0, "cgs", cgs, TRUE);
-        pr_blocka(debug, 0, "sblock", sblock, TRUE);
-        fflush(debug);
-    }
-
-    cgsnum = make_invblock(cgs, natoms+1);
-    shknum = make_invblocka(sblock, natoms+1);
-    border = mk_border(fp, natoms, cgsnum, shknum, &nbor);
-
-    snew(maxatom, nnodes);
-    tload  = capacity[0]*natoms;
-    tcap   = 1.0;
-    nodeid = 0;
-    /* Start at bor is 1, to force the first block on the first processor */
-    for (i = 0; (i < nbor) && (tload < natoms); i++)
-    {
-        if (i < (nbor-1))
-        {
-            b1 = border[i+1].atom;
-        }
-        else
-        {
-            b1 = natoms;
-        }
-
-        b0 = border[i].atom;
-
-        if ((fabs(b0-tload) < fabs(b1-tload)))
-        {
-            /* New nodeid time */
-            multinr_cgs[nodeid] = border[i].ic;
-            maxatom[nodeid]     = b0;
-            tcap               -= capacity[nodeid];
-            nodeid++;
-
-            /* Recompute target load */
-            tload = b0 + (natoms-b0)*capacity[nodeid]/tcap;
-
-            if (debug)
-            {
-                printf("tload: %g tcap: %g  nodeid: %d\n", tload, tcap, nodeid);
-            }
-        }
-    }
-    /* Now the last one... */
-    while (nodeid < nnodes)
-    {
-        multinr_cgs[nodeid] = cgs->nr;
-        /* Store atom number, see above */
-        maxatom[nodeid]     = natoms;
-        nodeid++;
-    }
-    if (nodeid != nnodes)
-    {
-        gmx_fatal(FARGS, "nodeid = %d, nnodes = %d, file %s, line %d",
-                  nodeid, nnodes, __FILE__, __LINE__);
-    }
-
-    for (i = nnodes-1; (i > 0); i--)
-    {
-        maxatom[i] -= maxatom[i-1];
-    }
-
-    if (fp)
-    {
-        fprintf(fp, "Division over nodes in atoms:\n");
-        for (i = 0; (i < nnodes); i++)
-        {
-            fprintf(fp, " %7d", maxatom[i]);
-        }
-        fprintf(fp, "\n");
-    }
-
-    sfree(maxatom);
-    sfree(shknum);
-    sfree(cgsnum);
-    sfree(border);
-}
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
 
 typedef struct {
     int atom, sid;
@@ -600,7 +72,7 @@ static int sid_comp(const void *a, const void *b)
     }
 }
 
-static int mk_grey(int nnodes, egCol egc[], t_graph *g, int *AtomI,
+static int mk_grey(egCol egc[], t_graph *g, int *AtomI,
                    int maxsid, t_sid sid[])
 {
     int  j, ng, ai, aj, g0;
@@ -735,7 +207,7 @@ static int mk_sblocks(FILE *fp, t_graph *g, int maxsid, t_sid sid[])
             /* Make all the neighbours of this black node grey
              * and set their block number
              */
-            ng = mk_grey(nnodes, egc, g, &fG, maxsid, sid);
+            ng = mk_grey(egc, g, &fG, maxsid, sid);
             /* ng is the number of white nodes made grey */
             nG += ng;
             nW -= ng;
@@ -773,7 +245,7 @@ static int ms_comp(const void *a, const void *b)
     }
 }
 
-static int merge_sid(int i0, int at_start, int at_end, int nsid, t_sid sid[],
+static int merge_sid(int at_start, int at_end, int nsid, t_sid sid[],
                      t_blocka *sblock)
 {
     int          i, j, k, n, isid, ndel;
@@ -941,7 +413,7 @@ void gen_sblocks(FILE *fp, int at_start, int at_end,
      * part of the shake block too. There may be cases where blocks overlap
      * and they will have to be merged.
      */
-    nsid = merge_sid(i0, at_start, at_end, nsid, sid, sblock);
+    nsid = merge_sid(at_start, at_end, nsid, sid, sblock);
     /* Now sort the shake blocks again... */
     /*qsort(sid,natoms,(size_t)sizeof(sid[0]),sid_comp);*/
 
@@ -980,180 +452,3 @@ void gen_sblocks(FILE *fp, int at_start, int at_end,
         fprintf(debug, "Done gen_sblocks\n");
     }
 }
-
-static t_blocka block2blocka(t_block *block)
-{
-    t_blocka blocka;
-    int      i;
-
-    blocka.nr           = block->nr;
-    blocka.nalloc_index = blocka.nr + 1;
-    snew(blocka.index, blocka.nalloc_index);
-    for (i = 0; i <= block->nr; i++)
-    {
-        blocka.index[i] = block->index[i];
-    }
-    blocka.nra      = block->index[block->nr];
-    blocka.nalloc_a = blocka.nra;
-    snew(blocka.a, blocka.nalloc_a);
-    for (i = 0; i < blocka.nra; i++)
-    {
-        blocka.a[i] = i;
-    }
-
-    return blocka;
-}
-
-typedef struct
-{
-    int   nconstr;
-    int   index[10];
-} pd_constraintlist_t;
-
-
-static void
-find_constraint_range_recursive(pd_constraintlist_t *  constraintlist,
-                                int                    thisatom,
-                                int                    depth,
-                                int *                  min_atomid,
-                                int *                  max_atomid)
-{
-    int i, j;
-    int nconstr;
-
-    for (i = 0; i < constraintlist[thisatom].nconstr; i++)
-    {
-        j = constraintlist[thisatom].index[i];
-
-        *min_atomid = (j < *min_atomid) ? j : *min_atomid;
-        *max_atomid = (j > *max_atomid) ? j : *max_atomid;
-
-        if (depth > 0)
-        {
-            find_constraint_range_recursive(constraintlist, j, depth-1, min_atomid, max_atomid);
-        }
-    }
-}
-
-static void
-pd_determine_constraints_range(t_inputrec *      ir,
-                               int               natoms,
-                               t_ilist *         ilist,
-                               int               hid[],
-                               int *             min_nodeid,
-                               int *             max_nodeid)
-{
-    int                  i, j, k;
-    int                  nratoms;
-    int                  depth;
-    int                  ai, aj;
-    int                  min_atomid, max_atomid;
-    pd_constraintlist_t *constraintlist;
-
-    nratoms = interaction_function[F_CONSTR].nratoms;
-    depth   = ir->nProjOrder;
-
-    snew(constraintlist, natoms);
-
-    /* Make a list of all the connections */
-    for (i = 0; i < ilist->nr; i += nratoms+1)
-    {
-        ai = ilist->iatoms[i+1];
-        aj = ilist->iatoms[i+2];
-        constraintlist[ai].index[constraintlist[ai].nconstr++] = aj;
-        constraintlist[aj].index[constraintlist[aj].nconstr++] = ai;
-    }
-
-    for (i = 0; i < natoms; i++)
-    {
-        min_atomid = i;
-        max_atomid = i;
-
-        find_constraint_range_recursive(constraintlist, i, depth, &min_atomid, &max_atomid);
-
-        min_nodeid[i] = hid[min_atomid];
-        max_nodeid[i] = hid[max_atomid];
-    }
-    sfree(constraintlist);
-}
-
-
-void split_top(FILE *fp, int nnodes, gmx_localtop_t *top, t_inputrec *ir, t_block *mols,
-               real *capacity, int *multinr_cgs, int **multinr_nre, int *left_range, int * right_range)
-{
-    int        natoms, i, j, k, mj, atom, maxatom, sstart, send, bstart, nodeid;
-    t_blocka   sblock;
-    int       *homeind;
-    int        ftype, nvsite_constr, nra, nrd;
-    t_iatom   *ia;
-    int        minhome, ihome, minidx;
-    int       *constr_min_nodeid;
-    int       *constr_max_nodeid;
-
-    if (nnodes <= 1)
-    {
-        return;
-    }
-
-    natoms = mols->index[mols->nr];
-
-    if (fp)
-    {
-        fprintf(fp, "splitting topology...\n");
-    }
-
-/*#define MOL_BORDER*/
-/*Removed the above to allow splitting molecules with h-bond constraints
-   over processors. The results in DP are the same. */
-    init_blocka(&sblock);
-    if (ir->eConstrAlg != econtLINCS)
-    {
-#ifndef MOL_BORDER
-        /* Make a special shake block that includes settles */
-        gen_sblocks(fp, 0, natoms, &top->idef, &sblock, TRUE);
-#else
-        sblock = block2blocka(mols);
-#endif
-    }
-
-    split_blocks(fp, ir, nnodes, &top->cgs, &sblock, capacity, multinr_cgs);
-
-    homeind = home_index(nnodes, &top->cgs, multinr_cgs);
-
-    snew(constr_min_nodeid, natoms);
-    snew(constr_max_nodeid, natoms);
-
-    if (top->idef.il[F_CONSTR].nr > 0)
-    {
-        pd_determine_constraints_range(ir, natoms, &top->idef.il[F_CONSTR], homeind, constr_min_nodeid, constr_max_nodeid);
-    }
-    else
-    {
-        /* Not 100% necessary, but it is a bad habit to have uninitialized arrays around... */
-        for (i = 0; i < natoms; i++)
-        {
-            constr_min_nodeid[i] = constr_max_nodeid[i] = homeind[i];
-        }
-    }
-
-    /* Default limits (no communication) for PD constraints */
-    left_range[0] = 0;
-    for (i = 1; i < nnodes; i++)
-    {
-        left_range[i]    = top->cgs.index[multinr_cgs[i-1]];
-        right_range[i-1] = left_range[i]-1;
-    }
-    right_range[nnodes-1] = top->cgs.index[multinr_cgs[nnodes-1]]-1;
-
-    for (j = 0; (j < F_NRE); j++)
-    {
-        split_force2(ir, nnodes, homeind, j, &top->idef.il[j], multinr_nre[j], constr_min_nodeid, constr_max_nodeid,
-                     left_range, right_range);
-    }
-
-    sfree(constr_min_nodeid);
-    sfree(constr_max_nodeid);
-
-    sfree(homeind);
-    done_blocka(&sblock);
-}