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.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
40 * Implements routines in pbc.h.
42 * Utility functions for handling periodic boundary conditions.
43 * Mainly used in analysis tools.
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/math/vecdump.h"
58 #include "gromacs/pbcutil/ishift.h"
59 #include "gromacs/pbcutil/mshift.h"
60 #include "gromacs/topology/topology.h"
61 #include "gromacs/utility/exceptions.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/gmxassert.h"
64 #include "gromacs/utility/smalloc.h"
66 const gmx::EnumerationArray<PbcType, std::string> c_pbcTypeNames = {
67 { "xyz", "no", "xy", "screw", "unset" }
70 /* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
73 epbcdxRECTANGULAR = 1,
85 //! Margin factor for error message
86 #define BOX_MARGIN 1.0010
87 //! Margin correction if the box is too skewed
88 #define BOX_MARGIN_CORRECT 1.0005
90 int numPbcDimensions(PbcType pbcType)
97 gmx_fatal(FARGS, "Number of PBC dimensions was requested before the PBC type set.");
98 case PbcType::Xyz: npbcdim = 3; break;
99 case PbcType::XY: npbcdim = 2; break;
100 case PbcType::Screw: npbcdim = 3; break;
101 case PbcType::No: npbcdim = 0; break;
103 gmx_fatal(FARGS, "Unknown pbcType=%s in numPbcDimensions", c_pbcTypeNames[pbcType].c_str());
109 void dump_pbc(FILE* fp, t_pbc* pbc)
113 fprintf(fp, "pbcTypeDX = %d\n", pbc->pbcTypeDX);
114 pr_rvecs(fp, 0, "box", pbc->box, DIM);
115 pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
116 pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
117 pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
118 rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
119 pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
120 fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
121 fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
122 if (pbc->ntric_vec > 0)
124 pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE);
125 pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
129 const char* check_box(PbcType pbcType, const matrix box)
133 if (pbcType == PbcType::Unset)
135 pbcType = guessPbcType(box);
138 if (pbcType == PbcType::No)
143 if ((box[XX][YY] != 0) || (box[XX][ZZ] != 0) || (box[YY][ZZ] != 0))
145 ptr = "Only triclinic boxes with the first vector parallel to the x-axis and the second "
146 "vector in the xy-plane are supported.";
148 else if (pbcType == PbcType::Screw && (box[YY][XX] != 0 || box[ZZ][XX] != 0))
150 ptr = "The unit cell can not have off-diagonal x-components with screw pbc";
152 else if (std::fabs(box[YY][XX]) > BOX_MARGIN * 0.5 * box[XX][XX]
153 || (pbcType != PbcType::XY
154 && (std::fabs(box[ZZ][XX]) > BOX_MARGIN * 0.5 * box[XX][XX]
155 || std::fabs(box[ZZ][YY]) > BOX_MARGIN * 0.5 * box[YY][YY])))
157 ptr = "Triclinic box is too skewed.";
167 void matrix_convert(matrix box, const rvec vec, const rvec angleInDegrees)
170 svmul(DEG2RAD, angleInDegrees, angle);
171 box[XX][XX] = vec[XX];
172 box[YY][XX] = vec[YY] * cos(angle[ZZ]);
173 box[YY][YY] = vec[YY] * sin(angle[ZZ]);
174 box[ZZ][XX] = vec[ZZ] * cos(angle[YY]);
175 box[ZZ][YY] = vec[ZZ] * (cos(angle[XX]) - cos(angle[YY]) * cos(angle[ZZ])) / sin(angle[ZZ]);
177 std::sqrt(gmx::square(vec[ZZ]) - box[ZZ][XX] * box[ZZ][XX] - box[ZZ][YY] * box[ZZ][YY]);
180 real max_cutoff2(PbcType pbcType, const matrix box)
182 real min_hv2, min_ss;
183 const real oneFourth = 0.25;
185 /* Physical limitation of the cut-off
186 * by half the length of the shortest box vector.
188 min_hv2 = oneFourth * std::min(norm2(box[XX]), norm2(box[YY]));
189 if (pbcType != PbcType::XY)
191 min_hv2 = std::min(min_hv2, oneFourth * norm2(box[ZZ]));
194 /* Limitation to the smallest diagonal element due to optimizations:
195 * checking only linear combinations of single box-vectors (2 in x)
196 * in the grid search and pbc_dx is a lot faster
197 * than checking all possible combinations.
199 if (pbcType == PbcType::XY)
201 min_ss = std::min(box[XX][XX], box[YY][YY]);
205 min_ss = std::min(box[XX][XX], std::min(box[YY][YY] - std::fabs(box[ZZ][YY]), box[ZZ][ZZ]));
208 return std::min(min_hv2, min_ss * min_ss);
211 //! Set to true if warning has been printed
212 static gmx_bool bWarnedGuess = FALSE;
214 PbcType guessPbcType(const matrix box)
218 if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] > 0)
220 pbcType = PbcType::Xyz;
222 else if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] == 0)
224 pbcType = PbcType::XY;
226 else if (box[XX][XX] == 0 && box[YY][YY] == 0 && box[ZZ][ZZ] == 0)
228 pbcType = PbcType::No;
235 "WARNING: Unsupported box diagonal %f %f %f, "
236 "will not use periodic boundary conditions\n\n",
242 pbcType = PbcType::No;
247 fprintf(debug, "Guessed pbc = %s from the box matrix\n", c_pbcTypeNames[pbcType].c_str());
253 //! Check if the box still obeys the restrictions, if not, correct it
254 static int correct_box_elem(FILE* fplog, int step, tensor box, int v, int d)
256 int shift, maxshift = 10;
260 /* correct elem d of vector v with vector d */
261 while (box[v][d] > BOX_MARGIN_CORRECT * 0.5 * box[d][d])
265 fprintf(fplog, "Step %d: correcting invalid box:\n", step);
266 pr_rvecs(fplog, 0, "old box", box, DIM);
268 rvec_dec(box[v], box[d]);
272 pr_rvecs(fplog, 0, "new box", box, DIM);
274 if (shift <= -maxshift)
276 gmx_fatal(FARGS, "Box was shifted at least %d times. Please see log-file.", maxshift);
279 while (box[v][d] < -BOX_MARGIN_CORRECT * 0.5 * box[d][d])
283 fprintf(fplog, "Step %d: correcting invalid box:\n", step);
284 pr_rvecs(fplog, 0, "old box", box, DIM);
286 rvec_inc(box[v], box[d]);
290 pr_rvecs(fplog, 0, "new box", box, DIM);
292 if (shift >= maxshift)
294 gmx_fatal(FARGS, "Box was shifted at least %d times. Please see log-file.", maxshift);
301 gmx_bool correct_box(FILE* fplog, int step, tensor box)
306 zy = correct_box_elem(fplog, step, box, ZZ, YY);
307 zx = correct_box_elem(fplog, step, box, ZZ, XX);
308 yx = correct_box_elem(fplog, step, box, YY, XX);
310 bCorrected = ((zy != 0) || (zx != 0) || (yx != 0));
315 //! Do the real arithmetic for filling the pbc struct
316 static void low_set_pbc(t_pbc* pbc, PbcType pbcType, const ivec dd_pbc, const matrix box)
318 int order[3] = { 0, -1, 1 };
322 pbc->pbcType = pbcType;
323 pbc->ndim_ePBC = numPbcDimensions(pbcType);
325 if (pbc->pbcType == PbcType::No)
327 pbc->pbcTypeDX = epbcdxNOPBC;
332 copy_mat(box, pbc->box);
333 pbc->max_cutoff2 = 0;
337 for (int i = 0; (i < DIM); i++)
339 pbc->fbox_diag[i] = box[i][i];
340 pbc->hbox_diag[i] = pbc->fbox_diag[i] * 0.5;
341 pbc->mhbox_diag[i] = -pbc->hbox_diag[i];
344 ptr = check_box(pbcType, box);
347 fprintf(stderr, "Warning: %s\n", ptr);
348 pr_rvecs(stderr, 0, " Box", box, DIM);
349 fprintf(stderr, " Can not fix pbc.\n\n");
350 pbc->pbcTypeDX = epbcdxUNSUPPORTED;
354 if (pbcType == PbcType::Screw && nullptr != dd_pbc)
356 /* This combinated should never appear here */
357 gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
361 for (int i = 0; i < DIM; i++)
363 if ((dd_pbc && dd_pbc[i] == 0) || (pbcType == PbcType::XY && i == ZZ))
376 /* 1D pbc is not an mdp option and it is therefore only used
377 * with single shifts.
379 pbc->pbcTypeDX = epbcdx1D_RECT;
380 for (int i = 0; i < DIM; i++)
387 GMX_ASSERT(pbc->dim < DIM, "Dimension for PBC incorrect");
388 for (int i = 0; i < pbc->dim; i++)
390 if (pbc->box[pbc->dim][i] != 0)
392 pbc->pbcTypeDX = epbcdx1D_TRIC;
397 pbc->pbcTypeDX = epbcdx2D_RECT;
398 for (int i = 0; i < DIM; i++)
405 for (int i = 0; i < DIM; i++)
409 for (int j = 0; j < i; j++)
411 if (pbc->box[i][j] != 0)
413 pbc->pbcTypeDX = epbcdx2D_TRIC;
420 if (pbcType != PbcType::Screw)
424 pbc->pbcTypeDX = epbcdxTRICLINIC;
428 pbc->pbcTypeDX = epbcdxRECTANGULAR;
433 pbc->pbcTypeDX = (box[ZZ][YY] == 0 ? epbcdxSCREW_RECT : epbcdxSCREW_TRIC);
434 if (pbc->pbcTypeDX == epbcdxSCREW_TRIC)
437 "Screw pbc is not yet implemented for triclinic boxes.\n"
438 "Can not fix pbc.\n");
439 pbc->pbcTypeDX = epbcdxUNSUPPORTED;
443 default: gmx_fatal(FARGS, "Incorrect number of pbc dimensions with DD: %d", npbcdim);
445 pbc->max_cutoff2 = max_cutoff2(pbcType, box);
447 if (pbc->pbcTypeDX == epbcdxTRICLINIC || pbc->pbcTypeDX == epbcdx2D_TRIC
448 || pbc->pbcTypeDX == epbcdxSCREW_TRIC)
452 pr_rvecs(debug, 0, "Box", box, DIM);
453 fprintf(debug, "max cutoff %.3f\n", sqrt(pbc->max_cutoff2));
455 /* We will only need single shifts here */
456 for (int kk = 0; kk < 3; kk++)
459 if (!bPBC[ZZ] && k != 0)
463 for (int jj = 0; jj < 3; jj++)
466 if (!bPBC[YY] && j != 0)
470 for (int ii = 0; ii < 3; ii++)
473 if (!bPBC[XX] && i != 0)
477 /* A shift is only useful when it is trilinic */
478 if (j != 0 || k != 0)
485 for (int d = 0; d < DIM; d++)
487 trial[d] = i * box[XX][d] + j * box[YY][d] + k * box[ZZ][d];
488 /* Choose the vector within the brick around 0,0,0 that
489 * will become the shortest due to shift try.
500 pos[d] = std::min(pbc->hbox_diag[d], -trial[d]);
504 pos[d] = std::max(-pbc->hbox_diag[d], -trial[d]);
507 d2old += gmx::square(pos[d]);
508 d2new += gmx::square(pos[d] + trial[d]);
510 if (BOX_MARGIN * d2new < d2old)
512 /* Check if shifts with one box vector less do better */
513 gmx_bool bUse = TRUE;
514 for (int dd = 0; dd < DIM; dd++)
516 int shift = (dd == 0 ? i : (dd == 1 ? j : k));
520 for (int d = 0; d < DIM; d++)
522 d2new_c += gmx::square(pos[d] + trial[d] - shift * box[dd][d]);
524 if (d2new_c <= BOX_MARGIN * d2new)
532 /* Accept this shift vector. */
533 if (pbc->ntric_vec >= MAX_NTRICVEC)
536 "\nWARNING: Found more than %d triclinic "
537 "correction vectors, ignoring some.\n"
538 " There is probably something wrong with your "
541 pr_rvecs(stderr, 0, " Box", box, DIM);
545 copy_rvec(trial, pbc->tric_vec[pbc->ntric_vec]);
546 pbc->tric_shift[pbc->ntric_vec][XX] = i;
547 pbc->tric_shift[pbc->ntric_vec][YY] = j;
548 pbc->tric_shift[pbc->ntric_vec][ZZ] = k;
554 " tricvec %2d = %2d %2d %2d %5.2f %5.2f "
555 "%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
580 void set_pbc(t_pbc* pbc, PbcType pbcType, const matrix box)
582 if (pbcType == PbcType::Unset)
584 pbcType = guessPbcType(box);
587 low_set_pbc(pbc, pbcType, nullptr, box);
590 t_pbc* set_pbc_dd(t_pbc* pbc, PbcType pbcType, const ivec domdecCells, gmx_bool bSingleDir, const matrix box)
592 if (pbcType == PbcType::No)
594 pbc->pbcType = pbcType;
599 if (nullptr == domdecCells)
601 low_set_pbc(pbc, pbcType, nullptr, box);
605 if (pbcType == PbcType::Screw && domdecCells[XX] > 1)
607 /* The rotation has been taken care of during coordinate communication */
608 pbcType = PbcType::Xyz;
613 for (int i = 0; i < DIM; i++)
616 if (domdecCells[i] <= (bSingleDir ? 1 : 2) && !(pbcType == PbcType::XY && i == ZZ))
625 low_set_pbc(pbc, pbcType, usePBC, box);
629 pbc->pbcType = PbcType::No;
633 return (pbc->pbcType != PbcType::No ? pbc : nullptr);
636 void pbc_dx(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx)
639 rvec dx_start, trial;
643 rvec_sub(x1, x2, dx);
645 switch (pbc->pbcTypeDX)
647 case epbcdxRECTANGULAR:
648 for (i = 0; i < DIM; i++)
650 while (dx[i] > pbc->hbox_diag[i])
652 dx[i] -= pbc->fbox_diag[i];
654 while (dx[i] <= pbc->mhbox_diag[i])
656 dx[i] += pbc->fbox_diag[i];
660 case epbcdxTRICLINIC:
661 for (i = DIM - 1; i >= 0; i--)
663 while (dx[i] > pbc->hbox_diag[i])
665 for (j = i; j >= 0; j--)
667 dx[j] -= pbc->box[i][j];
670 while (dx[i] <= pbc->mhbox_diag[i])
672 for (j = i; j >= 0; j--)
674 dx[j] += pbc->box[i][j];
678 /* dx is the distance in a rectangular box */
680 if (d2min > pbc->max_cutoff2)
682 copy_rvec(dx, dx_start);
684 /* Now try all possible shifts, when the distance is within max_cutoff
685 * it must be the shortest possible distance.
688 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
690 rvec_add(dx_start, pbc->tric_vec[i], trial);
691 d2trial = norm2(trial);
694 copy_rvec(trial, dx);
702 for (i = 0; i < DIM; i++)
706 while (dx[i] > pbc->hbox_diag[i])
708 dx[i] -= pbc->fbox_diag[i];
710 while (dx[i] <= pbc->mhbox_diag[i])
712 dx[i] += pbc->fbox_diag[i];
719 for (i = DIM - 1; i >= 0; i--)
723 while (dx[i] > pbc->hbox_diag[i])
725 for (j = i; j >= 0; j--)
727 dx[j] -= pbc->box[i][j];
730 while (dx[i] <= pbc->mhbox_diag[i])
732 for (j = i; j >= 0; j--)
734 dx[j] += pbc->box[i][j];
737 d2min += dx[i] * dx[i];
740 if (d2min > pbc->max_cutoff2)
742 copy_rvec(dx, dx_start);
744 /* Now try all possible shifts, when the distance is within max_cutoff
745 * it must be the shortest possible distance.
748 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
750 rvec_add(dx_start, pbc->tric_vec[i], trial);
752 for (j = 0; j < DIM; j++)
756 d2trial += trial[j] * trial[j];
761 copy_rvec(trial, dx);
768 case epbcdxSCREW_RECT:
769 /* The shift definition requires x first */
771 while (dx[XX] > pbc->hbox_diag[XX])
773 dx[XX] -= pbc->fbox_diag[XX];
776 while (dx[XX] <= pbc->mhbox_diag[XX])
778 dx[XX] += pbc->fbox_diag[YY];
783 /* Rotate around the x-axis in the middle of the box */
784 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
785 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
787 /* Normal pbc for y and z */
788 for (i = YY; i <= ZZ; i++)
790 while (dx[i] > pbc->hbox_diag[i])
792 dx[i] -= pbc->fbox_diag[i];
794 while (dx[i] <= pbc->mhbox_diag[i])
796 dx[i] += pbc->fbox_diag[i];
801 case epbcdxUNSUPPORTED: break;
802 default: gmx_fatal(FARGS, "Internal error in pbc_dx, set_pbc has not been called");
806 int pbc_dx_aiuc(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx)
809 rvec dx_start, trial;
811 ivec ishift, ishift_start;
813 rvec_sub(x1, x2, dx);
816 switch (pbc->pbcTypeDX)
818 case epbcdxRECTANGULAR:
819 for (i = 0; i < DIM; i++)
821 if (dx[i] > pbc->hbox_diag[i])
823 dx[i] -= pbc->fbox_diag[i];
826 else if (dx[i] <= pbc->mhbox_diag[i])
828 dx[i] += pbc->fbox_diag[i];
833 case epbcdxTRICLINIC:
834 /* For triclinic boxes the performance difference between
835 * if/else and two while loops is negligible.
836 * However, the while version can cause extreme delays
837 * before a simulation crashes due to large forces which
838 * can cause unlimited displacements.
839 * Also allowing multiple shifts would index fshift beyond bounds.
841 for (i = DIM - 1; i >= 1; i--)
843 if (dx[i] > pbc->hbox_diag[i])
845 for (j = i; j >= 0; j--)
847 dx[j] -= pbc->box[i][j];
851 else if (dx[i] <= pbc->mhbox_diag[i])
853 for (j = i; j >= 0; j--)
855 dx[j] += pbc->box[i][j];
860 /* Allow 2 shifts in x */
861 if (dx[XX] > pbc->hbox_diag[XX])
863 dx[XX] -= pbc->fbox_diag[XX];
865 if (dx[XX] > pbc->hbox_diag[XX])
867 dx[XX] -= pbc->fbox_diag[XX];
871 else if (dx[XX] <= pbc->mhbox_diag[XX])
873 dx[XX] += pbc->fbox_diag[XX];
875 if (dx[XX] <= pbc->mhbox_diag[XX])
877 dx[XX] += pbc->fbox_diag[XX];
881 /* dx is the distance in a rectangular box */
883 if (d2min > pbc->max_cutoff2)
885 copy_rvec(dx, dx_start);
886 copy_ivec(ishift, ishift_start);
888 /* Now try all possible shifts, when the distance is within max_cutoff
889 * it must be the shortest possible distance.
892 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
894 rvec_add(dx_start, pbc->tric_vec[i], trial);
895 d2trial = norm2(trial);
898 copy_rvec(trial, dx);
899 ivec_add(ishift_start, pbc->tric_shift[i], ishift);
907 for (i = 0; i < DIM; i++)
911 if (dx[i] > pbc->hbox_diag[i])
913 dx[i] -= pbc->fbox_diag[i];
916 else if (dx[i] <= pbc->mhbox_diag[i])
918 dx[i] += pbc->fbox_diag[i];
926 for (i = DIM - 1; i >= 1; i--)
930 if (dx[i] > pbc->hbox_diag[i])
932 for (j = i; j >= 0; j--)
934 dx[j] -= pbc->box[i][j];
938 else if (dx[i] <= pbc->mhbox_diag[i])
940 for (j = i; j >= 0; j--)
942 dx[j] += pbc->box[i][j];
946 d2min += dx[i] * dx[i];
951 /* Allow 2 shifts in x */
952 if (dx[XX] > pbc->hbox_diag[XX])
954 dx[XX] -= pbc->fbox_diag[XX];
956 if (dx[XX] > pbc->hbox_diag[XX])
958 dx[XX] -= pbc->fbox_diag[XX];
962 else if (dx[XX] <= pbc->mhbox_diag[XX])
964 dx[XX] += pbc->fbox_diag[XX];
966 if (dx[XX] <= pbc->mhbox_diag[XX])
968 dx[XX] += pbc->fbox_diag[XX];
972 d2min += dx[XX] * dx[XX];
974 if (d2min > pbc->max_cutoff2)
976 copy_rvec(dx, dx_start);
977 copy_ivec(ishift, ishift_start);
978 /* Now try all possible shifts, when the distance is within max_cutoff
979 * it must be the shortest possible distance.
982 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
984 rvec_add(dx_start, pbc->tric_vec[i], trial);
986 for (j = 0; j < DIM; j++)
990 d2trial += trial[j] * trial[j];
995 copy_rvec(trial, dx);
996 ivec_add(ishift_start, pbc->tric_shift[i], ishift);
1005 if (dx[i] > pbc->hbox_diag[i])
1007 dx[i] -= pbc->fbox_diag[i];
1010 else if (dx[i] <= pbc->mhbox_diag[i])
1012 dx[i] += pbc->fbox_diag[i];
1018 if (dx[i] > pbc->hbox_diag[i])
1020 rvec_dec(dx, pbc->box[i]);
1023 else if (dx[i] <= pbc->mhbox_diag[i])
1025 rvec_inc(dx, pbc->box[i]);
1029 case epbcdxSCREW_RECT:
1030 /* The shift definition requires x first */
1031 if (dx[XX] > pbc->hbox_diag[XX])
1033 dx[XX] -= pbc->fbox_diag[XX];
1036 else if (dx[XX] <= pbc->mhbox_diag[XX])
1038 dx[XX] += pbc->fbox_diag[XX];
1041 if (ishift[XX] == 1 || ishift[XX] == -1)
1043 /* Rotate around the x-axis in the middle of the box */
1044 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1045 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1047 /* Normal pbc for y and z */
1048 for (i = YY; i <= ZZ; i++)
1050 if (dx[i] > pbc->hbox_diag[i])
1052 dx[i] -= pbc->fbox_diag[i];
1055 else if (dx[i] <= pbc->mhbox_diag[i])
1057 dx[i] += pbc->fbox_diag[i];
1063 case epbcdxUNSUPPORTED: break;
1066 "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
1069 is = IVEC2IS(ishift);
1072 range_check_mesg(is, 0, SHIFTS, "PBC shift vector index range check.");
1078 //! Compute distance vector in double precision
1079 void pbc_dx_d(const t_pbc* pbc, const dvec x1, const dvec x2, dvec dx)
1082 dvec dx_start, trial;
1083 double d2min, d2trial;
1086 dvec_sub(x1, x2, dx);
1088 switch (pbc->pbcTypeDX)
1090 case epbcdxRECTANGULAR:
1092 for (i = 0; i < DIM; i++)
1096 while (dx[i] > pbc->hbox_diag[i])
1098 dx[i] -= pbc->fbox_diag[i];
1100 while (dx[i] <= pbc->mhbox_diag[i])
1102 dx[i] += pbc->fbox_diag[i];
1107 case epbcdxTRICLINIC:
1110 for (i = DIM - 1; i >= 0; i--)
1114 while (dx[i] > pbc->hbox_diag[i])
1116 for (j = i; j >= 0; j--)
1118 dx[j] -= pbc->box[i][j];
1121 while (dx[i] <= pbc->mhbox_diag[i])
1123 for (j = i; j >= 0; j--)
1125 dx[j] += pbc->box[i][j];
1128 d2min += dx[i] * dx[i];
1131 if (d2min > pbc->max_cutoff2)
1133 copy_dvec(dx, dx_start);
1134 /* Now try all possible shifts, when the distance is within max_cutoff
1135 * it must be the shortest possible distance.
1138 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
1140 for (j = 0; j < DIM; j++)
1142 trial[j] = dx_start[j] + pbc->tric_vec[i][j];
1145 for (j = 0; j < DIM; j++)
1149 d2trial += trial[j] * trial[j];
1152 if (d2trial < d2min)
1154 copy_dvec(trial, dx);
1161 case epbcdxSCREW_RECT:
1162 /* The shift definition requires x first */
1164 while (dx[XX] > pbc->hbox_diag[XX])
1166 dx[XX] -= pbc->fbox_diag[XX];
1169 while (dx[XX] <= pbc->mhbox_diag[XX])
1171 dx[XX] += pbc->fbox_diag[YY];
1176 /* Rotate around the x-axis in the middle of the box */
1177 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1178 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1180 /* Normal pbc for y and z */
1181 for (i = YY; i <= ZZ; i++)
1183 while (dx[i] > pbc->hbox_diag[i])
1185 dx[i] -= pbc->fbox_diag[i];
1187 while (dx[i] <= pbc->mhbox_diag[i])
1189 dx[i] += pbc->fbox_diag[i];
1194 case epbcdxUNSUPPORTED: break;
1195 default: gmx_fatal(FARGS, "Internal error in pbc_dx, set_pbc has not been called");
1199 void calc_shifts(const matrix box, rvec shift_vec[])
1201 for (int n = 0, m = -D_BOX_Z; m <= D_BOX_Z; m++)
1203 for (int l = -D_BOX_Y; l <= D_BOX_Y; l++)
1205 for (int k = -D_BOX_X; k <= D_BOX_X; k++, n++)
1207 for (int d = 0; d < DIM; d++)
1209 shift_vec[n][d] = k * box[XX][d] + l * box[YY][d] + m * box[ZZ][d];
1216 void calc_box_center(int ecenter, const matrix box, rvec box_center)
1220 clear_rvec(box_center);
1224 for (m = 0; (m < DIM); m++)
1226 for (d = 0; d < DIM; d++)
1228 box_center[d] += 0.5 * box[m][d];
1233 for (d = 0; d < DIM; d++)
1235 box_center[d] = 0.5 * box[d][d];
1238 case ecenterZERO: break;
1239 default: gmx_fatal(FARGS, "Unsupported value %d for ecenter", ecenter);
1243 void calc_triclinic_images(const matrix box, rvec img[])
1247 /* Calculate 3 adjacent images in the xy-plane */
1248 copy_rvec(box[0], img[0]);
1249 copy_rvec(box[1], img[1]);
1252 svmul(-1, img[1], img[1]);
1254 rvec_sub(img[1], img[0], img[2]);
1256 /* Get the next 3 in the xy-plane as mirror images */
1257 for (i = 0; i < 3; i++)
1259 svmul(-1, img[i], img[3 + i]);
1262 /* Calculate the first 4 out of xy-plane images */
1263 copy_rvec(box[2], img[6]);
1266 svmul(-1, img[6], img[6]);
1268 for (i = 0; i < 3; i++)
1270 rvec_add(img[6], img[i + 1], img[7 + i]);
1273 /* Mirror the last 4 from the previous in opposite rotation */
1274 for (i = 0; i < 4; i++)
1276 svmul(-1, img[6 + (2 + i) % 4], img[10 + i]);
1280 void calc_compact_unitcell_vertices(int ecenter, const matrix box, rvec vert[])
1282 rvec img[NTRICIMG], box_center;
1283 int n, i, j, tmp[4], d;
1284 const real oneFourth = 0.25;
1286 calc_triclinic_images(box, img);
1289 for (i = 2; i <= 5; i += 3)
1300 tmp[2] = (i + 1) % 6;
1301 tmp[3] = tmp[1] + 4;
1302 for (j = 0; j < 4; j++)
1304 for (d = 0; d < DIM; d++)
1306 vert[n][d] = img[i][d] + img[tmp[j]][d] + img[tmp[(j + 1) % 4]][d];
1311 for (i = 7; i <= 13; i += 6)
1313 tmp[0] = (i - 7) / 2;
1314 tmp[1] = tmp[0] + 1;
1324 for (j = 0; j < 4; j++)
1326 for (d = 0; d < DIM; d++)
1328 vert[n][d] = img[i][d] + img[tmp[j]][d] + img[tmp[(j + 1) % 4]][d];
1333 for (i = 9; i <= 11; i += 2)
1343 tmp[1] = tmp[0] + 1;
1353 for (j = 0; j < 4; j++)
1355 for (d = 0; d < DIM; d++)
1357 vert[n][d] = img[i][d] + img[tmp[j]][d] + img[tmp[(j + 1) % 4]][d];
1363 calc_box_center(ecenter, box, box_center);
1364 for (i = 0; i < NCUCVERT; i++)
1366 for (d = 0; d < DIM; d++)
1368 vert[i][d] = vert[i][d] * oneFourth + box_center[d];
1373 int* compact_unitcell_edges()
1375 /* this is an index in vert[] (see calc_box_vertices) */
1376 /*static int edge[NCUCEDGE*2];*/
1378 static const int hexcon[24] = { 0, 9, 1, 19, 2, 15, 3, 21, 4, 17, 5, 11,
1379 6, 23, 7, 13, 8, 20, 10, 18, 12, 16, 14, 22 };
1382 snew(edge, NCUCEDGE * 2);
1385 for (i = 0; i < 6; i++)
1387 for (j = 0; j < 4; j++)
1389 edge[e++] = 4 * i + j;
1390 edge[e++] = 4 * i + (j + 1) % 4;
1393 for (i = 0; i < 12 * 2; i++)
1395 edge[e++] = hexcon[i];
1401 void put_atoms_in_box(PbcType pbcType, const matrix box, gmx::ArrayRef<gmx::RVec> x)
1405 if (pbcType == PbcType::Screw)
1407 gmx_fatal(FARGS, "Sorry, %s pbc is not yet supported", c_pbcTypeNames[pbcType].c_str());
1410 if (pbcType == PbcType::XY)
1421 for (gmx::index i = 0; (i < x.ssize()); ++i)
1423 for (m = npbcdim - 1; m >= 0; m--)
1427 for (d = 0; d <= m; d++)
1429 x[i][d] += box[m][d];
1432 while (x[i][m] >= box[m][m])
1434 for (d = 0; d <= m; d++)
1436 x[i][d] -= box[m][d];
1444 for (gmx::index i = 0; (i < x.ssize()); ++i)
1446 for (d = 0; d < npbcdim; d++)
1450 x[i][d] += box[d][d];
1452 while (x[i][d] >= box[d][d])
1454 x[i][d] -= box[d][d];
1461 void put_atoms_in_box_omp(PbcType pbcType, const matrix box, gmx::ArrayRef<gmx::RVec> x, gmx_unused int nth)
1463 #pragma omp parallel for num_threads(nth) schedule(static)
1464 for (int t = 0; t < nth; t++)
1468 size_t natoms = x.size();
1469 size_t offset = (natoms * t) / nth;
1470 size_t len = (natoms * (t + 1)) / nth - offset;
1471 put_atoms_in_box(pbcType, box, x.subArray(offset, len));
1473 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1477 void put_atoms_in_triclinic_unitcell(int ecenter, const matrix box, gmx::ArrayRef<gmx::RVec> x)
1479 rvec box_center, shift_center;
1480 real shm01, shm02, shm12, shift;
1483 calc_box_center(ecenter, box, box_center);
1485 /* The product of matrix shm with a coordinate gives the shift vector
1486 which is required determine the periodic cell position */
1487 shm01 = box[1][0] / box[1][1];
1488 shm02 = (box[1][1] * box[2][0] - box[2][1] * box[1][0]) / (box[1][1] * box[2][2]);
1489 shm12 = box[2][1] / box[2][2];
1491 clear_rvec(shift_center);
1492 for (d = 0; d < DIM; d++)
1494 rvec_inc(shift_center, box[d]);
1496 svmul(0.5, shift_center, shift_center);
1497 rvec_sub(box_center, shift_center, shift_center);
1499 shift_center[0] = shm01 * shift_center[1] + shm02 * shift_center[2];
1500 shift_center[1] = shm12 * shift_center[2];
1501 shift_center[2] = 0;
1503 for (gmx::index i = 0; (i < x.ssize()); ++i)
1505 for (m = DIM - 1; m >= 0; m--)
1507 shift = shift_center[m];
1510 shift += shm01 * x[i][1] + shm02 * x[i][2];
1514 shift += shm12 * x[i][2];
1516 while (x[i][m] - shift < 0)
1518 for (d = 0; d <= m; d++)
1520 x[i][d] += box[m][d];
1523 while (x[i][m] - shift >= box[m][m])
1525 for (d = 0; d <= m; d++)
1527 x[i][d] -= box[m][d];
1534 void put_atoms_in_compact_unitcell(PbcType pbcType, int ecenter, const matrix box, gmx::ArrayRef<gmx::RVec> x)
1537 rvec box_center, dx;
1539 set_pbc(&pbc, pbcType, box);
1541 if (pbc.pbcTypeDX == epbcdxUNSUPPORTED)
1543 gmx_fatal(FARGS, "Can not put atoms in compact unitcell with unsupported PBC");
1546 calc_box_center(ecenter, box, box_center);
1547 for (gmx::index i = 0; (i < x.ssize()); ++i)
1549 pbc_dx(&pbc, x[i], box_center, dx);
1550 rvec_add(box_center, dx, x[i]);
1554 /*! \brief Make molecules whole by shifting positions
1556 * \param[in] fplog Log file
1557 * \param[in] pbcType The PBC type
1558 * \param[in] box The simulation box
1559 * \param[in] mtop System topology definition
1560 * \param[in,out] x The coordinates of the atoms
1561 * \param[in] bFirst Specifier for first-time PBC removal
1564 low_do_pbc_mtop(FILE* fplog, PbcType pbcType, const matrix box, const gmx_mtop_t* mtop, rvec x[], gmx_bool bFirst)
1568 if (bFirst && fplog)
1570 fprintf(fplog, "Removing pbc first time\n");
1574 for (const gmx_molblock_t& molb : mtop->molblock)
1576 const gmx_moltype_t& moltype = mtop->moltype[molb.type];
1577 if (moltype.atoms.nr == 1 || (!bFirst && moltype.atoms.nr == 1))
1579 /* Just one atom or charge group in the molecule, no PBC required */
1580 as += molb.nmol * moltype.atoms.nr;
1584 t_graph graph = mk_graph_moltype(moltype);
1586 for (mol = 0; mol < molb.nmol; mol++)
1588 mk_mshift(fplog, &graph, pbcType, box, x + as);
1590 shift_self(graph, box, x + as);
1591 /* The molecule is whole now.
1592 * We don't need the second mk_mshift call as in do_pbc_first,
1593 * since we no longer need this graph.
1596 as += moltype.atoms.nr;
1602 void do_pbc_first_mtop(FILE* fplog, PbcType pbcType, const matrix box, const gmx_mtop_t* mtop, rvec x[])
1604 low_do_pbc_mtop(fplog, pbcType, box, mtop, x, TRUE);
1607 void do_pbc_mtop(PbcType pbcType, const matrix box, const gmx_mtop_t* mtop, rvec x[])
1609 low_do_pbc_mtop(nullptr, pbcType, box, mtop, x, FALSE);