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 char* epbc_names
[epbcNR
+ 1] = { "xyz", "no", "xy", "screw", nullptr };
68 /* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
71 epbcdxRECTANGULAR
= 1,
83 //! Margin factor for error message
84 #define BOX_MARGIN 1.0010
85 //! Margin correction if the box is too skewed
86 #define BOX_MARGIN_CORRECT 1.0005
88 int ePBC2npbcdim(int ePBC
)
94 case epbcXYZ
: npbcdim
= 3; break;
95 case epbcXY
: npbcdim
= 2; break;
96 case epbcSCREW
: npbcdim
= 3; break;
97 case epbcNONE
: npbcdim
= 0; break;
98 default: gmx_fatal(FARGS
, "Unknown ePBC=%d in ePBC2npbcdim", ePBC
);
104 void dump_pbc(FILE* fp
, t_pbc
* pbc
)
108 fprintf(fp
, "ePBCDX = %d\n", pbc
->ePBCDX
);
109 pr_rvecs(fp
, 0, "box", pbc
->box
, DIM
);
110 pr_rvecs(fp
, 0, "fbox_diag", &pbc
->fbox_diag
, 1);
111 pr_rvecs(fp
, 0, "hbox_diag", &pbc
->hbox_diag
, 1);
112 pr_rvecs(fp
, 0, "mhbox_diag", &pbc
->mhbox_diag
, 1);
113 rvec_add(pbc
->hbox_diag
, pbc
->mhbox_diag
, sum_box
);
114 pr_rvecs(fp
, 0, "sum of the above two", &sum_box
, 1);
115 fprintf(fp
, "max_cutoff2 = %g\n", pbc
->max_cutoff2
);
116 fprintf(fp
, "ntric_vec = %d\n", pbc
->ntric_vec
);
117 if (pbc
->ntric_vec
> 0)
119 pr_ivecs(fp
, 0, "tric_shift", pbc
->tric_shift
, pbc
->ntric_vec
, FALSE
);
120 pr_rvecs(fp
, 0, "tric_vec", pbc
->tric_vec
, pbc
->ntric_vec
);
124 const char* check_box(int ePBC
, const matrix box
)
130 ePBC
= guess_ePBC(box
);
133 if (ePBC
== epbcNONE
)
138 if ((box
[XX
][YY
] != 0) || (box
[XX
][ZZ
] != 0) || (box
[YY
][ZZ
] != 0))
140 ptr
= "Only triclinic boxes with the first vector parallel to the x-axis and the second "
141 "vector in the xy-plane are supported.";
143 else if (ePBC
== epbcSCREW
&& (box
[YY
][XX
] != 0 || box
[ZZ
][XX
] != 0))
145 ptr
= "The unit cell can not have off-diagonal x-components with screw pbc";
147 else if (std::fabs(box
[YY
][XX
]) > BOX_MARGIN
* 0.5 * box
[XX
][XX
]
149 && (std::fabs(box
[ZZ
][XX
]) > BOX_MARGIN
* 0.5 * box
[XX
][XX
]
150 || std::fabs(box
[ZZ
][YY
]) > BOX_MARGIN
* 0.5 * box
[YY
][YY
])))
152 ptr
= "Triclinic box is too skewed.";
162 void matrix_convert(matrix box
, const rvec vec
, const rvec angleInDegrees
)
165 svmul(DEG2RAD
, angleInDegrees
, angle
);
166 box
[XX
][XX
] = vec
[XX
];
167 box
[YY
][XX
] = vec
[YY
] * cos(angle
[ZZ
]);
168 box
[YY
][YY
] = vec
[YY
] * sin(angle
[ZZ
]);
169 box
[ZZ
][XX
] = vec
[ZZ
] * cos(angle
[YY
]);
170 box
[ZZ
][YY
] = vec
[ZZ
] * (cos(angle
[XX
]) - cos(angle
[YY
]) * cos(angle
[ZZ
])) / sin(angle
[ZZ
]);
172 std::sqrt(gmx::square(vec
[ZZ
]) - box
[ZZ
][XX
] * box
[ZZ
][XX
] - box
[ZZ
][YY
] * box
[ZZ
][YY
]);
175 real
max_cutoff2(int ePBC
, const matrix box
)
177 real min_hv2
, min_ss
;
178 const real oneFourth
= 0.25;
180 /* Physical limitation of the cut-off
181 * by half the length of the shortest box vector.
183 min_hv2
= oneFourth
* std::min(norm2(box
[XX
]), norm2(box
[YY
]));
186 min_hv2
= std::min(min_hv2
, oneFourth
* norm2(box
[ZZ
]));
189 /* Limitation to the smallest diagonal element due to optimizations:
190 * checking only linear combinations of single box-vectors (2 in x)
191 * in the grid search and pbc_dx is a lot faster
192 * than checking all possible combinations.
196 min_ss
= std::min(box
[XX
][XX
], box
[YY
][YY
]);
200 min_ss
= std::min(box
[XX
][XX
], std::min(box
[YY
][YY
] - std::fabs(box
[ZZ
][YY
]), box
[ZZ
][ZZ
]));
203 return std::min(min_hv2
, min_ss
* min_ss
);
206 //! Set to true if warning has been printed
207 static gmx_bool bWarnedGuess
= FALSE
;
209 int guess_ePBC(const matrix box
)
213 if (box
[XX
][XX
] > 0 && box
[YY
][YY
] > 0 && box
[ZZ
][ZZ
] > 0)
217 else if (box
[XX
][XX
] > 0 && box
[YY
][YY
] > 0 && box
[ZZ
][ZZ
] == 0)
221 else if (box
[XX
][XX
] == 0 && box
[YY
][YY
] == 0 && box
[ZZ
][ZZ
] == 0)
230 "WARNING: Unsupported box diagonal %f %f %f, "
231 "will not use periodic boundary conditions\n\n",
232 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
240 fprintf(debug
, "Guessed pbc = %s from the box matrix\n", epbc_names
[ePBC
]);
246 //! Check if the box still obeys the restrictions, if not, correct it
247 static int correct_box_elem(FILE* fplog
, int step
, tensor box
, int v
, int d
)
249 int shift
, maxshift
= 10;
253 /* correct elem d of vector v with vector d */
254 while (box
[v
][d
] > BOX_MARGIN_CORRECT
* 0.5 * box
[d
][d
])
258 fprintf(fplog
, "Step %d: correcting invalid box:\n", step
);
259 pr_rvecs(fplog
, 0, "old box", box
, DIM
);
261 rvec_dec(box
[v
], box
[d
]);
265 pr_rvecs(fplog
, 0, "new box", box
, DIM
);
267 if (shift
<= -maxshift
)
269 gmx_fatal(FARGS
, "Box was shifted at least %d times. Please see log-file.", maxshift
);
272 while (box
[v
][d
] < -BOX_MARGIN_CORRECT
* 0.5 * box
[d
][d
])
276 fprintf(fplog
, "Step %d: correcting invalid box:\n", step
);
277 pr_rvecs(fplog
, 0, "old box", box
, DIM
);
279 rvec_inc(box
[v
], box
[d
]);
283 pr_rvecs(fplog
, 0, "new box", box
, DIM
);
285 if (shift
>= maxshift
)
287 gmx_fatal(FARGS
, "Box was shifted at least %d times. Please see log-file.", maxshift
);
294 gmx_bool
correct_box(FILE* fplog
, int step
, tensor box
, t_graph
* graph
)
299 zy
= correct_box_elem(fplog
, step
, box
, ZZ
, YY
);
300 zx
= correct_box_elem(fplog
, step
, box
, ZZ
, XX
);
301 yx
= correct_box_elem(fplog
, step
, box
, YY
, XX
);
303 bCorrected
= ((zy
!= 0) || (zx
!= 0) || (yx
!= 0));
305 if (bCorrected
&& graph
)
307 /* correct the graph */
308 for (i
= graph
->at_start
; i
< graph
->at_end
; i
++)
310 graph
->ishift
[i
][YY
] -= graph
->ishift
[i
][ZZ
] * zy
;
311 graph
->ishift
[i
][XX
] -= graph
->ishift
[i
][ZZ
] * zx
;
312 graph
->ishift
[i
][XX
] -= graph
->ishift
[i
][YY
] * yx
;
319 //! Do the real arithmetic for filling the pbc struct
320 static void low_set_pbc(t_pbc
* pbc
, int ePBC
, const ivec dd_pbc
, const matrix box
)
322 int order
[3] = { 0, -1, 1 };
327 pbc
->ndim_ePBC
= ePBC2npbcdim(ePBC
);
329 if (pbc
->ePBC
== epbcNONE
)
331 pbc
->ePBCDX
= epbcdxNOPBC
;
336 copy_mat(box
, pbc
->box
);
337 pbc
->max_cutoff2
= 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
);
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
;
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");
365 for (int i
= 0; i
< DIM
; i
++)
367 if ((dd_pbc
&& dd_pbc
[i
] == 0) || (ePBC
== epbcXY
&& i
== ZZ
))
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
++)
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
;
401 pbc
->ePBCDX
= epbcdx2D_RECT
;
402 for (int i
= 0; i
< DIM
; i
++)
409 for (int i
= 0; i
< DIM
; i
++)
413 for (int j
= 0; j
< i
; j
++)
415 if (pbc
->box
[i
][j
] != 0)
417 pbc
->ePBCDX
= epbcdx2D_TRIC
;
424 if (ePBC
!= epbcSCREW
)
428 pbc
->ePBCDX
= epbcdxTRICLINIC
;
432 pbc
->ePBCDX
= epbcdxRECTANGULAR
;
437 pbc
->ePBCDX
= (box
[ZZ
][YY
] == 0 ? epbcdxSCREW_RECT
: epbcdxSCREW_TRIC
);
438 if (pbc
->ePBCDX
== epbcdxSCREW_TRIC
)
441 "Screw pbc is not yet implemented for triclinic boxes.\n"
442 "Can not fix pbc.\n");
443 pbc
->ePBCDX
= epbcdxUNSUPPORTED
;
447 default: gmx_fatal(FARGS
, "Incorrect number of pbc dimensions with DD: %d", npbcdim
);
449 pbc
->max_cutoff2
= max_cutoff2(ePBC
, box
);
451 if (pbc
->ePBCDX
== epbcdxTRICLINIC
|| pbc
->ePBCDX
== epbcdx2D_TRIC
|| pbc
->ePBCDX
== epbcdxSCREW_TRIC
)
455 pr_rvecs(debug
, 0, "Box", box
, DIM
);
456 fprintf(debug
, "max cutoff %.3f\n", sqrt(pbc
->max_cutoff2
));
458 /* We will only need single shifts here */
459 for (int kk
= 0; kk
< 3; kk
++)
462 if (!bPBC
[ZZ
] && k
!= 0)
466 for (int jj
= 0; jj
< 3; jj
++)
469 if (!bPBC
[YY
] && j
!= 0)
473 for (int ii
= 0; ii
< 3; ii
++)
476 if (!bPBC
[XX
] && i
!= 0)
480 /* A shift is only useful when it is trilinic */
481 if (j
!= 0 || k
!= 0)
488 for (int d
= 0; d
< DIM
; d
++)
490 trial
[d
] = i
* box
[XX
][d
] + j
* box
[YY
][d
] + k
* box
[ZZ
][d
];
491 /* Choose the vector within the brick around 0,0,0 that
492 * will become the shortest due to shift try.
503 pos
[d
] = std::min(pbc
->hbox_diag
[d
], -trial
[d
]);
507 pos
[d
] = std::max(-pbc
->hbox_diag
[d
], -trial
[d
]);
510 d2old
+= gmx::square(pos
[d
]);
511 d2new
+= gmx::square(pos
[d
] + trial
[d
]);
513 if (BOX_MARGIN
* d2new
< d2old
)
515 /* Check if shifts with one box vector less do better */
516 gmx_bool bUse
= TRUE
;
517 for (int dd
= 0; dd
< DIM
; dd
++)
519 int shift
= (dd
== 0 ? i
: (dd
== 1 ? j
: k
));
523 for (int d
= 0; d
< DIM
; d
++)
525 d2new_c
+= gmx::square(pos
[d
] + trial
[d
] - shift
* box
[dd
][d
]);
527 if (d2new_c
<= BOX_MARGIN
* d2new
)
535 /* Accept this shift vector. */
536 if (pbc
->ntric_vec
>= MAX_NTRICVEC
)
539 "\nWARNING: Found more than %d triclinic "
540 "correction vectors, ignoring some.\n"
541 " There is probably something wrong with your "
544 pr_rvecs(stderr
, 0, " Box", box
, DIM
);
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
;
557 " tricvec %2d = %2d %2d %2d %5.2f %5.2f "
558 "%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
559 pbc
->ntric_vec
, i
, j
, k
, sqrt(d2old
),
560 sqrt(d2new
), trial
[XX
], trial
[YY
], trial
[ZZ
],
561 pos
[XX
], pos
[YY
], pos
[ZZ
]);
574 void set_pbc(t_pbc
* pbc
, int ePBC
, const matrix box
)
578 ePBC
= guess_ePBC(box
);
581 low_set_pbc(pbc
, ePBC
, nullptr, box
);
584 t_pbc
* set_pbc_dd(t_pbc
* pbc
, int ePBC
, const ivec domdecCells
, gmx_bool bSingleDir
, const matrix box
)
586 if (ePBC
== epbcNONE
)
593 if (nullptr == domdecCells
)
595 low_set_pbc(pbc
, ePBC
, nullptr, box
);
599 if (ePBC
== epbcSCREW
&& domdecCells
[XX
] > 1)
601 /* The rotation has been taken care of during coordinate communication */
607 for (int i
= 0; i
< DIM
; i
++)
610 if (domdecCells
[i
] <= (bSingleDir
? 1 : 2) && !(ePBC
== epbcXY
&& i
== ZZ
))
619 low_set_pbc(pbc
, ePBC
, usePBC
, box
);
623 pbc
->ePBC
= epbcNONE
;
627 return (pbc
->ePBC
!= epbcNONE
? pbc
: nullptr);
630 void pbc_dx(const t_pbc
* pbc
, const rvec x1
, const rvec x2
, rvec dx
)
633 rvec dx_start
, trial
;
637 rvec_sub(x1
, x2
, dx
);
641 case epbcdxRECTANGULAR
:
642 for (i
= 0; i
< DIM
; i
++)
644 while (dx
[i
] > pbc
->hbox_diag
[i
])
646 dx
[i
] -= pbc
->fbox_diag
[i
];
648 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
650 dx
[i
] += pbc
->fbox_diag
[i
];
654 case epbcdxTRICLINIC
:
655 for (i
= DIM
- 1; i
>= 0; i
--)
657 while (dx
[i
] > pbc
->hbox_diag
[i
])
659 for (j
= i
; j
>= 0; j
--)
661 dx
[j
] -= pbc
->box
[i
][j
];
664 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
666 for (j
= i
; j
>= 0; j
--)
668 dx
[j
] += pbc
->box
[i
][j
];
672 /* dx is the distance in a rectangular box */
674 if (d2min
> pbc
->max_cutoff2
)
676 copy_rvec(dx
, dx_start
);
678 /* Now try all possible shifts, when the distance is within max_cutoff
679 * it must be the shortest possible distance.
682 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
684 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
685 d2trial
= norm2(trial
);
688 copy_rvec(trial
, dx
);
696 for (i
= 0; i
< DIM
; i
++)
700 while (dx
[i
] > pbc
->hbox_diag
[i
])
702 dx
[i
] -= pbc
->fbox_diag
[i
];
704 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
706 dx
[i
] += pbc
->fbox_diag
[i
];
713 for (i
= DIM
- 1; i
>= 0; i
--)
717 while (dx
[i
] > pbc
->hbox_diag
[i
])
719 for (j
= i
; j
>= 0; j
--)
721 dx
[j
] -= pbc
->box
[i
][j
];
724 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
726 for (j
= i
; j
>= 0; j
--)
728 dx
[j
] += pbc
->box
[i
][j
];
731 d2min
+= dx
[i
] * dx
[i
];
734 if (d2min
> pbc
->max_cutoff2
)
736 copy_rvec(dx
, dx_start
);
738 /* Now try all possible shifts, when the distance is within max_cutoff
739 * it must be the shortest possible distance.
742 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
744 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
746 for (j
= 0; j
< DIM
; j
++)
750 d2trial
+= trial
[j
] * trial
[j
];
755 copy_rvec(trial
, dx
);
762 case epbcdxSCREW_RECT
:
763 /* The shift definition requires x first */
765 while (dx
[XX
] > pbc
->hbox_diag
[XX
])
767 dx
[XX
] -= pbc
->fbox_diag
[XX
];
770 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
772 dx
[XX
] += pbc
->fbox_diag
[YY
];
777 /* Rotate around the x-axis in the middle of the box */
778 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
779 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
781 /* Normal pbc for y and z */
782 for (i
= YY
; i
<= ZZ
; i
++)
784 while (dx
[i
] > pbc
->hbox_diag
[i
])
786 dx
[i
] -= pbc
->fbox_diag
[i
];
788 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
790 dx
[i
] += pbc
->fbox_diag
[i
];
795 case epbcdxUNSUPPORTED
: break;
796 default: gmx_fatal(FARGS
, "Internal error in pbc_dx, set_pbc has not been called");
800 int pbc_dx_aiuc(const t_pbc
* pbc
, const rvec x1
, const rvec x2
, rvec dx
)
803 rvec dx_start
, trial
;
805 ivec ishift
, ishift_start
;
807 rvec_sub(x1
, x2
, dx
);
812 case epbcdxRECTANGULAR
:
813 for (i
= 0; i
< DIM
; i
++)
815 if (dx
[i
] > pbc
->hbox_diag
[i
])
817 dx
[i
] -= pbc
->fbox_diag
[i
];
820 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
822 dx
[i
] += pbc
->fbox_diag
[i
];
827 case epbcdxTRICLINIC
:
828 /* For triclinic boxes the performance difference between
829 * if/else and two while loops is negligible.
830 * However, the while version can cause extreme delays
831 * before a simulation crashes due to large forces which
832 * can cause unlimited displacements.
833 * Also allowing multiple shifts would index fshift beyond bounds.
835 for (i
= DIM
- 1; i
>= 1; i
--)
837 if (dx
[i
] > pbc
->hbox_diag
[i
])
839 for (j
= i
; j
>= 0; j
--)
841 dx
[j
] -= pbc
->box
[i
][j
];
845 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
847 for (j
= i
; j
>= 0; j
--)
849 dx
[j
] += pbc
->box
[i
][j
];
854 /* Allow 2 shifts in x */
855 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
857 dx
[XX
] -= pbc
->fbox_diag
[XX
];
859 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
861 dx
[XX
] -= pbc
->fbox_diag
[XX
];
865 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
867 dx
[XX
] += pbc
->fbox_diag
[XX
];
869 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
871 dx
[XX
] += pbc
->fbox_diag
[XX
];
875 /* dx is the distance in a rectangular box */
877 if (d2min
> pbc
->max_cutoff2
)
879 copy_rvec(dx
, dx_start
);
880 copy_ivec(ishift
, ishift_start
);
882 /* Now try all possible shifts, when the distance is within max_cutoff
883 * it must be the shortest possible distance.
886 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
888 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
889 d2trial
= norm2(trial
);
892 copy_rvec(trial
, dx
);
893 ivec_add(ishift_start
, pbc
->tric_shift
[i
], ishift
);
901 for (i
= 0; i
< DIM
; i
++)
905 if (dx
[i
] > pbc
->hbox_diag
[i
])
907 dx
[i
] -= pbc
->fbox_diag
[i
];
910 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
912 dx
[i
] += pbc
->fbox_diag
[i
];
920 for (i
= DIM
- 1; i
>= 1; i
--)
924 if (dx
[i
] > pbc
->hbox_diag
[i
])
926 for (j
= i
; j
>= 0; j
--)
928 dx
[j
] -= pbc
->box
[i
][j
];
932 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
934 for (j
= i
; j
>= 0; j
--)
936 dx
[j
] += pbc
->box
[i
][j
];
940 d2min
+= dx
[i
] * dx
[i
];
945 /* Allow 2 shifts in x */
946 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
948 dx
[XX
] -= pbc
->fbox_diag
[XX
];
950 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
952 dx
[XX
] -= pbc
->fbox_diag
[XX
];
956 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
958 dx
[XX
] += pbc
->fbox_diag
[XX
];
960 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
962 dx
[XX
] += pbc
->fbox_diag
[XX
];
966 d2min
+= dx
[XX
] * dx
[XX
];
968 if (d2min
> pbc
->max_cutoff2
)
970 copy_rvec(dx
, dx_start
);
971 copy_ivec(ishift
, ishift_start
);
972 /* Now try all possible shifts, when the distance is within max_cutoff
973 * it must be the shortest possible distance.
976 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
978 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
980 for (j
= 0; j
< DIM
; j
++)
984 d2trial
+= trial
[j
] * trial
[j
];
989 copy_rvec(trial
, dx
);
990 ivec_add(ishift_start
, pbc
->tric_shift
[i
], ishift
);
999 if (dx
[i
] > pbc
->hbox_diag
[i
])
1001 dx
[i
] -= pbc
->fbox_diag
[i
];
1004 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1006 dx
[i
] += pbc
->fbox_diag
[i
];
1012 if (dx
[i
] > pbc
->hbox_diag
[i
])
1014 rvec_dec(dx
, pbc
->box
[i
]);
1017 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1019 rvec_inc(dx
, pbc
->box
[i
]);
1023 case epbcdxSCREW_RECT
:
1024 /* The shift definition requires x first */
1025 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
1027 dx
[XX
] -= pbc
->fbox_diag
[XX
];
1030 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
1032 dx
[XX
] += pbc
->fbox_diag
[XX
];
1035 if (ishift
[XX
] == 1 || ishift
[XX
] == -1)
1037 /* Rotate around the x-axis in the middle of the box */
1038 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
1039 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
1041 /* Normal pbc for y and z */
1042 for (i
= YY
; i
<= ZZ
; i
++)
1044 if (dx
[i
] > pbc
->hbox_diag
[i
])
1046 dx
[i
] -= pbc
->fbox_diag
[i
];
1049 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1051 dx
[i
] += pbc
->fbox_diag
[i
];
1057 case epbcdxUNSUPPORTED
: break;
1060 "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
1063 is
= IVEC2IS(ishift
);
1066 range_check_mesg(is
, 0, SHIFTS
, "PBC shift vector index range check.");
1072 //! Compute distance vector in double precision
1073 void pbc_dx_d(const t_pbc
* pbc
, const dvec x1
, const dvec x2
, dvec dx
)
1076 dvec dx_start
, trial
;
1077 double d2min
, d2trial
;
1080 dvec_sub(x1
, x2
, dx
);
1082 switch (pbc
->ePBCDX
)
1084 case epbcdxRECTANGULAR
:
1086 for (i
= 0; i
< DIM
; i
++)
1090 while (dx
[i
] > pbc
->hbox_diag
[i
])
1092 dx
[i
] -= pbc
->fbox_diag
[i
];
1094 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
1096 dx
[i
] += pbc
->fbox_diag
[i
];
1101 case epbcdxTRICLINIC
:
1104 for (i
= DIM
- 1; i
>= 0; i
--)
1108 while (dx
[i
] > pbc
->hbox_diag
[i
])
1110 for (j
= i
; j
>= 0; j
--)
1112 dx
[j
] -= pbc
->box
[i
][j
];
1115 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
1117 for (j
= i
; j
>= 0; j
--)
1119 dx
[j
] += pbc
->box
[i
][j
];
1122 d2min
+= dx
[i
] * dx
[i
];
1125 if (d2min
> pbc
->max_cutoff2
)
1127 copy_dvec(dx
, dx_start
);
1128 /* Now try all possible shifts, when the distance is within max_cutoff
1129 * it must be the shortest possible distance.
1132 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
1134 for (j
= 0; j
< DIM
; j
++)
1136 trial
[j
] = dx_start
[j
] + pbc
->tric_vec
[i
][j
];
1139 for (j
= 0; j
< DIM
; j
++)
1143 d2trial
+= trial
[j
] * trial
[j
];
1146 if (d2trial
< d2min
)
1148 copy_dvec(trial
, dx
);
1155 case epbcdxSCREW_RECT
:
1156 /* The shift definition requires x first */
1158 while (dx
[XX
] > pbc
->hbox_diag
[XX
])
1160 dx
[XX
] -= pbc
->fbox_diag
[XX
];
1163 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
1165 dx
[XX
] += pbc
->fbox_diag
[YY
];
1170 /* Rotate around the x-axis in the middle of the box */
1171 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
1172 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
1174 /* Normal pbc for y and z */
1175 for (i
= YY
; i
<= ZZ
; i
++)
1177 while (dx
[i
] > pbc
->hbox_diag
[i
])
1179 dx
[i
] -= pbc
->fbox_diag
[i
];
1181 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
1183 dx
[i
] += pbc
->fbox_diag
[i
];
1188 case epbcdxUNSUPPORTED
: break;
1189 default: gmx_fatal(FARGS
, "Internal error in pbc_dx, set_pbc has not been called");
1193 void calc_shifts(const matrix box
, rvec shift_vec
[])
1195 for (int n
= 0, m
= -D_BOX_Z
; m
<= D_BOX_Z
; m
++)
1197 for (int l
= -D_BOX_Y
; l
<= D_BOX_Y
; l
++)
1199 for (int k
= -D_BOX_X
; k
<= D_BOX_X
; k
++, n
++)
1201 for (int d
= 0; d
< DIM
; d
++)
1203 shift_vec
[n
][d
] = k
* box
[XX
][d
] + l
* box
[YY
][d
] + m
* box
[ZZ
][d
];
1210 void calc_box_center(int ecenter
, const matrix box
, rvec box_center
)
1214 clear_rvec(box_center
);
1218 for (m
= 0; (m
< DIM
); m
++)
1220 for (d
= 0; d
< DIM
; d
++)
1222 box_center
[d
] += 0.5 * box
[m
][d
];
1227 for (d
= 0; d
< DIM
; d
++)
1229 box_center
[d
] = 0.5 * box
[d
][d
];
1232 case ecenterZERO
: break;
1233 default: gmx_fatal(FARGS
, "Unsupported value %d for ecenter", ecenter
);
1237 void calc_triclinic_images(const matrix box
, rvec img
[])
1241 /* Calculate 3 adjacent images in the xy-plane */
1242 copy_rvec(box
[0], img
[0]);
1243 copy_rvec(box
[1], img
[1]);
1246 svmul(-1, img
[1], img
[1]);
1248 rvec_sub(img
[1], img
[0], img
[2]);
1250 /* Get the next 3 in the xy-plane as mirror images */
1251 for (i
= 0; i
< 3; i
++)
1253 svmul(-1, img
[i
], img
[3 + i
]);
1256 /* Calculate the first 4 out of xy-plane images */
1257 copy_rvec(box
[2], img
[6]);
1260 svmul(-1, img
[6], img
[6]);
1262 for (i
= 0; i
< 3; i
++)
1264 rvec_add(img
[6], img
[i
+ 1], img
[7 + i
]);
1267 /* Mirror the last 4 from the previous in opposite rotation */
1268 for (i
= 0; i
< 4; i
++)
1270 svmul(-1, img
[6 + (2 + i
) % 4], img
[10 + i
]);
1274 void calc_compact_unitcell_vertices(int ecenter
, const matrix box
, rvec vert
[])
1276 rvec img
[NTRICIMG
], box_center
;
1277 int n
, i
, j
, tmp
[4], d
;
1278 const real oneFourth
= 0.25;
1280 calc_triclinic_images(box
, img
);
1283 for (i
= 2; i
<= 5; i
+= 3)
1294 tmp
[2] = (i
+ 1) % 6;
1295 tmp
[3] = tmp
[1] + 4;
1296 for (j
= 0; j
< 4; j
++)
1298 for (d
= 0; d
< DIM
; d
++)
1300 vert
[n
][d
] = img
[i
][d
] + img
[tmp
[j
]][d
] + img
[tmp
[(j
+ 1) % 4]][d
];
1305 for (i
= 7; i
<= 13; i
+= 6)
1307 tmp
[0] = (i
- 7) / 2;
1308 tmp
[1] = tmp
[0] + 1;
1318 for (j
= 0; j
< 4; j
++)
1320 for (d
= 0; d
< DIM
; d
++)
1322 vert
[n
][d
] = img
[i
][d
] + img
[tmp
[j
]][d
] + img
[tmp
[(j
+ 1) % 4]][d
];
1327 for (i
= 9; i
<= 11; i
+= 2)
1337 tmp
[1] = tmp
[0] + 1;
1347 for (j
= 0; j
< 4; j
++)
1349 for (d
= 0; d
< DIM
; d
++)
1351 vert
[n
][d
] = img
[i
][d
] + img
[tmp
[j
]][d
] + img
[tmp
[(j
+ 1) % 4]][d
];
1357 calc_box_center(ecenter
, box
, box_center
);
1358 for (i
= 0; i
< NCUCVERT
; i
++)
1360 for (d
= 0; d
< DIM
; d
++)
1362 vert
[i
][d
] = vert
[i
][d
] * oneFourth
+ box_center
[d
];
1367 int* compact_unitcell_edges()
1369 /* this is an index in vert[] (see calc_box_vertices) */
1370 /*static int edge[NCUCEDGE*2];*/
1372 static const int hexcon
[24] = { 0, 9, 1, 19, 2, 15, 3, 21, 4, 17, 5, 11,
1373 6, 23, 7, 13, 8, 20, 10, 18, 12, 16, 14, 22 };
1376 snew(edge
, NCUCEDGE
* 2);
1379 for (i
= 0; i
< 6; i
++)
1381 for (j
= 0; j
< 4; j
++)
1383 edge
[e
++] = 4 * i
+ j
;
1384 edge
[e
++] = 4 * i
+ (j
+ 1) % 4;
1387 for (i
= 0; i
< 12 * 2; i
++)
1389 edge
[e
++] = hexcon
[i
];
1395 void put_atoms_in_box(int ePBC
, const matrix box
, gmx::ArrayRef
<gmx::RVec
> x
)
1399 if (ePBC
== epbcSCREW
)
1401 gmx_fatal(FARGS
, "Sorry, %s pbc is not yet supported", epbc_names
[ePBC
]);
1415 for (gmx::index i
= 0; (i
< x
.ssize()); ++i
)
1417 for (m
= npbcdim
- 1; m
>= 0; m
--)
1421 for (d
= 0; d
<= m
; d
++)
1423 x
[i
][d
] += box
[m
][d
];
1426 while (x
[i
][m
] >= box
[m
][m
])
1428 for (d
= 0; d
<= m
; d
++)
1430 x
[i
][d
] -= box
[m
][d
];
1438 for (gmx::index i
= 0; (i
< x
.ssize()); ++i
)
1440 for (d
= 0; d
< npbcdim
; d
++)
1444 x
[i
][d
] += box
[d
][d
];
1446 while (x
[i
][d
] >= box
[d
][d
])
1448 x
[i
][d
] -= box
[d
][d
];
1455 void put_atoms_in_box_omp(int ePBC
, const matrix box
, gmx::ArrayRef
<gmx::RVec
> x
, gmx_unused
int nth
)
1457 #pragma omp parallel for num_threads(nth) schedule(static)
1458 for (int t
= 0; t
< nth
; t
++)
1462 size_t natoms
= x
.size();
1463 size_t offset
= (natoms
* t
) / nth
;
1464 size_t len
= (natoms
* (t
+ 1)) / nth
- offset
;
1465 put_atoms_in_box(ePBC
, box
, x
.subArray(offset
, len
));
1467 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1471 void put_atoms_in_triclinic_unitcell(int ecenter
, const matrix box
, gmx::ArrayRef
<gmx::RVec
> x
)
1473 rvec box_center
, shift_center
;
1474 real shm01
, shm02
, shm12
, shift
;
1477 calc_box_center(ecenter
, box
, box_center
);
1479 /* The product of matrix shm with a coordinate gives the shift vector
1480 which is required determine the periodic cell position */
1481 shm01
= box
[1][0] / box
[1][1];
1482 shm02
= (box
[1][1] * box
[2][0] - box
[2][1] * box
[1][0]) / (box
[1][1] * box
[2][2]);
1483 shm12
= box
[2][1] / box
[2][2];
1485 clear_rvec(shift_center
);
1486 for (d
= 0; d
< DIM
; d
++)
1488 rvec_inc(shift_center
, box
[d
]);
1490 svmul(0.5, shift_center
, shift_center
);
1491 rvec_sub(box_center
, shift_center
, shift_center
);
1493 shift_center
[0] = shm01
* shift_center
[1] + shm02
* shift_center
[2];
1494 shift_center
[1] = shm12
* shift_center
[2];
1495 shift_center
[2] = 0;
1497 for (gmx::index i
= 0; (i
< x
.ssize()); ++i
)
1499 for (m
= DIM
- 1; m
>= 0; m
--)
1501 shift
= shift_center
[m
];
1504 shift
+= shm01
* x
[i
][1] + shm02
* x
[i
][2];
1508 shift
+= shm12
* x
[i
][2];
1510 while (x
[i
][m
] - shift
< 0)
1512 for (d
= 0; d
<= m
; d
++)
1514 x
[i
][d
] += box
[m
][d
];
1517 while (x
[i
][m
] - shift
>= box
[m
][m
])
1519 for (d
= 0; d
<= m
; d
++)
1521 x
[i
][d
] -= box
[m
][d
];
1528 void put_atoms_in_compact_unitcell(int ePBC
, int ecenter
, const matrix box
, gmx::ArrayRef
<gmx::RVec
> x
)
1531 rvec box_center
, dx
;
1533 set_pbc(&pbc
, ePBC
, box
);
1535 if (pbc
.ePBCDX
== epbcdxUNSUPPORTED
)
1537 gmx_fatal(FARGS
, "Can not put atoms in compact unitcell with unsupported PBC");
1540 calc_box_center(ecenter
, box
, box_center
);
1541 for (gmx::index i
= 0; (i
< x
.ssize()); ++i
)
1543 pbc_dx(&pbc
, x
[i
], box_center
, dx
);
1544 rvec_add(box_center
, dx
, x
[i
]);
1548 /*! \brief Make molecules whole by shifting positions
1550 * \param[in] fplog Log file
1551 * \param[in] ePBC The PBC type
1552 * \param[in] box The simulation box
1553 * \param[in] mtop System topology definition
1554 * \param[in,out] x The coordinates of the atoms
1555 * \param[in] bFirst Specifier for first-time PBC removal
1557 static void low_do_pbc_mtop(FILE* fplog
, int ePBC
, const matrix box
, const gmx_mtop_t
* mtop
, rvec x
[], gmx_bool bFirst
)
1562 if (bFirst
&& fplog
)
1564 fprintf(fplog
, "Removing pbc first time\n");
1569 for (const gmx_molblock_t
& molb
: mtop
->molblock
)
1571 const gmx_moltype_t
& moltype
= mtop
->moltype
[molb
.type
];
1572 if (moltype
.atoms
.nr
== 1 || (!bFirst
&& moltype
.atoms
.nr
== 1))
1574 /* Just one atom or charge group in the molecule, no PBC required */
1575 as
+= molb
.nmol
* moltype
.atoms
.nr
;
1579 mk_graph_moltype(moltype
, graph
);
1581 for (mol
= 0; mol
< molb
.nmol
; mol
++)
1583 mk_mshift(fplog
, graph
, ePBC
, box
, x
+ as
);
1585 shift_self(graph
, box
, x
+ as
);
1586 /* The molecule is whole now.
1587 * We don't need the second mk_mshift call as in do_pbc_first,
1588 * since we no longer need this graph.
1591 as
+= moltype
.atoms
.nr
;
1599 void do_pbc_first_mtop(FILE* fplog
, int ePBC
, const matrix box
, const gmx_mtop_t
* mtop
, rvec x
[])
1601 low_do_pbc_mtop(fplog
, ePBC
, box
, mtop
, x
, TRUE
);
1604 void do_pbc_mtop(int ePBC
, const matrix box
, const gmx_mtop_t
* mtop
, rvec x
[])
1606 low_do_pbc_mtop(nullptr, ePBC
, box
, mtop
, x
, FALSE
);