2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
53 #include "gmx_fatal.h"
61 static t_sf *init_sf(int nr)
67 for (i = 0; (i < nr); i++)
76 static void done_sf(int nr, t_sf *sf)
80 for (i = 0; (i < nr); i++)
89 static void push_sf(t_sf *sf, int nr, t_iatom ia[])
93 srenew(sf->ia, sf->nr+nr);
94 for (i = 0; (i < nr); i++)
96 sf->ia[sf->nr+i] = ia[i];
101 static int min_nodeid(int nr, atom_id list[], int hid[])
103 int i, nodeid, minnodeid;
107 gmx_incons("Invalid node number");
109 minnodeid = hid[list[0]];
110 for (i = 1; (i < nr); i++)
112 if ((nodeid = hid[list[i]]) < minnodeid)
122 static void split_force2(t_inputrec *ir, int nnodes, int hid[], int ftype, t_ilist *ilist,
124 int *constr_min_nodeid, int * constr_max_nodeid,
125 int *left_range, int *right_range)
127 int i, j, k, type, nodeid, nratoms, tnr;
130 int node_low_ai, node_low_aj, node_high_ai, node_high_aj;
131 int node_low, node_high;
136 sf = init_sf(nnodes);
138 node_high = node_low = 0;
141 /* Walk along all the bonded forces, find the appropriate node
142 * to calc it on, and add it to that nodes list.
144 for (i = 0; i < ilist->nr; i += (1+nratoms))
146 type = ilist->iatoms[i];
147 nratoms = interaction_function[ftype].nratoms;
149 if (ftype == F_CONSTR)
151 ai = ilist->iatoms[i+1];
152 aj = ilist->iatoms[i+2];
158 if (ir->eConstrAlg == econtLINCS)
160 node_low_ai = constr_min_nodeid[ai];
161 node_low_aj = constr_min_nodeid[aj];
162 node_high_ai = constr_max_nodeid[ai];
163 node_high_aj = constr_max_nodeid[aj];
165 node_low = min(node_low_ai, node_low_aj);
166 node_high = max(node_high_ai, node_high_aj);
168 if (node_high-nodei > 1 || nodei-node_low > 1 ||
169 node_high-nodej > 1 || nodej-node_low > 1)
171 gmx_fatal(FARGS, "Constraint dependencies further away than next-neighbor\n"
172 "in particle decomposition. Constraint between atoms %d--%d evaluated\n"
173 "on node %d and %d, but atom %d has connections within %d bonds (lincs_order)\n"
174 "of node %d, and atom %d has connections within %d bonds of node %d.\n"
175 "Reduce the # nodes, lincs_order, or\n"
176 "try domain decomposition.", ai, aj, nodei, nodej, ai, ir->nProjOrder, node_low, aj, ir->nProjOrder, node_high);
179 if (node_low < nodei || node_low < nodej)
181 right_range[node_low] = max(right_range[node_low], aj);
183 if (node_high > nodei || node_high > nodej)
185 left_range[node_high] = min(left_range[node_high], ai);
191 if (hid[ilist->iatoms[i+2]] != nodei)
193 gmx_fatal(FARGS, "Shake block crossing node boundaries\n"
194 "constraint between atoms (%d,%d) (try LINCS instead!)",
195 ilist->iatoms[i+1]+1, ilist->iatoms[i+2]+1);
199 else if (ftype == F_SETTLE)
201 /* Only the first particle is stored for settles ... */
202 ai = ilist->iatoms[i+1];
204 if (nodeid != hid[ilist->iatoms[i+2]] ||
205 nodeid != hid[ilist->iatoms[i+3]])
207 gmx_fatal(FARGS, "Settle block crossing node boundaries\n"
208 "constraint between atoms %d, %d, %d)",
209 ai, ilist->iatoms[i+2], ilist->iatoms[i+3]);
212 else if (interaction_function[ftype].flags & IF_VSITE)
214 /* Virtual sites are special, since we need to pre-communicate
215 * their coordinates to construct vsites before then main
216 * coordinate communication.
217 * Vsites can have constructing atoms both larger and smaller than themselves.
218 * To minimize communication and book-keeping, each vsite is constructed on
219 * the home node of the atomnr of the vsite.
220 * Since the vsite coordinates too have to be communicated to the next node,
223 * 1. Pre-communicate coordinates of constructing atoms
224 * 2. Construct the vsite
225 * 3. Perform main coordinate communication
227 * Note that this has change from gromacs 4.0 and earlier, where the vsite
228 * was constructed on the home node of the lowest index of any of the constructing
229 * atoms and the vsite itself.
232 if (ftype == F_VSITE2)
236 else if (ftype == F_VSITE4FD || ftype == F_VSITE4FDN)
245 /* Vsites are constructed on the home node of the actual site to save communication
246 * and simplify the book-keeping.
248 nodeid = hid[ilist->iatoms[i+1]];
250 for (k = 2; k < nvsite_constr+2; k++)
252 if (hid[ilist->iatoms[i+k]] < (nodeid-1) ||
253 hid[ilist->iatoms[i+k]] > (nodeid+1))
255 gmx_fatal(FARGS, "Virtual site %d and its constructing"
256 " atoms are not on the same or adjacent\n"
257 " nodes. This is necessary to avoid a lot\n"
258 " of extra communication. The easiest way"
259 " to ensure this is to place virtual sites\n"
260 " close to the constructing atoms.\n"
261 " Sorry, but you will have to rework your topology!\n",
268 nodeid = min_nodeid(nratoms, &ilist->iatoms[i+1], hid);
271 if (ftype == F_CONSTR && ir->eConstrAlg == econtLINCS)
273 push_sf(&(sf[nodeid]), nratoms+1, &(ilist->iatoms[i]));
275 if (node_low < nodeid)
277 push_sf(&(sf[node_low]), nratoms+1, &(ilist->iatoms[i]));
280 if (node_high > nodeid)
282 push_sf(&(sf[node_high]), nratoms+1, &(ilist->iatoms[i]));
288 push_sf(&(sf[nodeid]), nratoms+1, &(ilist->iatoms[i]));
295 srenew(ilist->iatoms, ilist->nr);
299 for (nodeid = 0; (nodeid < nnodes); nodeid++)
301 for (i = 0; (i < sf[nodeid].nr); i++)
303 ilist->iatoms[tnr++] = sf[nodeid].ia[i];
306 multinr[nodeid] = (nodeid == 0) ? 0 : multinr[nodeid-1];
307 multinr[nodeid] += sf[nodeid].nr;
310 if (tnr != ilist->nr)
312 gmx_incons("Splitting forces over processors");
318 static int *home_index(int nnodes, t_block *cgs, int *multinr)
320 /* This routine determines the node id for each particle */
322 int nodeid, j0, j1, j, k;
324 snew(hid, cgs->index[cgs->nr]);
325 /* Initiate to -1 to make it possible to check afterwards,
326 * all hid's should be set in the loop below
328 for (k = 0; (k < cgs->index[cgs->nr]); k++)
333 /* loop over nodes */
334 for (nodeid = 0; (nodeid < nnodes); nodeid++)
336 j0 = (nodeid == 0) ? 0 : multinr[nodeid-1];
337 j1 = multinr[nodeid];
339 /* j0 and j1 are the boundariesin the index array */
340 for (j = j0; (j < j1); j++)
342 for (k = cgs->index[j]; (k < cgs->index[j+1]); k++)
348 /* Now verify that all hid's are not -1 */
349 for (k = 0; (k < cgs->index[cgs->nr]); k++)
353 gmx_fatal(FARGS, "hid[%d] = -1, cgs->nr = %d, natoms = %d",
354 k, cgs->nr, cgs->index[cgs->nr]);
365 void set_bor(t_border *b, int atom, int ic, int is)
369 fprintf(debug, "border @ atom %5d [ ic = %5d, is = %5d ]\n", atom, ic, is);
376 static gmx_bool is_bor(atom_id ai[], int i)
378 return ((ai[i] != ai[i-1]) || ((ai[i] == NO_ATID) && (ai[i-1] == NO_ATID)));
381 static t_border *mk_border(FILE *fp, int natom, atom_id *invcgs,
382 atom_id *invshk, int *nb)
385 atom_id *sbor, *cbor;
386 int i, j, is, ic, ns, nc, nbor;
390 for (i = 0; (i < natom); i++)
392 fprintf(debug, "atom: %6d cgindex: %6d shkindex: %6d\n",
393 i, invcgs[i], invshk[i]);
400 for (i = 1; (i < natom); i++)
402 if (is_bor(invcgs, i))
406 if (is_bor(invshk, i))
415 fprintf(fp, "There are %d charge group borders", nc);
418 fprintf(fp, " and %d shake borders", ns);
422 snew(bor, max(nc, ns));
424 while ((ic < nc) || (is < ns))
426 if (sbor[is] == cbor[ic])
428 set_bor(&(bor[nbor]), cbor[ic], ic, is);
439 else if (cbor[ic] > sbor[is])
443 set_bor(&(bor[nbor]), cbor[ic], ic, is);
461 is++; /*gmx_fatal(FARGS,"Can't happen is=%d, ic=%d (%s, %d)",
462 is,ic,__FILE__,__LINE__);*/
467 fprintf(fp, "There are %d total borders\n", nbor);
472 fprintf(debug, "There are %d actual bor entries\n", nbor);
473 for (i = 0; (i < nbor); i++)
475 fprintf(debug, "bor[%5d] = atom: %d ic: %d is: %d\n", i,
476 bor[i].atom, bor[i].ic, bor[i].is);
485 static void split_blocks(FILE *fp, int nnodes,
486 t_block *cgs, t_blocka *sblock, real capacity[],
489 int natoms, *maxatom;
490 int i, ii, ai, b0, b1;
491 int nodeid, last_shk, nbor;
496 atom_id *shknum, *cgsnum;
498 natoms = cgs->index[cgs->nr];
502 pr_block(debug, 0, "cgs", cgs, TRUE);
503 pr_blocka(debug, 0, "sblock", sblock, TRUE);
507 cgsnum = make_invblock(cgs, natoms+1);
508 shknum = make_invblocka(sblock, natoms+1);
509 border = mk_border(fp, natoms, cgsnum, shknum, &nbor);
511 snew(maxatom, nnodes);
512 tload = capacity[0]*natoms;
515 /* Start at bor is 1, to force the first block on the first processor */
516 for (i = 0; (i < nbor) && (tload < natoms); i++)
520 b1 = border[i+1].atom;
529 if ((fabs(b0-tload) < fabs(b1-tload)))
531 /* New nodeid time */
532 multinr_cgs[nodeid] = border[i].ic;
533 maxatom[nodeid] = b0;
534 tcap -= capacity[nodeid];
537 /* Recompute target load */
538 tload = b0 + (natoms-b0)*capacity[nodeid]/tcap;
542 printf("tload: %g tcap: %g nodeid: %d\n", tload, tcap, nodeid);
546 /* Now the last one... */
547 while (nodeid < nnodes)
549 multinr_cgs[nodeid] = cgs->nr;
550 /* Store atom number, see above */
551 maxatom[nodeid] = natoms;
554 if (nodeid != nnodes)
556 gmx_fatal(FARGS, "nodeid = %d, nnodes = %d, file %s, line %d",
557 nodeid, nnodes, __FILE__, __LINE__);
560 for (i = nnodes-1; (i > 0); i--)
562 maxatom[i] -= maxatom[i-1];
567 fprintf(fp, "Division over nodes in atoms:\n");
568 for (i = 0; (i < nnodes); i++)
570 fprintf(fp, " %7d", maxatom[i]);
585 static int sid_comp(const void *a, const void *b)
593 dd = sa->sid-sb->sid;
596 return (sa->atom-sb->atom);
604 static int mk_grey(egCol egc[], t_graph *g, int *AtomI,
605 int maxsid, t_sid sid[])
607 int j, ng, ai, aj, g0;
613 /* Loop over all the bonds */
614 for (j = 0; (j < g->nedge[ai]); j++)
616 aj = g->edge[ai][j]-g0;
617 /* If there is a white one, make it gray and set pbc */
618 if (egc[aj] == egcolWhite)
626 /* Check whether this one has been set before... */
627 range_check(aj+g0, 0, maxsid);
628 range_check(ai+g0, 0, maxsid);
629 if (sid[aj+g0].sid != -1)
631 gmx_fatal(FARGS, "sid[%d]=%d, sid[%d]=%d, file %s, line %d",
632 ai, sid[ai+g0].sid, aj, sid[aj+g0].sid, __FILE__, __LINE__);
636 sid[aj+g0].sid = sid[ai+g0].sid;
637 sid[aj+g0].atom = aj+g0;
645 static int first_colour(int fC, egCol Col, t_graph *g, egCol egc[])
646 /* Return the first node with colour Col starting at fC.
647 * return -1 if none found.
652 for (i = fC; (i < g->nnodes); i++)
654 if ((g->nedge[i] > 0) && (egc[i] == Col))
663 static int mk_sblocks(FILE *fp, t_graph *g, int maxsid, t_sid sid[])
666 int nW, nG, nB; /* Number of Grey, Black, White */
667 int fW, fG; /* First of each category */
668 egCol *egc = NULL; /* The colour of each node */
688 /* We even have a loop invariant:
689 * nW+nG+nB == g->nbound
694 fprintf(fp, "Walking down the molecule graph to make constraint-blocks\n");
699 /* Find the first white, this will allways be a larger
700 * number than before, because no nodes are made white
703 if ((fW = first_colour(fW, egcolWhite, g, egc)) == -1)
705 gmx_fatal(FARGS, "No WHITE nodes found while nW=%d\n", nW);
708 /* Make the first white node grey, and set the block number */
710 range_check(fW+g0, 0, maxsid);
711 sid[fW+g0].sid = nblock++;
715 /* Initial value for the first grey */
720 fprintf(debug, "Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
721 nW, nG, nB, nW+nG+nB);
726 if ((fG = first_colour(fG, egcolGrey, g, egc)) == -1)
728 gmx_fatal(FARGS, "No GREY nodes found while nG=%d\n", nG);
731 /* Make the first grey node black */
732 egc[fG] = egcolBlack;
736 /* Make all the neighbours of this black node grey
737 * and set their block number
739 ng = mk_grey(egc, g, &fG, maxsid, sid);
740 /* ng is the number of white nodes made grey */
749 fprintf(debug, "Found %d shake blocks\n", nblock);
757 int first, last, sid;
760 static int ms_comp(const void *a, const void *b)
762 t_merge_sid *ma = (t_merge_sid *)a;
763 t_merge_sid *mb = (t_merge_sid *)b;
766 d = ma->first-mb->first;
769 return ma->last-mb->last;
777 static int merge_sid(int at_start, int at_end, int nsid, t_sid sid[],
780 int i, j, k, n, isid, ndel;
784 /* We try to remdy the following problem:
785 * Atom: 1 2 3 4 5 6 7 8 9 10
786 * Sid: 0 -1 1 0 -1 1 1 1 1 1
789 /* Determine first and last atom in each shake ID */
792 for (k = 0; (k < nsid); k++)
794 ms[k].first = at_end+1;
798 for (i = at_start; (i < at_end); i++)
801 range_check(isid, -1, nsid);
804 ms[isid].first = min(ms[isid].first, sid[i].atom);
805 ms[isid].last = max(ms[isid].last, sid[i].atom);
808 qsort(ms, nsid, sizeof(ms[0]), ms_comp);
810 /* Now merge the overlapping ones */
812 for (k = 0; (k < nsid); )
814 for (j = k+1; (j < nsid); )
816 if (ms[j].first <= ms[k].last)
818 ms[k].last = max(ms[k].last, ms[j].last);
819 ms[k].first = min(ms[k].first, ms[j].first);
835 for (k = 0; (k < nsid); k++)
837 while ((k < nsid-1) && (ms[k].sid == -1))
839 for (j = k+1; (j < nsid); j++)
841 memcpy(&(ms[j-1]), &(ms[j]), sizeof(ms[0]));
847 for (k = at_start; (k < at_end); k++)
853 sblock->nalloc_index = sblock->nr+1;
854 snew(sblock->index, sblock->nalloc_index);
855 sblock->nra = at_end - at_start;
856 sblock->nalloc_a = sblock->nra;
857 snew(sblock->a, sblock->nalloc_a);
858 sblock->index[0] = 0;
859 for (k = n = 0; (k < nsid); k++)
861 sblock->index[k+1] = sblock->index[k] + ms[k].last - ms[k].first+1;
862 for (j = ms[k].first; (j <= ms[k].last); j++)
864 range_check(n, 0, sblock->nra);
866 range_check(j, 0, at_end);
867 if (sid[j].sid == -1)
873 fprintf(stderr, "Double sids (%d, %d) for atom %d\n", sid[j].sid, k, j);
878 /* Removed 2007-09-04
879 sblock->index[k+1] = natoms;
880 for(k=0; (k<natoms); k++)
881 if (sid[k].sid == -1)
886 assert(sblock->index[k] == sblock->nra);
892 void gen_sblocks(FILE *fp, int at_start, int at_end,
893 t_idef *idef, t_blocka *sblock,
897 int i, i0, j, k, istart, n;
901 g = mk_graph(NULL, idef, at_start, at_end, TRUE, bSettle);
904 p_graph(debug, "Graaf Dracula", g);
907 for (i = at_start; (i < at_end); i++)
912 nsid = mk_sblocks(fp, g, at_end, sid);
919 /* Now sort the shake blocks... */
920 qsort(sid+at_start, at_end-at_start, (size_t)sizeof(sid[0]), sid_comp);
924 fprintf(debug, "Sorted shake block\n");
925 for (i = at_start; (i < at_end); i++)
927 fprintf(debug, "sid[%5d] = atom:%5d sid:%5d\n", i, sid[i].atom, sid[i].sid);
930 /* Now check how many are NOT -1, i.e. how many have to be shaken */
931 for (i0 = at_start; (i0 < at_end); i0++)
933 if (sid[i0].sid > -1)
939 /* Now we have the sids that have to be shaken. We'll check the min and
940 * max atom numbers and this determines the shake block. DvdS 2007-07-19.
941 * For the purpose of making boundaries all atoms in between need to be
942 * part of the shake block too. There may be cases where blocks overlap
943 * and they will have to be merged.
945 nsid = merge_sid(at_start, at_end, nsid, sid, sblock);
946 /* Now sort the shake blocks again... */
947 /*qsort(sid,natoms,(size_t)sizeof(sid[0]),sid_comp);*/
949 /* Fill the sblock struct */
950 /* sblock->nr = nsid;
951 sblock->nra = natoms;
952 srenew(sblock->a,sblock->nra);
953 srenew(sblock->index,sblock->nr+1);
958 sblock->index[n++]=k;
960 istart = sid[i].atom;
961 while ((i<natoms-1) && (sid[i+1].sid == isid))
963 /* After while: we found a new block, or are thru with the atoms */
964 /* for(j=istart; (j<=sid[i].atom); j++,k++)
966 sblock->index[n] = k;
970 gmx_fatal(FARGS,"Death Horror: nsid = %d, n= %d",nsid,n);
976 /* Due to unknown reason this free generates a problem sometimes */
981 fprintf(debug, "Done gen_sblocks\n");
985 static t_blocka block2blocka(t_block *block)
990 blocka.nr = block->nr;
991 blocka.nalloc_index = blocka.nr + 1;
992 snew(blocka.index, blocka.nalloc_index);
993 for (i = 0; i <= block->nr; i++)
995 blocka.index[i] = block->index[i];
997 blocka.nra = block->index[block->nr];
998 blocka.nalloc_a = blocka.nra;
999 snew(blocka.a, blocka.nalloc_a);
1000 for (i = 0; i < blocka.nra; i++)
1012 } pd_constraintlist_t;
1016 find_constraint_range_recursive(pd_constraintlist_t * constraintlist,
1025 for (i = 0; i < constraintlist[thisatom].nconstr; i++)
1027 j = constraintlist[thisatom].index[i];
1029 *min_atomid = (j < *min_atomid) ? j : *min_atomid;
1030 *max_atomid = (j > *max_atomid) ? j : *max_atomid;
1034 find_constraint_range_recursive(constraintlist, j, depth-1, min_atomid, max_atomid);
1040 pd_determine_constraints_range(t_inputrec * ir,
1051 int min_atomid, max_atomid;
1052 pd_constraintlist_t *constraintlist;
1054 nratoms = interaction_function[F_CONSTR].nratoms;
1055 depth = ir->nProjOrder;
1057 snew(constraintlist, natoms);
1059 /* Make a list of all the connections */
1060 for (i = 0; i < ilist->nr; i += nratoms+1)
1062 ai = ilist->iatoms[i+1];
1063 aj = ilist->iatoms[i+2];
1064 constraintlist[ai].index[constraintlist[ai].nconstr++] = aj;
1065 constraintlist[aj].index[constraintlist[aj].nconstr++] = ai;
1068 for (i = 0; i < natoms; i++)
1073 find_constraint_range_recursive(constraintlist, i, depth, &min_atomid, &max_atomid);
1075 min_nodeid[i] = hid[min_atomid];
1076 max_nodeid[i] = hid[max_atomid];
1078 sfree(constraintlist);
1082 void split_top(FILE *fp, int nnodes, gmx_localtop_t *top, t_inputrec *ir, t_block *mols,
1083 real *capacity, int *multinr_cgs, int **multinr_nre, int *left_range, int * right_range)
1085 int natoms, i, j, k, mj, atom, maxatom, sstart, send, bstart, nodeid;
1088 int ftype, nvsite_constr, nra, nrd;
1090 int minhome, ihome, minidx;
1091 int *constr_min_nodeid;
1092 int *constr_max_nodeid;
1099 natoms = mols->index[mols->nr];
1103 fprintf(fp, "splitting topology...\n");
1106 /*#define MOL_BORDER*/
1107 /*Removed the above to allow splitting molecules with h-bond constraints
1108 over processors. The results in DP are the same. */
1109 init_blocka(&sblock);
1110 if (ir->eConstrAlg != econtLINCS)
1113 /* Make a special shake block that includes settles */
1114 gen_sblocks(fp, 0, natoms, &top->idef, &sblock, TRUE);
1116 sblock = block2blocka(mols);
1120 split_blocks(fp, nnodes, &top->cgs, &sblock, capacity, multinr_cgs);
1122 homeind = home_index(nnodes, &top->cgs, multinr_cgs);
1124 snew(constr_min_nodeid, natoms);
1125 snew(constr_max_nodeid, natoms);
1127 if (top->idef.il[F_CONSTR].nr > 0)
1129 pd_determine_constraints_range(ir, natoms, &top->idef.il[F_CONSTR], homeind, constr_min_nodeid, constr_max_nodeid);
1133 /* Not 100% necessary, but it is a bad habit to have uninitialized arrays around... */
1134 for (i = 0; i < natoms; i++)
1136 constr_min_nodeid[i] = constr_max_nodeid[i] = homeind[i];
1140 /* Default limits (no communication) for PD constraints */
1142 for (i = 1; i < nnodes; i++)
1144 left_range[i] = top->cgs.index[multinr_cgs[i-1]];
1145 right_range[i-1] = left_range[i]-1;
1147 right_range[nnodes-1] = top->cgs.index[multinr_cgs[nnodes-1]]-1;
1149 for (j = 0; (j < F_NRE); j++)
1151 split_force2(ir, nnodes, homeind, j, &top->idef.il[j], multinr_nre[j], constr_min_nodeid, constr_max_nodeid,
1152 left_range, right_range);
1155 sfree(constr_min_nodeid);
1156 sfree(constr_max_nodeid);
1159 done_blocka(&sblock);