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, 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/fileio/txtdump.h"
53 #include "gromacs/gmxlib/gmx_omp_nthreads.h"
54 #include "gromacs/legacyheaders/types/commrec.h"
55 #include "gromacs/math/utilities.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdtypes/inputrec.h"
58 #include "gromacs/pbcutil/ishift.h"
59 #include "gromacs/pbcutil/mshift.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", NULL
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 int inputrec2nboundeddim(const t_inputrec
*ir
)
102 if (ir
->nwall
== 2 && ir
->ePBC
== epbcXY
)
108 return ePBC2npbcdim(ir
->ePBC
);
112 void dump_pbc(FILE *fp
, t_pbc
*pbc
)
116 fprintf(fp
, "ePBCDX = %d\n", pbc
->ePBCDX
);
117 pr_rvecs(fp
, 0, "box", pbc
->box
, DIM
);
118 pr_rvecs(fp
, 0, "fbox_diag", &pbc
->fbox_diag
, 1);
119 pr_rvecs(fp
, 0, "hbox_diag", &pbc
->hbox_diag
, 1);
120 pr_rvecs(fp
, 0, "mhbox_diag", &pbc
->mhbox_diag
, 1);
121 rvec_add(pbc
->hbox_diag
, pbc
->mhbox_diag
, sum_box
);
122 pr_rvecs(fp
, 0, "sum of the above two", &sum_box
, 1);
123 fprintf(fp
, "max_cutoff2 = %g\n", pbc
->max_cutoff2
);
124 fprintf(fp
, "ntric_vec = %d\n", pbc
->ntric_vec
);
125 if (pbc
->ntric_vec
> 0)
127 pr_ivecs(fp
, 0, "tric_shift", pbc
->tric_shift
, pbc
->ntric_vec
, FALSE
);
128 pr_rvecs(fp
, 0, "tric_vec", pbc
->tric_vec
, pbc
->ntric_vec
);
132 const char *check_box(int ePBC
, matrix box
)
138 ePBC
= guess_ePBC(box
);
141 if (ePBC
== epbcNONE
)
146 if ((box
[XX
][YY
] != 0) || (box
[XX
][ZZ
] != 0) || (box
[YY
][ZZ
] != 0))
148 ptr
= "Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.";
150 else if (ePBC
== epbcSCREW
&& (box
[YY
][XX
] != 0 || box
[ZZ
][XX
] != 0))
152 ptr
= "The unit cell can not have off-diagonal x-components with screw pbc";
154 else if (std::fabs(box
[YY
][XX
]) > BOX_MARGIN
*0.5*box
[XX
][XX
] ||
156 (std::fabs(box
[ZZ
][XX
]) > BOX_MARGIN
*0.5*box
[XX
][XX
] ||
157 std::fabs(box
[ZZ
][YY
]) > BOX_MARGIN
*0.5*box
[YY
][YY
])))
159 ptr
= "Triclinic box is too skewed.";
169 real
max_cutoff2(int ePBC
, matrix box
)
171 real min_hv2
, min_ss
;
172 const real oneFourth
= 0.25;
174 /* Physical limitation of the cut-off
175 * by half the length of the shortest box vector.
177 min_hv2
= oneFourth
* std::min(norm2(box
[XX
]), norm2(box
[YY
]));
180 min_hv2
= std::min(min_hv2
, oneFourth
* norm2(box
[ZZ
]));
183 /* Limitation to the smallest diagonal element due to optimizations:
184 * checking only linear combinations of single box-vectors (2 in x)
185 * in the grid search and pbc_dx is a lot faster
186 * than checking all possible combinations.
190 min_ss
= std::min(box
[XX
][XX
], box
[YY
][YY
]);
194 min_ss
= std::min(box
[XX
][XX
], std::min(box
[YY
][YY
] - std::fabs(box
[ZZ
][YY
]), box
[ZZ
][ZZ
]));
197 return std::min(min_hv2
, min_ss
*min_ss
);
200 //! Set to true if warning has been printed
201 static gmx_bool bWarnedGuess
= FALSE
;
203 int guess_ePBC(matrix box
)
207 if (box
[XX
][XX
] > 0 && box
[YY
][YY
] > 0 && box
[ZZ
][ZZ
] > 0)
211 else if (box
[XX
][XX
] > 0 && box
[YY
][YY
] > 0 && box
[ZZ
][ZZ
] == 0)
215 else if (box
[XX
][XX
] == 0 && box
[YY
][YY
] == 0 && box
[ZZ
][ZZ
] == 0)
223 fprintf(stderr
, "WARNING: Unsupported box diagonal %f %f %f, "
224 "will not use periodic boundary conditions\n\n",
225 box
[XX
][XX
], box
[YY
][YY
], box
[ZZ
][ZZ
]);
233 fprintf(debug
, "Guessed pbc = %s from the box matrix\n", epbc_names
[ePBC
]);
239 //! Check if the box still obeys the restrictions, if not, correct it
240 static int correct_box_elem(FILE *fplog
, int step
, tensor box
, int v
, int d
)
242 int shift
, maxshift
= 10;
246 /* correct elem d of vector v with vector d */
247 while (box
[v
][d
] > BOX_MARGIN_CORRECT
*0.5*box
[d
][d
])
251 fprintf(fplog
, "Step %d: correcting invalid box:\n", step
);
252 pr_rvecs(fplog
, 0, "old box", box
, DIM
);
254 rvec_dec(box
[v
], box
[d
]);
258 pr_rvecs(fplog
, 0, "new box", box
, DIM
);
260 if (shift
<= -maxshift
)
263 "Box was shifted at least %d times. Please see log-file.",
267 while (box
[v
][d
] < -BOX_MARGIN_CORRECT
*0.5*box
[d
][d
])
271 fprintf(fplog
, "Step %d: correcting invalid box:\n", step
);
272 pr_rvecs(fplog
, 0, "old box", box
, DIM
);
274 rvec_inc(box
[v
], box
[d
]);
278 pr_rvecs(fplog
, 0, "new box", box
, DIM
);
280 if (shift
>= maxshift
)
283 "Box was shifted at least %d times. Please see log-file.",
291 gmx_bool
correct_box(FILE *fplog
, int step
, tensor box
, t_graph
*graph
)
296 zy
= correct_box_elem(fplog
, step
, box
, ZZ
, YY
);
297 zx
= correct_box_elem(fplog
, step
, box
, ZZ
, XX
);
298 yx
= correct_box_elem(fplog
, step
, box
, YY
, XX
);
300 bCorrected
= (zy
|| zx
|| yx
);
302 if (bCorrected
&& graph
)
304 /* correct the graph */
305 for (i
= graph
->at_start
; i
< graph
->at_end
; i
++)
307 graph
->ishift
[i
][YY
] -= graph
->ishift
[i
][ZZ
]*zy
;
308 graph
->ishift
[i
][XX
] -= graph
->ishift
[i
][ZZ
]*zx
;
309 graph
->ishift
[i
][XX
] -= graph
->ishift
[i
][YY
]*yx
;
316 int ndof_com(t_inputrec
*ir
)
327 n
= (ir
->nwall
== 0 ? 3 : 2);
333 gmx_incons("Unknown pbc in calc_nrdf");
339 //! Do the real arithmetic for filling the pbc struct
340 static void low_set_pbc(t_pbc
*pbc
, int ePBC
, ivec
*dd_nc
, matrix box
)
342 int order
[3] = { 0, -1, 1 };
347 pbc
->ndim_ePBC
= ePBC2npbcdim(ePBC
);
349 copy_mat(box
, pbc
->box
);
350 pbc
->max_cutoff2
= 0;
354 for (int i
= 0; (i
< DIM
); i
++)
356 pbc
->fbox_diag
[i
] = box
[i
][i
];
357 pbc
->hbox_diag
[i
] = pbc
->fbox_diag
[i
]*0.5;
358 pbc
->mhbox_diag
[i
] = -pbc
->hbox_diag
[i
];
361 ptr
= check_box(ePBC
, box
);
362 if (ePBC
== epbcNONE
)
364 pbc
->ePBCDX
= epbcdxNOPBC
;
368 fprintf(stderr
, "Warning: %s\n", ptr
);
369 pr_rvecs(stderr
, 0, " Box", box
, DIM
);
370 fprintf(stderr
, " Can not fix pbc.\n\n");
371 pbc
->ePBCDX
= epbcdxUNSUPPORTED
;
375 if (ePBC
== epbcSCREW
&& NULL
!= dd_nc
)
377 /* This combinated should never appear here */
378 gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
382 for (int i
= 0; i
< DIM
; i
++)
384 if ((dd_nc
&& (*dd_nc
)[i
] > 1) || (ePBC
== epbcXY
&& i
== ZZ
))
397 /* 1D pbc is not an mdp option and it is therefore only used
398 * with single shifts.
400 pbc
->ePBCDX
= epbcdx1D_RECT
;
401 for (int i
= 0; i
< DIM
; i
++)
408 GMX_ASSERT(pbc
->dim
< DIM
, "Dimension for PBC incorrect");
409 for (int i
= 0; i
< pbc
->dim
; i
++)
411 if (pbc
->box
[pbc
->dim
][i
] != 0)
413 pbc
->ePBCDX
= epbcdx1D_TRIC
;
418 pbc
->ePBCDX
= epbcdx2D_RECT
;
419 for (int i
= 0; i
< DIM
; i
++)
426 for (int i
= 0; i
< DIM
; i
++)
430 for (int j
= 0; j
< i
; j
++)
432 if (pbc
->box
[i
][j
] != 0)
434 pbc
->ePBCDX
= epbcdx2D_TRIC
;
441 if (ePBC
!= epbcSCREW
)
445 pbc
->ePBCDX
= epbcdxTRICLINIC
;
449 pbc
->ePBCDX
= epbcdxRECTANGULAR
;
454 pbc
->ePBCDX
= (box
[ZZ
][YY
] == 0 ? epbcdxSCREW_RECT
: epbcdxSCREW_TRIC
);
455 if (pbc
->ePBCDX
== epbcdxSCREW_TRIC
)
458 "Screw pbc is not yet implemented for triclinic boxes.\n"
459 "Can not fix pbc.\n");
460 pbc
->ePBCDX
= epbcdxUNSUPPORTED
;
465 gmx_fatal(FARGS
, "Incorrect number of pbc dimensions with DD: %d",
468 pbc
->max_cutoff2
= max_cutoff2(ePBC
, box
);
470 if (pbc
->ePBCDX
== epbcdxTRICLINIC
||
471 pbc
->ePBCDX
== epbcdx2D_TRIC
||
472 pbc
->ePBCDX
== epbcdxSCREW_TRIC
)
476 pr_rvecs(debug
, 0, "Box", box
, DIM
);
477 fprintf(debug
, "max cutoff %.3f\n", sqrt(pbc
->max_cutoff2
));
479 /* We will only need single shifts here */
480 for (int kk
= 0; kk
< 3; kk
++)
483 if (!bPBC
[ZZ
] && k
!= 0)
487 for (int jj
= 0; jj
< 3; jj
++)
490 if (!bPBC
[YY
] && j
!= 0)
494 for (int ii
= 0; ii
< 3; ii
++)
497 if (!bPBC
[XX
] && i
!= 0)
501 /* A shift is only useful when it is trilinic */
502 if (j
!= 0 || k
!= 0)
509 for (int d
= 0; d
< DIM
; d
++)
511 trial
[d
] = i
*box
[XX
][d
] + j
*box
[YY
][d
] + k
*box
[ZZ
][d
];
512 /* Choose the vector within the brick around 0,0,0 that
513 * will become the shortest due to shift try.
524 pos
[d
] = std::min( pbc
->hbox_diag
[d
], -trial
[d
]);
528 pos
[d
] = std::max(-pbc
->hbox_diag
[d
], -trial
[d
]);
531 d2old
+= sqr(pos
[d
]);
532 d2new
+= sqr(pos
[d
] + trial
[d
]);
534 if (BOX_MARGIN
*d2new
< d2old
)
536 /* Check if shifts with one box vector less do better */
537 gmx_bool bUse
= TRUE
;
538 for (int dd
= 0; dd
< DIM
; dd
++)
540 int shift
= (dd
== 0 ? i
: (dd
== 1 ? j
: k
));
544 for (int d
= 0; d
< DIM
; d
++)
546 d2new_c
+= sqr(pos
[d
] + trial
[d
] - shift
*box
[dd
][d
]);
548 if (d2new_c
<= BOX_MARGIN
*d2new
)
556 /* Accept this shift vector. */
557 if (pbc
->ntric_vec
>= MAX_NTRICVEC
)
559 fprintf(stderr
, "\nWARNING: Found more than %d triclinic correction vectors, ignoring some.\n"
560 " There is probably something wrong with your box.\n", MAX_NTRICVEC
);
561 pr_rvecs(stderr
, 0, " Box", box
, DIM
);
565 copy_rvec(trial
, pbc
->tric_vec
[pbc
->ntric_vec
]);
566 pbc
->tric_shift
[pbc
->ntric_vec
][XX
] = i
;
567 pbc
->tric_shift
[pbc
->ntric_vec
][YY
] = j
;
568 pbc
->tric_shift
[pbc
->ntric_vec
][ZZ
] = k
;
573 fprintf(debug
, " tricvec %2d = %2d %2d %2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
574 pbc
->ntric_vec
, i
, j
, k
,
575 sqrt(d2old
), sqrt(d2new
),
576 trial
[XX
], trial
[YY
], trial
[ZZ
],
577 pos
[XX
], pos
[YY
], pos
[ZZ
]);
590 void set_pbc(t_pbc
*pbc
, int ePBC
, matrix box
)
594 ePBC
= guess_ePBC(box
);
597 low_set_pbc(pbc
, ePBC
, NULL
, box
);
600 t_pbc
*set_pbc_dd(t_pbc
*pbc
, int ePBC
,
601 struct gmx_domdec_t
*dd
, gmx_bool bSingleDir
, matrix box
)
612 if (ePBC
== epbcSCREW
&& dd
->nc
[XX
] > 1)
614 /* The rotation has been taken care of during coordinate communication */
618 for (i
= 0; i
< DIM
; i
++)
620 if (dd
->nc
[i
] <= (bSingleDir
? 1 : 2))
623 if (!(ePBC
== epbcXY
&& i
== ZZ
))
637 low_set_pbc(pbc
, ePBC
, npbcdim
< DIM
? &nc2
: NULL
, box
);
640 return (npbcdim
> 0 ? pbc
: NULL
);
643 void pbc_dx(const t_pbc
*pbc
, const rvec x1
, const rvec x2
, rvec dx
)
646 rvec dx_start
, trial
;
650 rvec_sub(x1
, x2
, dx
);
654 case epbcdxRECTANGULAR
:
655 for (i
= 0; i
< DIM
; i
++)
657 while (dx
[i
] > pbc
->hbox_diag
[i
])
659 dx
[i
] -= pbc
->fbox_diag
[i
];
661 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
663 dx
[i
] += pbc
->fbox_diag
[i
];
667 case epbcdxTRICLINIC
:
668 for (i
= DIM
-1; i
>= 0; i
--)
670 while (dx
[i
] > pbc
->hbox_diag
[i
])
672 for (j
= i
; j
>= 0; j
--)
674 dx
[j
] -= pbc
->box
[i
][j
];
677 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
679 for (j
= i
; j
>= 0; j
--)
681 dx
[j
] += pbc
->box
[i
][j
];
685 /* dx is the distance in a rectangular box */
687 if (d2min
> pbc
->max_cutoff2
)
689 copy_rvec(dx
, dx_start
);
691 /* Now try all possible shifts, when the distance is within max_cutoff
692 * it must be the shortest possible distance.
695 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
697 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
698 d2trial
= norm2(trial
);
701 copy_rvec(trial
, dx
);
709 for (i
= 0; i
< DIM
; i
++)
713 while (dx
[i
] > pbc
->hbox_diag
[i
])
715 dx
[i
] -= pbc
->fbox_diag
[i
];
717 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
719 dx
[i
] += pbc
->fbox_diag
[i
];
726 for (i
= DIM
-1; i
>= 0; i
--)
730 while (dx
[i
] > pbc
->hbox_diag
[i
])
732 for (j
= i
; j
>= 0; j
--)
734 dx
[j
] -= pbc
->box
[i
][j
];
737 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
739 for (j
= i
; j
>= 0; j
--)
741 dx
[j
] += pbc
->box
[i
][j
];
744 d2min
+= dx
[i
]*dx
[i
];
747 if (d2min
> pbc
->max_cutoff2
)
749 copy_rvec(dx
, dx_start
);
751 /* Now try all possible shifts, when the distance is within max_cutoff
752 * it must be the shortest possible distance.
755 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
757 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
759 for (j
= 0; j
< DIM
; j
++)
763 d2trial
+= trial
[j
]*trial
[j
];
768 copy_rvec(trial
, dx
);
775 case epbcdxSCREW_RECT
:
776 /* The shift definition requires x first */
778 while (dx
[XX
] > pbc
->hbox_diag
[XX
])
780 dx
[XX
] -= pbc
->fbox_diag
[XX
];
783 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
785 dx
[XX
] += pbc
->fbox_diag
[YY
];
790 /* Rotate around the x-axis in the middle of the box */
791 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
792 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
794 /* Normal pbc for y and z */
795 for (i
= YY
; i
<= ZZ
; i
++)
797 while (dx
[i
] > pbc
->hbox_diag
[i
])
799 dx
[i
] -= pbc
->fbox_diag
[i
];
801 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
803 dx
[i
] += pbc
->fbox_diag
[i
];
808 case epbcdxUNSUPPORTED
:
811 gmx_fatal(FARGS
, "Internal error in pbc_dx, set_pbc has not been called");
816 int pbc_dx_aiuc(const t_pbc
*pbc
, const rvec x1
, const rvec x2
, rvec dx
)
819 rvec dx_start
, trial
;
821 ivec ishift
, ishift_start
;
823 rvec_sub(x1
, x2
, dx
);
828 case epbcdxRECTANGULAR
:
829 for (i
= 0; i
< DIM
; i
++)
831 if (dx
[i
] > pbc
->hbox_diag
[i
])
833 dx
[i
] -= pbc
->fbox_diag
[i
];
836 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
838 dx
[i
] += pbc
->fbox_diag
[i
];
843 case epbcdxTRICLINIC
:
844 /* For triclinic boxes the performance difference between
845 * if/else and two while loops is negligible.
846 * However, the while version can cause extreme delays
847 * before a simulation crashes due to large forces which
848 * can cause unlimited displacements.
849 * Also allowing multiple shifts would index fshift beyond bounds.
851 for (i
= DIM
-1; i
>= 1; i
--)
853 if (dx
[i
] > pbc
->hbox_diag
[i
])
855 for (j
= i
; j
>= 0; j
--)
857 dx
[j
] -= pbc
->box
[i
][j
];
861 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
863 for (j
= i
; j
>= 0; j
--)
865 dx
[j
] += pbc
->box
[i
][j
];
870 /* Allow 2 shifts in x */
871 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
873 dx
[XX
] -= pbc
->fbox_diag
[XX
];
875 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
877 dx
[XX
] -= pbc
->fbox_diag
[XX
];
881 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
883 dx
[XX
] += pbc
->fbox_diag
[XX
];
885 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
887 dx
[XX
] += pbc
->fbox_diag
[XX
];
891 /* dx is the distance in a rectangular box */
893 if (d2min
> pbc
->max_cutoff2
)
895 copy_rvec(dx
, dx_start
);
896 copy_ivec(ishift
, ishift_start
);
898 /* Now try all possible shifts, when the distance is within max_cutoff
899 * it must be the shortest possible distance.
902 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
904 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
905 d2trial
= norm2(trial
);
908 copy_rvec(trial
, dx
);
909 ivec_add(ishift_start
, pbc
->tric_shift
[i
], ishift
);
917 for (i
= 0; i
< DIM
; i
++)
921 if (dx
[i
] > pbc
->hbox_diag
[i
])
923 dx
[i
] -= pbc
->fbox_diag
[i
];
926 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
928 dx
[i
] += pbc
->fbox_diag
[i
];
936 for (i
= DIM
-1; i
>= 1; i
--)
940 if (dx
[i
] > pbc
->hbox_diag
[i
])
942 for (j
= i
; j
>= 0; j
--)
944 dx
[j
] -= pbc
->box
[i
][j
];
948 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
950 for (j
= i
; j
>= 0; j
--)
952 dx
[j
] += pbc
->box
[i
][j
];
956 d2min
+= dx
[i
]*dx
[i
];
961 /* Allow 2 shifts in x */
962 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
964 dx
[XX
] -= pbc
->fbox_diag
[XX
];
966 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
968 dx
[XX
] -= pbc
->fbox_diag
[XX
];
972 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
974 dx
[XX
] += pbc
->fbox_diag
[XX
];
976 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
978 dx
[XX
] += pbc
->fbox_diag
[XX
];
982 d2min
+= dx
[XX
]*dx
[XX
];
984 if (d2min
> pbc
->max_cutoff2
)
986 copy_rvec(dx
, dx_start
);
987 copy_ivec(ishift
, ishift_start
);
988 /* Now try all possible shifts, when the distance is within max_cutoff
989 * it must be the shortest possible distance.
992 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
994 rvec_add(dx_start
, pbc
->tric_vec
[i
], trial
);
996 for (j
= 0; j
< DIM
; j
++)
1000 d2trial
+= trial
[j
]*trial
[j
];
1003 if (d2trial
< d2min
)
1005 copy_rvec(trial
, dx
);
1006 ivec_add(ishift_start
, pbc
->tric_shift
[i
], ishift
);
1015 if (dx
[i
] > pbc
->hbox_diag
[i
])
1017 dx
[i
] -= pbc
->fbox_diag
[i
];
1020 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1022 dx
[i
] += pbc
->fbox_diag
[i
];
1028 if (dx
[i
] > pbc
->hbox_diag
[i
])
1030 rvec_dec(dx
, pbc
->box
[i
]);
1033 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1035 rvec_inc(dx
, pbc
->box
[i
]);
1039 case epbcdxSCREW_RECT
:
1040 /* The shift definition requires x first */
1041 if (dx
[XX
] > pbc
->hbox_diag
[XX
])
1043 dx
[XX
] -= pbc
->fbox_diag
[XX
];
1046 else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
1048 dx
[XX
] += pbc
->fbox_diag
[XX
];
1051 if (ishift
[XX
] == 1 || ishift
[XX
] == -1)
1053 /* Rotate around the x-axis in the middle of the box */
1054 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
1055 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
1057 /* Normal pbc for y and z */
1058 for (i
= YY
; i
<= ZZ
; i
++)
1060 if (dx
[i
] > pbc
->hbox_diag
[i
])
1062 dx
[i
] -= pbc
->fbox_diag
[i
];
1065 else if (dx
[i
] <= pbc
->mhbox_diag
[i
])
1067 dx
[i
] += pbc
->fbox_diag
[i
];
1073 case epbcdxUNSUPPORTED
:
1076 gmx_fatal(FARGS
, "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
1080 is
= IVEC2IS(ishift
);
1083 range_check_mesg(is
, 0, SHIFTS
, "PBC shift vector index range check.");
1089 //! Compute distance vector in double precision
1090 void pbc_dx_d(const t_pbc
*pbc
, const dvec x1
, const dvec x2
, dvec dx
)
1093 dvec dx_start
, trial
;
1094 double d2min
, d2trial
;
1097 dvec_sub(x1
, x2
, dx
);
1099 switch (pbc
->ePBCDX
)
1101 case epbcdxRECTANGULAR
:
1103 for (i
= 0; i
< DIM
; i
++)
1107 while (dx
[i
] > pbc
->hbox_diag
[i
])
1109 dx
[i
] -= pbc
->fbox_diag
[i
];
1111 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
1113 dx
[i
] += pbc
->fbox_diag
[i
];
1118 case epbcdxTRICLINIC
:
1121 for (i
= DIM
-1; i
>= 0; i
--)
1125 while (dx
[i
] > pbc
->hbox_diag
[i
])
1127 for (j
= i
; j
>= 0; j
--)
1129 dx
[j
] -= pbc
->box
[i
][j
];
1132 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
1134 for (j
= i
; j
>= 0; j
--)
1136 dx
[j
] += pbc
->box
[i
][j
];
1139 d2min
+= dx
[i
]*dx
[i
];
1142 if (d2min
> pbc
->max_cutoff2
)
1144 copy_dvec(dx
, dx_start
);
1145 /* Now try all possible shifts, when the distance is within max_cutoff
1146 * it must be the shortest possible distance.
1149 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
))
1151 for (j
= 0; j
< DIM
; j
++)
1153 trial
[j
] = dx_start
[j
] + pbc
->tric_vec
[i
][j
];
1156 for (j
= 0; j
< DIM
; j
++)
1160 d2trial
+= trial
[j
]*trial
[j
];
1163 if (d2trial
< d2min
)
1165 copy_dvec(trial
, dx
);
1172 case epbcdxSCREW_RECT
:
1173 /* The shift definition requires x first */
1175 while (dx
[XX
] > pbc
->hbox_diag
[XX
])
1177 dx
[XX
] -= pbc
->fbox_diag
[XX
];
1180 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
])
1182 dx
[XX
] += pbc
->fbox_diag
[YY
];
1187 /* Rotate around the x-axis in the middle of the box */
1188 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
1189 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
1191 /* Normal pbc for y and z */
1192 for (i
= YY
; i
<= ZZ
; i
++)
1194 while (dx
[i
] > pbc
->hbox_diag
[i
])
1196 dx
[i
] -= pbc
->fbox_diag
[i
];
1198 while (dx
[i
] <= pbc
->mhbox_diag
[i
])
1200 dx
[i
] += pbc
->fbox_diag
[i
];
1205 case epbcdxUNSUPPORTED
:
1208 gmx_fatal(FARGS
, "Internal error in pbc_dx, set_pbc has not been called");
1213 //! Compute the box image corresponding to input vectors
1214 gmx_bool
image_rect(ivec xi
, ivec xj
, ivec box_size
, real rlong2
, int *shift
, real
*r2
)
1222 for (m
= 0; (m
< DIM
); m
++)
1254 gmx_bool
image_cylindric(ivec xi
, ivec xj
, ivec box_size
, real rlong2
,
1255 int *shift
, real
*r2
)
1263 for (m
= 0; (m
< DIM
); m
++)
1299 void calc_shifts(matrix box
, rvec shift_vec
[])
1301 int k
, l
, m
, d
, n
, test
;
1304 for (m
= -D_BOX_Z
; m
<= D_BOX_Z
; m
++)
1306 for (l
= -D_BOX_Y
; l
<= D_BOX_Y
; l
++)
1308 for (k
= -D_BOX_X
; k
<= D_BOX_X
; k
++, n
++)
1310 test
= XYZ2IS(k
, l
, m
);
1313 gmx_incons("inconsistent shift numbering");
1315 for (d
= 0; d
< DIM
; d
++)
1317 shift_vec
[n
][d
] = k
*box
[XX
][d
] + l
*box
[YY
][d
] + m
*box
[ZZ
][d
];
1324 void calc_box_center(int ecenter
, matrix box
, rvec box_center
)
1328 clear_rvec(box_center
);
1332 for (m
= 0; (m
< DIM
); m
++)
1334 for (d
= 0; d
< DIM
; d
++)
1336 box_center
[d
] += 0.5*box
[m
][d
];
1341 for (d
= 0; d
< DIM
; d
++)
1343 box_center
[d
] = 0.5*box
[d
][d
];
1349 gmx_fatal(FARGS
, "Unsupported value %d for ecenter", ecenter
);
1353 void calc_triclinic_images(matrix box
, rvec img
[])
1357 /* Calculate 3 adjacent images in the xy-plane */
1358 copy_rvec(box
[0], img
[0]);
1359 copy_rvec(box
[1], img
[1]);
1362 svmul(-1, img
[1], img
[1]);
1364 rvec_sub(img
[1], img
[0], img
[2]);
1366 /* Get the next 3 in the xy-plane as mirror images */
1367 for (i
= 0; i
< 3; i
++)
1369 svmul(-1, img
[i
], img
[3+i
]);
1372 /* Calculate the first 4 out of xy-plane images */
1373 copy_rvec(box
[2], img
[6]);
1376 svmul(-1, img
[6], img
[6]);
1378 for (i
= 0; i
< 3; i
++)
1380 rvec_add(img
[6], img
[i
+1], img
[7+i
]);
1383 /* Mirror the last 4 from the previous in opposite rotation */
1384 for (i
= 0; i
< 4; i
++)
1386 svmul(-1, img
[6 + (2+i
) % 4], img
[10+i
]);
1390 void calc_compact_unitcell_vertices(int ecenter
, matrix box
, rvec vert
[])
1392 rvec img
[NTRICIMG
], box_center
;
1393 int n
, i
, j
, tmp
[4], d
;
1394 const real oneFourth
= 0.25;
1396 calc_triclinic_images(box
, img
);
1399 for (i
= 2; i
<= 5; i
+= 3)
1412 for (j
= 0; j
< 4; j
++)
1414 for (d
= 0; d
< DIM
; d
++)
1416 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1421 for (i
= 7; i
<= 13; i
+= 6)
1434 for (j
= 0; j
< 4; j
++)
1436 for (d
= 0; d
< DIM
; d
++)
1438 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1443 for (i
= 9; i
<= 11; i
+= 2)
1463 for (j
= 0; j
< 4; j
++)
1465 for (d
= 0; d
< DIM
; d
++)
1467 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1473 calc_box_center(ecenter
, box
, box_center
);
1474 for (i
= 0; i
< NCUCVERT
; i
++)
1476 for (d
= 0; d
< DIM
; d
++)
1478 vert
[i
][d
] = vert
[i
][d
]*oneFourth
+box_center
[d
];
1483 int *compact_unitcell_edges()
1485 /* this is an index in vert[] (see calc_box_vertices) */
1486 /*static int edge[NCUCEDGE*2];*/
1488 static const int hexcon
[24] = {
1489 0, 9, 1, 19, 2, 15, 3, 21,
1490 4, 17, 5, 11, 6, 23, 7, 13,
1491 8, 20, 10, 18, 12, 16, 14, 22
1495 snew(edge
, NCUCEDGE
*2);
1498 for (i
= 0; i
< 6; i
++)
1500 for (j
= 0; j
< 4; j
++)
1502 edge
[e
++] = 4*i
+ j
;
1503 edge
[e
++] = 4*i
+ (j
+1) % 4;
1506 for (i
= 0; i
< 12*2; i
++)
1508 edge
[e
++] = hexcon
[i
];
1514 void put_atoms_in_box_omp(int ePBC
, matrix box
, int natoms
, rvec x
[])
1517 nth
= gmx_omp_nthreads_get(emntDefault
);
1519 #pragma omp parallel for num_threads(nth) schedule(static)
1520 for (t
= 0; t
< nth
; t
++)
1526 offset
= (natoms
*t
)/nth
;
1527 len
= (natoms
*(t
+ 1))/nth
- offset
;
1528 put_atoms_in_box(ePBC
, box
, len
, x
+ offset
);
1530 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
1534 void put_atoms_in_box(int ePBC
, matrix box
, int natoms
, rvec x
[])
1536 int npbcdim
, i
, m
, d
;
1538 if (ePBC
== epbcSCREW
)
1540 gmx_fatal(FARGS
, "Sorry, %s pbc is not yet supported", epbc_names
[ePBC
]);
1554 for (i
= 0; (i
< natoms
); i
++)
1556 for (m
= npbcdim
-1; m
>= 0; m
--)
1560 for (d
= 0; d
<= m
; d
++)
1562 x
[i
][d
] += box
[m
][d
];
1565 while (x
[i
][m
] >= box
[m
][m
])
1567 for (d
= 0; d
<= m
; d
++)
1569 x
[i
][d
] -= box
[m
][d
];
1577 for (i
= 0; i
< natoms
; i
++)
1579 for (d
= 0; d
< npbcdim
; d
++)
1583 x
[i
][d
] += box
[d
][d
];
1585 while (x
[i
][d
] >= box
[d
][d
])
1587 x
[i
][d
] -= box
[d
][d
];
1594 void put_atoms_in_triclinic_unitcell(int ecenter
, matrix box
,
1595 int natoms
, rvec x
[])
1597 rvec box_center
, shift_center
;
1598 real shm01
, shm02
, shm12
, shift
;
1601 calc_box_center(ecenter
, box
, box_center
);
1603 /* The product of matrix shm with a coordinate gives the shift vector
1604 which is required determine the periodic cell position */
1605 shm01
= box
[1][0]/box
[1][1];
1606 shm02
= (box
[1][1]*box
[2][0] - box
[2][1]*box
[1][0])/(box
[1][1]*box
[2][2]);
1607 shm12
= box
[2][1]/box
[2][2];
1609 clear_rvec(shift_center
);
1610 for (d
= 0; d
< DIM
; d
++)
1612 rvec_inc(shift_center
, box
[d
]);
1614 svmul(0.5, shift_center
, shift_center
);
1615 rvec_sub(box_center
, shift_center
, shift_center
);
1617 shift_center
[0] = shm01
*shift_center
[1] + shm02
*shift_center
[2];
1618 shift_center
[1] = shm12
*shift_center
[2];
1619 shift_center
[2] = 0;
1621 for (i
= 0; (i
< natoms
); i
++)
1623 for (m
= DIM
-1; m
>= 0; m
--)
1625 shift
= shift_center
[m
];
1628 shift
+= shm01
*x
[i
][1] + shm02
*x
[i
][2];
1632 shift
+= shm12
*x
[i
][2];
1634 while (x
[i
][m
]-shift
< 0)
1636 for (d
= 0; d
<= m
; d
++)
1638 x
[i
][d
] += box
[m
][d
];
1641 while (x
[i
][m
]-shift
>= box
[m
][m
])
1643 for (d
= 0; d
<= m
; d
++)
1645 x
[i
][d
] -= box
[m
][d
];
1652 void put_atoms_in_compact_unitcell(int ePBC
, int ecenter
, matrix box
,
1653 int natoms
, rvec x
[])
1656 rvec box_center
, dx
;
1659 set_pbc(&pbc
, ePBC
, box
);
1661 if (pbc
.ePBCDX
== epbcdxUNSUPPORTED
)
1663 gmx_fatal(FARGS
, "Can not put atoms in compact unitcell with unsupported PBC");
1666 calc_box_center(ecenter
, box
, box_center
);
1667 for (i
= 0; i
< natoms
; i
++)
1669 pbc_dx(&pbc
, x
[i
], box_center
, dx
);
1670 rvec_add(box_center
, dx
, x
[i
]);