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,2014,2015,2016,2017,2018, 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.
45 #include "gromacs/math/vec.h"
46 #include "gromacs/pbcutil/pbc.h"
47 #include "gromacs/topology/idef.h"
48 #include "gromacs/topology/ifunc.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/utility/smalloc.h"
51 #include "gromacs/utility/stringutil.h"
53 /************************************************************
55 * S H I F T U T I L I T I E S
57 ************************************************************/
60 /************************************************************
62 * G R A P H G E N E R A T I O N C O D E
64 ************************************************************/
66 static void add_gbond(t_graph *g, int a0, int a1)
72 inda0 = a0 - g->at_start;
73 inda1 = a1 - g->at_start;
75 /* Search for a direct edge between a0 and a1.
76 * All egdes are bidirectional, so we only need to search one way.
78 for (i = 0; (i < g->nedge[inda0] && !bFound); i++)
80 bFound = (g->edge[inda0][i] == a1);
85 g->edge[inda0][g->nedge[inda0]++] = a1;
86 g->edge[inda1][g->nedge[inda1]++] = a0;
90 /* Make the actual graph for an ilist, returns whether an edge was added.
92 * When a non-null part array is supplied with part indices for each atom,
93 * edges are only added when atoms have a different part index.
95 static bool mk_igraph(t_graph *g, int ftype, const t_ilist *il,
96 int at_start, int at_end,
102 bool addedEdge = false;
110 np = interaction_function[ftype].nratoms;
112 if (np > 1 && ia[1] >= at_start && ia[1] < at_end)
114 if (ia[np] >= at_end)
117 "Molecule in topology has atom numbers below and "
118 "above natoms (%d).\n"
119 "You are probably trying to use a trajectory which does "
120 "not match the first %d atoms of the run input file.\n"
121 "You can make a matching run input file with gmx convert-tpr.",
124 if (ftype == F_SETTLE)
126 /* Bond all the atoms in the settle */
127 add_gbond(g, ia[1], ia[2]);
128 add_gbond(g, ia[1], ia[3]);
131 else if (part == nullptr)
133 /* Simply add this bond */
134 for (j = 1; j < np; j++)
136 add_gbond(g, ia[j], ia[j+1]);
142 /* Add this bond when it connects two unlinked parts of the graph */
143 for (j = 1; j < np; j++)
145 if (part[ia[j]] != part[ia[j+1]])
147 add_gbond(g, ia[j], ia[j+1]);
160 [[noreturn]] static void g_error(int line, const char *file)
162 gmx_fatal(FARGS, "Trying to print nonexistent graph (file %s, line %d)",
165 #define GCHECK(g) if ((g) == NULL) g_error(__LINE__, __FILE__)
167 void p_graph(FILE *log, const char *title, t_graph *g)
170 const char *cc[egcolNR] = { "W", "G", "B" };
173 fprintf(log, "graph: %s\n", title);
174 fprintf(log, "nnodes: %d\n", g->nnodes);
175 fprintf(log, "nbound: %d\n", g->nbound);
176 fprintf(log, "start: %d\n", g->at_start);
177 fprintf(log, "end: %d\n", g->at_end);
178 fprintf(log, " atom shiftx shifty shiftz C nedg e1 e2 etc.\n");
179 for (i = 0; (i < g->nnodes); i++)
183 fprintf(log, "%5d%7d%7d%7d %1s%5d", g->at_start+i+1,
184 g->ishift[g->at_start+i][XX],
185 g->ishift[g->at_start+i][YY],
186 g->ishift[g->at_start+i][ZZ],
187 (g->negc > 0) ? cc[g->egc[i]] : " ",
189 for (j = 0; (j < g->nedge[i]); j++)
191 fprintf(log, " %5d", g->edge[i][j]+1);
199 static void calc_1se(t_graph *g, int ftype, const t_ilist *il,
200 int nbond[], int at_start, int at_end)
202 int k, nratoms, end, j;
208 for (j = 0; (j < end); j += nratoms+1, ia += nratoms+1)
210 nratoms = interaction_function[ftype].nratoms;
212 if (ftype == F_SETTLE)
215 if (iaa >= at_start && iaa < at_end)
220 g->at_start = std::min(g->at_start, iaa);
221 g->at_end = std::max(g->at_end, iaa+2+1);
226 for (k = 1; (k <= nratoms); k++)
229 if (iaa >= at_start && iaa < at_end)
231 g->at_start = std::min(g->at_start, iaa);
232 g->at_end = std::max(g->at_end, iaa+1);
233 /* When making the graph we (might) link all atoms in an interaction
234 * sequentially. Therefore the end atoms add 1 to the count,
235 * the middle atoms 2.
237 if (k == 1 || k == nratoms)
251 static int calc_start_end(FILE *fplog, t_graph *g, const t_ilist il[],
252 int at_start, int at_end,
257 g->at_start = at_end;
260 /* First add all the real bonds: they should determine the molecular
263 for (i = 0; (i < F_NRE); i++)
265 if (interaction_function[i].flags & IF_CHEMBOND)
267 calc_1se(g, i, &il[i], nbond, at_start, at_end);
270 /* Then add all the other interactions in fixed lists, but first
271 * check to see what's there already.
273 for (i = 0; (i < F_NRE); i++)
275 if (!(interaction_function[i].flags & IF_CHEMBOND))
277 calc_1se(g, i, &il[i], nbond, at_start, at_end);
283 for (i = g->at_start; (i < g->at_end); i++)
286 nnb = std::max(nnb, nbond[i]);
290 fprintf(fplog, "Max number of connections per atom is %d\n", nnb);
291 fprintf(fplog, "Total number of connections is %d\n", nbtot);
298 static void compact_graph(FILE *fplog, t_graph *g)
300 int i, j, n, max_nedge;
304 for (i = 0; i < g->nnodes; i++)
306 for (j = 0; j < g->nedge[i]; j++)
308 g->edge[0][n++] = g->edge[i][j];
310 max_nedge = std::max(max_nedge, g->nedge[i]);
312 srenew(g->edge[0], n);
313 /* set pointers after srenew because edge[0] might move */
314 for (i = 1; i < g->nnodes; i++)
316 g->edge[i] = g->edge[i-1] + g->nedge[i-1];
321 fprintf(fplog, "Max number of graph edges per atom is %d\n",
323 fprintf(fplog, "Total number of graph edges is %d\n", n);
327 static gmx_bool determine_graph_parts(t_graph *g, int *part)
334 /* Initialize the part array with all entries different */
335 for (at_i = g->at_start; at_i < g->at_end; at_i++)
340 /* Loop over the graph until the part array is fixed */
345 for (i = 0; (i < g->nnodes); i++)
347 at_i = g->at_start + i;
349 for (e = 0; e < g->nedge[i]; e++)
351 /* Set part for both nodes to the minimum */
352 if (part[at_i2[e]] > part[at_i])
354 part[at_i2[e]] = part[at_i];
357 else if (part[at_i2[e]] < part[at_i])
359 part[at_i] = part[at_i2[e]];
363 if (part[at_i] != part[g->at_start])
370 fprintf(debug, "graph part[] nchanged=%d, bMultiPart=%d\n",
371 nchanged, bMultiPart);
374 while (nchanged > 0);
379 void mk_graph_ilist(FILE *fplog,
380 const t_ilist *ilist, int at_start, int at_end,
381 gmx_bool bShakeOnly, gmx_bool bSettle,
388 /* The naming is somewhat confusing, but we need g->at0 and g->at1
389 * for shifthing coordinates to a new array (not in place) when
390 * some atoms are not connected by the graph, which runs from
391 * g->at_start (>= g->at0) to g->at_end (<= g->at1).
395 g->parts = t_graph::BondedParts::Single;
398 nbtot = calc_start_end(fplog, g, ilist, at_start, at_end, nbond);
400 if (g->at_start >= g->at_end)
402 g->at_start = at_start;
409 g->nnodes = g->at_end - g->at_start;
410 snew(g->nedge, g->nnodes);
411 snew(g->edge, g->nnodes);
412 /* Allocate a single array and set pointers into it */
413 snew(g->edge[0], nbtot);
414 for (i = 1; (i < g->nnodes); i++)
416 g->edge[i] = g->edge[i-1] + nbond[g->at_start+i-1];
421 /* First add all the real bonds: they should determine the molecular
424 for (i = 0; (i < F_NRE); i++)
426 if (interaction_function[i].flags & IF_CHEMBOND)
428 mk_igraph(g, i, &(ilist[i]), at_start, at_end, nullptr);
432 /* Determine of which separated parts the IF_CHEMBOND graph consists.
433 * Store the parts in the nbond array.
435 bMultiPart = determine_graph_parts(g, nbond);
439 /* Then add all the other interactions in fixed lists,
440 * but only when they connect parts of the graph
441 * that are not connected through IF_CHEMBOND interactions.
443 bool addedEdge = false;
444 for (i = 0; (i < F_NRE); i++)
446 if (!(interaction_function[i].flags & IF_CHEMBOND))
448 bool addedEdgeForType =
449 mk_igraph(g, i, &(ilist[i]), at_start, at_end, nbond);
450 addedEdge = (addedEdge || addedEdgeForType);
456 g->parts = t_graph::BondedParts::MultipleConnected;
460 g->parts = t_graph::BondedParts::MultipleDisconnected;
464 /* Removed all the unused space from the edge array */
465 compact_graph(fplog, g);
469 /* This is a special thing used in splitter.c to generate shake-blocks */
470 mk_igraph(g, F_CONSTR, &(ilist[F_CONSTR]), at_start, at_end, nullptr);
473 mk_igraph(g, F_SETTLE, &(ilist[F_SETTLE]), at_start, at_end, nullptr);
477 for (i = 0; (i < g->nnodes); i++)
491 snew(g->ishift, g->at1);
495 p_graph(debug, "graph", g);
499 t_graph *mk_graph(FILE *fplog,
500 const t_idef *idef, int at_start, int at_end,
501 gmx_bool bShakeOnly, gmx_bool bSettle)
507 mk_graph_ilist(fplog, idef->il, at_start, at_end, bShakeOnly, bSettle, g);
512 void done_graph(t_graph *g)
525 /************************************************************
527 * S H I F T C A L C U L A T I O N C O D E
529 ************************************************************/
531 static void mk_1shift_tric(int npbcdim, const matrix box, const rvec hbox,
532 const rvec xi, const rvec xj, const int *mi, int *mj)
534 /* Calculate periodicity for triclinic box... */
538 rvec_sub(xi, xj, dx);
541 for (m = npbcdim-1; (m >= 0); m--)
543 /* If dx < hbox, then xj will be reduced by box, so that
544 * xi - xj will be bigger
546 if (dx[m] < -hbox[m])
549 for (d = m-1; d >= 0; d--)
554 else if (dx[m] >= hbox[m])
557 for (d = m-1; d >= 0; d--)
569 static void mk_1shift(int npbcdim, const rvec hbox, const rvec xi, const rvec xj,
570 const int *mi, int *mj)
572 /* Calculate periodicity for rectangular box... */
576 rvec_sub(xi, xj, dx);
579 for (m = 0; (m < npbcdim); m++)
581 /* If dx < hbox, then xj will be reduced by box, so that
582 * xi - xj will be bigger
584 if (dx[m] < -hbox[m])
588 else if (dx[m] >= hbox[m])
599 static void mk_1shift_screw(const matrix box, const rvec hbox,
600 const rvec xi, const rvec xj, const int *mi, int *mj)
602 /* Calculate periodicity for rectangular box... */
606 if ((mi[XX] > 0 && mi[XX] % 2 == 1) ||
607 (mi[XX] < 0 && -mi[XX] % 2 == 1))
616 rvec_sub(xi, xj, dx);
618 if (dx[XX] < -hbox[XX])
622 else if (dx[XX] >= hbox[XX])
630 if (mj[XX] != mi[XX])
633 dx[YY] = xi[YY] - (box[YY][YY] + box[ZZ][YY] - xj[YY]);
634 dx[ZZ] = xi[ZZ] - (box[ZZ][ZZ] - xj[ZZ]);
636 for (m = 1; (m < DIM); m++)
638 /* The signs are taken such that we can first shift x and rotate
639 * and then shift y and z.
641 if (dx[m] < -hbox[m])
643 mj[m] = mi[m] - signi;
645 else if (dx[m] >= hbox[m])
647 mj[m] = mi[m] + signi;
656 static int mk_grey(egCol egc[], t_graph *g, int *AtomI,
657 int npbcdim, const matrix box, const rvec x[], int *nerror)
659 int m, j, ng, ai, aj, g0;
665 for (m = 0; (m < DIM); m++)
667 hbox[m] = box[m][m]*0.5;
669 bTriclinic = TRICLINIC(box);
675 /* Loop over all the bonds */
676 for (j = 0; (j < g->nedge[ai-g0]); j++)
678 aj = g->edge[ai-g0][j];
679 /* If there is a white one, make it grey and set pbc */
682 mk_1shift_screw(box, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
686 mk_1shift_tric(npbcdim, box, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
690 mk_1shift(npbcdim, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
693 if (egc[aj-g0] == egcolWhite)
695 if (aj - g0 < *AtomI)
699 egc[aj-g0] = egcolGrey;
701 copy_ivec(is_aj, g->ishift[aj]);
705 else if ((is_aj[XX] != g->ishift[aj][XX]) ||
706 (is_aj[YY] != g->ishift[aj][YY]) ||
707 (is_aj[ZZ] != g->ishift[aj][ZZ]))
711 set_pbc(&pbc, -1, box);
712 pbc_dx(&pbc, x[ai], x[aj], dx);
713 fprintf(debug, "mk_grey: shifts for atom %d due to atom %d\n"
714 "are (%d,%d,%d), should be (%d,%d,%d)\n"
716 aj+1, ai+1, is_aj[XX], is_aj[YY], is_aj[ZZ],
717 g->ishift[aj][XX], g->ishift[aj][YY], g->ishift[aj][ZZ],
718 dx[XX], dx[YY], dx[ZZ]);
726 static int first_colour(int fC, egCol Col, t_graph *g, const egCol egc[])
727 /* Return the first node with colour Col starting at fC.
728 * return -1 if none found.
733 for (i = fC; (i < g->nnodes); i++)
735 if ((g->nedge[i] > 0) && (egc[i] == Col))
744 /* Returns the maximum length of the graph edges for coordinates x */
745 static real maxEdgeLength(const t_graph g,
752 set_pbc(&pbc, ePBC, box);
754 real maxEdgeLength2 = 0;
756 for (int node = 0; node < g.nnodes; node++)
758 for (int edge = 0; edge < g.nedge[node]; edge++)
760 int nodeJ = g.edge[node][edge];
762 pbc_dx(&pbc, x[g.at0 + node], x[g.at0 + nodeJ], dx);
763 maxEdgeLength2 = std::max(maxEdgeLength2, norm2(dx));
767 return std::sqrt(maxEdgeLength2);
770 void mk_mshift(FILE *log, t_graph *g, int ePBC,
771 const matrix box, const rvec x[])
773 static int nerror_tot = 0;
776 int nW, nG, nB; /* Number of Grey, Black, White */
777 int fW, fG; /* First of each category */
780 g->bScrewPBC = (ePBC == epbcSCREW);
792 /* This puts everything in the central box, that is does not move it
793 * at all. If we return without doing this for a system without bonds
794 * (i.e. only settles) all water molecules are moved to the opposite octant
796 for (i = g->at0; (i < g->at1); i++)
798 g->ishift[i][XX] = g->ishift[i][YY] = g->ishift[i][ZZ] = 0;
807 if (nnodes > g->negc)
810 srenew(g->egc, g->negc);
812 memset(g->egc, 0, (size_t)(nnodes*sizeof(g->egc[0])));
820 /* We even have a loop invariant:
821 * nW+nG+nB == g->nbound
824 fprintf(log, "Starting W loop\n");
828 /* Find the first white, this will allways be a larger
829 * number than before, because no nodes are made white
832 if ((fW = first_colour(fW, egcolWhite, g, g->egc)) == -1)
834 gmx_fatal(FARGS, "No WHITE nodes found while nW=%d\n", nW);
837 /* Make the first white node grey */
838 g->egc[fW] = egcolGrey;
842 /* Initial value for the first grey */
845 fprintf(log, "Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
846 nW, nG, nB, nW+nG+nB);
850 if ((fG = first_colour(fG, egcolGrey, g, g->egc)) == -1)
852 gmx_fatal(FARGS, "No GREY nodes found while nG=%d\n", nG);
855 /* Make the first grey node black */
856 g->egc[fG] = egcolBlack;
860 /* Make all the neighbours of this black node grey
861 * and set their periodicity
863 ng = mk_grey(g->egc, g, &fG, npbcdim, box, x, &nerror);
864 /* ng is the number of white nodes made grey */
871 /* We use a threshold of 0.25*boxSize for generating a fatal error
872 * to allow removing PBC for systems with periodic molecules.
874 * TODO: Consider a better solution for systems with periodic
875 * molecules. Ideally analysis tools should only ask to make
876 * non-periodic molecules whole in a system with periodic mols.
878 constexpr real c_relativeDistanceThreshold = 0.25;
880 int numPbcDimensions = ePBC2npbcdim(ePBC);
881 GMX_RELEASE_ASSERT(numPbcDimensions > 0, "Expect PBC with graph");
882 real minBoxSize = norm(box[XX]);
883 for (int d = 1; d < numPbcDimensions; d++)
885 minBoxSize = std::min(minBoxSize, norm(box[d]));
887 real maxDistance = maxEdgeLength(*g, ePBC, box, x);
888 if (maxDistance >= c_relativeDistanceThreshold*minBoxSize)
890 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.",
891 g->at1 - g->at0, maxDistance,
892 maxDistance >= 0.5*minBoxSize ? "above" : "close to");
896 case t_graph::BondedParts::MultipleConnected:
897 /* Ideally we should check if the long distances are
898 * actually between the parts, but that would require
899 * a lot of extra code.
901 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.";
904 mesg += " Either you have excessively large distances between atoms in bonded interactions or your system is exploding.";
906 gmx_fatal(FARGS, "%s", mesg.c_str());
909 /* The most likely reason for arriving here is a periodic molecule. */
912 if (nerror_tot <= 100)
914 fprintf(stderr, "There were %d inconsistent shifts. Check your topology\n",
918 fprintf(log, "There were %d inconsistent shifts. Check your topology\n",
922 if (nerror_tot == 100)
924 fprintf(stderr, "Will stop reporting inconsistent shifts\n");
927 fprintf(log, "Will stop reporting inconsistent shifts\n");
933 /************************************************************
935 * A C T U A L S H I F T C O D E
937 ************************************************************/
939 void shift_x(const t_graph *g, const matrix box, const rvec x[], rvec x_s[])
950 for (j = g->at0; j < g0; j++)
952 copy_rvec(x[j], x_s[j]);
957 for (j = g0; (j < g1); j++)
963 if ((tx > 0 && tx % 2 == 1) ||
964 (tx < 0 && -tx %2 == 1))
966 x_s[j][XX] = x[j][XX] + tx*box[XX][XX];
967 x_s[j][YY] = box[YY][YY] + box[ZZ][YY] - x[j][YY];
968 x_s[j][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
972 x_s[j][XX] = x[j][XX];
974 x_s[j][YY] = x[j][YY] + ty*box[YY][YY] + tz*box[ZZ][YY];
975 x_s[j][ZZ] = x[j][ZZ] + tz*box[ZZ][ZZ];
978 else if (TRICLINIC(box))
980 for (j = g0; (j < g1); j++)
986 x_s[j][XX] = x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
987 x_s[j][YY] = x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
988 x_s[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
993 for (j = g0; (j < g1); j++)
999 x_s[j][XX] = x[j][XX]+tx*box[XX][XX];
1000 x_s[j][YY] = x[j][YY]+ty*box[YY][YY];
1001 x_s[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
1005 for (j = g1; j < g->at1; j++)
1007 copy_rvec(x[j], x_s[j]);
1011 void shift_self(const t_graph *g, const matrix box, rvec x[])
1019 gmx_incons("screw pbc not implemented for shift_self");
1027 fprintf(stderr, "Shifting atoms %d to %d\n", g0, g0+gn);
1031 for (j = g0; (j < g1); j++)
1037 x[j][XX] = x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
1038 x[j][YY] = x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
1039 x[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
1044 for (j = g0; (j < g1); j++)
1050 x[j][XX] = x[j][XX]+tx*box[XX][XX];
1051 x[j][YY] = x[j][YY]+ty*box[YY][YY];
1052 x[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
1057 void unshift_x(const t_graph *g, const matrix box, rvec x[], const rvec x_s[])
1065 gmx_incons("screw pbc not implemented (yet) for unshift_x");
1072 for (j = g->at0; j < g0; j++)
1074 copy_rvec(x_s[j], x[j]);
1079 for (j = g0; (j < g1); j++)
1085 x[j][XX] = x_s[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1086 x[j][YY] = x_s[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1087 x[j][ZZ] = x_s[j][ZZ]-tz*box[ZZ][ZZ];
1092 for (j = g0; (j < g1); j++)
1098 x[j][XX] = x_s[j][XX]-tx*box[XX][XX];
1099 x[j][YY] = x_s[j][YY]-ty*box[YY][YY];
1100 x[j][ZZ] = x_s[j][ZZ]-tz*box[ZZ][ZZ];
1104 for (j = g1; j < g->at1; j++)
1106 copy_rvec(x_s[j], x[j]);
1110 void unshift_self(const t_graph *g, const matrix box, rvec x[])
1118 gmx_incons("screw pbc not implemented for unshift_self");
1127 for (j = g0; (j < g1); j++)
1133 x[j][XX] = x[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1134 x[j][YY] = x[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1135 x[j][ZZ] = x[j][ZZ]-tz*box[ZZ][ZZ];
1140 for (j = g0; (j < g1); j++)
1146 x[j][XX] = x[j][XX]-tx*box[XX][XX];
1147 x[j][YY] = x[j][YY]-ty*box[YY][YY];
1148 x[j][ZZ] = x[j][ZZ]-tz*box[ZZ][ZZ];