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, 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"
52 /************************************************************
54 * S H I F T U T I L I T I E S
56 ************************************************************/
59 /************************************************************
61 * G R A P H G E N E R A T I O N C O D E
63 ************************************************************/
65 static void add_gbond(t_graph *g, int a0, int a1)
71 inda0 = a0 - g->at_start;
72 inda1 = a1 - g->at_start;
74 /* Search for a direct edge between a0 and a1.
75 * All egdes are bidirectional, so we only need to search one way.
77 for (i = 0; (i < g->nedge[inda0] && !bFound); i++)
79 bFound = (g->edge[inda0][i] == a1);
84 g->edge[inda0][g->nedge[inda0]++] = a1;
85 g->edge[inda1][g->nedge[inda1]++] = a0;
89 static void mk_igraph(t_graph *g, int ftype, const t_ilist *il,
90 int at_start, int at_end,
103 np = interaction_function[ftype].nratoms;
105 if (ia[1] >= at_start && ia[1] < at_end)
107 if (ia[np] >= at_end)
110 "Molecule in topology has atom numbers below and "
111 "above natoms (%d).\n"
112 "You are probably trying to use a trajectory which does "
113 "not match the first %d atoms of the run input file.\n"
114 "You can make a matching run input file with gmx convert-tpr.",
117 if (ftype == F_SETTLE)
119 /* Bond all the atoms in the settle */
120 add_gbond(g, ia[1], ia[2]);
121 add_gbond(g, ia[1], ia[3]);
123 else if (part == nullptr)
125 /* Simply add this bond */
126 for (j = 1; j < np; j++)
128 add_gbond(g, ia[j], ia[j+1]);
133 /* Add this bond when it connects two unlinked parts of the graph */
134 for (j = 1; j < np; j++)
136 if (part[ia[j]] != part[ia[j+1]])
138 add_gbond(g, ia[j], ia[j+1]);
148 gmx_noreturn static void g_error(int line, const char *file)
150 gmx_fatal(FARGS, "Trying to print nonexistent graph (file %s, line %d)",
153 #define GCHECK(g) if (g == NULL) g_error(__LINE__, __FILE__)
155 void p_graph(FILE *log, const char *title, t_graph *g)
158 const char *cc[egcolNR] = { "W", "G", "B" };
161 fprintf(log, "graph: %s\n", title);
162 fprintf(log, "nnodes: %d\n", g->nnodes);
163 fprintf(log, "nbound: %d\n", g->nbound);
164 fprintf(log, "start: %d\n", g->at_start);
165 fprintf(log, "end: %d\n", g->at_end);
166 fprintf(log, " atom shiftx shifty shiftz C nedg e1 e2 etc.\n");
167 for (i = 0; (i < g->nnodes); i++)
171 fprintf(log, "%5d%7d%7d%7d %1s%5d", g->at_start+i+1,
172 g->ishift[g->at_start+i][XX],
173 g->ishift[g->at_start+i][YY],
174 g->ishift[g->at_start+i][ZZ],
175 (g->negc > 0) ? cc[g->egc[i]] : " ",
177 for (j = 0; (j < g->nedge[i]); j++)
179 fprintf(log, " %5d", g->edge[i][j]+1);
187 static void calc_1se(t_graph *g, int ftype, const t_ilist *il,
188 int nbond[], int at_start, int at_end)
190 int k, nratoms, end, j;
196 for (j = 0; (j < end); j += nratoms+1, ia += nratoms+1)
198 nratoms = interaction_function[ftype].nratoms;
200 if (ftype == F_SETTLE)
203 if (iaa >= at_start && iaa < at_end)
208 g->at_start = std::min(g->at_start, iaa);
209 g->at_end = std::max(g->at_end, iaa+2+1);
214 for (k = 1; (k <= nratoms); k++)
217 if (iaa >= at_start && iaa < at_end)
219 g->at_start = std::min(g->at_start, iaa);
220 g->at_end = std::max(g->at_end, iaa+1);
221 /* When making the graph we (might) link all atoms in an interaction
222 * sequentially. Therefore the end atoms add 1 to the count,
223 * the middle atoms 2.
225 if (k == 1 || k == nratoms)
239 static int calc_start_end(FILE *fplog, t_graph *g, const t_ilist il[],
240 int at_start, int at_end,
245 g->at_start = at_end;
248 /* First add all the real bonds: they should determine the molecular
251 for (i = 0; (i < F_NRE); i++)
253 if (interaction_function[i].flags & IF_CHEMBOND)
255 calc_1se(g, i, &il[i], nbond, at_start, at_end);
258 /* Then add all the other interactions in fixed lists, but first
259 * check to see what's there already.
261 for (i = 0; (i < F_NRE); i++)
263 if (!(interaction_function[i].flags & IF_CHEMBOND))
265 calc_1se(g, i, &il[i], nbond, at_start, at_end);
271 for (i = g->at_start; (i < g->at_end); i++)
274 nnb = std::max(nnb, nbond[i]);
278 fprintf(fplog, "Max number of connections per atom is %d\n", nnb);
279 fprintf(fplog, "Total number of connections is %d\n", nbtot);
286 static void compact_graph(FILE *fplog, t_graph *g)
288 int i, j, n, max_nedge;
292 for (i = 0; i < g->nnodes; i++)
294 for (j = 0; j < g->nedge[i]; j++)
296 g->edge[0][n++] = g->edge[i][j];
298 max_nedge = std::max(max_nedge, g->nedge[i]);
300 srenew(g->edge[0], n);
301 /* set pointers after srenew because edge[0] might move */
302 for (i = 1; i < g->nnodes; i++)
304 g->edge[i] = g->edge[i-1] + g->nedge[i-1];
309 fprintf(fplog, "Max number of graph edges per atom is %d\n",
311 fprintf(fplog, "Total number of graph edges is %d\n", n);
315 static gmx_bool determine_graph_parts(t_graph *g, int *part)
322 /* Initialize the part array with all entries different */
323 for (at_i = g->at_start; at_i < g->at_end; at_i++)
328 /* Loop over the graph until the part array is fixed */
333 for (i = 0; (i < g->nnodes); i++)
335 at_i = g->at_start + i;
337 for (e = 0; e < g->nedge[i]; e++)
339 /* Set part for both nodes to the minimum */
340 if (part[at_i2[e]] > part[at_i])
342 part[at_i2[e]] = part[at_i];
345 else if (part[at_i2[e]] < part[at_i])
347 part[at_i] = part[at_i2[e]];
351 if (part[at_i] != part[g->at_start])
358 fprintf(debug, "graph part[] nchanged=%d, bMultiPart=%d\n",
359 nchanged, bMultiPart);
362 while (nchanged > 0);
367 void mk_graph_ilist(FILE *fplog,
368 const t_ilist *ilist, int at_start, int at_end,
369 gmx_bool bShakeOnly, gmx_bool bSettle,
376 /* The naming is somewhat confusing, but we need g->at0 and g->at1
377 * for shifthing coordinates to a new array (not in place) when
378 * some atoms are not connected by the graph, which runs from
379 * g->at_start (>= g->at0) to g->at_end (<= g->at1).
385 nbtot = calc_start_end(fplog, g, ilist, at_start, at_end, nbond);
387 if (g->at_start >= g->at_end)
389 g->at_start = at_start;
396 g->nnodes = g->at_end - g->at_start;
397 snew(g->nedge, g->nnodes);
398 snew(g->edge, g->nnodes);
399 /* Allocate a single array and set pointers into it */
400 snew(g->edge[0], nbtot);
401 for (i = 1; (i < g->nnodes); i++)
403 g->edge[i] = g->edge[i-1] + nbond[g->at_start+i-1];
408 /* First add all the real bonds: they should determine the molecular
411 for (i = 0; (i < F_NRE); i++)
413 if (interaction_function[i].flags & IF_CHEMBOND)
415 mk_igraph(g, i, &(ilist[i]), at_start, at_end, nullptr);
419 /* Determine of which separated parts the IF_CHEMBOND graph consists.
420 * Store the parts in the nbond array.
422 bMultiPart = determine_graph_parts(g, nbond);
426 /* Then add all the other interactions in fixed lists,
427 * but only when they connect parts of the graph
428 * that are not connected through IF_CHEMBOND interactions.
430 for (i = 0; (i < F_NRE); i++)
432 if (!(interaction_function[i].flags & IF_CHEMBOND))
434 mk_igraph(g, i, &(ilist[i]), at_start, at_end, nbond);
439 /* Removed all the unused space from the edge array */
440 compact_graph(fplog, g);
444 /* This is a special thing used in splitter.c to generate shake-blocks */
445 mk_igraph(g, F_CONSTR, &(ilist[F_CONSTR]), at_start, at_end, nullptr);
448 mk_igraph(g, F_SETTLE, &(ilist[F_SETTLE]), at_start, at_end, nullptr);
452 for (i = 0; (i < g->nnodes); i++)
466 snew(g->ishift, g->at1);
470 p_graph(debug, "graph", g);
474 t_graph *mk_graph(FILE *fplog,
475 const t_idef *idef, int at_start, int at_end,
476 gmx_bool bShakeOnly, gmx_bool bSettle)
482 mk_graph_ilist(fplog, idef->il, at_start, at_end, bShakeOnly, bSettle, g);
487 void done_graph(t_graph *g)
500 /************************************************************
502 * S H I F T C A L C U L A T I O N C O D E
504 ************************************************************/
506 static void mk_1shift_tric(int npbcdim, const matrix box, const rvec hbox,
507 const rvec xi, const rvec xj, int *mi, int *mj)
509 /* Calculate periodicity for triclinic box... */
513 rvec_sub(xi, xj, dx);
516 for (m = npbcdim-1; (m >= 0); m--)
518 /* If dx < hbox, then xj will be reduced by box, so that
519 * xi - xj will be bigger
521 if (dx[m] < -hbox[m])
524 for (d = m-1; d >= 0; d--)
529 else if (dx[m] >= hbox[m])
532 for (d = m-1; d >= 0; d--)
544 static void mk_1shift(int npbcdim, const rvec hbox, const rvec xi, const rvec xj,
547 /* Calculate periodicity for rectangular box... */
551 rvec_sub(xi, xj, dx);
554 for (m = 0; (m < npbcdim); m++)
556 /* If dx < hbox, then xj will be reduced by box, so that
557 * xi - xj will be bigger
559 if (dx[m] < -hbox[m])
563 else if (dx[m] >= hbox[m])
574 static void mk_1shift_screw(const matrix box, const rvec hbox,
575 const rvec xi, const rvec xj, int *mi, int *mj)
577 /* Calculate periodicity for rectangular box... */
581 if ((mi[XX] > 0 && mi[XX] % 2 == 1) ||
582 (mi[XX] < 0 && -mi[XX] % 2 == 1))
591 rvec_sub(xi, xj, dx);
593 if (dx[XX] < -hbox[XX])
597 else if (dx[XX] >= hbox[XX])
605 if (mj[XX] != mi[XX])
608 dx[YY] = xi[YY] - (box[YY][YY] + box[ZZ][YY] - xj[YY]);
609 dx[ZZ] = xi[ZZ] - (box[ZZ][ZZ] - xj[ZZ]);
611 for (m = 1; (m < DIM); m++)
613 /* The signs are taken such that we can first shift x and rotate
614 * and then shift y and z.
616 if (dx[m] < -hbox[m])
618 mj[m] = mi[m] - signi;
620 else if (dx[m] >= hbox[m])
622 mj[m] = mi[m] + signi;
631 static int mk_grey(egCol egc[], t_graph *g, int *AtomI,
632 int npbcdim, const matrix box, const rvec x[], int *nerror)
634 int m, j, ng, ai, aj, g0;
640 for (m = 0; (m < DIM); m++)
642 hbox[m] = box[m][m]*0.5;
644 bTriclinic = TRICLINIC(box);
650 /* Loop over all the bonds */
651 for (j = 0; (j < g->nedge[ai-g0]); j++)
653 aj = g->edge[ai-g0][j];
654 /* If there is a white one, make it grey and set pbc */
657 mk_1shift_screw(box, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
661 mk_1shift_tric(npbcdim, box, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
665 mk_1shift(npbcdim, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
668 if (egc[aj-g0] == egcolWhite)
670 if (aj - g0 < *AtomI)
674 egc[aj-g0] = egcolGrey;
676 copy_ivec(is_aj, g->ishift[aj]);
680 else if ((is_aj[XX] != g->ishift[aj][XX]) ||
681 (is_aj[YY] != g->ishift[aj][YY]) ||
682 (is_aj[ZZ] != g->ishift[aj][ZZ]))
686 set_pbc(&pbc, -1, box);
687 pbc_dx(&pbc, x[ai], x[aj], dx);
688 fprintf(debug, "mk_grey: shifts for atom %d due to atom %d\n"
689 "are (%d,%d,%d), should be (%d,%d,%d)\n"
691 aj+1, ai+1, is_aj[XX], is_aj[YY], is_aj[ZZ],
692 g->ishift[aj][XX], g->ishift[aj][YY], g->ishift[aj][ZZ],
693 dx[XX], dx[YY], dx[ZZ]);
701 static int first_colour(int fC, egCol Col, t_graph *g, egCol egc[])
702 /* Return the first node with colour Col starting at fC.
703 * return -1 if none found.
708 for (i = fC; (i < g->nnodes); i++)
710 if ((g->nedge[i] > 0) && (egc[i] == Col))
719 void mk_mshift(FILE *log, t_graph *g, int ePBC,
720 const matrix box, const rvec x[])
722 static int nerror_tot = 0;
725 int nW, nG, nB; /* Number of Grey, Black, White */
726 int fW, fG; /* First of each category */
729 g->bScrewPBC = (ePBC == epbcSCREW);
741 /* This puts everything in the central box, that is does not move it
742 * at all. If we return without doing this for a system without bonds
743 * (i.e. only settles) all water molecules are moved to the opposite octant
745 for (i = g->at0; (i < g->at1); i++)
747 g->ishift[i][XX] = g->ishift[i][YY] = g->ishift[i][ZZ] = 0;
756 if (nnodes > g->negc)
759 srenew(g->egc, g->negc);
761 memset(g->egc, 0, (size_t)(nnodes*sizeof(g->egc[0])));
769 /* We even have a loop invariant:
770 * nW+nG+nB == g->nbound
773 fprintf(log, "Starting W loop\n");
777 /* Find the first white, this will allways be a larger
778 * number than before, because no nodes are made white
781 if ((fW = first_colour(fW, egcolWhite, g, g->egc)) == -1)
783 gmx_fatal(FARGS, "No WHITE nodes found while nW=%d\n", nW);
786 /* Make the first white node grey */
787 g->egc[fW] = egcolGrey;
791 /* Initial value for the first grey */
794 fprintf(log, "Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
795 nW, nG, nB, nW+nG+nB);
799 if ((fG = first_colour(fG, egcolGrey, g, g->egc)) == -1)
801 gmx_fatal(FARGS, "No GREY nodes found while nG=%d\n", nG);
804 /* Make the first grey node black */
805 g->egc[fG] = egcolBlack;
809 /* Make all the neighbours of this black node grey
810 * and set their periodicity
812 ng = mk_grey(g->egc, g, &fG, npbcdim, box, x, &nerror);
813 /* ng is the number of white nodes made grey */
821 if (nerror_tot <= 100)
823 fprintf(stderr, "There were %d inconsistent shifts. Check your topology\n",
827 fprintf(log, "There were %d inconsistent shifts. Check your topology\n",
831 if (nerror_tot == 100)
833 fprintf(stderr, "Will stop reporting inconsistent shifts\n");
836 fprintf(log, "Will stop reporting inconsistent shifts\n");
842 /************************************************************
844 * A C T U A L S H I F T C O D E
846 ************************************************************/
848 void shift_x(const t_graph *g, const matrix box, const rvec x[], rvec x_s[])
859 for (j = g->at0; j < g0; j++)
861 copy_rvec(x[j], x_s[j]);
866 for (j = g0; (j < g1); j++)
872 if ((tx > 0 && tx % 2 == 1) ||
873 (tx < 0 && -tx %2 == 1))
875 x_s[j][XX] = x[j][XX] + tx*box[XX][XX];
876 x_s[j][YY] = box[YY][YY] + box[ZZ][YY] - x[j][YY];
877 x_s[j][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
881 x_s[j][XX] = x[j][XX];
883 x_s[j][YY] = x[j][YY] + ty*box[YY][YY] + tz*box[ZZ][YY];
884 x_s[j][ZZ] = x[j][ZZ] + tz*box[ZZ][ZZ];
887 else if (TRICLINIC(box))
889 for (j = g0; (j < g1); j++)
895 x_s[j][XX] = x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
896 x_s[j][YY] = x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
897 x_s[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
902 for (j = g0; (j < g1); j++)
908 x_s[j][XX] = x[j][XX]+tx*box[XX][XX];
909 x_s[j][YY] = x[j][YY]+ty*box[YY][YY];
910 x_s[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
914 for (j = g1; j < g->at1; j++)
916 copy_rvec(x[j], x_s[j]);
920 void shift_self(const t_graph *g, const matrix box, rvec x[])
928 gmx_incons("screw pbc not implemented for shift_self");
936 fprintf(stderr, "Shifting atoms %d to %d\n", g0, g0+gn);
940 for (j = g0; (j < g1); j++)
946 x[j][XX] = x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
947 x[j][YY] = x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
948 x[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
953 for (j = g0; (j < g1); j++)
959 x[j][XX] = x[j][XX]+tx*box[XX][XX];
960 x[j][YY] = x[j][YY]+ty*box[YY][YY];
961 x[j][ZZ] = x[j][ZZ]+tz*box[ZZ][ZZ];
966 void unshift_x(const t_graph *g, const matrix box, rvec x[], const rvec x_s[])
974 gmx_incons("screw pbc not implemented (yet) for unshift_x");
981 for (j = g->at0; j < g0; j++)
983 copy_rvec(x_s[j], x[j]);
988 for (j = g0; (j < g1); j++)
994 x[j][XX] = x_s[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
995 x[j][YY] = x_s[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
996 x[j][ZZ] = x_s[j][ZZ]-tz*box[ZZ][ZZ];
1001 for (j = g0; (j < g1); j++)
1007 x[j][XX] = x_s[j][XX]-tx*box[XX][XX];
1008 x[j][YY] = x_s[j][YY]-ty*box[YY][YY];
1009 x[j][ZZ] = x_s[j][ZZ]-tz*box[ZZ][ZZ];
1013 for (j = g1; j < g->at1; j++)
1015 copy_rvec(x_s[j], x[j]);
1019 void unshift_self(const t_graph *g, const matrix box, rvec x[])
1027 gmx_incons("screw pbc not implemented for unshift_self");
1036 for (j = g0; (j < g1); j++)
1042 x[j][XX] = x[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1043 x[j][YY] = x[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1044 x[j][ZZ] = x[j][ZZ]-tz*box[ZZ][ZZ];
1049 for (j = g0; (j < g1); j++)
1055 x[j][XX] = x[j][XX]-tx*box[XX][XX];
1056 x[j][YY] = x[j][YY]-ty*box[YY][YY];
1057 x[j][ZZ] = x[j][ZZ]-tz*box[ZZ][ZZ];