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.
39 * Implements routines in pbc.h.
41 * Utility functions for handling periodic boundary conditions.
42 * Mainly used in analysis tools.
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. */
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
)
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
);
100 void dump_pbc(FILE *fp
, t_pbc
*pbc
)
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
)
126 ePBC
= guess_ePBC(box
);
129 if (ePBC
== epbcNONE
)
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
] ||
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.";
157 void matrix_convert(matrix box
, const rvec vec
, const rvec angleInDegrees
)
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
]));
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.
192 min_ss
= std::min(box
[XX
][XX
], box
[YY
][YY
]);
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
)
209 if (box
[XX
][XX
] > 0 && box
[YY
][YY
] > 0 && box
[ZZ
][ZZ
] > 0)
213 else 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)
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
]);
235 fprintf(debug
, "Guessed pbc = %s from the box matrix\n", epbc_names
[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;
248 /* correct elem d of vector v with vector d */
249 while (box
[v
][d
] > BOX_MARGIN_CORRECT
*0.5*box
[d
][d
])
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
]);
260 pr_rvecs(fplog
, 0, "new box", box
, DIM
);
262 if (shift
<= -maxshift
)
265 "Box was shifted at least %d times. Please see log-file.",
269 while (box
[v
][d
] < -BOX_MARGIN_CORRECT
*0.5*box
[d
][d
])
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
]);
280 pr_rvecs(fplog
, 0, "new box", box
, DIM
);
282 if (shift
>= maxshift
)
285 "Box was shifted at least %d times. Please see log-file.",
293 gmx_bool
correct_box(FILE *fplog
, int step
, tensor box
, t_graph
*graph
)
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
;
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 };
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
;
448 gmx_fatal(FARGS
, "Incorrect number of pbc dimensions with DD: %d",
451 pbc
->max_cutoff2
= max_cutoff2(ePBC
, box
);
453 if (pbc
->ePBCDX
== epbcdxTRICLINIC
||
454 pbc
->ePBCDX
== epbcdx2D_TRIC
||
455 pbc
->ePBCDX
== epbcdxSCREW_TRIC
)
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
++)
466 if (!bPBC
[ZZ
] && k
!= 0)
470 for (int jj
= 0; jj
< 3; jj
++)
473 if (!bPBC
[YY
] && j
!= 0)
477 for (int ii
= 0; ii
< 3; ii
++)
480 if (!bPBC
[XX
] && i
!= 0)
484 /* A shift is only useful when it is trilinic */
485 if (j
!= 0 || k
!= 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.
507 pos
[d
] = std::min( pbc
->hbox_diag
[d
], -trial
[d
]);
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
));
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
)
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
);
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
;
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
)
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
)
594 if (nullptr == domdecCells
)
596 low_set_pbc(pbc
, ePBC
, nullptr, box
);
600 if (ePBC
== epbcSCREW
&& domdecCells
[XX
] > 1)
602 /* The rotation has been taken care of during coordinate communication */
608 for (int i
= 0; i
< DIM
; i
++)
611 if (domdecCells
[i
] <= (bSingleDir
? 1 : 2) &&
612 !(ePBC
== epbcXY
&& i
== ZZ
))
621 low_set_pbc(pbc
, ePBC
, usePBC
, box
);
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
)
635 rvec dx_start
, trial
;
639 rvec_sub(x1
, x2
, dx
);
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
];
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 */
676 if (d2min
> pbc
->max_cutoff2
)
678 copy_rvec(dx
, dx_start
);
680 /* Now try all possible shifts, when the distance is within max_cutoff
681 * it must be the shortest possible distance.
684 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
686 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
687 d2trial
= norm2(trial
);
690 copy_rvec(trial
, dx
);
698 for (i
= 0; i
< DIM
; i
++)
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
];
715 for (i
= DIM
-1; i
>= 0; i
--)
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
);
740 /* Now try all possible shifts, when the distance is within max_cutoff
741 * it must be the shortest possible distance.
744 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
746 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
748 for (j
= 0; j
< DIM
; j
++)
752 d2trial
+= trial
[j
]*trial
[j
];
757 copy_rvec(trial
, dx
);
764 case epbcdxSCREW_RECT
:
765 /* The shift definition requires x first */
767 while (dx
[XX
] > pbc
->hbox_diag
[XX
])
769 dx
[XX
] -= pbc
->fbox_diag
[XX
];
772 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
774 dx
[XX
] += pbc
->fbox_diag
[YY
];
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
];
797 case epbcdxUNSUPPORTED
:
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
)
807 rvec dx_start
, trial
;
809 ivec ishift
, ishift_start
;
811 rvec_sub(x1
, x2
, dx
);
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
];
824 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
826 dx
[i
] += pbc
->fbox_diag
[i
];
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
];
849 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
851 for (j
= i
; j
>= 0; j
--)
853 dx
[j
] += pbc
->box
[i
][j
];
858 /* Allow 2 shifts in x */
859 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
861 dx
[XX
] -= pbc
->fbox_diag
[XX
];
863 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
865 dx
[XX
] -= pbc
->fbox_diag
[XX
];
869 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
871 dx
[XX
] += pbc
->fbox_diag
[XX
];
873 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
875 dx
[XX
] += pbc
->fbox_diag
[XX
];
879 /* dx is the distance in a rectangular box */
881 if (d2min
> pbc
->max_cutoff2
)
883 copy_rvec(dx
, dx_start
);
884 copy_ivec(ishift
, ishift_start
);
886 /* Now try all possible shifts, when the distance is within max_cutoff
887 * it must be the shortest possible distance.
890 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
892 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
893 d2trial
= norm2(trial
);
896 copy_rvec(trial
, dx
);
897 ivec_add(ishift_start
, pbc
->tric_shift
[i
], ishift
);
905 for (i
= 0; i
< DIM
; i
++)
909 if (dx
[i
] > pbc
->hbox_diag
[i
])
911 dx
[i
] -= pbc
->fbox_diag
[i
];
914 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
916 dx
[i
] += pbc
->fbox_diag
[i
];
924 for (i
= DIM
-1; i
>= 1; i
--)
928 if (dx
[i
] > pbc
->hbox_diag
[i
])
930 for (j
= i
; j
>= 0; j
--)
932 dx
[j
] -= pbc
->box
[i
][j
];
936 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
938 for (j
= i
; j
>= 0; j
--)
940 dx
[j
] += pbc
->box
[i
][j
];
944 d2min
+= dx
[i
]*dx
[i
];
949 /* Allow 2 shifts in x */
950 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
952 dx
[XX
] -= pbc
->fbox_diag
[XX
];
954 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
956 dx
[XX
] -= pbc
->fbox_diag
[XX
];
960 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
962 dx
[XX
] += pbc
->fbox_diag
[XX
];
964 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
966 dx
[XX
] += pbc
->fbox_diag
[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.
980 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
982 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
984 for (j
= 0; j
< DIM
; j
++)
988 d2trial
+= trial
[j
]*trial
[j
];
993 copy_rvec(trial
, dx
);
994 ivec_add(ishift_start
, pbc
->tric_shift
[i
], ishift
);
1003 if (dx
[i
] > pbc
->hbox_diag
[i
])
1005 dx
[i
] -= pbc
->fbox_diag
[i
];
1008 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1010 dx
[i
] += pbc
->fbox_diag
[i
];
1016 if (dx
[i
] > pbc
->hbox_diag
[i
])
1018 rvec_dec(dx
, pbc
->box
[i
]);
1021 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1023 rvec_inc(dx
, pbc
->box
[i
]);
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
];
1034 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
1036 dx
[XX
] += pbc
->fbox_diag
[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
];
1053 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1055 dx
[i
] += pbc
->fbox_diag
[i
];
1061 case epbcdxUNSUPPORTED
:
1064 gmx_fatal(FARGS
, "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
1067 is
= IVEC2IS(ishift
);
1070 range_check_mesg(is
, 0, SHIFTS
, "PBC shift vector index range check.");
1076 //! Compute distance vector in double precision
1077 void pbc_dx_d(const t_pbc
*pbc
, const dvec x1
, const dvec x2
, dvec dx
)
1080 dvec dx_start
, trial
;
1081 double d2min
, d2trial
;
1084 dvec_sub(x1
, x2
, dx
);
1086 switch (pbc
->ePBCDX
)
1088 case epbcdxRECTANGULAR
:
1090 for (i
= 0; i
< DIM
; i
++)
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
];
1105 case epbcdxTRICLINIC
:
1108 for (i
= DIM
-1; i
>= 0; i
--)
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.
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
];
1143 for (j
= 0; j
< DIM
; j
++)
1147 d2trial
+= trial
[j
]*trial
[j
];
1150 if (d2trial
< d2min
)
1152 copy_dvec(trial
, dx
);
1159 case epbcdxSCREW_RECT
:
1160 /* The shift definition requires x first */
1162 while (dx
[XX
] > pbc
->hbox_diag
[XX
])
1164 dx
[XX
] -= pbc
->fbox_diag
[XX
];
1167 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
1169 dx
[XX
] += pbc
->fbox_diag
[YY
];
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
];
1192 case epbcdxUNSUPPORTED
:
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
;
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
);
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
)
1228 clear_rvec(box_center
);
1232 for (m
= 0; (m
< DIM
); m
++)
1234 for (d
= 0; d
< DIM
; d
++)
1236 box_center
[d
] += 0.5*box
[m
][d
];
1241 for (d
= 0; d
< DIM
; d
++)
1243 box_center
[d
] = 0.5*box
[d
][d
];
1249 gmx_fatal(FARGS
, "Unsupported value %d for ecenter", ecenter
);
1253 void calc_triclinic_images(const matrix box
, rvec img
[])
1257 /* Calculate 3 adjacent images in the xy-plane */
1258 copy_rvec(box
[0], img
[0]);
1259 copy_rvec(box
[1], img
[1]);
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]);
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
);
1299 for (i
= 2; i
<= 5; i
+= 3)
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
];
1321 for (i
= 7; i
<= 13; i
+= 6)
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
];
1343 for (i
= 9; i
<= 11; i
+= 2)
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
];
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];*/
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
1395 snew(edge
, NCUCEDGE
*2);
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
];
1414 void put_atoms_in_box(int ePBC
, const matrix box
, gmx::ArrayRef
<gmx::RVec
> x
)
1418 if (ePBC
== epbcSCREW
)
1420 gmx_fatal(FARGS
, "Sorry, %s pbc is not yet supported", epbc_names
[ePBC
]);
1434 for (gmx::index i
= 0; (i
< x
.ssize()); ++i
)
1436 for (m
= npbcdim
-1; m
>= 0; m
--)
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
];
1457 for (gmx::index i
= 0; (i
< x
.ssize()); ++i
)
1459 for (d
= 0; d
< npbcdim
; d
++)
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
)
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
;
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
];
1526 shift
+= shm01
*x
[i
][1] + shm02
*x
[i
][2];
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
)
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
[],
1587 if (bFirst
&& fplog
)
1589 fprintf(fplog
, "Removing pbc first time\n");
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
;
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
;
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
);