1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROningen Mixture of Alchemy and Childrens' Stories
49 #include "gmx_fatal.h"
52 /* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
53 enum { epbcdxRECTANGULAR
=1, epbcdxTRICLINIC
,
54 epbcdx2D_RECT
, epbcdx2D_TRIC
,
55 epbcdx1D_RECT
, epbcdx1D_TRIC
,
56 epbcdxSCREW_RECT
, epbcdxSCREW_TRIC
,
57 epbcdxNOPBC
, epbcdxUNSUPPORTED
};
59 /* Margin factor for error message and correction if the box is too skewed */
60 #define BOX_MARGIN 1.0010
61 #define BOX_MARGIN_CORRECT 1.0005
63 int ePBC2npbcdim(int ePBC
)
68 case epbcXYZ
: npbcdim
= 3; break;
69 case epbcXY
: npbcdim
= 2; break;
70 case epbcSCREW
: npbcdim
= 3; break;
71 case epbcNONE
: npbcdim
= 0; break;
72 default: gmx_fatal(FARGS
,"Unknown ePBC=%d in ePBC2npbcdim",ePBC
);
78 int inputrec2nboundeddim(t_inputrec
*ir
)
80 if (ir
->nwall
== 2 && ir
->ePBC
== epbcXY
)
86 return ePBC2npbcdim(ir
->ePBC
);
90 void dump_pbc(FILE *fp
,t_pbc
*pbc
)
94 fprintf(fp
,"ePBCDX = %d\n",pbc
->ePBCDX
);
95 pr_rvecs(fp
,0,"box",pbc
->box
,DIM
);
96 pr_rvecs(fp
,0,"fbox_diag",&pbc
->fbox_diag
,1);
97 pr_rvecs(fp
,0,"hbox_diag",&pbc
->hbox_diag
,1);
98 pr_rvecs(fp
,0,"mhbox_diag",&pbc
->mhbox_diag
,1);
99 rvec_add(pbc
->hbox_diag
,pbc
->mhbox_diag
,sum_box
);
100 pr_rvecs(fp
,0,"sum of the above two",&sum_box
,1);
101 fprintf(fp
,"max_cutoff2 = %g\n",pbc
->max_cutoff2
);
102 fprintf(fp
,"bLimitDistance = %s\n",BOOL(pbc
->bLimitDistance
));
103 fprintf(fp
,"limit_distance2 = %g\n",pbc
->limit_distance2
);
104 fprintf(fp
,"ntric_vec = %d\n",pbc
->ntric_vec
);
105 if (pbc
->ntric_vec
> 0) {
106 pr_ivecs(fp
,0,"tric_shift",pbc
->tric_shift
,pbc
->ntric_vec
,FALSE
);
107 pr_rvecs(fp
,0,"tric_vec",pbc
->tric_vec
,pbc
->ntric_vec
);
111 const char *check_box(int ePBC
,matrix box
)
116 ePBC
= guess_ePBC(box
);
118 if (ePBC
== epbcNONE
)
121 if ((box
[XX
][YY
] != 0) || (box
[XX
][ZZ
] != 0) || (box
[YY
][ZZ
] != 0)) {
122 ptr
= "Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.";
123 } else if (ePBC
== epbcSCREW
&& (box
[YY
][XX
] != 0 || box
[ZZ
][XX
] != 0)) {
124 ptr
= "The unit cell can not have off-diagonal x-components with screw pbc";
125 } else if (fabs(box
[YY
][XX
]) > BOX_MARGIN
*0.5*box
[XX
][XX
] ||
127 (fabs(box
[ZZ
][XX
]) > BOX_MARGIN
*0.5*box
[XX
][XX
] ||
128 fabs(box
[ZZ
][YY
]) > BOX_MARGIN
*0.5*box
[YY
][YY
]))) {
129 ptr
= "Triclinic box is too skewed.";
137 real
max_cutoff2(int ePBC
,matrix box
)
141 /* Physical limitation of the cut-off
142 * by half the length of the shortest box vector.
144 min_hv2
= min(0.25*norm2(box
[XX
]),0.25*norm2(box
[YY
]));
146 min_hv2
= min(min_hv2
,0.25*norm2(box
[ZZ
]));
148 /* Limitation to the smallest diagonal element due to optimizations:
149 * checking only linear combinations of single box-vectors (2 in x)
150 * in the grid search and pbc_dx is a lot faster
151 * than checking all possible combinations.
153 if (ePBC
== epbcXY
) {
154 min_ss
= min(box
[XX
][XX
],box
[YY
][YY
]);
156 min_ss
= min(box
[XX
][XX
],min(box
[YY
][YY
]-fabs(box
[ZZ
][YY
]),box
[ZZ
][ZZ
]));
159 return min(min_hv2
,min_ss
*min_ss
);
162 /* this one is mostly harmless... */
163 static gmx_bool bWarnedGuess
=FALSE
;
165 int guess_ePBC(matrix box
)
169 if (box
[XX
][XX
]>0 && box
[YY
][YY
]>0 && box
[ZZ
][ZZ
]>0) {
171 } else if (box
[XX
][XX
]>0 && box
[YY
][YY
]>0 && box
[ZZ
][ZZ
]==0) {
173 } else if (box
[XX
][XX
]==0 && box
[YY
][YY
]==0 && box
[ZZ
][ZZ
]==0) {
177 fprintf(stderr
,"WARNING: Unsupported box diagonal %f %f %f, "
178 "will not use periodic boundary conditions\n\n",
179 box
[XX
][XX
],box
[YY
][YY
],box
[ZZ
][ZZ
]);
186 fprintf(debug
,"Guessed pbc = %s from the box matrix\n",epbc_names
[ePBC
]);
191 static int correct_box_elem(FILE *fplog
,int step
,tensor box
,int v
,int d
)
193 int shift
,maxshift
=10;
197 /* correct elem d of vector v with vector d */
198 while (box
[v
][d
] > BOX_MARGIN_CORRECT
*0.5*box
[d
][d
]) {
200 fprintf(fplog
,"Step %d: correcting invalid box:\n",step
);
201 pr_rvecs(fplog
,0,"old box",box
,DIM
);
203 rvec_dec(box
[v
],box
[d
]);
206 pr_rvecs(fplog
,0,"new box",box
,DIM
);
208 if (shift
<= -maxshift
)
210 "Box was shifted at least %d times. Please see log-file.",
213 while (box
[v
][d
] < -BOX_MARGIN_CORRECT
*0.5*box
[d
][d
]) {
215 fprintf(fplog
,"Step %d: correcting invalid box:\n",step
);
216 pr_rvecs(fplog
,0,"old box",box
,DIM
);
218 rvec_inc(box
[v
],box
[d
]);
221 pr_rvecs(fplog
,0,"new box",box
,DIM
);
223 if (shift
>= maxshift
)
225 "Box was shifted at least %d times. Please see log-file.",
232 gmx_bool
correct_box(FILE *fplog
,int step
,tensor box
,t_graph
*graph
)
237 /* check if the box still obeys the restrictions, if not, correct it */
238 zy
= correct_box_elem(fplog
,step
,box
,ZZ
,YY
);
239 zx
= correct_box_elem(fplog
,step
,box
,ZZ
,XX
);
240 yx
= correct_box_elem(fplog
,step
,box
,YY
,XX
);
242 bCorrected
= (zy
|| zx
|| yx
);
244 if (bCorrected
&& graph
) {
245 /* correct the graph */
246 for(i
=0; i
<graph
->nnodes
; i
++) {
247 graph
->ishift
[i
][YY
] -= graph
->ishift
[i
][ZZ
]*zy
;
248 graph
->ishift
[i
][XX
] -= graph
->ishift
[i
][ZZ
]*zx
;
249 graph
->ishift
[i
][XX
] -= graph
->ishift
[i
][YY
]*yx
;
256 int ndof_com(t_inputrec
*ir
)
266 n
= (ir
->nwall
== 0 ? 3 : 2);
272 gmx_incons("Unknown pbc in calc_nrdf");
278 static void low_set_pbc(t_pbc
*pbc
,int ePBC
,ivec
*dd_nc
,matrix box
)
280 int order
[5]={0,-1,1,-2,2};
281 int ii
,jj
,kk
,i
,j
,k
,d
,dd
,jc
,kc
,npbcdim
,shift
;
283 real d2old
,d2new
,d2new_c
;
288 pbc
->ndim_ePBC
= ePBC2npbcdim(ePBC
);
290 copy_mat(box
,pbc
->box
);
291 pbc
->bLimitDistance
= FALSE
;
292 pbc
->max_cutoff2
= 0;
295 for(i
=0; (i
<DIM
); i
++) {
296 pbc
->fbox_diag
[i
] = box
[i
][i
];
297 pbc
->hbox_diag
[i
] = pbc
->fbox_diag
[i
]*0.5;
298 pbc
->mhbox_diag
[i
] = -pbc
->hbox_diag
[i
];
301 ptr
= check_box(ePBC
,box
);
302 if (ePBC
== epbcNONE
) {
303 pbc
->ePBCDX
= epbcdxNOPBC
;
305 fprintf(stderr
, "Warning: %s\n",ptr
);
306 pr_rvecs(stderr
,0," Box",box
,DIM
);
307 fprintf(stderr
, " Can not fix pbc.\n");
308 pbc
->ePBCDX
= epbcdxUNSUPPORTED
;
309 pbc
->bLimitDistance
= TRUE
;
310 pbc
->limit_distance2
= 0;
312 if (ePBC
== epbcSCREW
&& dd_nc
) {
313 /* This combinated should never appear here */
314 gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
318 for(i
=0; i
<DIM
; i
++) {
319 if ((dd_nc
&& (*dd_nc
)[i
] > 1) || (ePBC
== epbcXY
&& i
== ZZ
)) {
328 /* 1D pbc is not an mdp option and it is therefore only used
329 * with single shifts.
331 pbc
->ePBCDX
= epbcdx1D_RECT
;
335 for(i
=0; i
<pbc
->dim
; i
++)
336 if (pbc
->box
[pbc
->dim
][i
] != 0)
337 pbc
->ePBCDX
= epbcdx1D_TRIC
;
340 pbc
->ePBCDX
= epbcdx2D_RECT
;
347 if (pbc
->box
[i
][j
] != 0)
348 pbc
->ePBCDX
= epbcdx2D_TRIC
;
351 if (ePBC
!= epbcSCREW
) {
352 if (TRICLINIC(box
)) {
353 pbc
->ePBCDX
= epbcdxTRICLINIC
;
355 pbc
->ePBCDX
= epbcdxRECTANGULAR
;
358 pbc
->ePBCDX
= (box
[ZZ
][YY
]==0 ? epbcdxSCREW_RECT
: epbcdxSCREW_TRIC
);
359 if (pbc
->ePBCDX
== epbcdxSCREW_TRIC
) {
361 "Screw pbc is not yet implemented for triclinic boxes.\n"
362 "Can not fix pbc.\n");
363 pbc
->ePBCDX
= epbcdxUNSUPPORTED
;
368 gmx_fatal(FARGS
,"Incorrect number of pbc dimensions with DD: %d",
371 pbc
->max_cutoff2
= max_cutoff2(ePBC
,box
);
373 if (pbc
->ePBCDX
== epbcdxTRICLINIC
||
374 pbc
->ePBCDX
== epbcdx2D_TRIC
||
375 pbc
->ePBCDX
== epbcdxSCREW_TRIC
) {
377 pr_rvecs(debug
,0,"Box",box
,DIM
);
378 fprintf(debug
,"max cutoff %.3f\n",sqrt(pbc
->max_cutoff2
));
381 /* We will only use single shifts, but we will check a few
382 * more shifts to see if there is a limiting distance
383 * above which we can not be sure of the correct distance.
385 for(kk
=0; kk
<5; kk
++) {
387 if (!bPBC
[ZZ
] && k
!= 0)
389 for(jj
=0; jj
<5; jj
++) {
391 if (!bPBC
[YY
] && j
!= 0)
393 for(ii
=0; ii
<3; ii
++) {
395 if (!bPBC
[XX
] && i
!= 0)
397 /* A shift is only useful when it is trilinic */
398 if (j
!= 0 || k
!= 0) {
401 for(d
=0; d
<DIM
; d
++) {
402 trial
[d
] = i
*box
[XX
][d
] + j
*box
[YY
][d
] + k
*box
[ZZ
][d
];
403 /* Choose the vector within the brick around 0,0,0 that
404 * will become the shortest due to shift try.
411 pos
[d
] = min( pbc
->hbox_diag
[d
],-trial
[d
]);
413 pos
[d
] = max(-pbc
->hbox_diag
[d
],-trial
[d
]);
415 d2old
+= sqr(pos
[d
]);
416 d2new
+= sqr(pos
[d
] + trial
[d
]);
418 if (BOX_MARGIN
*d2new
< d2old
) {
419 if (j
< -1 || j
> 1 || k
< -1 || k
> 1) {
420 /* Check if there is a single shift vector
421 * that decreases this distance even more.
431 d2new_c
+= sqr(pos
[d
] + trial
[d
]
432 - jc
*box
[YY
][d
] - kc
*box
[ZZ
][d
]);
433 if (d2new_c
> BOX_MARGIN
*d2new
) {
434 /* Reject this shift vector, as there is no a priori limit
435 * to the number of shifts that decrease distances.
437 if (!pbc
->bLimitDistance
|| d2new
< pbc
->limit_distance2
)
438 pbc
->limit_distance2
= d2new
;
439 pbc
->bLimitDistance
= TRUE
;
442 /* Check if shifts with one box vector less do better */
444 for(dd
=0; dd
<DIM
; dd
++) {
445 shift
= (dd
==0 ? i
: (dd
==1 ? j
: k
));
449 d2new_c
+= sqr(pos
[d
] + trial
[d
] - shift
*box
[dd
][d
]);
450 if (d2new_c
<= BOX_MARGIN
*d2new
)
455 /* Accept this shift vector. */
456 if (pbc
->ntric_vec
>= MAX_NTRICVEC
) {
457 fprintf(stderr
,"\nWARNING: Found more than %d triclinic correction vectors, ignoring some.\n"
458 " There is probably something wrong with your box.\n",MAX_NTRICVEC
);
459 pr_rvecs(stderr
,0," Box",box
,DIM
);
461 copy_rvec(trial
,pbc
->tric_vec
[pbc
->ntric_vec
]);
462 pbc
->tric_shift
[pbc
->ntric_vec
][XX
] = i
;
463 pbc
->tric_shift
[pbc
->ntric_vec
][YY
] = j
;
464 pbc
->tric_shift
[pbc
->ntric_vec
][ZZ
] = k
;
470 fprintf(debug
," tricvec %2d = %2d %2d %2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
471 pbc
->ntric_vec
,i
,j
,k
,
472 sqrt(d2old
),sqrt(d2new
),
473 trial
[XX
],trial
[YY
],trial
[ZZ
],
474 pos
[XX
],pos
[YY
],pos
[ZZ
]);
485 void set_pbc(t_pbc
*pbc
,int ePBC
,matrix box
)
488 ePBC
= guess_ePBC(box
);
490 low_set_pbc(pbc
,ePBC
,NULL
,box
);
493 t_pbc
*set_pbc_dd(t_pbc
*pbc
,int ePBC
,
494 gmx_domdec_t
*dd
,gmx_bool bSingleDir
,matrix box
)
502 if (ePBC
== epbcSCREW
&& dd
->nc
[XX
] > 1) {
503 /* The rotation has been taken care of during coordinate communication */
507 for(i
=0; i
<DIM
; i
++) {
508 if (dd
->nc
[i
] <= (bSingleDir
? 1 : 2)) {
510 if (!(ePBC
== epbcXY
&& i
== ZZ
))
519 low_set_pbc(pbc
,ePBC
,npbcdim
<DIM
? &nc2
: NULL
,box
);
521 return (npbcdim
> 0 ? pbc
: NULL
);
524 void pbc_dx(const t_pbc
*pbc
,const rvec x1
, const rvec x2
, rvec dx
)
533 switch (pbc
->ePBCDX
) {
534 case epbcdxRECTANGULAR
:
535 for(i
=0; i
<DIM
; i
++) {
536 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
537 dx
[i
] -= pbc
->fbox_diag
[i
];
539 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
540 dx
[i
] += pbc
->fbox_diag
[i
];
544 case epbcdxTRICLINIC
:
545 for(i
=DIM
-1; i
>=0; i
--) {
546 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
548 dx
[j
] -= pbc
->box
[i
][j
];
550 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
552 dx
[j
] += pbc
->box
[i
][j
];
555 /* dx is the distance in a rectangular box */
557 if (d2min
> pbc
->max_cutoff2
) {
558 copy_rvec(dx
,dx_start
);
560 /* Now try all possible shifts, when the distance is within max_cutoff
561 * it must be the shortest possible distance.
564 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
565 rvec_add(dx_start
,pbc
->tric_vec
[i
],trial
);
566 d2trial
= norm2(trial
);
567 if (d2trial
< d2min
) {
576 for(i
=0; i
<DIM
; i
++) {
578 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
579 dx
[i
] -= pbc
->fbox_diag
[i
];
581 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
582 dx
[i
] += pbc
->fbox_diag
[i
];
589 for(i
=DIM
-1; i
>=0; i
--) {
591 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
593 dx
[j
] -= pbc
->box
[i
][j
];
595 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
597 dx
[j
] += pbc
->box
[i
][j
];
599 d2min
+= dx
[i
]*dx
[i
];
602 if (d2min
> pbc
->max_cutoff2
) {
603 copy_rvec(dx
,dx_start
);
605 /* Now try all possible shifts, when the distance is within max_cutoff
606 * it must be the shortest possible distance.
609 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
610 rvec_add(dx_start
,pbc
->tric_vec
[i
],trial
);
612 for(j
=0; j
<DIM
; j
++) {
614 d2trial
+= trial
[j
]*trial
[j
];
617 if (d2trial
< d2min
) {
625 case epbcdxSCREW_RECT
:
626 /* The shift definition requires x first */
628 while (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
629 dx
[XX
] -= pbc
->fbox_diag
[XX
];
632 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
633 dx
[XX
] += pbc
->fbox_diag
[YY
];
637 /* Rotate around the x-axis in the middle of the box */
638 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
639 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
641 /* Normal pbc for y and z */
642 for(i
=YY
; i
<=ZZ
; i
++) {
643 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
644 dx
[i
] -= pbc
->fbox_diag
[i
];
646 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
647 dx
[i
] += pbc
->fbox_diag
[i
];
652 case epbcdxUNSUPPORTED
:
655 gmx_fatal(FARGS
,"Internal error in pbc_dx, set_pbc has not been called");
660 int pbc_dx_aiuc(const t_pbc
*pbc
,const rvec x1
, const rvec x2
, rvec dx
)
665 ivec ishift
,ishift_start
;
670 switch (pbc
->ePBCDX
) {
671 case epbcdxRECTANGULAR
:
672 for(i
=0; i
<DIM
; i
++) {
673 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
674 dx
[i
] -= pbc
->fbox_diag
[i
];
676 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
677 dx
[i
] += pbc
->fbox_diag
[i
];
682 case epbcdxTRICLINIC
:
683 /* For triclinic boxes the performance difference between
684 * if/else and two while loops is negligible.
685 * However, the while version can cause extreme delays
686 * before a simulation crashes due to large forces which
687 * can cause unlimited displacements.
688 * Also allowing multiple shifts would index fshift beyond bounds.
690 for(i
=DIM
-1; i
>=1; i
--) {
691 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
693 dx
[j
] -= pbc
->box
[i
][j
];
695 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
697 dx
[j
] += pbc
->box
[i
][j
];
701 /* Allow 2 shifts in x */
702 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
703 dx
[XX
] -= pbc
->fbox_diag
[XX
];
705 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
706 dx
[XX
] -= pbc
->fbox_diag
[XX
];
709 } else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
710 dx
[XX
] += pbc
->fbox_diag
[XX
];
712 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
713 dx
[XX
] += pbc
->fbox_diag
[XX
];
717 /* dx is the distance in a rectangular box */
719 if (d2min
> pbc
->max_cutoff2
) {
720 copy_rvec(dx
,dx_start
);
721 copy_ivec(ishift
,ishift_start
);
723 /* Now try all possible shifts, when the distance is within max_cutoff
724 * it must be the shortest possible distance.
727 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
728 rvec_add(dx_start
,pbc
->tric_vec
[i
],trial
);
729 d2trial
= norm2(trial
);
730 if (d2trial
< d2min
) {
732 ivec_add(ishift_start
,pbc
->tric_shift
[i
],ishift
);
740 for(i
=0; i
<DIM
; i
++) {
742 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
743 dx
[i
] -= pbc
->fbox_diag
[i
];
745 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
746 dx
[i
] += pbc
->fbox_diag
[i
];
754 for(i
=DIM
-1; i
>=1; i
--) {
756 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
758 dx
[j
] -= pbc
->box
[i
][j
];
760 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
762 dx
[j
] += pbc
->box
[i
][j
];
765 d2min
+= dx
[i
]*dx
[i
];
768 if (pbc
->dim
!= XX
) {
769 /* Allow 2 shifts in x */
770 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
771 dx
[XX
] -= pbc
->fbox_diag
[XX
];
773 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
774 dx
[XX
] -= pbc
->fbox_diag
[XX
];
777 } else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
778 dx
[XX
] += pbc
->fbox_diag
[XX
];
780 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
781 dx
[XX
] += pbc
->fbox_diag
[XX
];
785 d2min
+= dx
[XX
]*dx
[XX
];
787 if (d2min
> pbc
->max_cutoff2
) {
788 copy_rvec(dx
,dx_start
);
789 copy_ivec(ishift
,ishift_start
);
790 /* Now try all possible shifts, when the distance is within max_cutoff
791 * it must be the shortest possible distance.
794 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
795 rvec_add(dx_start
,pbc
->tric_vec
[i
],trial
);
797 for(j
=0; j
<DIM
; j
++) {
799 d2trial
+= trial
[j
]*trial
[j
];
802 if (d2trial
< d2min
) {
804 ivec_add(ishift_start
,pbc
->tric_shift
[i
],ishift
);
813 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
814 dx
[i
] -= pbc
->fbox_diag
[i
];
816 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
817 dx
[i
] += pbc
->fbox_diag
[i
];
823 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
824 rvec_dec(dx
,pbc
->box
[i
]);
826 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
827 rvec_inc(dx
,pbc
->box
[i
]);
831 case epbcdxSCREW_RECT
:
832 /* The shift definition requires x first */
833 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
834 dx
[XX
] -= pbc
->fbox_diag
[XX
];
836 } else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
837 dx
[XX
] += pbc
->fbox_diag
[XX
];
840 if (ishift
[XX
] == 1 || ishift
[XX
] == -1) {
841 /* Rotate around the x-axis in the middle of the box */
842 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
843 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
845 /* Normal pbc for y and z */
846 for(i
=YY
; i
<=ZZ
; i
++) {
847 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
848 dx
[i
] -= pbc
->fbox_diag
[i
];
850 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
851 dx
[i
] += pbc
->fbox_diag
[i
];
857 case epbcdxUNSUPPORTED
:
860 gmx_fatal(FARGS
,"Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
864 is
= IVEC2IS(ishift
);
867 range_check_mesg(is
,0,SHIFTS
,"PBC shift vector index range check.");
873 void pbc_dx_d(const t_pbc
*pbc
,const dvec x1
, const dvec x2
, dvec dx
)
877 double d2min
,d2trial
;
882 switch (pbc
->ePBCDX
) {
883 case epbcdxRECTANGULAR
:
885 for(i
=0; i
<DIM
; i
++) {
887 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
888 dx
[i
] -= pbc
->fbox_diag
[i
];
890 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
891 dx
[i
] += pbc
->fbox_diag
[i
];
896 case epbcdxTRICLINIC
:
899 for(i
=DIM
-1; i
>=0; i
--) {
901 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
903 dx
[j
] -= pbc
->box
[i
][j
];
905 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
907 dx
[j
] += pbc
->box
[i
][j
];
909 d2min
+= dx
[i
]*dx
[i
];
912 if (d2min
> pbc
->max_cutoff2
) {
913 copy_dvec(dx
,dx_start
);
914 /* Now try all possible shifts, when the distance is within max_cutoff
915 * it must be the shortest possible distance.
918 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
919 for(j
=0; j
<DIM
; j
++) {
920 trial
[j
] = dx_start
[j
] + pbc
->tric_vec
[i
][j
];
923 for(j
=0; j
<DIM
; j
++) {
925 d2trial
+= trial
[j
]*trial
[j
];
928 if (d2trial
< d2min
) {
936 case epbcdxSCREW_RECT
:
937 /* The shift definition requires x first */
939 while (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
940 dx
[XX
] -= pbc
->fbox_diag
[XX
];
943 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
944 dx
[XX
] += pbc
->fbox_diag
[YY
];
948 /* Rotate around the x-axis in the middle of the box */
949 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
950 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
952 /* Normal pbc for y and z */
953 for(i
=YY
; i
<=ZZ
; i
++) {
954 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
955 dx
[i
] -= pbc
->fbox_diag
[i
];
957 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
958 dx
[i
] += pbc
->fbox_diag
[i
];
963 case epbcdxUNSUPPORTED
:
966 gmx_fatal(FARGS
,"Internal error in pbc_dx, set_pbc has not been called");
971 gmx_bool
image_rect(ivec xi
,ivec xj
,ivec box_size
,real rlong2
,int *shift
,real
*r2
)
979 for(m
=0; (m
<DIM
); m
++) {
1003 gmx_bool
image_cylindric(ivec xi
,ivec xj
,ivec box_size
,real rlong2
,
1004 int *shift
,real
*r2
)
1012 for(m
=0; (m
<DIM
); m
++) {
1039 void calc_shifts(matrix box
,rvec shift_vec
[])
1044 for(m
= -D_BOX_Z
; m
<= D_BOX_Z
; m
++)
1045 for(l
= -D_BOX_Y
; l
<= D_BOX_Y
; l
++)
1046 for(k
= -D_BOX_X
; k
<= D_BOX_X
; k
++,n
++) {
1047 test
= XYZ2IS(k
,l
,m
);
1049 gmx_incons("inconsistent shift numbering");
1050 for(d
= 0; d
< DIM
; d
++)
1051 shift_vec
[n
][d
] = k
*box
[XX
][d
] + l
*box
[YY
][d
] + m
*box
[ZZ
][d
];
1055 void calc_box_center(int ecenter
,matrix box
,rvec box_center
)
1059 clear_rvec(box_center
);
1062 for(m
=0; (m
<DIM
); m
++)
1063 for(d
=0; d
<DIM
; d
++)
1064 box_center
[d
] += 0.5*box
[m
][d
];
1067 for(d
=0; d
<DIM
; d
++)
1068 box_center
[d
] = 0.5*box
[d
][d
];
1073 gmx_fatal(FARGS
,"Unsupported value %d for ecenter",ecenter
);
1077 void calc_triclinic_images(matrix box
,rvec img
[])
1081 /* Calculate 3 adjacent images in the xy-plane */
1082 copy_rvec(box
[0],img
[0]);
1083 copy_rvec(box
[1],img
[1]);
1085 svmul(-1,img
[1],img
[1]);
1086 rvec_sub(img
[1],img
[0],img
[2]);
1088 /* Get the next 3 in the xy-plane as mirror images */
1090 svmul(-1,img
[i
],img
[3+i
]);
1092 /* Calculate the first 4 out of xy-plane images */
1093 copy_rvec(box
[2],img
[6]);
1095 svmul(-1,img
[6],img
[6]);
1097 rvec_add(img
[6],img
[i
+1],img
[7+i
]);
1099 /* Mirror the last 4 from the previous in opposite rotation */
1101 svmul(-1,img
[6 + (2+i
) % 4],img
[10+i
]);
1104 void calc_compact_unitcell_vertices(int ecenter
,matrix box
,rvec vert
[])
1106 rvec img
[NTRICIMG
],box_center
;
1109 calc_triclinic_images(box
,img
);
1112 for(i
=2; i
<=5; i
+=3) {
1120 for(j
=0; j
<4; j
++) {
1121 for(d
=0; d
<DIM
; d
++)
1122 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1126 for(i
=7; i
<=13; i
+=6) {
1134 for(j
=0; j
<4; j
++) {
1135 for(d
=0; d
<DIM
; d
++)
1136 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1140 for(i
=9; i
<=11; i
+=2) {
1151 for(j
=0; j
<4; j
++) {
1152 for(d
=0; d
<DIM
; d
++)
1153 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1158 calc_box_center(ecenter
,box
,box_center
);
1159 for(i
=0; i
<NCUCVERT
; i
++)
1160 for(d
=0; d
<DIM
; d
++)
1161 vert
[i
][d
] = vert
[i
][d
]*0.25+box_center
[d
];
1164 int *compact_unitcell_edges()
1166 /* this is an index in vert[] (see calc_box_vertices) */
1167 /*static int edge[NCUCEDGE*2];*/
1169 static const int hexcon
[24] = { 0,9, 1,19, 2,15, 3,21,
1170 4,17, 5,11, 6,23, 7,13,
1171 8,20, 10,18, 12,16, 14,22 };
1173 gmx_bool bFirst
= TRUE
;
1175 snew(edge
,NCUCEDGE
*2);
1180 for(j
=0; j
<4; j
++) {
1181 edge
[e
++] = 4*i
+ j
;
1182 edge
[e
++] = 4*i
+ (j
+1) % 4;
1184 for(i
=0; i
<12*2; i
++)
1185 edge
[e
++] = hexcon
[i
];
1193 void put_atom_in_box(matrix box
,rvec x
)
1197 for(m
=DIM
-1; m
>=0; m
--) {
1201 while (x
[m
] >= box
[m
][m
])
1207 void put_atoms_in_box(matrix box
,int natoms
,rvec x
[])
1211 for(i
=0; (i
<natoms
); i
++)
1212 put_atom_in_box(box
,x
[i
]);
1215 void put_atoms_in_triclinic_unitcell(int ecenter
,matrix box
,
1216 int natoms
,rvec x
[])
1218 rvec box_center
,shift_center
;
1219 real shm01
,shm02
,shm12
,shift
;
1222 calc_box_center(ecenter
,box
,box_center
);
1224 /* The product of matrix shm with a coordinate gives the shift vector
1225 which is required determine the periodic cell position */
1226 shm01
= box
[1][0]/box
[1][1];
1227 shm02
= (box
[1][1]*box
[2][0] - box
[2][1]*box
[1][0])/(box
[1][1]*box
[2][2]);
1228 shm12
= box
[2][1]/box
[2][2];
1230 clear_rvec(shift_center
);
1231 for(d
=0; d
<DIM
; d
++)
1232 rvec_inc(shift_center
,box
[d
]);
1233 svmul(0.5,shift_center
,shift_center
);
1234 rvec_sub(box_center
,shift_center
,shift_center
);
1236 shift_center
[0] = shm01
*shift_center
[1] + shm02
*shift_center
[2];
1237 shift_center
[1] = shm12
*shift_center
[2];
1238 shift_center
[2] = 0;
1240 for(i
=0; (i
<natoms
); i
++)
1241 for(m
=DIM
-1; m
>=0; m
--) {
1242 shift
= shift_center
[m
];
1244 shift
+= shm01
*x
[i
][1] + shm02
*x
[i
][2];
1245 } else if (m
== 1) {
1246 shift
+= shm12
*x
[i
][2];
1248 while (x
[i
][m
]-shift
< 0)
1250 x
[i
][d
] += box
[m
][d
];
1251 while (x
[i
][m
]-shift
>= box
[m
][m
])
1253 x
[i
][d
] -= box
[m
][d
];
1258 put_atoms_in_compact_unitcell(int ePBC
,int ecenter
,matrix box
,
1259 int natoms
,rvec x
[])
1265 set_pbc(&pbc
,ePBC
,box
);
1266 calc_box_center(ecenter
,box
,box_center
);
1267 for(i
=0; i
<natoms
; i
++) {
1268 pbc_dx(&pbc
,x
[i
],box_center
,dx
);
1269 rvec_add(box_center
,dx
,x
[i
]);
1272 return pbc
.bLimitDistance
?
1273 "WARNING: Could not put all atoms in the compact unitcell\n"