Bug Summary

File:gromacs/gmxlib/pbc.c
Location:line 1538, column 9
Description:Value stored to 'bFirst' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, 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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include <assert.h>
43
44#include "typedefs.h"
45#include "types/commrec.h"
46#include "gromacs/math/vec.h"
47#include "gromacs/math/utilities.h"
48#include "pbc.h"
49#include "gromacs/utility/smalloc.h"
50#include "txtdump.h"
51#include "gromacs/utility/fatalerror.h"
52#include "names.h"
53#include "macros.h"
54#include "gmx_omp_nthreads.h"
55
56/* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
57enum {
58 epbcdxRECTANGULAR = 1, epbcdxTRICLINIC,
59 epbcdx2D_RECT, epbcdx2D_TRIC,
60 epbcdx1D_RECT, epbcdx1D_TRIC,
61 epbcdxSCREW_RECT, epbcdxSCREW_TRIC,
62 epbcdxNOPBC, epbcdxUNSUPPORTED
63};
64
65/* Margin factor for error message and correction if the box is too skewed */
66#define BOX_MARGIN1.0010 1.0010
67#define BOX_MARGIN_CORRECT1.0005 1.0005
68
69int ePBC2npbcdim(int ePBC)
70{
71 int npbcdim = 0;
72
73 switch (ePBC)
74 {
75 case epbcXYZ: npbcdim = 3; break;
76 case epbcXY: npbcdim = 2; break;
77 case epbcSCREW: npbcdim = 3; break;
78 case epbcNONE: npbcdim = 0; break;
79 default: gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 79, "Unknown ePBC=%d in ePBC2npbcdim", ePBC);
80 }
81
82 return npbcdim;
83}
84
85int inputrec2nboundeddim(t_inputrec *ir)
86{
87 if (ir->nwall == 2 && ir->ePBC == epbcXY)
88 {
89 return 3;
90 }
91 else
92 {
93 return ePBC2npbcdim(ir->ePBC);
94 }
95}
96
97void dump_pbc(FILE *fp, t_pbc *pbc)
98{
99 rvec sum_box;
100
101 fprintf(fp, "ePBCDX = %d\n", pbc->ePBCDX);
102 pr_rvecs(fp, 0, "box", pbc->box, DIM3);
103 pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
104 pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
105 pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
106 rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
107 pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
108 fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
109 fprintf(fp, "bLimitDistance = %s\n", EBOOL(pbc->bLimitDistance)((((pbc->bLimitDistance) < 0) || ((pbc->bLimitDistance
) >= (2))) ? "UNDEFINED" : (bool_names)[pbc->bLimitDistance
]));
110 fprintf(fp, "limit_distance2 = %g\n", pbc->limit_distance2);
111 fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
112 if (pbc->ntric_vec > 0)
113 {
114 pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE0);
115 pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
116 }
117}
118
119const char *check_box(int ePBC, matrix box)
120{
121 const char *ptr;
122
123 if (ePBC == -1)
124 {
125 ePBC = guess_ePBC(box);
126 }
127
128 if (ePBC == epbcNONE)
129 {
130 return NULL((void*)0);
131 }
132
133 if ((box[XX0][YY1] != 0) || (box[XX0][ZZ2] != 0) || (box[YY1][ZZ2] != 0))
134 {
135 ptr = "Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.";
136 }
137 else if (ePBC == epbcSCREW && (box[YY1][XX0] != 0 || box[ZZ2][XX0] != 0))
138 {
139 ptr = "The unit cell can not have off-diagonal x-components with screw pbc";
140 }
141 else if (fabs(box[YY1][XX0]) > BOX_MARGIN1.0010*0.5*box[XX0][XX0] ||
142 (ePBC != epbcXY &&
143 (fabs(box[ZZ2][XX0]) > BOX_MARGIN1.0010*0.5*box[XX0][XX0] ||
144 fabs(box[ZZ2][YY1]) > BOX_MARGIN1.0010*0.5*box[YY1][YY1])))
145 {
146 ptr = "Triclinic box is too skewed.";
147 }
148 else
149 {
150 ptr = NULL((void*)0);
151 }
152
153 return ptr;
154}
155
156real max_cutoff2(int ePBC, matrix box)
157{
158 real min_hv2, min_ss;
159
160 /* Physical limitation of the cut-off
161 * by half the length of the shortest box vector.
162 */
163 min_hv2 = min(0.25*norm2(box[XX]), 0.25*norm2(box[YY]))(((0.25*norm2(box[0])) < (0.25*norm2(box[1]))) ? (0.25*norm2
(box[0])) : (0.25*norm2(box[1])) );
164 if (ePBC != epbcXY)
165 {
166 min_hv2 = min(min_hv2, 0.25*norm2(box[ZZ]))(((min_hv2) < (0.25*norm2(box[2]))) ? (min_hv2) : (0.25*norm2
(box[2])) );
167 }
168
169 /* Limitation to the smallest diagonal element due to optimizations:
170 * checking only linear combinations of single box-vectors (2 in x)
171 * in the grid search and pbc_dx is a lot faster
172 * than checking all possible combinations.
173 */
174 if (ePBC == epbcXY)
175 {
176 min_ss = min(box[XX][XX], box[YY][YY])(((box[0][0]) < (box[1][1])) ? (box[0][0]) : (box[1][1]) );
177 }
178 else
179 {
180 min_ss = min(box[XX][XX], min(box[YY][YY]-fabs(box[ZZ][YY]), box[ZZ][ZZ]))(((box[0][0]) < ((((box[1][1]-fabs(box[2][1])) < (box[2
][2])) ? (box[1][1]-fabs(box[2][1])) : (box[2][2]) ))) ? (box
[0][0]) : ((((box[1][1]-fabs(box[2][1])) < (box[2][2])) ? (
box[1][1]-fabs(box[2][1])) : (box[2][2]) )) );
181 }
182
183 return min(min_hv2, min_ss*min_ss)(((min_hv2) < (min_ss*min_ss)) ? (min_hv2) : (min_ss*min_ss
) );
184}
185
186/* this one is mostly harmless... */
187static gmx_bool bWarnedGuess = FALSE0;
188
189int guess_ePBC(matrix box)
190{
191 int ePBC;
192
193 if (box[XX0][XX0] > 0 && box[YY1][YY1] > 0 && box[ZZ2][ZZ2] > 0)
194 {
195 ePBC = epbcXYZ;
196 }
197 else if (box[XX0][XX0] > 0 && box[YY1][YY1] > 0 && box[ZZ2][ZZ2] == 0)
198 {
199 ePBC = epbcXY;
200 }
201 else if (box[XX0][XX0] == 0 && box[YY1][YY1] == 0 && box[ZZ2][ZZ2] == 0)
202 {
203 ePBC = epbcNONE;
204 }
205 else
206 {
207 if (!bWarnedGuess)
208 {
209 fprintf(stderrstderr, "WARNING: Unsupported box diagonal %f %f %f, "
210 "will not use periodic boundary conditions\n\n",
211 box[XX0][XX0], box[YY1][YY1], box[ZZ2][ZZ2]);
212 bWarnedGuess = TRUE1;
213 }
214 ePBC = epbcNONE;
215 }
216
217 if (debug)
218 {
219 fprintf(debug, "Guessed pbc = %s from the box matrix\n", epbc_names[ePBC]);
220 }
221
222 return ePBC;
223}
224
225static int correct_box_elem(FILE *fplog, int step, tensor box, int v, int d)
226{
227 int shift, maxshift = 10;
228
229 shift = 0;
230
231 /* correct elem d of vector v with vector d */
232 while (box[v][d] > BOX_MARGIN_CORRECT1.0005*0.5*box[d][d])
233 {
234 if (fplog)
235 {
236 fprintf(fplog, "Step %d: correcting invalid box:\n", step);
237 pr_rvecs(fplog, 0, "old box", box, DIM3);
238 }
239 rvec_dec(box[v], box[d]);
240 shift--;
241 if (fplog)
242 {
243 pr_rvecs(fplog, 0, "new box", box, DIM3);
244 }
245 if (shift <= -maxshift)
246 {
247 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 247,
248 "Box was shifted at least %d times. Please see log-file.",
249 maxshift);
250 }
251 }
252 while (box[v][d] < -BOX_MARGIN_CORRECT1.0005*0.5*box[d][d])
253 {
254 if (fplog)
255 {
256 fprintf(fplog, "Step %d: correcting invalid box:\n", step);
257 pr_rvecs(fplog, 0, "old box", box, DIM3);
258 }
259 rvec_inc(box[v], box[d]);
260 shift++;
261 if (fplog)
262 {
263 pr_rvecs(fplog, 0, "new box", box, DIM3);
264 }
265 if (shift >= maxshift)
266 {
267 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 267,
268 "Box was shifted at least %d times. Please see log-file.",
269 maxshift);
270 }
271 }
272
273 return shift;
274}
275
276gmx_bool correct_box(FILE *fplog, int step, tensor box, t_graph *graph)
277{
278 int zy, zx, yx, i;
279 gmx_bool bCorrected;
280
281 /* check if the box still obeys the restrictions, if not, correct it */
282 zy = correct_box_elem(fplog, step, box, ZZ2, YY1);
283 zx = correct_box_elem(fplog, step, box, ZZ2, XX0);
284 yx = correct_box_elem(fplog, step, box, YY1, XX0);
285
286 bCorrected = (zy || zx || yx);
287
288 if (bCorrected && graph)
289 {
290 /* correct the graph */
291 for (i = graph->at_start; i < graph->at_end; i++)
292 {
293 graph->ishift[i][YY1] -= graph->ishift[i][ZZ2]*zy;
294 graph->ishift[i][XX0] -= graph->ishift[i][ZZ2]*zx;
295 graph->ishift[i][XX0] -= graph->ishift[i][YY1]*yx;
296 }
297 }
298
299 return bCorrected;
300}
301
302int ndof_com(t_inputrec *ir)
303{
304 int n = 0;
305
306 switch (ir->ePBC)
307 {
308 case epbcXYZ:
309 case epbcNONE:
310 n = 3;
311 break;
312 case epbcXY:
313 n = (ir->nwall == 0 ? 3 : 2);
314 break;
315 case epbcSCREW:
316 n = 1;
317 break;
318 default:
319 gmx_incons("Unknown pbc in calc_nrdf")_gmx_error("incons", "Unknown pbc in calc_nrdf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c"
, 319);
320 }
321
322 return n;
323}
324
325static void low_set_pbc(t_pbc *pbc, int ePBC, ivec *dd_nc, matrix box)
326{
327 int order[5] = {0, -1, 1, -2, 2};
328 int ii, jj, kk, i, j, k, d, dd, jc, kc, npbcdim, shift;
329 ivec bPBC;
330 real d2old, d2new, d2new_c;
331 rvec trial, pos;
332 gmx_bool bXY, bUse;
333 const char *ptr;
334
335 pbc->ePBC = ePBC;
336 pbc->ndim_ePBC = ePBC2npbcdim(ePBC);
337
338 copy_mat(box, pbc->box);
339 pbc->bLimitDistance = FALSE0;
340 pbc->max_cutoff2 = 0;
341 pbc->dim = -1;
342
343 for (i = 0; (i < DIM3); i++)
344 {
345 pbc->fbox_diag[i] = box[i][i];
346 pbc->hbox_diag[i] = pbc->fbox_diag[i]*0.5;
347 pbc->mhbox_diag[i] = -pbc->hbox_diag[i];
348 }
349
350 ptr = check_box(ePBC, box);
351 if (ePBC == epbcNONE)
352 {
353 pbc->ePBCDX = epbcdxNOPBC;
354 }
355 else if (ptr)
356 {
357 fprintf(stderrstderr, "Warning: %s\n", ptr);
358 pr_rvecs(stderrstderr, 0, " Box", box, DIM3);
359 fprintf(stderrstderr, " Can not fix pbc.\n");
360 pbc->ePBCDX = epbcdxUNSUPPORTED;
361 pbc->bLimitDistance = TRUE1;
362 pbc->limit_distance2 = 0;
363 }
364 else
365 {
366 if (ePBC == epbcSCREW && dd_nc)
367 {
368 /* This combinated should never appear here */
369 gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL")_gmx_error("incons", "low_set_pbc called with screw pbc and dd_nc != NULL"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 369
);
370 }
371
372 npbcdim = 0;
373 for (i = 0; i < DIM3; i++)
374 {
375 if ((dd_nc && (*dd_nc)[i] > 1) || (ePBC == epbcXY && i == ZZ2))
376 {
377 bPBC[i] = 0;
378 }
379 else
380 {
381 bPBC[i] = 1;
382 npbcdim++;
383 }
384 }
385 switch (npbcdim)
386 {
387 case 1:
388 /* 1D pbc is not an mdp option and it is therefore only used
389 * with single shifts.
390 */
391 pbc->ePBCDX = epbcdx1D_RECT;
392 for (i = 0; i < DIM3; i++)
393 {
394 if (bPBC[i])
395 {
396 pbc->dim = i;
397 }
398 }
399 assert(pbc->dim < DIM)((void) (0));
400 for (i = 0; i < pbc->dim; i++)
401 {
402 if (pbc->box[pbc->dim][i] != 0)
403 {
404 pbc->ePBCDX = epbcdx1D_TRIC;
405 }
406 }
407 break;
408 case 2:
409 pbc->ePBCDX = epbcdx2D_RECT;
410 for (i = 0; i < DIM3; i++)
411 {
412 if (!bPBC[i])
413 {
414 pbc->dim = i;
415 }
416 }
417 for (i = 0; i < DIM3; i++)
418 {
419 if (bPBC[i])
420 {
421 for (j = 0; j < i; j++)
422 {
423 if (pbc->box[i][j] != 0)
424 {
425 pbc->ePBCDX = epbcdx2D_TRIC;
426 }
427 }
428 }
429 }
430 break;
431 case 3:
432 if (ePBC != epbcSCREW)
433 {
434 if (TRICLINIC(box)(box[1][0] != 0 || box[2][0] != 0 || box[2][1] != 0))
435 {
436 pbc->ePBCDX = epbcdxTRICLINIC;
437 }
438 else
439 {
440 pbc->ePBCDX = epbcdxRECTANGULAR;
441 }
442 }
443 else
444 {
445 pbc->ePBCDX = (box[ZZ2][YY1] == 0 ? epbcdxSCREW_RECT : epbcdxSCREW_TRIC);
446 if (pbc->ePBCDX == epbcdxSCREW_TRIC)
447 {
448 fprintf(stderrstderr,
449 "Screw pbc is not yet implemented for triclinic boxes.\n"
450 "Can not fix pbc.\n");
451 pbc->ePBCDX = epbcdxUNSUPPORTED;
452 }
453 }
454 break;
455 default:
456 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 456, "Incorrect number of pbc dimensions with DD: %d",
457 npbcdim);
458 }
459 pbc->max_cutoff2 = max_cutoff2(ePBC, box);
460
461 if (pbc->ePBCDX == epbcdxTRICLINIC ||
462 pbc->ePBCDX == epbcdx2D_TRIC ||
463 pbc->ePBCDX == epbcdxSCREW_TRIC)
464 {
465 if (debug)
466 {
467 pr_rvecs(debug, 0, "Box", box, DIM3);
468 fprintf(debug, "max cutoff %.3f\n", sqrt(pbc->max_cutoff2));
469 }
470 pbc->ntric_vec = 0;
471 /* We will only use single shifts, but we will check a few
472 * more shifts to see if there is a limiting distance
473 * above which we can not be sure of the correct distance.
474 */
475 for (kk = 0; kk < 5; kk++)
476 {
477 k = order[kk];
478 if (!bPBC[ZZ2] && k != 0)
479 {
480 continue;
481 }
482 for (jj = 0; jj < 5; jj++)
483 {
484 j = order[jj];
485 if (!bPBC[YY1] && j != 0)
486 {
487 continue;
488 }
489 for (ii = 0; ii < 3; ii++)
490 {
491 i = order[ii];
492 if (!bPBC[XX0] && i != 0)
493 {
494 continue;
495 }
496 /* A shift is only useful when it is trilinic */
497 if (j != 0 || k != 0)
498 {
499 d2old = 0;
500 d2new = 0;
501 for (d = 0; d < DIM3; d++)
502 {
503 trial[d] = i*box[XX0][d] + j*box[YY1][d] + k*box[ZZ2][d];
504 /* Choose the vector within the brick around 0,0,0 that
505 * will become the shortest due to shift try.
506 */
507 if (d == pbc->dim)
508 {
509 trial[d] = 0;
510 pos[d] = 0;
511 }
512 else
513 {
514 if (trial[d] < 0)
515 {
516 pos[d] = min( pbc->hbox_diag[d], -trial[d])(((pbc->hbox_diag[d]) < (-trial[d])) ? (pbc->hbox_diag
[d]) : (-trial[d]) );
517 }
518 else
519 {
520 pos[d] = max(-pbc->hbox_diag[d], -trial[d])(((-pbc->hbox_diag[d]) > (-trial[d])) ? (-pbc->hbox_diag
[d]) : (-trial[d]) );
521 }
522 }
523 d2old += sqr(pos[d]);
524 d2new += sqr(pos[d] + trial[d]);
525 }
526 if (BOX_MARGIN1.0010*d2new < d2old)
527 {
528 if (j < -1 || j > 1 || k < -1 || k > 1)
529 {
530 /* Check if there is a single shift vector
531 * that decreases this distance even more.
532 */
533 jc = 0;
534 kc = 0;
535 if (j < -1 || j > 1)
536 {
537 jc = j/2;
538 }
539 if (k < -1 || k > 1)
540 {
541 kc = k/2;
542 }
543 d2new_c = 0;
544 for (d = 0; d < DIM3; d++)
545 {
546 d2new_c += sqr(pos[d] + trial[d]
547 - jc*box[YY1][d] - kc*box[ZZ2][d]);
548 }
549 if (d2new_c > BOX_MARGIN1.0010*d2new)
550 {
551 /* Reject this shift vector, as there is no a priori limit
552 * to the number of shifts that decrease distances.
553 */
554 if (!pbc->bLimitDistance || d2new < pbc->limit_distance2)
555 {
556 pbc->limit_distance2 = d2new;
557 }
558 pbc->bLimitDistance = TRUE1;
559 }
560 }
561 else
562 {
563 /* Check if shifts with one box vector less do better */
564 bUse = TRUE1;
565 for (dd = 0; dd < DIM3; dd++)
566 {
567 shift = (dd == 0 ? i : (dd == 1 ? j : k));
568 if (shift)
569 {
570 d2new_c = 0;
571 for (d = 0; d < DIM3; d++)
572 {
573 d2new_c += sqr(pos[d] + trial[d] - shift*box[dd][d]);
574 }
575 if (d2new_c <= BOX_MARGIN1.0010*d2new)
576 {
577 bUse = FALSE0;
578 }
579 }
580 }
581 if (bUse)
582 {
583 /* Accept this shift vector. */
584 if (pbc->ntric_vec >= MAX_NTRICVEC12)
585 {
586 fprintf(stderrstderr, "\nWARNING: Found more than %d triclinic correction vectors, ignoring some.\n"
587 " There is probably something wrong with your box.\n", MAX_NTRICVEC12);
588 pr_rvecs(stderrstderr, 0, " Box", box, DIM3);
589 }
590 else
591 {
592 copy_rvec(trial, pbc->tric_vec[pbc->ntric_vec]);
593 pbc->tric_shift[pbc->ntric_vec][XX0] = i;
594 pbc->tric_shift[pbc->ntric_vec][YY1] = j;
595 pbc->tric_shift[pbc->ntric_vec][ZZ2] = k;
596 pbc->ntric_vec++;
597 }
598 }
599 }
600 if (debug)
601 {
602 fprintf(debug, " tricvec %2d = %2d %2d %2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
603 pbc->ntric_vec, i, j, k,
604 sqrt(d2old), sqrt(d2new),
605 trial[XX0], trial[YY1], trial[ZZ2],
606 pos[XX0], pos[YY1], pos[ZZ2]);
607 }
608 }
609 }
610 }
611 }
612 }
613 }
614 }
615}
616
617void set_pbc(t_pbc *pbc, int ePBC, matrix box)
618{
619 if (ePBC == -1)
620 {
621 ePBC = guess_ePBC(box);
622 }
623
624 low_set_pbc(pbc, ePBC, NULL((void*)0), box);
625}
626
627t_pbc *set_pbc_dd(t_pbc *pbc, int ePBC,
628 gmx_domdec_t *dd, gmx_bool bSingleDir, matrix box)
629{
630 ivec nc2;
631 int npbcdim, i;
632
633 if (dd == NULL((void*)0))
634 {
635 npbcdim = DIM3;
636 }
637 else
638 {
639 if (ePBC == epbcSCREW && dd->nc[XX0] > 1)
640 {
641 /* The rotation has been taken care of during coordinate communication */
642 ePBC = epbcXYZ;
643 }
644 npbcdim = 0;
645 for (i = 0; i < DIM3; i++)
646 {
647 if (dd->nc[i] <= (bSingleDir ? 1 : 2))
648 {
649 nc2[i] = 1;
650 if (!(ePBC == epbcXY && i == ZZ2))
651 {
652 npbcdim++;
653 }
654 }
655 else
656 {
657 nc2[i] = dd->nc[i];
658 }
659 }
660 }
661
662 if (npbcdim > 0)
663 {
664 low_set_pbc(pbc, ePBC, npbcdim < DIM3 ? &nc2 : NULL((void*)0), box);
665 }
666
667 return (npbcdim > 0 ? pbc : NULL((void*)0));
668}
669
670void pbc_dx(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
671{
672 int i, j;
673 rvec dx_start, trial;
674 real d2min, d2trial;
675 gmx_bool bRot;
676
677 rvec_sub(x1, x2, dx);
678
679 switch (pbc->ePBCDX)
680 {
681 case epbcdxRECTANGULAR:
682 for (i = 0; i < DIM3; i++)
683 {
684 while (dx[i] > pbc->hbox_diag[i])
685 {
686 dx[i] -= pbc->fbox_diag[i];
687 }
688 while (dx[i] <= pbc->mhbox_diag[i])
689 {
690 dx[i] += pbc->fbox_diag[i];
691 }
692 }
693 break;
694 case epbcdxTRICLINIC:
695 for (i = DIM3-1; i >= 0; i--)
696 {
697 while (dx[i] > pbc->hbox_diag[i])
698 {
699 for (j = i; j >= 0; j--)
700 {
701 dx[j] -= pbc->box[i][j];
702 }
703 }
704 while (dx[i] <= pbc->mhbox_diag[i])
705 {
706 for (j = i; j >= 0; j--)
707 {
708 dx[j] += pbc->box[i][j];
709 }
710 }
711 }
712 /* dx is the distance in a rectangular box */
713 d2min = norm2(dx);
714 if (d2min > pbc->max_cutoff2)
715 {
716 copy_rvec(dx, dx_start);
717 d2min = norm2(dx);
718 /* Now try all possible shifts, when the distance is within max_cutoff
719 * it must be the shortest possible distance.
720 */
721 i = 0;
722 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
723 {
724 rvec_add(dx_start, pbc->tric_vec[i], trial);
725 d2trial = norm2(trial);
726 if (d2trial < d2min)
727 {
728 copy_rvec(trial, dx);
729 d2min = d2trial;
730 }
731 i++;
732 }
733 }
734 break;
735 case epbcdx2D_RECT:
736 for (i = 0; i < DIM3; i++)
737 {
738 if (i != pbc->dim)
739 {
740 while (dx[i] > pbc->hbox_diag[i])
741 {
742 dx[i] -= pbc->fbox_diag[i];
743 }
744 while (dx[i] <= pbc->mhbox_diag[i])
745 {
746 dx[i] += pbc->fbox_diag[i];
747 }
748 }
749 }
750 break;
751 case epbcdx2D_TRIC:
752 d2min = 0;
753 for (i = DIM3-1; i >= 0; i--)
754 {
755 if (i != pbc->dim)
756 {
757 while (dx[i] > pbc->hbox_diag[i])
758 {
759 for (j = i; j >= 0; j--)
760 {
761 dx[j] -= pbc->box[i][j];
762 }
763 }
764 while (dx[i] <= pbc->mhbox_diag[i])
765 {
766 for (j = i; j >= 0; j--)
767 {
768 dx[j] += pbc->box[i][j];
769 }
770 }
771 d2min += dx[i]*dx[i];
772 }
773 }
774 if (d2min > pbc->max_cutoff2)
775 {
776 copy_rvec(dx, dx_start);
777 d2min = norm2(dx);
778 /* Now try all possible shifts, when the distance is within max_cutoff
779 * it must be the shortest possible distance.
780 */
781 i = 0;
782 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
783 {
784 rvec_add(dx_start, pbc->tric_vec[i], trial);
785 d2trial = 0;
786 for (j = 0; j < DIM3; j++)
787 {
788 if (j != pbc->dim)
789 {
790 d2trial += trial[j]*trial[j];
791 }
792 }
793 if (d2trial < d2min)
794 {
795 copy_rvec(trial, dx);
796 d2min = d2trial;
797 }
798 i++;
799 }
800 }
801 break;
802 case epbcdxSCREW_RECT:
803 /* The shift definition requires x first */
804 bRot = FALSE0;
805 while (dx[XX0] > pbc->hbox_diag[XX0])
806 {
807 dx[XX0] -= pbc->fbox_diag[XX0];
808 bRot = !bRot;
809 }
810 while (dx[XX0] <= pbc->mhbox_diag[XX0])
811 {
812 dx[XX0] += pbc->fbox_diag[YY1];
813 bRot = !bRot;
814 }
815 if (bRot)
816 {
817 /* Rotate around the x-axis in the middle of the box */
818 dx[YY1] = pbc->box[YY1][YY1] - x1[YY1] - x2[YY1];
819 dx[ZZ2] = pbc->box[ZZ2][ZZ2] - x1[ZZ2] - x2[ZZ2];
820 }
821 /* Normal pbc for y and z */
822 for (i = YY1; i <= ZZ2; i++)
823 {
824 while (dx[i] > pbc->hbox_diag[i])
825 {
826 dx[i] -= pbc->fbox_diag[i];
827 }
828 while (dx[i] <= pbc->mhbox_diag[i])
829 {
830 dx[i] += pbc->fbox_diag[i];
831 }
832 }
833 break;
834 case epbcdxNOPBC:
835 case epbcdxUNSUPPORTED:
836 break;
837 default:
838 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 838, "Internal error in pbc_dx, set_pbc has not been called");
839 break;
840 }
841}
842
843int pbc_dx_aiuc(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
844{
845 int i, j, is;
846 rvec dx_start, trial;
847 real d2min, d2trial;
848 ivec ishift, ishift_start;
849
850 rvec_sub(x1, x2, dx);
851 clear_ivec(ishift);
852
853 switch (pbc->ePBCDX)
854 {
855 case epbcdxRECTANGULAR:
856 for (i = 0; i < DIM3; i++)
857 {
858 if (dx[i] > pbc->hbox_diag[i])
859 {
860 dx[i] -= pbc->fbox_diag[i];
861 ishift[i]--;
862 }
863 else if (dx[i] <= pbc->mhbox_diag[i])
864 {
865 dx[i] += pbc->fbox_diag[i];
866 ishift[i]++;
867 }
868 }
869 break;
870 case epbcdxTRICLINIC:
871 /* For triclinic boxes the performance difference between
872 * if/else and two while loops is negligible.
873 * However, the while version can cause extreme delays
874 * before a simulation crashes due to large forces which
875 * can cause unlimited displacements.
876 * Also allowing multiple shifts would index fshift beyond bounds.
877 */
878 for (i = DIM3-1; i >= 1; i--)
879 {
880 if (dx[i] > pbc->hbox_diag[i])
881 {
882 for (j = i; j >= 0; j--)
883 {
884 dx[j] -= pbc->box[i][j];
885 }
886 ishift[i]--;
887 }
888 else if (dx[i] <= pbc->mhbox_diag[i])
889 {
890 for (j = i; j >= 0; j--)
891 {
892 dx[j] += pbc->box[i][j];
893 }
894 ishift[i]++;
895 }
896 }
897 /* Allow 2 shifts in x */
898 if (dx[XX0] > pbc->hbox_diag[XX0])
899 {
900 dx[XX0] -= pbc->fbox_diag[XX0];
901 ishift[XX0]--;
902 if (dx[XX0] > pbc->hbox_diag[XX0])
903 {
904 dx[XX0] -= pbc->fbox_diag[XX0];
905 ishift[XX0]--;
906 }
907 }
908 else if (dx[XX0] <= pbc->mhbox_diag[XX0])
909 {
910 dx[XX0] += pbc->fbox_diag[XX0];
911 ishift[XX0]++;
912 if (dx[XX0] <= pbc->mhbox_diag[XX0])
913 {
914 dx[XX0] += pbc->fbox_diag[XX0];
915 ishift[XX0]++;
916 }
917 }
918 /* dx is the distance in a rectangular box */
919 d2min = norm2(dx);
920 if (d2min > pbc->max_cutoff2)
921 {
922 copy_rvec(dx, dx_start);
923 copy_ivec(ishift, ishift_start);
924 d2min = norm2(dx);
925 /* Now try all possible shifts, when the distance is within max_cutoff
926 * it must be the shortest possible distance.
927 */
928 i = 0;
929 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
930 {
931 rvec_add(dx_start, pbc->tric_vec[i], trial);
932 d2trial = norm2(trial);
933 if (d2trial < d2min)
934 {
935 copy_rvec(trial, dx);
936 ivec_add(ishift_start, pbc->tric_shift[i], ishift);
937 d2min = d2trial;
938 }
939 i++;
940 }
941 }
942 break;
943 case epbcdx2D_RECT:
944 for (i = 0; i < DIM3; i++)
945 {
946 if (i != pbc->dim)
947 {
948 if (dx[i] > pbc->hbox_diag[i])
949 {
950 dx[i] -= pbc->fbox_diag[i];
951 ishift[i]--;
952 }
953 else if (dx[i] <= pbc->mhbox_diag[i])
954 {
955 dx[i] += pbc->fbox_diag[i];
956 ishift[i]++;
957 }
958 }
959 }
960 break;
961 case epbcdx2D_TRIC:
962 d2min = 0;
963 for (i = DIM3-1; i >= 1; i--)
964 {
965 if (i != pbc->dim)
966 {
967 if (dx[i] > pbc->hbox_diag[i])
968 {
969 for (j = i; j >= 0; j--)
970 {
971 dx[j] -= pbc->box[i][j];
972 }
973 ishift[i]--;
974 }
975 else if (dx[i] <= pbc->mhbox_diag[i])
976 {
977 for (j = i; j >= 0; j--)
978 {
979 dx[j] += pbc->box[i][j];
980 }
981 ishift[i]++;
982 }
983 d2min += dx[i]*dx[i];
984 }
985 }
986 if (pbc->dim != XX0)
987 {
988 /* Allow 2 shifts in x */
989 if (dx[XX0] > pbc->hbox_diag[XX0])
990 {
991 dx[XX0] -= pbc->fbox_diag[XX0];
992 ishift[XX0]--;
993 if (dx[XX0] > pbc->hbox_diag[XX0])
994 {
995 dx[XX0] -= pbc->fbox_diag[XX0];
996 ishift[XX0]--;
997 }
998 }
999 else if (dx[XX0] <= pbc->mhbox_diag[XX0])
1000 {
1001 dx[XX0] += pbc->fbox_diag[XX0];
1002 ishift[XX0]++;
1003 if (dx[XX0] <= pbc->mhbox_diag[XX0])
1004 {
1005 dx[XX0] += pbc->fbox_diag[XX0];
1006 ishift[XX0]++;
1007 }
1008 }
1009 d2min += dx[XX0]*dx[XX0];
1010 }
1011 if (d2min > pbc->max_cutoff2)
1012 {
1013 copy_rvec(dx, dx_start);
1014 copy_ivec(ishift, ishift_start);
1015 /* Now try all possible shifts, when the distance is within max_cutoff
1016 * it must be the shortest possible distance.
1017 */
1018 i = 0;
1019 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
1020 {
1021 rvec_add(dx_start, pbc->tric_vec[i], trial);
1022 d2trial = 0;
1023 for (j = 0; j < DIM3; j++)
1024 {
1025 if (j != pbc->dim)
1026 {
1027 d2trial += trial[j]*trial[j];
1028 }
1029 }
1030 if (d2trial < d2min)
1031 {
1032 copy_rvec(trial, dx);
1033 ivec_add(ishift_start, pbc->tric_shift[i], ishift);
1034 d2min = d2trial;
1035 }
1036 i++;
1037 }
1038 }
1039 break;
1040 case epbcdx1D_RECT:
1041 i = pbc->dim;
1042 if (dx[i] > pbc->hbox_diag[i])
1043 {
1044 dx[i] -= pbc->fbox_diag[i];
1045 ishift[i]--;
1046 }
1047 else if (dx[i] <= pbc->mhbox_diag[i])
1048 {
1049 dx[i] += pbc->fbox_diag[i];
1050 ishift[i]++;
1051 }
1052 break;
1053 case epbcdx1D_TRIC:
1054 i = pbc->dim;
1055 if (dx[i] > pbc->hbox_diag[i])
1056 {
1057 rvec_dec(dx, pbc->box[i]);
1058 ishift[i]--;
1059 }
1060 else if (dx[i] <= pbc->mhbox_diag[i])
1061 {
1062 rvec_inc(dx, pbc->box[i]);
1063 ishift[i]++;
1064 }
1065 break;
1066 case epbcdxSCREW_RECT:
1067 /* The shift definition requires x first */
1068 if (dx[XX0] > pbc->hbox_diag[XX0])
1069 {
1070 dx[XX0] -= pbc->fbox_diag[XX0];
1071 ishift[XX0]--;
1072 }
1073 else if (dx[XX0] <= pbc->mhbox_diag[XX0])
1074 {
1075 dx[XX0] += pbc->fbox_diag[XX0];
1076 ishift[XX0]++;
1077 }
1078 if (ishift[XX0] == 1 || ishift[XX0] == -1)
1079 {
1080 /* Rotate around the x-axis in the middle of the box */
1081 dx[YY1] = pbc->box[YY1][YY1] - x1[YY1] - x2[YY1];
1082 dx[ZZ2] = pbc->box[ZZ2][ZZ2] - x1[ZZ2] - x2[ZZ2];
1083 }
1084 /* Normal pbc for y and z */
1085 for (i = YY1; i <= ZZ2; i++)
1086 {
1087 if (dx[i] > pbc->hbox_diag[i])
1088 {
1089 dx[i] -= pbc->fbox_diag[i];
1090 ishift[i]--;
1091 }
1092 else if (dx[i] <= pbc->mhbox_diag[i])
1093 {
1094 dx[i] += pbc->fbox_diag[i];
1095 ishift[i]++;
1096 }
1097 }
1098 break;
1099 case epbcdxNOPBC:
1100 case epbcdxUNSUPPORTED:
1101 break;
1102 default:
1103 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 1103, "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
1104 break;
1105 }
1106
1107 is = IVEC2IS(ishift)(((2*2 +1)*((2*1 +1)*(((ishift)[2])+1)+((ishift)[1])+1)+((ishift
)[0])+2));
1108 if (debug)
1109 {
1110 range_check_mesg(is, 0, SHIFTS, "PBC shift vector index range check.")_range_check(is, 0, ((2*1 +1)*(2*1 +1)*(2*2 +1)), "PBC shift vector index range check."
,"is", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c"
, 1110);
1111 }
1112
1113 return is;
1114}
1115
1116void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx)
1117{
1118 int i, j;
1119 dvec dx_start, trial;
1120 double d2min, d2trial;
1121 gmx_bool bRot;
1122
1123 dvec_sub(x1, x2, dx);
1124
1125 switch (pbc->ePBCDX)
1126 {
1127 case epbcdxRECTANGULAR:
1128 case epbcdx2D_RECT:
1129 for (i = 0; i < DIM3; i++)
1130 {
1131 if (i != pbc->dim)
1132 {
1133 while (dx[i] > pbc->hbox_diag[i])
1134 {
1135 dx[i] -= pbc->fbox_diag[i];
1136 }
1137 while (dx[i] <= pbc->mhbox_diag[i])
1138 {
1139 dx[i] += pbc->fbox_diag[i];
1140 }
1141 }
1142 }
1143 break;
1144 case epbcdxTRICLINIC:
1145 case epbcdx2D_TRIC:
1146 d2min = 0;
1147 for (i = DIM3-1; i >= 0; i--)
1148 {
1149 if (i != pbc->dim)
1150 {
1151 while (dx[i] > pbc->hbox_diag[i])
1152 {
1153 for (j = i; j >= 0; j--)
1154 {
1155 dx[j] -= pbc->box[i][j];
1156 }
1157 }
1158 while (dx[i] <= pbc->mhbox_diag[i])
1159 {
1160 for (j = i; j >= 0; j--)
1161 {
1162 dx[j] += pbc->box[i][j];
1163 }
1164 }
1165 d2min += dx[i]*dx[i];
1166 }
1167 }
1168 if (d2min > pbc->max_cutoff2)
1169 {
1170 copy_dvec(dx, dx_start);
1171 /* Now try all possible shifts, when the distance is within max_cutoff
1172 * it must be the shortest possible distance.
1173 */
1174 i = 0;
1175 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
1176 {
1177 for (j = 0; j < DIM3; j++)
1178 {
1179 trial[j] = dx_start[j] + pbc->tric_vec[i][j];
1180 }
1181 d2trial = 0;
1182 for (j = 0; j < DIM3; j++)
1183 {
1184 if (j != pbc->dim)
1185 {
1186 d2trial += trial[j]*trial[j];
1187 }
1188 }
1189 if (d2trial < d2min)
1190 {
1191 copy_dvec(trial, dx);
1192 d2min = d2trial;
1193 }
1194 i++;
1195 }
1196 }
1197 break;
1198 case epbcdxSCREW_RECT:
1199 /* The shift definition requires x first */
1200 bRot = FALSE0;
1201 while (dx[XX0] > pbc->hbox_diag[XX0])
1202 {
1203 dx[XX0] -= pbc->fbox_diag[XX0];
1204 bRot = !bRot;
1205 }
1206 while (dx[XX0] <= pbc->mhbox_diag[XX0])
1207 {
1208 dx[XX0] += pbc->fbox_diag[YY1];
1209 bRot = !bRot;
1210 }
1211 if (bRot)
1212 {
1213 /* Rotate around the x-axis in the middle of the box */
1214 dx[YY1] = pbc->box[YY1][YY1] - x1[YY1] - x2[YY1];
1215 dx[ZZ2] = pbc->box[ZZ2][ZZ2] - x1[ZZ2] - x2[ZZ2];
1216 }
1217 /* Normal pbc for y and z */
1218 for (i = YY1; i <= ZZ2; i++)
1219 {
1220 while (dx[i] > pbc->hbox_diag[i])
1221 {
1222 dx[i] -= pbc->fbox_diag[i];
1223 }
1224 while (dx[i] <= pbc->mhbox_diag[i])
1225 {
1226 dx[i] += pbc->fbox_diag[i];
1227 }
1228 }
1229 break;
1230 case epbcdxNOPBC:
1231 case epbcdxUNSUPPORTED:
1232 break;
1233 default:
1234 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 1234, "Internal error in pbc_dx, set_pbc has not been called");
1235 break;
1236 }
1237}
1238
1239gmx_bool image_rect(ivec xi, ivec xj, ivec box_size, real rlong2, int *shift, real *r2)
1240{
1241 int m, t;
1242 int dx, b, b_2;
1243 real dxr, rij2;
1244
1245 rij2 = 0.0;
1246 t = 0;
1247 for (m = 0; (m < DIM3); m++)
1248 {
1249 dx = xi[m]-xj[m];
1250 t *= DIM3;
1251 b = box_size[m];
1252 b_2 = b/2;
1253 if (dx < -b_2)
1254 {
1255 t += 2;
1256 dx += b;
1257 }
1258 else if (dx > b_2)
1259 {
1260 dx -= b;
1261 }
1262 else
1263 {
1264 t += 1;
1265 }
1266 dxr = dx;
1267 rij2 += dxr*dxr;
1268 if (rij2 >= rlong2)
1269 {
1270 return FALSE0;
1271 }
1272 }
1273
1274 *shift = t;
1275 *r2 = rij2;
1276 return TRUE1;
1277}
1278
1279gmx_bool image_cylindric(ivec xi, ivec xj, ivec box_size, real rlong2,
1280 int *shift, real *r2)
1281{
1282 int m, t;
1283 int dx, b, b_2;
1284 real dxr, rij2;
1285
1286 rij2 = 0.0;
1287 t = 0;
1288 for (m = 0; (m < DIM3); m++)
1289 {
1290 dx = xi[m]-xj[m];
1291 t *= DIM3;
1292 b = box_size[m];
1293 b_2 = b/2;
1294 if (dx < -b_2)
1295 {
1296 t += 2;
1297 dx += b;
1298 }
1299 else if (dx > b_2)
1300 {
1301 dx -= b;
1302 }
1303 else
1304 {
1305 t += 1;
1306 }
1307
1308 dxr = dx;
1309 rij2 += dxr*dxr;
1310 if (m < ZZ2)
1311 {
1312 if (rij2 >= rlong2)
1313 {
1314 return FALSE0;
1315 }
1316 }
1317 }
1318
1319 *shift = t;
1320 *r2 = rij2;
1321 return TRUE1;
1322}
1323
1324void calc_shifts(matrix box, rvec shift_vec[])
1325{
1326 int k, l, m, d, n, test;
1327
1328 n = 0;
1329 for (m = -D_BOX_Z1; m <= D_BOX_Z1; m++)
1330 {
1331 for (l = -D_BOX_Y1; l <= D_BOX_Y1; l++)
1332 {
1333 for (k = -D_BOX_X2; k <= D_BOX_X2; k++, n++)
1334 {
1335 test = XYZ2IS(k, l, m)((2*2 +1)*((2*1 +1)*((m)+1)+(l)+1)+(k)+2);
1336 if (n != test)
1337 {
1338 gmx_incons("inconsistent shift numbering")_gmx_error("incons", "inconsistent shift numbering", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c"
, 1338);
1339 }
1340 for (d = 0; d < DIM3; d++)
1341 {
1342 shift_vec[n][d] = k*box[XX0][d] + l*box[YY1][d] + m*box[ZZ2][d];
1343 }
1344 }
1345 }
1346 }
1347}
1348
1349void calc_box_center(int ecenter, matrix box, rvec box_center)
1350{
1351 int d, m;
1352
1353 clear_rvec(box_center);
1354 switch (ecenter)
1355 {
1356 case ecenterTRIC:
1357 for (m = 0; (m < DIM3); m++)
1358 {
1359 for (d = 0; d < DIM3; d++)
1360 {
1361 box_center[d] += 0.5*box[m][d];
1362 }
1363 }
1364 break;
1365 case ecenterRECT:
1366 for (d = 0; d < DIM3; d++)
1367 {
1368 box_center[d] = 0.5*box[d][d];
1369 }
1370 break;
1371 case ecenterZERO:
1372 break;
1373 default:
1374 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 1374, "Unsupported value %d for ecenter", ecenter);
1375 }
1376}
1377
1378void calc_triclinic_images(matrix box, rvec img[])
1379{
1380 int i;
1381
1382 /* Calculate 3 adjacent images in the xy-plane */
1383 copy_rvec(box[0], img[0]);
1384 copy_rvec(box[1], img[1]);
1385 if (img[1][XX0] < 0)
1386 {
1387 svmul(-1, img[1], img[1]);
1388 }
1389 rvec_sub(img[1], img[0], img[2]);
1390
1391 /* Get the next 3 in the xy-plane as mirror images */
1392 for (i = 0; i < 3; i++)
1393 {
1394 svmul(-1, img[i], img[3+i]);
1395 }
1396
1397 /* Calculate the first 4 out of xy-plane images */
1398 copy_rvec(box[2], img[6]);
1399 if (img[6][XX0] < 0)
1400 {
1401 svmul(-1, img[6], img[6]);
1402 }
1403 for (i = 0; i < 3; i++)
1404 {
1405 rvec_add(img[6], img[i+1], img[7+i]);
1406 }
1407
1408 /* Mirror the last 4 from the previous in opposite rotation */
1409 for (i = 0; i < 4; i++)
1410 {
1411 svmul(-1, img[6 + (2+i) % 4], img[10+i]);
1412 }
1413}
1414
1415void calc_compact_unitcell_vertices(int ecenter, matrix box, rvec vert[])
1416{
1417 rvec img[NTRICIMG14], box_center;
1418 int n, i, j, tmp[4], d;
1419
1420 calc_triclinic_images(box, img);
1421
1422 n = 0;
1423 for (i = 2; i <= 5; i += 3)
1424 {
1425 tmp[0] = i-1;
1426 if (i == 2)
1427 {
1428 tmp[1] = 8;
1429 }
1430 else
1431 {
1432 tmp[1] = 6;
1433 }
1434 tmp[2] = (i+1) % 6;
1435 tmp[3] = tmp[1]+4;
1436 for (j = 0; j < 4; j++)
1437 {
1438 for (d = 0; d < DIM3; d++)
1439 {
1440 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1441 }
1442 n++;
1443 }
1444 }
1445 for (i = 7; i <= 13; i += 6)
1446 {
1447 tmp[0] = (i-7)/2;
1448 tmp[1] = tmp[0]+1;
1449 if (i == 7)
1450 {
1451 tmp[2] = 8;
1452 }
1453 else
1454 {
1455 tmp[2] = 10;
1456 }
1457 tmp[3] = i-1;
1458 for (j = 0; j < 4; j++)
1459 {
1460 for (d = 0; d < DIM3; d++)
1461 {
1462 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1463 }
1464 n++;
1465 }
1466 }
1467 for (i = 9; i <= 11; i += 2)
1468 {
1469 if (i == 9)
1470 {
1471 tmp[0] = 3;
1472 }
1473 else
1474 {
1475 tmp[0] = 0;
1476 }
1477 tmp[1] = tmp[0]+1;
1478 if (i == 9)
1479 {
1480 tmp[2] = 6;
1481 }
1482 else
1483 {
1484 tmp[2] = 12;
1485 }
1486 tmp[3] = i-1;
1487 for (j = 0; j < 4; j++)
1488 {
1489 for (d = 0; d < DIM3; d++)
1490 {
1491 vert[n][d] = img[i][d]+img[tmp[j]][d]+img[tmp[(j+1)%4]][d];
1492 }
1493 n++;
1494 }
1495 }
1496
1497 calc_box_center(ecenter, box, box_center);
1498 for (i = 0; i < NCUCVERT24; i++)
1499 {
1500 for (d = 0; d < DIM3; d++)
1501 {
1502 vert[i][d] = vert[i][d]*0.25+box_center[d];
1503 }
1504 }
1505}
1506
1507int *compact_unitcell_edges()
1508{
1509 /* this is an index in vert[] (see calc_box_vertices) */
1510 /*static int edge[NCUCEDGE*2];*/
1511 int *edge;
1512 static const int hexcon[24] = {
1513 0, 9, 1, 19, 2, 15, 3, 21,
1514 4, 17, 5, 11, 6, 23, 7, 13,
1515 8, 20, 10, 18, 12, 16, 14, 22
1516 };
1517 int e, i, j;
1518 gmx_bool bFirst = TRUE1;
1519
1520 snew(edge, NCUCEDGE*2)(edge) = save_calloc("edge", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c"
, 1520, (36*2), sizeof(*(edge)));
1521
1522 if (bFirst)
1523 {
1524 e = 0;
1525 for (i = 0; i < 6; i++)
1526 {
1527 for (j = 0; j < 4; j++)
1528 {
1529 edge[e++] = 4*i + j;
1530 edge[e++] = 4*i + (j+1) % 4;
1531 }
1532 }
1533 for (i = 0; i < 12*2; i++)
1534 {
1535 edge[e++] = hexcon[i];
1536 }
1537
1538 bFirst = FALSE0;
Value stored to 'bFirst' is never read
1539 }
1540
1541 return edge;
1542}
1543
1544void put_atoms_in_box_omp(int ePBC, matrix box, int natoms, rvec x[])
1545{
1546 int t, nth;
1547 nth = gmx_omp_nthreads_get(emntDefault);
1548
1549#pragma omp parallel for num_threads(nth) schedule(static)
1550 for (t = 0; t < nth; t++)
1551 {
1552 int offset, len;
1553
1554 offset = (natoms*t )/nth;
1555 len = (natoms*(t + 1))/nth - offset;
1556 put_atoms_in_box(ePBC, box, len, x + offset);
1557 }
1558}
1559
1560void put_atoms_in_box(int ePBC, matrix box, int natoms, rvec x[])
1561{
1562 int npbcdim, i, m, d;
1563
1564 if (ePBC == epbcSCREW)
1565 {
1566 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/pbc.c", 1566, "Sorry, %s pbc is not yet supported", epbc_names[ePBC]);
1567 }
1568
1569 if (ePBC == epbcXY)
1570 {
1571 npbcdim = 2;
1572 }
1573 else
1574 {
1575 npbcdim = 3;
1576 }
1577
1578 if (TRICLINIC(box)(box[1][0] != 0 || box[2][0] != 0 || box[2][1] != 0))
1579 {
1580 for (i = 0; (i < natoms); i++)
1581 {
1582 for (m = npbcdim-1; m >= 0; m--)
1583 {
1584 while (x[i][m] < 0)
1585 {
1586 for (d = 0; d <= m; d++)
1587 {
1588 x[i][d] += box[m][d];
1589 }
1590 }
1591 while (x[i][m] >= box[m][m])
1592 {
1593 for (d = 0; d <= m; d++)
1594 {
1595 x[i][d] -= box[m][d];
1596 }
1597 }
1598 }
1599 }
1600 }
1601 else
1602 {
1603 for (i = 0; i < natoms; i++)
1604 {
1605 for (d = 0; d < npbcdim; d++)
1606 {
1607 while (x[i][d] < 0)
1608 {
1609 x[i][d] += box[d][d];
1610 }
1611 while (x[i][d] >= box[d][d])
1612 {
1613 x[i][d] -= box[d][d];
1614 }
1615 }
1616 }
1617 }
1618}
1619
1620void put_atoms_in_triclinic_unitcell(int ecenter, matrix box,
1621 int natoms, rvec x[])
1622{
1623 rvec box_center, shift_center;
1624 real shm01, shm02, shm12, shift;
1625 int i, m, d;
1626
1627 calc_box_center(ecenter, box, box_center);
1628
1629 /* The product of matrix shm with a coordinate gives the shift vector
1630 which is required determine the periodic cell position */
1631 shm01 = box[1][0]/box[1][1];
1632 shm02 = (box[1][1]*box[2][0] - box[2][1]*box[1][0])/(box[1][1]*box[2][2]);
1633 shm12 = box[2][1]/box[2][2];
1634
1635 clear_rvec(shift_center);
1636 for (d = 0; d < DIM3; d++)
1637 {
1638 rvec_inc(shift_center, box[d]);
1639 }
1640 svmul(0.5, shift_center, shift_center);
1641 rvec_sub(box_center, shift_center, shift_center);
1642
1643 shift_center[0] = shm01*shift_center[1] + shm02*shift_center[2];
1644 shift_center[1] = shm12*shift_center[2];
1645 shift_center[2] = 0;
1646
1647 for (i = 0; (i < natoms); i++)
1648 {
1649 for (m = DIM3-1; m >= 0; m--)
1650 {
1651 shift = shift_center[m];
1652 if (m == 0)
1653 {
1654 shift += shm01*x[i][1] + shm02*x[i][2];
1655 }
1656 else if (m == 1)
1657 {
1658 shift += shm12*x[i][2];
1659 }
1660 while (x[i][m]-shift < 0)
1661 {
1662 for (d = 0; d <= m; d++)
1663 {
1664 x[i][d] += box[m][d];
1665 }
1666 }
1667 while (x[i][m]-shift >= box[m][m])
1668 {
1669 for (d = 0; d <= m; d++)
1670 {
1671 x[i][d] -= box[m][d];
1672 }
1673 }
1674 }
1675 }
1676}
1677
1678const char *
1679put_atoms_in_compact_unitcell(int ePBC, int ecenter, matrix box,
1680 int natoms, rvec x[])
1681{
1682 t_pbc pbc;
1683 rvec box_center, dx;
1684 int i;
1685
1686 set_pbc(&pbc, ePBC, box);
1687 calc_box_center(ecenter, box, box_center);
1688 for (i = 0; i < natoms; i++)
1689 {
1690 pbc_dx(&pbc, x[i], box_center, dx);
1691 rvec_add(box_center, dx, x[i]);
1692 }
1693
1694 return pbc.bLimitDistance ?
1695 "WARNING: Could not put all atoms in the compact unitcell\n"
1696 : NULL((void*)0);
1697}