Add replacements for pbc enumerations
[gromacs.git] / src / gromacs / pbcutil / pbc.cpp
blob0397683415186f86a1bc23231bff768680814190
1 /*
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,2019, 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.
37 /*! \internal \file
38 * \brief
39 * Implements routines in pbc.h.
41 * Utility functions for handling periodic boundary conditions.
42 * Mainly used in analysis tools.
44 #include "gmxpre.h"
46 #include "pbc.h"
48 #include <cmath>
50 #include <algorithm>
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/math/vecdump.h"
57 #include "gromacs/pbcutil/ishift.h"
58 #include "gromacs/pbcutil/mshift.h"
59 #include "gromacs/topology/topology.h"
60 #include "gromacs/utility/exceptions.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/gmxassert.h"
63 #include "gromacs/utility/smalloc.h"
65 const char *epbc_names[epbcNR+1] =
67 "xyz", "no", "xy", "screw", nullptr
70 /* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
71 enum {
72 epbcdxRECTANGULAR = 1, epbcdxTRICLINIC,
73 epbcdx2D_RECT, epbcdx2D_TRIC,
74 epbcdx1D_RECT, epbcdx1D_TRIC,
75 epbcdxSCREW_RECT, epbcdxSCREW_TRIC,
76 epbcdxNOPBC, epbcdxUNSUPPORTED
79 //! Margin factor for error message
80 #define BOX_MARGIN 1.0010
81 //! Margin correction if the box is too skewed
82 #define BOX_MARGIN_CORRECT 1.0005
84 int ePBC2npbcdim(int ePBC)
86 int npbcdim = 0;
88 switch (ePBC)
90 case epbcXYZ: npbcdim = 3; break;
91 case epbcXY: npbcdim = 2; break;
92 case epbcSCREW: npbcdim = 3; break;
93 case epbcNONE: npbcdim = 0; break;
94 default: gmx_fatal(FARGS, "Unknown ePBC=%d in ePBC2npbcdim", ePBC);
97 return npbcdim;
100 void dump_pbc(FILE *fp, t_pbc *pbc)
102 rvec sum_box;
104 fprintf(fp, "ePBCDX = %d\n", pbc->ePBCDX);
105 pr_rvecs(fp, 0, "box", pbc->box, DIM);
106 pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
107 pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
108 pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
109 rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
110 pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
111 fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
112 fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
113 if (pbc->ntric_vec > 0)
115 pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE);
116 pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
120 const char *check_box(int ePBC, const matrix box)
122 const char *ptr;
124 if (ePBC == -1)
126 ePBC = guess_ePBC(box);
129 if (ePBC == epbcNONE)
131 return nullptr;
134 if ((box[XX][YY] != 0) || (box[XX][ZZ] != 0) || (box[YY][ZZ] != 0))
136 ptr = "Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.";
138 else if (ePBC == epbcSCREW && (box[YY][XX] != 0 || box[ZZ][XX] != 0))
140 ptr = "The unit cell can not have off-diagonal x-components with screw pbc";
142 else if (std::fabs(box[YY][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
143 (ePBC != epbcXY &&
144 (std::fabs(box[ZZ][XX]) > BOX_MARGIN*0.5*box[XX][XX] ||
145 std::fabs(box[ZZ][YY]) > BOX_MARGIN*0.5*box[YY][YY])))
147 ptr = "Triclinic box is too skewed.";
149 else
151 ptr = nullptr;
154 return ptr;
157 void matrix_convert(matrix box, const rvec vec, const rvec angleInDegrees)
159 rvec angle;
160 svmul(DEG2RAD, angleInDegrees, angle);
161 box[XX][XX] = vec[XX];
162 box[YY][XX] = vec[YY]*cos(angle[ZZ]);
163 box[YY][YY] = vec[YY]*sin(angle[ZZ]);
164 box[ZZ][XX] = vec[ZZ]*cos(angle[YY]);
165 box[ZZ][YY] = vec[ZZ]
166 *(cos(angle[XX])-cos(angle[YY])*cos(angle[ZZ]))/sin(angle[ZZ]);
167 box[ZZ][ZZ] = std::sqrt(gmx::square(vec[ZZ])
168 -box[ZZ][XX]*box[ZZ][XX]-box[ZZ][YY]*box[ZZ][YY]);
171 real max_cutoff2(int ePBC, const matrix box)
173 real min_hv2, min_ss;
174 const real oneFourth = 0.25;
176 /* Physical limitation of the cut-off
177 * by half the length of the shortest box vector.
179 min_hv2 = oneFourth * std::min(norm2(box[XX]), norm2(box[YY]));
180 if (ePBC != epbcXY)
182 min_hv2 = std::min(min_hv2, oneFourth * norm2(box[ZZ]));
185 /* Limitation to the smallest diagonal element due to optimizations:
186 * checking only linear combinations of single box-vectors (2 in x)
187 * in the grid search and pbc_dx is a lot faster
188 * than checking all possible combinations.
190 if (ePBC == epbcXY)
192 min_ss = std::min(box[XX][XX], box[YY][YY]);
194 else
196 min_ss = std::min(box[XX][XX], std::min(box[YY][YY] - std::fabs(box[ZZ][YY]), box[ZZ][ZZ]));
199 return std::min(min_hv2, min_ss*min_ss);
202 //! Set to true if warning has been printed
203 static gmx_bool bWarnedGuess = FALSE;
205 int guess_ePBC(const matrix box)
207 int ePBC;
209 if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] > 0)
211 ePBC = epbcXYZ;
213 else if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] == 0)
215 ePBC = epbcXY;
217 else if (box[XX][XX] == 0 && box[YY][YY] == 0 && box[ZZ][ZZ] == 0)
219 ePBC = epbcNONE;
221 else
223 if (!bWarnedGuess)
225 fprintf(stderr, "WARNING: Unsupported box diagonal %f %f %f, "
226 "will not use periodic boundary conditions\n\n",
227 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
228 bWarnedGuess = TRUE;
230 ePBC = epbcNONE;
233 if (debug)
235 fprintf(debug, "Guessed pbc = %s from the box matrix\n", epbc_names[ePBC]);
238 return ePBC;
241 //! Check if the box still obeys the restrictions, if not, correct it
242 static int correct_box_elem(FILE *fplog, int step, tensor box, int v, int d)
244 int shift, maxshift = 10;
246 shift = 0;
248 /* correct elem d of vector v with vector d */
249 while (box[v][d] > BOX_MARGIN_CORRECT*0.5*box[d][d])
251 if (fplog)
253 fprintf(fplog, "Step %d: correcting invalid box:\n", step);
254 pr_rvecs(fplog, 0, "old box", box, DIM);
256 rvec_dec(box[v], box[d]);
257 shift--;
258 if (fplog)
260 pr_rvecs(fplog, 0, "new box", box, DIM);
262 if (shift <= -maxshift)
264 gmx_fatal(FARGS,
265 "Box was shifted at least %d times. Please see log-file.",
266 maxshift);
269 while (box[v][d] < -BOX_MARGIN_CORRECT*0.5*box[d][d])
271 if (fplog)
273 fprintf(fplog, "Step %d: correcting invalid box:\n", step);
274 pr_rvecs(fplog, 0, "old box", box, DIM);
276 rvec_inc(box[v], box[d]);
277 shift++;
278 if (fplog)
280 pr_rvecs(fplog, 0, "new box", box, DIM);
282 if (shift >= maxshift)
284 gmx_fatal(FARGS,
285 "Box was shifted at least %d times. Please see log-file.",
286 maxshift);
290 return shift;
293 gmx_bool correct_box(FILE *fplog, int step, tensor box, t_graph *graph)
295 int zy, zx, yx, i;
296 gmx_bool bCorrected;
298 zy = correct_box_elem(fplog, step, box, ZZ, YY);
299 zx = correct_box_elem(fplog, step, box, ZZ, XX);
300 yx = correct_box_elem(fplog, step, box, YY, XX);
302 bCorrected = ((zy != 0) || (zx != 0) || (yx != 0));
304 if (bCorrected && graph)
306 /* correct the graph */
307 for (i = graph->at_start; i < graph->at_end; i++)
309 graph->ishift[i][YY] -= graph->ishift[i][ZZ]*zy;
310 graph->ishift[i][XX] -= graph->ishift[i][ZZ]*zx;
311 graph->ishift[i][XX] -= graph->ishift[i][YY]*yx;
315 return bCorrected;
318 //! Do the real arithmetic for filling the pbc struct
319 static void low_set_pbc(t_pbc *pbc, int ePBC,
320 const ivec dd_pbc, const matrix box)
322 int order[3] = { 0, -1, 1 };
323 ivec bPBC;
324 const char *ptr;
326 pbc->ePBC = ePBC;
327 pbc->ndim_ePBC = ePBC2npbcdim(ePBC);
329 if (pbc->ePBC == epbcNONE)
331 pbc->ePBCDX = epbcdxNOPBC;
333 return;
336 copy_mat(box, pbc->box);
337 pbc->max_cutoff2 = 0;
338 pbc->dim = -1;
339 pbc->ntric_vec = 0;
341 for (int i = 0; (i < DIM); i++)
343 pbc->fbox_diag[i] = box[i][i];
344 pbc->hbox_diag[i] = pbc->fbox_diag[i]*0.5;
345 pbc->mhbox_diag[i] = -pbc->hbox_diag[i];
348 ptr = check_box(ePBC, box);
349 if (ptr)
351 fprintf(stderr, "Warning: %s\n", ptr);
352 pr_rvecs(stderr, 0, " Box", box, DIM);
353 fprintf(stderr, " Can not fix pbc.\n\n");
354 pbc->ePBCDX = epbcdxUNSUPPORTED;
356 else
358 if (ePBC == epbcSCREW && nullptr != dd_pbc)
360 /* This combinated should never appear here */
361 gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
364 int npbcdim = 0;
365 for (int i = 0; i < DIM; i++)
367 if ((dd_pbc && dd_pbc[i] == 0) || (ePBC == epbcXY && i == ZZ))
369 bPBC[i] = 0;
371 else
373 bPBC[i] = 1;
374 npbcdim++;
377 switch (npbcdim)
379 case 1:
380 /* 1D pbc is not an mdp option and it is therefore only used
381 * with single shifts.
383 pbc->ePBCDX = epbcdx1D_RECT;
384 for (int i = 0; i < DIM; i++)
386 if (bPBC[i])
388 pbc->dim = i;
391 GMX_ASSERT(pbc->dim < DIM, "Dimension for PBC incorrect");
392 for (int i = 0; i < pbc->dim; i++)
394 if (pbc->box[pbc->dim][i] != 0)
396 pbc->ePBCDX = epbcdx1D_TRIC;
399 break;
400 case 2:
401 pbc->ePBCDX = epbcdx2D_RECT;
402 for (int i = 0; i < DIM; i++)
404 if (!bPBC[i])
406 pbc->dim = i;
409 for (int i = 0; i < DIM; i++)
411 if (bPBC[i])
413 for (int j = 0; j < i; j++)
415 if (pbc->box[i][j] != 0)
417 pbc->ePBCDX = epbcdx2D_TRIC;
422 break;
423 case 3:
424 if (ePBC != epbcSCREW)
426 if (TRICLINIC(box))
428 pbc->ePBCDX = epbcdxTRICLINIC;
430 else
432 pbc->ePBCDX = epbcdxRECTANGULAR;
435 else
437 pbc->ePBCDX = (box[ZZ][YY] == 0 ? epbcdxSCREW_RECT : epbcdxSCREW_TRIC);
438 if (pbc->ePBCDX == epbcdxSCREW_TRIC)
440 fprintf(stderr,
441 "Screw pbc is not yet implemented for triclinic boxes.\n"
442 "Can not fix pbc.\n");
443 pbc->ePBCDX = epbcdxUNSUPPORTED;
446 break;
447 default:
448 gmx_fatal(FARGS, "Incorrect number of pbc dimensions with DD: %d",
449 npbcdim);
451 pbc->max_cutoff2 = max_cutoff2(ePBC, box);
453 if (pbc->ePBCDX == epbcdxTRICLINIC ||
454 pbc->ePBCDX == epbcdx2D_TRIC ||
455 pbc->ePBCDX == epbcdxSCREW_TRIC)
457 if (debug)
459 pr_rvecs(debug, 0, "Box", box, DIM);
460 fprintf(debug, "max cutoff %.3f\n", sqrt(pbc->max_cutoff2));
462 /* We will only need single shifts here */
463 for (int kk = 0; kk < 3; kk++)
465 int k = order[kk];
466 if (!bPBC[ZZ] && k != 0)
468 continue;
470 for (int jj = 0; jj < 3; jj++)
472 int j = order[jj];
473 if (!bPBC[YY] && j != 0)
475 continue;
477 for (int ii = 0; ii < 3; ii++)
479 int i = order[ii];
480 if (!bPBC[XX] && i != 0)
482 continue;
484 /* A shift is only useful when it is trilinic */
485 if (j != 0 || k != 0)
487 rvec trial;
488 rvec pos;
489 real d2old = 0;
490 real d2new = 0;
492 for (int d = 0; d < DIM; d++)
494 trial[d] = i*box[XX][d] + j*box[YY][d] + k*box[ZZ][d];
495 /* Choose the vector within the brick around 0,0,0 that
496 * will become the shortest due to shift try.
498 if (d == pbc->dim)
500 trial[d] = 0;
501 pos[d] = 0;
503 else
505 if (trial[d] < 0)
507 pos[d] = std::min( pbc->hbox_diag[d], -trial[d]);
509 else
511 pos[d] = std::max(-pbc->hbox_diag[d], -trial[d]);
514 d2old += gmx::square(pos[d]);
515 d2new += gmx::square(pos[d] + trial[d]);
517 if (BOX_MARGIN*d2new < d2old)
519 /* Check if shifts with one box vector less do better */
520 gmx_bool bUse = TRUE;
521 for (int dd = 0; dd < DIM; dd++)
523 int shift = (dd == 0 ? i : (dd == 1 ? j : k));
524 if (shift)
526 real d2new_c = 0;
527 for (int d = 0; d < DIM; d++)
529 d2new_c += gmx::square(pos[d] + trial[d] - shift*box[dd][d]);
531 if (d2new_c <= BOX_MARGIN*d2new)
533 bUse = FALSE;
537 if (bUse)
539 /* Accept this shift vector. */
540 if (pbc->ntric_vec >= MAX_NTRICVEC)
542 fprintf(stderr, "\nWARNING: Found more than %d triclinic correction vectors, ignoring some.\n"
543 " There is probably something wrong with your box.\n", MAX_NTRICVEC);
544 pr_rvecs(stderr, 0, " Box", box, DIM);
546 else
548 copy_rvec(trial, pbc->tric_vec[pbc->ntric_vec]);
549 pbc->tric_shift[pbc->ntric_vec][XX] = i;
550 pbc->tric_shift[pbc->ntric_vec][YY] = j;
551 pbc->tric_shift[pbc->ntric_vec][ZZ] = k;
552 pbc->ntric_vec++;
554 if (debug)
556 fprintf(debug, " tricvec %2d = %2d %2d %2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
557 pbc->ntric_vec, i, j, k,
558 sqrt(d2old), sqrt(d2new),
559 trial[XX], trial[YY], trial[ZZ],
560 pos[XX], pos[YY], pos[ZZ]);
573 void set_pbc(t_pbc *pbc, int ePBC, const matrix box)
575 if (ePBC == -1)
577 ePBC = guess_ePBC(box);
580 low_set_pbc(pbc, ePBC, nullptr, box);
583 t_pbc *set_pbc_dd(t_pbc *pbc, int ePBC,
584 const ivec domdecCells,
585 gmx_bool bSingleDir, const matrix box)
587 if (ePBC == epbcNONE)
589 pbc->ePBC = ePBC;
591 return nullptr;
594 if (nullptr == domdecCells)
596 low_set_pbc(pbc, ePBC, nullptr, box);
598 else
600 if (ePBC == epbcSCREW && domdecCells[XX] > 1)
602 /* The rotation has been taken care of during coordinate communication */
603 ePBC = epbcXYZ;
606 ivec usePBC;
607 int npbcdim = 0;
608 for (int i = 0; i < DIM; i++)
610 usePBC[i] = 0;
611 if (domdecCells[i] <= (bSingleDir ? 1 : 2) &&
612 !(ePBC == epbcXY && i == ZZ))
614 usePBC[i] = 1;
615 npbcdim++;
619 if (npbcdim > 0)
621 low_set_pbc(pbc, ePBC, usePBC, box);
623 else
625 pbc->ePBC = epbcNONE;
629 return (pbc->ePBC != epbcNONE ? pbc : nullptr);
632 void pbc_dx(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
634 int i, j;
635 rvec dx_start, trial;
636 real d2min, d2trial;
637 gmx_bool bRot;
639 rvec_sub(x1, x2, dx);
641 switch (pbc->ePBCDX)
643 case epbcdxRECTANGULAR:
644 for (i = 0; i < DIM; i++)
646 while (dx[i] > pbc->hbox_diag[i])
648 dx[i] -= pbc->fbox_diag[i];
650 while (dx[i] <= pbc->mhbox_diag[i])
652 dx[i] += pbc->fbox_diag[i];
655 break;
656 case epbcdxTRICLINIC:
657 for (i = DIM-1; i >= 0; i--)
659 while (dx[i] > pbc->hbox_diag[i])
661 for (j = i; j >= 0; j--)
663 dx[j] -= pbc->box[i][j];
666 while (dx[i] <= pbc->mhbox_diag[i])
668 for (j = i; j >= 0; j--)
670 dx[j] += pbc->box[i][j];
674 /* dx is the distance in a rectangular box */
675 d2min = norm2(dx);
676 if (d2min > pbc->max_cutoff2)
678 copy_rvec(dx, dx_start);
679 d2min = norm2(dx);
680 /* Now try all possible shifts, when the distance is within max_cutoff
681 * it must be the shortest possible distance.
683 i = 0;
684 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
686 rvec_add(dx_start, pbc->tric_vec[i], trial);
687 d2trial = norm2(trial);
688 if (d2trial < d2min)
690 copy_rvec(trial, dx);
691 d2min = d2trial;
693 i++;
696 break;
697 case epbcdx2D_RECT:
698 for (i = 0; i < DIM; i++)
700 if (i != pbc->dim)
702 while (dx[i] > pbc->hbox_diag[i])
704 dx[i] -= pbc->fbox_diag[i];
706 while (dx[i] <= pbc->mhbox_diag[i])
708 dx[i] += pbc->fbox_diag[i];
712 break;
713 case epbcdx2D_TRIC:
714 d2min = 0;
715 for (i = DIM-1; i >= 0; i--)
717 if (i != pbc->dim)
719 while (dx[i] > pbc->hbox_diag[i])
721 for (j = i; j >= 0; j--)
723 dx[j] -= pbc->box[i][j];
726 while (dx[i] <= pbc->mhbox_diag[i])
728 for (j = i; j >= 0; j--)
730 dx[j] += pbc->box[i][j];
733 d2min += dx[i]*dx[i];
736 if (d2min > pbc->max_cutoff2)
738 copy_rvec(dx, dx_start);
739 d2min = norm2(dx);
740 /* Now try all possible shifts, when the distance is within max_cutoff
741 * it must be the shortest possible distance.
743 i = 0;
744 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
746 rvec_add(dx_start, pbc->tric_vec[i], trial);
747 d2trial = 0;
748 for (j = 0; j < DIM; j++)
750 if (j != pbc->dim)
752 d2trial += trial[j]*trial[j];
755 if (d2trial < d2min)
757 copy_rvec(trial, dx);
758 d2min = d2trial;
760 i++;
763 break;
764 case epbcdxSCREW_RECT:
765 /* The shift definition requires x first */
766 bRot = FALSE;
767 while (dx[XX] > pbc->hbox_diag[XX])
769 dx[XX] -= pbc->fbox_diag[XX];
770 bRot = !bRot;
772 while (dx[XX] <= pbc->mhbox_diag[XX])
774 dx[XX] += pbc->fbox_diag[YY];
775 bRot = !bRot;
777 if (bRot)
779 /* Rotate around the x-axis in the middle of the box */
780 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
781 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
783 /* Normal pbc for y and z */
784 for (i = YY; i <= ZZ; i++)
786 while (dx[i] > pbc->hbox_diag[i])
788 dx[i] -= pbc->fbox_diag[i];
790 while (dx[i] <= pbc->mhbox_diag[i])
792 dx[i] += pbc->fbox_diag[i];
795 break;
796 case epbcdxNOPBC:
797 case epbcdxUNSUPPORTED:
798 break;
799 default:
800 gmx_fatal(FARGS, "Internal error in pbc_dx, set_pbc has not been called");
804 int pbc_dx_aiuc(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
806 int i, j, is;
807 rvec dx_start, trial;
808 real d2min, d2trial;
809 ivec ishift, ishift_start;
811 rvec_sub(x1, x2, dx);
812 clear_ivec(ishift);
814 switch (pbc->ePBCDX)
816 case epbcdxRECTANGULAR:
817 for (i = 0; i < DIM; i++)
819 if (dx[i] > pbc->hbox_diag[i])
821 dx[i] -= pbc->fbox_diag[i];
822 ishift[i]--;
824 else if (dx[i] <= pbc->mhbox_diag[i])
826 dx[i] += pbc->fbox_diag[i];
827 ishift[i]++;
830 break;
831 case epbcdxTRICLINIC:
832 /* For triclinic boxes the performance difference between
833 * if/else and two while loops is negligible.
834 * However, the while version can cause extreme delays
835 * before a simulation crashes due to large forces which
836 * can cause unlimited displacements.
837 * Also allowing multiple shifts would index fshift beyond bounds.
839 for (i = DIM-1; i >= 1; i--)
841 if (dx[i] > pbc->hbox_diag[i])
843 for (j = i; j >= 0; j--)
845 dx[j] -= pbc->box[i][j];
847 ishift[i]--;
849 else if (dx[i] <= pbc->mhbox_diag[i])
851 for (j = i; j >= 0; j--)
853 dx[j] += pbc->box[i][j];
855 ishift[i]++;
858 /* Allow 2 shifts in x */
859 if (dx[XX] > pbc->hbox_diag[XX])
861 dx[XX] -= pbc->fbox_diag[XX];
862 ishift[XX]--;
863 if (dx[XX] > pbc->hbox_diag[XX])
865 dx[XX] -= pbc->fbox_diag[XX];
866 ishift[XX]--;
869 else if (dx[XX] <= pbc->mhbox_diag[XX])
871 dx[XX] += pbc->fbox_diag[XX];
872 ishift[XX]++;
873 if (dx[XX] <= pbc->mhbox_diag[XX])
875 dx[XX] += pbc->fbox_diag[XX];
876 ishift[XX]++;
879 /* dx is the distance in a rectangular box */
880 d2min = norm2(dx);
881 if (d2min > pbc->max_cutoff2)
883 copy_rvec(dx, dx_start);
884 copy_ivec(ishift, ishift_start);
885 d2min = norm2(dx);
886 /* Now try all possible shifts, when the distance is within max_cutoff
887 * it must be the shortest possible distance.
889 i = 0;
890 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
892 rvec_add(dx_start, pbc->tric_vec[i], trial);
893 d2trial = norm2(trial);
894 if (d2trial < d2min)
896 copy_rvec(trial, dx);
897 ivec_add(ishift_start, pbc->tric_shift[i], ishift);
898 d2min = d2trial;
900 i++;
903 break;
904 case epbcdx2D_RECT:
905 for (i = 0; i < DIM; i++)
907 if (i != pbc->dim)
909 if (dx[i] > pbc->hbox_diag[i])
911 dx[i] -= pbc->fbox_diag[i];
912 ishift[i]--;
914 else if (dx[i] <= pbc->mhbox_diag[i])
916 dx[i] += pbc->fbox_diag[i];
917 ishift[i]++;
921 break;
922 case epbcdx2D_TRIC:
923 d2min = 0;
924 for (i = DIM-1; i >= 1; i--)
926 if (i != pbc->dim)
928 if (dx[i] > pbc->hbox_diag[i])
930 for (j = i; j >= 0; j--)
932 dx[j] -= pbc->box[i][j];
934 ishift[i]--;
936 else if (dx[i] <= pbc->mhbox_diag[i])
938 for (j = i; j >= 0; j--)
940 dx[j] += pbc->box[i][j];
942 ishift[i]++;
944 d2min += dx[i]*dx[i];
947 if (pbc->dim != XX)
949 /* Allow 2 shifts in x */
950 if (dx[XX] > pbc->hbox_diag[XX])
952 dx[XX] -= pbc->fbox_diag[XX];
953 ishift[XX]--;
954 if (dx[XX] > pbc->hbox_diag[XX])
956 dx[XX] -= pbc->fbox_diag[XX];
957 ishift[XX]--;
960 else if (dx[XX] <= pbc->mhbox_diag[XX])
962 dx[XX] += pbc->fbox_diag[XX];
963 ishift[XX]++;
964 if (dx[XX] <= pbc->mhbox_diag[XX])
966 dx[XX] += pbc->fbox_diag[XX];
967 ishift[XX]++;
970 d2min += dx[XX]*dx[XX];
972 if (d2min > pbc->max_cutoff2)
974 copy_rvec(dx, dx_start);
975 copy_ivec(ishift, ishift_start);
976 /* Now try all possible shifts, when the distance is within max_cutoff
977 * it must be the shortest possible distance.
979 i = 0;
980 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
982 rvec_add(dx_start, pbc->tric_vec[i], trial);
983 d2trial = 0;
984 for (j = 0; j < DIM; j++)
986 if (j != pbc->dim)
988 d2trial += trial[j]*trial[j];
991 if (d2trial < d2min)
993 copy_rvec(trial, dx);
994 ivec_add(ishift_start, pbc->tric_shift[i], ishift);
995 d2min = d2trial;
997 i++;
1000 break;
1001 case epbcdx1D_RECT:
1002 i = pbc->dim;
1003 if (dx[i] > pbc->hbox_diag[i])
1005 dx[i] -= pbc->fbox_diag[i];
1006 ishift[i]--;
1008 else if (dx[i] <= pbc->mhbox_diag[i])
1010 dx[i] += pbc->fbox_diag[i];
1011 ishift[i]++;
1013 break;
1014 case epbcdx1D_TRIC:
1015 i = pbc->dim;
1016 if (dx[i] > pbc->hbox_diag[i])
1018 rvec_dec(dx, pbc->box[i]);
1019 ishift[i]--;
1021 else if (dx[i] <= pbc->mhbox_diag[i])
1023 rvec_inc(dx, pbc->box[i]);
1024 ishift[i]++;
1026 break;
1027 case epbcdxSCREW_RECT:
1028 /* The shift definition requires x first */
1029 if (dx[XX] > pbc->hbox_diag[XX])
1031 dx[XX] -= pbc->fbox_diag[XX];
1032 ishift[XX]--;
1034 else if (dx[XX] <= pbc->mhbox_diag[XX])
1036 dx[XX] += pbc->fbox_diag[XX];
1037 ishift[XX]++;
1039 if (ishift[XX] == 1 || ishift[XX] == -1)
1041 /* Rotate around the x-axis in the middle of the box */
1042 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1043 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1045 /* Normal pbc for y and z */
1046 for (i = YY; i <= ZZ; i++)
1048 if (dx[i] > pbc->hbox_diag[i])
1050 dx[i] -= pbc->fbox_diag[i];
1051 ishift[i]--;
1053 else if (dx[i] <= pbc->mhbox_diag[i])
1055 dx[i] += pbc->fbox_diag[i];
1056 ishift[i]++;
1059 break;
1060 case epbcdxNOPBC:
1061 case epbcdxUNSUPPORTED:
1062 break;
1063 default:
1064 gmx_fatal(FARGS, "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
1067 is = IVEC2IS(ishift);
1068 if (debug)
1070 range_check_mesg(is, 0, SHIFTS, "PBC shift vector index range check.");
1073 return is;
1076 //! Compute distance vector in double precision
1077 void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx)
1079 int i, j;
1080 dvec dx_start, trial;
1081 double d2min, d2trial;
1082 gmx_bool bRot;
1084 dvec_sub(x1, x2, dx);
1086 switch (pbc->ePBCDX)
1088 case epbcdxRECTANGULAR:
1089 case epbcdx2D_RECT:
1090 for (i = 0; i < DIM; i++)
1092 if (i != pbc->dim)
1094 while (dx[i] > pbc->hbox_diag[i])
1096 dx[i] -= pbc->fbox_diag[i];
1098 while (dx[i] <= pbc->mhbox_diag[i])
1100 dx[i] += pbc->fbox_diag[i];
1104 break;
1105 case epbcdxTRICLINIC:
1106 case epbcdx2D_TRIC:
1107 d2min = 0;
1108 for (i = DIM-1; i >= 0; i--)
1110 if (i != pbc->dim)
1112 while (dx[i] > pbc->hbox_diag[i])
1114 for (j = i; j >= 0; j--)
1116 dx[j] -= pbc->box[i][j];
1119 while (dx[i] <= pbc->mhbox_diag[i])
1121 for (j = i; j >= 0; j--)
1123 dx[j] += pbc->box[i][j];
1126 d2min += dx[i]*dx[i];
1129 if (d2min > pbc->max_cutoff2)
1131 copy_dvec(dx, dx_start);
1132 /* Now try all possible shifts, when the distance is within max_cutoff
1133 * it must be the shortest possible distance.
1135 i = 0;
1136 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
1138 for (j = 0; j < DIM; j++)
1140 trial[j] = dx_start[j] + pbc->tric_vec[i][j];
1142 d2trial = 0;
1143 for (j = 0; j < DIM; j++)
1145 if (j != pbc->dim)
1147 d2trial += trial[j]*trial[j];
1150 if (d2trial < d2min)
1152 copy_dvec(trial, dx);
1153 d2min = d2trial;
1155 i++;
1158 break;
1159 case epbcdxSCREW_RECT:
1160 /* The shift definition requires x first */
1161 bRot = FALSE;
1162 while (dx[XX] > pbc->hbox_diag[XX])
1164 dx[XX] -= pbc->fbox_diag[XX];
1165 bRot = !bRot;
1167 while (dx[XX] <= pbc->mhbox_diag[XX])
1169 dx[XX] += pbc->fbox_diag[YY];
1170 bRot = !bRot;
1172 if (bRot)
1174 /* Rotate around the x-axis in the middle of the box */
1175 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1176 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1178 /* Normal pbc for y and z */
1179 for (i = YY; i <= ZZ; i++)
1181 while (dx[i] > pbc->hbox_diag[i])
1183 dx[i] -= pbc->fbox_diag[i];
1185 while (dx[i] <= pbc->mhbox_diag[i])
1187 dx[i] += pbc->fbox_diag[i];
1190 break;
1191 case epbcdxNOPBC:
1192 case epbcdxUNSUPPORTED:
1193 break;
1194 default:
1195 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 int k, l, m, d, n, test;
1203 n = 0;
1204 for (m = -D_BOX_Z; m <= D_BOX_Z; m++)
1206 for (l = -D_BOX_Y; l <= D_BOX_Y; l++)
1208 for (k = -D_BOX_X; k <= D_BOX_X; k++, n++)
1210 test = XYZ2IS(k, l, m);
1211 if (n != test)
1213 gmx_incons("inconsistent shift numbering");
1215 for (d = 0; d < DIM; d++)
1217 shift_vec[n][d] = k*box[XX][d] + l*box[YY][d] + m*box[ZZ][d];
1224 void calc_box_center(int ecenter, const matrix box, rvec box_center)
1226 int d, m;
1228 clear_rvec(box_center);
1229 switch (ecenter)
1231 case ecenterTRIC:
1232 for (m = 0; (m < DIM); m++)
1234 for (d = 0; d < DIM; d++)
1236 box_center[d] += 0.5*box[m][d];
1239 break;
1240 case ecenterRECT:
1241 for (d = 0; d < DIM; d++)
1243 box_center[d] = 0.5*box[d][d];
1245 break;
1246 case ecenterZERO:
1247 break;
1248 default:
1249 gmx_fatal(FARGS, "Unsupported value %d for ecenter", ecenter);
1253 void calc_triclinic_images(const matrix box, rvec img[])
1255 int i;
1257 /* Calculate 3 adjacent images in the xy-plane */
1258 copy_rvec(box[0], img[0]);
1259 copy_rvec(box[1], img[1]);
1260 if (img[1][XX] < 0)
1262 svmul(-1, img[1], img[1]);
1264 rvec_sub(img[1], img[0], img[2]);
1266 /* Get the next 3 in the xy-plane as mirror images */
1267 for (i = 0; i < 3; i++)
1269 svmul(-1, img[i], img[3+i]);
1272 /* Calculate the first 4 out of xy-plane images */
1273 copy_rvec(box[2], img[6]);
1274 if (img[6][XX] < 0)
1276 svmul(-1, img[6], img[6]);
1278 for (i = 0; i < 3; i++)
1280 rvec_add(img[6], img[i+1], img[7+i]);
1283 /* Mirror the last 4 from the previous in opposite rotation */
1284 for (i = 0; i < 4; i++)
1286 svmul(-1, img[6 + (2+i) % 4], img[10+i]);
1290 void calc_compact_unitcell_vertices(int ecenter, const matrix box, rvec vert[])
1292 rvec img[NTRICIMG], box_center;
1293 int n, i, j, tmp[4], d;
1294 const real oneFourth = 0.25;
1296 calc_triclinic_images(box, img);
1298 n = 0;
1299 for (i = 2; i <= 5; i += 3)
1301 tmp[0] = i-1;
1302 if (i == 2)
1304 tmp[1] = 8;
1306 else
1308 tmp[1] = 6;
1310 tmp[2] = (i+1) % 6;
1311 tmp[3] = tmp[1]+4;
1312 for (j = 0; j < 4; j++)
1314 for (d = 0; d < DIM; d++)
1316 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1318 n++;
1321 for (i = 7; i <= 13; i += 6)
1323 tmp[0] = (i-7)/2;
1324 tmp[1] = tmp[0]+1;
1325 if (i == 7)
1327 tmp[2] = 8;
1329 else
1331 tmp[2] = 10;
1333 tmp[3] = i-1;
1334 for (j = 0; j < 4; j++)
1336 for (d = 0; d < DIM; d++)
1338 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1340 n++;
1343 for (i = 9; i <= 11; i += 2)
1345 if (i == 9)
1347 tmp[0] = 3;
1349 else
1351 tmp[0] = 0;
1353 tmp[1] = tmp[0]+1;
1354 if (i == 9)
1356 tmp[2] = 6;
1358 else
1360 tmp[2] = 12;
1362 tmp[3] = i-1;
1363 for (j = 0; j < 4; j++)
1365 for (d = 0; d < DIM; d++)
1367 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1369 n++;
1373 calc_box_center(ecenter, box, box_center);
1374 for (i = 0; i < NCUCVERT; i++)
1376 for (d = 0; d < DIM; d++)
1378 vert[i][d] = vert[i][d]*oneFourth+box_center[d];
1383 int *compact_unitcell_edges()
1385 /* this is an index in vert[] (see calc_box_vertices) */
1386 /*static int edge[NCUCEDGE*2];*/
1387 int *edge;
1388 static const int hexcon[24] = {
1389 0, 9, 1, 19, 2, 15, 3, 21,
1390 4, 17, 5, 11, 6, 23, 7, 13,
1391 8, 20, 10, 18, 12, 16, 14, 22
1393 int e, i, j;
1395 snew(edge, NCUCEDGE*2);
1397 e = 0;
1398 for (i = 0; i < 6; i++)
1400 for (j = 0; j < 4; j++)
1402 edge[e++] = 4*i + j;
1403 edge[e++] = 4*i + (j+1) % 4;
1406 for (i = 0; i < 12*2; i++)
1408 edge[e++] = hexcon[i];
1411 return edge;
1414 void put_atoms_in_box(int ePBC, const matrix box, gmx::ArrayRef<gmx::RVec> x)
1416 int npbcdim, m, d;
1418 if (ePBC == epbcSCREW)
1420 gmx_fatal(FARGS, "Sorry, %s pbc is not yet supported", epbc_names[ePBC]);
1423 if (ePBC == epbcXY)
1425 npbcdim = 2;
1427 else
1429 npbcdim = 3;
1432 if (TRICLINIC(box))
1434 for (gmx::index i = 0; (i < x.ssize()); ++i)
1436 for (m = npbcdim-1; m >= 0; m--)
1438 while (x[i][m] < 0)
1440 for (d = 0; d <= m; d++)
1442 x[i][d] += box[m][d];
1445 while (x[i][m] >= box[m][m])
1447 for (d = 0; d <= m; d++)
1449 x[i][d] -= box[m][d];
1455 else
1457 for (gmx::index i = 0; (i < x.ssize()); ++i)
1459 for (d = 0; d < npbcdim; d++)
1461 while (x[i][d] < 0)
1463 x[i][d] += box[d][d];
1465 while (x[i][d] >= box[d][d])
1467 x[i][d] -= box[d][d];
1474 void put_atoms_in_box_omp(int ePBC, const matrix box, gmx::ArrayRef<gmx::RVec> x, gmx_unused int nth)
1476 int t;
1478 #pragma omp parallel for num_threads(nth) schedule(static)
1479 for (t = 0; t < nth; t++)
1483 size_t natoms = x.size();
1484 size_t offset = (natoms*t )/nth;
1485 size_t len = (natoms*(t + 1))/nth - offset;
1486 put_atoms_in_box(ePBC, box, x.subArray(offset, len));
1488 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1492 void put_atoms_in_triclinic_unitcell(int ecenter, const matrix box,
1493 gmx::ArrayRef<gmx::RVec> x)
1495 rvec box_center, shift_center;
1496 real shm01, shm02, shm12, shift;
1497 int m, d;
1499 calc_box_center(ecenter, box, box_center);
1501 /* The product of matrix shm with a coordinate gives the shift vector
1502 which is required determine the periodic cell position */
1503 shm01 = box[1][0]/box[1][1];
1504 shm02 = (box[1][1]*box[2][0] - box[2][1]*box[1][0])/(box[1][1]*box[2][2]);
1505 shm12 = box[2][1]/box[2][2];
1507 clear_rvec(shift_center);
1508 for (d = 0; d < DIM; d++)
1510 rvec_inc(shift_center, box[d]);
1512 svmul(0.5, shift_center, shift_center);
1513 rvec_sub(box_center, shift_center, shift_center);
1515 shift_center[0] = shm01*shift_center[1] + shm02*shift_center[2];
1516 shift_center[1] = shm12*shift_center[2];
1517 shift_center[2] = 0;
1519 for (gmx::index i = 0; (i < x.ssize()); ++i)
1521 for (m = DIM-1; m >= 0; m--)
1523 shift = shift_center[m];
1524 if (m == 0)
1526 shift += shm01*x[i][1] + shm02*x[i][2];
1528 else if (m == 1)
1530 shift += shm12*x[i][2];
1532 while (x[i][m]-shift < 0)
1534 for (d = 0; d <= m; d++)
1536 x[i][d] += box[m][d];
1539 while (x[i][m]-shift >= box[m][m])
1541 for (d = 0; d <= m; d++)
1543 x[i][d] -= box[m][d];
1550 void put_atoms_in_compact_unitcell(int ePBC, int ecenter, const matrix box,
1551 gmx::ArrayRef<gmx::RVec> x)
1553 t_pbc pbc;
1554 rvec box_center, dx;
1556 set_pbc(&pbc, ePBC, box);
1558 if (pbc.ePBCDX == epbcdxUNSUPPORTED)
1560 gmx_fatal(FARGS, "Can not put atoms in compact unitcell with unsupported PBC");
1563 calc_box_center(ecenter, box, box_center);
1564 for (gmx::index i = 0; (i < x.ssize()); ++i)
1566 pbc_dx(&pbc, x[i], box_center, dx);
1567 rvec_add(box_center, dx, x[i]);
1571 /*! \brief Make molecules whole by shifting positions
1573 * \param[in] fplog Log file
1574 * \param[in] ePBC The PBC type
1575 * \param[in] box The simulation box
1576 * \param[in] mtop System topology definition
1577 * \param[in,out] x The coordinates of the atoms
1578 * \param[in] bFirst Specifier for first-time PBC removal
1580 static void low_do_pbc_mtop(FILE *fplog, int ePBC, const matrix box,
1581 const gmx_mtop_t *mtop, rvec x[],
1582 gmx_bool bFirst)
1584 t_graph *graph;
1585 int as, mol;
1587 if (bFirst && fplog)
1589 fprintf(fplog, "Removing pbc first time\n");
1592 snew(graph, 1);
1593 as = 0;
1594 for (const gmx_molblock_t &molb : mtop->molblock)
1596 const gmx_moltype_t &moltype = mtop->moltype[molb.type];
1597 if (moltype.atoms.nr == 1 ||
1598 (!bFirst && moltype.cgs.nr == 1))
1600 /* Just one atom or charge group in the molecule, no PBC required */
1601 as += molb.nmol*moltype.atoms.nr;
1603 else
1605 mk_graph_moltype(moltype, graph);
1607 for (mol = 0; mol < molb.nmol; mol++)
1609 mk_mshift(fplog, graph, ePBC, box, x+as);
1611 shift_self(graph, box, x+as);
1612 /* The molecule is whole now.
1613 * We don't need the second mk_mshift call as in do_pbc_first,
1614 * since we no longer need this graph.
1617 as += moltype.atoms.nr;
1619 done_graph(graph);
1622 sfree(graph);
1625 void do_pbc_first_mtop(FILE *fplog, int ePBC, const matrix box,
1626 const gmx_mtop_t *mtop, rvec x[])
1628 low_do_pbc_mtop(fplog, ePBC, box, mtop, x, TRUE);
1631 void do_pbc_mtop(int ePBC, const matrix box,
1632 const gmx_mtop_t *mtop, rvec x[])
1634 low_do_pbc_mtop(nullptr, ePBC, box, mtop, x, FALSE);