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 static bool bWarnedGuess
=FALSE
;
164 int guess_ePBC(matrix box
)
168 if (box
[XX
][XX
]>0 && box
[YY
][YY
]>0 && box
[ZZ
][ZZ
]>0) {
170 } else if (box
[XX
][XX
]>0 && box
[YY
][YY
]>0 && box
[ZZ
][ZZ
]==0) {
172 } else if (box
[XX
][XX
]==0 && box
[YY
][YY
]==0 && box
[ZZ
][ZZ
]==0) {
176 fprintf(stderr
,"WARNING: Unsupported box diagonal %f %f %f, "
177 "will not use periodic boundary conditions\n\n",
178 box
[XX
][XX
],box
[YY
][YY
],box
[ZZ
][ZZ
]);
185 fprintf(debug
,"Guessed pbc = %s from the box matrix\n",epbc_names
[ePBC
]);
190 static int correct_box_elem(FILE *fplog
,int step
,tensor box
,int v
,int d
)
192 int shift
,maxshift
=10;
196 /* correct elem d of vector v with vector d */
197 while (box
[v
][d
] > BOX_MARGIN_CORRECT
*0.5*box
[d
][d
]) {
199 fprintf(fplog
,"Step %d: correcting invalid box:\n",step
);
200 pr_rvecs(fplog
,0,"old box",box
,DIM
);
202 rvec_dec(box
[v
],box
[d
]);
205 pr_rvecs(fplog
,0,"new box",box
,DIM
);
207 if (shift
<= -maxshift
)
209 "Box was shifted at least %d times. Please see log-file.",
212 while (box
[v
][d
] < -BOX_MARGIN_CORRECT
*0.5*box
[d
][d
]) {
214 fprintf(fplog
,"Step %d: correcting invalid box:\n",step
);
215 pr_rvecs(fplog
,0,"old box",box
,DIM
);
217 rvec_inc(box
[v
],box
[d
]);
220 pr_rvecs(fplog
,0,"new box",box
,DIM
);
222 if (shift
>= maxshift
)
224 "Box was shifted at least %d times. Please see log-file.",
231 bool correct_box(FILE *fplog
,int step
,tensor box
,t_graph
*graph
)
236 /* check if the box still obeys the restrictions, if not, correct it */
237 zy
= correct_box_elem(fplog
,step
,box
,ZZ
,YY
);
238 zx
= correct_box_elem(fplog
,step
,box
,ZZ
,XX
);
239 yx
= correct_box_elem(fplog
,step
,box
,YY
,XX
);
241 bCorrected
= (zy
|| zx
|| yx
);
243 if (bCorrected
&& graph
) {
244 /* correct the graph */
245 for(i
=0; i
<graph
->nnodes
; i
++) {
246 graph
->ishift
[i
][YY
] -= graph
->ishift
[i
][ZZ
]*zy
;
247 graph
->ishift
[i
][XX
] -= graph
->ishift
[i
][ZZ
]*zx
;
248 graph
->ishift
[i
][XX
] -= graph
->ishift
[i
][YY
]*yx
;
255 int ndof_com(t_inputrec
*ir
)
265 n
= (ir
->nwall
== 0 ? 3 : 2);
271 gmx_incons("Unknown pbc in calc_nrdf");
277 static void low_set_pbc(t_pbc
*pbc
,int ePBC
,ivec
*dd_nc
,matrix box
)
279 int order
[5]={0,-1,1,-2,2};
280 int ii
,jj
,kk
,i
,j
,k
,d
,dd
,jc
,kc
,npbcdim
,shift
;
282 real d2old
,d2new
,d2new_c
;
287 pbc
->ndim_ePBC
= ePBC2npbcdim(ePBC
);
289 copy_mat(box
,pbc
->box
);
290 pbc
->bLimitDistance
= FALSE
;
291 pbc
->max_cutoff2
= 0;
294 for(i
=0; (i
<DIM
); i
++) {
295 pbc
->fbox_diag
[i
] = box
[i
][i
];
296 pbc
->hbox_diag
[i
] = pbc
->fbox_diag
[i
]*0.5;
297 pbc
->mhbox_diag
[i
] = -pbc
->hbox_diag
[i
];
300 ptr
= check_box(ePBC
,box
);
301 if (ePBC
== epbcNONE
) {
302 pbc
->ePBCDX
= epbcdxNOPBC
;
304 fprintf(stderr
, "Warning: %s\n",ptr
);
305 pr_rvecs(stderr
,0," Box",box
,DIM
);
306 fprintf(stderr
, " Can not fix pbc.\n");
307 pbc
->ePBCDX
= epbcdxUNSUPPORTED
;
308 pbc
->bLimitDistance
= TRUE
;
309 pbc
->limit_distance2
= 0;
311 if (ePBC
== epbcSCREW
&& dd_nc
) {
312 /* This combinated should never appear here */
313 gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
317 for(i
=0; i
<DIM
; i
++) {
318 if ((dd_nc
&& (*dd_nc
)[i
] > 1) || (ePBC
== epbcXY
&& i
== ZZ
)) {
327 /* 1D pbc is not an mdp option and it is therefore only used
328 * with single shifts.
330 pbc
->ePBCDX
= epbcdx1D_RECT
;
334 for(i
=0; i
<pbc
->dim
; i
++)
335 if (pbc
->box
[pbc
->dim
][i
] != 0)
336 pbc
->ePBCDX
= epbcdx1D_TRIC
;
339 pbc
->ePBCDX
= epbcdx2D_RECT
;
346 if (pbc
->box
[i
][j
] != 0)
347 pbc
->ePBCDX
= epbcdx2D_TRIC
;
350 if (ePBC
!= epbcSCREW
) {
351 if (TRICLINIC(box
)) {
352 pbc
->ePBCDX
= epbcdxTRICLINIC
;
354 pbc
->ePBCDX
= epbcdxRECTANGULAR
;
357 pbc
->ePBCDX
= (box
[ZZ
][YY
]==0 ? epbcdxSCREW_RECT
: epbcdxSCREW_TRIC
);
358 if (pbc
->ePBCDX
== epbcdxSCREW_TRIC
) {
360 "Screw pbc is not yet implemented for triclinic boxes.\n"
361 "Can not fix pbc.\n");
362 pbc
->ePBCDX
= epbcdxUNSUPPORTED
;
367 gmx_fatal(FARGS
,"Incorrect number of pbc dimensions with DD: %d",
370 pbc
->max_cutoff2
= max_cutoff2(ePBC
,box
);
372 if (pbc
->ePBCDX
== epbcdxTRICLINIC
||
373 pbc
->ePBCDX
== epbcdx2D_TRIC
||
374 pbc
->ePBCDX
== epbcdxSCREW_TRIC
) {
376 pr_rvecs(debug
,0,"Box",box
,DIM
);
377 fprintf(debug
,"max cutoff %.3f\n",sqrt(pbc
->max_cutoff2
));
380 /* We will only use single shifts, but we will check a few
381 * more shifts to see if there is a limiting distance
382 * above which we can not be sure of the correct distance.
384 for(kk
=0; kk
<5; kk
++) {
386 if (!bPBC
[ZZ
] && k
!= 0)
388 for(jj
=0; jj
<5; jj
++) {
390 if (!bPBC
[YY
] && j
!= 0)
392 for(ii
=0; ii
<3; ii
++) {
394 if (!bPBC
[XX
] && i
!= 0)
396 /* A shift is only useful when it is trilinic */
397 if (j
!= 0 || k
!= 0) {
400 for(d
=0; d
<DIM
; d
++) {
401 trial
[d
] = i
*box
[XX
][d
] + j
*box
[YY
][d
] + k
*box
[ZZ
][d
];
402 /* Choose the vector within the brick around 0,0,0 that
403 * will become the shortest due to shift try.
410 pos
[d
] = min( pbc
->hbox_diag
[d
],-trial
[d
]);
412 pos
[d
] = max(-pbc
->hbox_diag
[d
],-trial
[d
]);
414 d2old
+= sqr(pos
[d
]);
415 d2new
+= sqr(pos
[d
] + trial
[d
]);
417 if (BOX_MARGIN
*d2new
< d2old
) {
418 if (j
< -1 || j
> 1 || k
< -1 || k
> 1) {
419 /* Check if there is a single shift vector
420 * that decreases this distance even more.
430 d2new_c
+= sqr(pos
[d
] + trial
[d
]
431 - jc
*box
[YY
][d
] - kc
*box
[ZZ
][d
]);
432 if (d2new_c
> BOX_MARGIN
*d2new
) {
433 /* Reject this shift vector, as there is no a priori limit
434 * to the number of shifts that decrease distances.
436 if (!pbc
->bLimitDistance
|| d2new
< pbc
->limit_distance2
)
437 pbc
->limit_distance2
= d2new
;
438 pbc
->bLimitDistance
= TRUE
;
441 /* Check if shifts with one box vector less do better */
443 for(dd
=0; dd
<DIM
; dd
++) {
444 shift
= (dd
==0 ? i
: (dd
==1 ? j
: k
));
448 d2new_c
+= sqr(pos
[d
] + trial
[d
] - shift
*box
[dd
][d
]);
449 if (d2new_c
<= BOX_MARGIN
*d2new
)
454 /* Accept this shift vector. */
455 if (pbc
->ntric_vec
>= MAX_NTRICVEC
) {
456 fprintf(stderr
,"\nWARNING: Found more than %d triclinic correction vectors, ignoring some.\n"
457 " There is probably something wrong with your box.\n",MAX_NTRICVEC
);
458 pr_rvecs(stderr
,0," Box",box
,DIM
);
460 copy_rvec(trial
,pbc
->tric_vec
[pbc
->ntric_vec
]);
461 pbc
->tric_shift
[pbc
->ntric_vec
][XX
] = i
;
462 pbc
->tric_shift
[pbc
->ntric_vec
][YY
] = j
;
463 pbc
->tric_shift
[pbc
->ntric_vec
][ZZ
] = k
;
469 fprintf(debug
," tricvec %2d = %2d %2d %2d %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
470 pbc
->ntric_vec
,i
,j
,k
,
471 sqrt(d2old
),sqrt(d2new
),
472 trial
[XX
],trial
[YY
],trial
[ZZ
],
473 pos
[XX
],pos
[YY
],pos
[ZZ
]);
484 void set_pbc(t_pbc
*pbc
,int ePBC
,matrix box
)
487 ePBC
= guess_ePBC(box
);
489 low_set_pbc(pbc
,ePBC
,NULL
,box
);
492 t_pbc
*set_pbc_dd(t_pbc
*pbc
,int ePBC
,
493 gmx_domdec_t
*dd
,bool bSingleDir
,matrix box
)
501 if (ePBC
== epbcSCREW
&& dd
->nc
[XX
] > 1) {
502 /* The rotation has been taken care of during coordinate communication */
506 for(i
=0; i
<DIM
; i
++) {
507 if (dd
->nc
[i
] <= (bSingleDir
? 1 : 2)) {
509 if (!(ePBC
== epbcXY
&& i
== ZZ
))
518 low_set_pbc(pbc
,ePBC
,npbcdim
<DIM
? &nc2
: NULL
,box
);
520 return (npbcdim
> 0 ? pbc
: NULL
);
523 void pbc_dx(const t_pbc
*pbc
,const rvec x1
, const rvec x2
, rvec dx
)
532 switch (pbc
->ePBCDX
) {
533 case epbcdxRECTANGULAR
:
534 for(i
=0; i
<DIM
; i
++) {
535 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
536 dx
[i
] -= pbc
->fbox_diag
[i
];
538 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
539 dx
[i
] += pbc
->fbox_diag
[i
];
543 case epbcdxTRICLINIC
:
544 for(i
=DIM
-1; i
>=0; i
--) {
545 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
547 dx
[j
] -= pbc
->box
[i
][j
];
549 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
551 dx
[j
] += pbc
->box
[i
][j
];
554 /* dx is the distance in a rectangular box */
556 if (d2min
> pbc
->max_cutoff2
) {
557 copy_rvec(dx
,dx_start
);
559 /* Now try all possible shifts, when the distance is within max_cutoff
560 * it must be the shortest possible distance.
563 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
564 rvec_add(dx_start
,pbc
->tric_vec
[i
],trial
);
565 d2trial
= norm2(trial
);
566 if (d2trial
< d2min
) {
575 for(i
=0; i
<DIM
; i
++) {
577 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
578 dx
[i
] -= pbc
->fbox_diag
[i
];
580 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
581 dx
[i
] += pbc
->fbox_diag
[i
];
588 for(i
=DIM
-1; i
>=0; i
--) {
590 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
592 dx
[j
] -= pbc
->box
[i
][j
];
594 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
596 dx
[j
] += pbc
->box
[i
][j
];
598 d2min
+= dx
[i
]*dx
[i
];
601 if (d2min
> pbc
->max_cutoff2
) {
602 copy_rvec(dx
,dx_start
);
604 /* Now try all possible shifts, when the distance is within max_cutoff
605 * it must be the shortest possible distance.
608 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
609 rvec_add(dx_start
,pbc
->tric_vec
[i
],trial
);
611 for(j
=0; j
<DIM
; j
++) {
613 d2trial
+= trial
[j
]*trial
[j
];
616 if (d2trial
< d2min
) {
624 case epbcdxSCREW_RECT
:
625 /* The shift definition requires x first */
627 while (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
628 dx
[XX
] -= pbc
->fbox_diag
[XX
];
631 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
632 dx
[XX
] += pbc
->fbox_diag
[YY
];
636 /* Rotate around the x-axis in the middle of the box */
637 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
638 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
640 /* Normal pbc for y and z */
641 for(i
=YY
; i
<=ZZ
; i
++) {
642 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
643 dx
[i
] -= pbc
->fbox_diag
[i
];
645 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
646 dx
[i
] += pbc
->fbox_diag
[i
];
651 case epbcdxUNSUPPORTED
:
654 gmx_fatal(FARGS
,"Internal error in pbc_dx, set_pbc has not been called");
659 int pbc_dx_aiuc(const t_pbc
*pbc
,const rvec x1
, const rvec x2
, rvec dx
)
664 ivec ishift
,ishift_start
;
669 switch (pbc
->ePBCDX
) {
670 case epbcdxRECTANGULAR
:
671 for(i
=0; i
<DIM
; i
++) {
672 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
673 dx
[i
] -= pbc
->fbox_diag
[i
];
675 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
676 dx
[i
] += pbc
->fbox_diag
[i
];
681 case epbcdxTRICLINIC
:
682 /* For triclinic boxes the performance difference between
683 * if/else and two while loops is negligible.
684 * However, the while version can cause extreme delays
685 * before a simulation crashes due to large forces which
686 * can cause unlimited displacements.
687 * Also allowing multiple shifts would index fshift beyond bounds.
689 for(i
=DIM
-1; i
>=1; i
--) {
690 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
692 dx
[j
] -= pbc
->box
[i
][j
];
694 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
696 dx
[j
] += pbc
->box
[i
][j
];
700 /* Allow 2 shifts in x */
701 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
702 dx
[XX
] -= pbc
->fbox_diag
[XX
];
704 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
705 dx
[XX
] -= pbc
->fbox_diag
[XX
];
708 } else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
709 dx
[XX
] += pbc
->fbox_diag
[XX
];
711 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
712 dx
[XX
] += pbc
->fbox_diag
[XX
];
716 /* dx is the distance in a rectangular box */
718 if (d2min
> pbc
->max_cutoff2
) {
719 copy_rvec(dx
,dx_start
);
720 copy_ivec(ishift
,ishift_start
);
722 /* Now try all possible shifts, when the distance is within max_cutoff
723 * it must be the shortest possible distance.
726 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
727 rvec_add(dx_start
,pbc
->tric_vec
[i
],trial
);
728 d2trial
= norm2(trial
);
729 if (d2trial
< d2min
) {
731 ivec_add(ishift_start
,pbc
->tric_shift
[i
],ishift
);
739 for(i
=0; i
<DIM
; i
++) {
741 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
742 dx
[i
] -= pbc
->fbox_diag
[i
];
744 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
745 dx
[i
] += pbc
->fbox_diag
[i
];
753 for(i
=DIM
-1; i
>=1; i
--) {
755 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
757 dx
[j
] -= pbc
->box
[i
][j
];
759 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
761 dx
[j
] += pbc
->box
[i
][j
];
764 d2min
+= dx
[i
]*dx
[i
];
767 if (pbc
->dim
!= XX
) {
768 /* Allow 2 shifts in x */
769 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
770 dx
[XX
] -= pbc
->fbox_diag
[XX
];
772 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
773 dx
[XX
] -= pbc
->fbox_diag
[XX
];
776 } else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
777 dx
[XX
] += pbc
->fbox_diag
[XX
];
779 if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
780 dx
[XX
] += pbc
->fbox_diag
[XX
];
784 d2min
+= dx
[XX
]*dx
[XX
];
786 if (d2min
> pbc
->max_cutoff2
) {
787 copy_rvec(dx
,dx_start
);
788 copy_ivec(ishift
,ishift_start
);
789 /* Now try all possible shifts, when the distance is within max_cutoff
790 * it must be the shortest possible distance.
793 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
794 rvec_add(dx_start
,pbc
->tric_vec
[i
],trial
);
796 for(j
=0; j
<DIM
; j
++) {
798 d2trial
+= trial
[j
]*trial
[j
];
801 if (d2trial
< d2min
) {
803 ivec_add(ishift_start
,pbc
->tric_shift
[i
],ishift
);
812 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
813 dx
[i
] -= pbc
->fbox_diag
[i
];
815 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
816 dx
[i
] += pbc
->fbox_diag
[i
];
822 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
823 rvec_dec(dx
,pbc
->box
[i
]);
825 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
826 rvec_inc(dx
,pbc
->box
[i
]);
830 case epbcdxSCREW_RECT
:
831 /* The shift definition requires x first */
832 if (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
833 dx
[XX
] -= pbc
->fbox_diag
[XX
];
835 } else if (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
836 dx
[XX
] += pbc
->fbox_diag
[XX
];
839 if (ishift
[XX
] == 1 || ishift
[XX
] == -1) {
840 /* Rotate around the x-axis in the middle of the box */
841 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
842 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
844 /* Normal pbc for y and z */
845 for(i
=YY
; i
<=ZZ
; i
++) {
846 if (dx
[i
] > pbc
->hbox_diag
[i
]) {
847 dx
[i
] -= pbc
->fbox_diag
[i
];
849 } else if (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
850 dx
[i
] += pbc
->fbox_diag
[i
];
856 case epbcdxUNSUPPORTED
:
859 gmx_fatal(FARGS
,"Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
863 is
= IVEC2IS(ishift
);
866 range_check_mesg(is
,0,SHIFTS
,"PBC shift vector index range check.");
872 void pbc_dx_d(const t_pbc
*pbc
,const dvec x1
, const dvec x2
, dvec dx
)
876 double d2min
,d2trial
;
881 switch (pbc
->ePBCDX
) {
882 case epbcdxRECTANGULAR
:
884 for(i
=0; i
<DIM
; i
++) {
886 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
887 dx
[i
] -= pbc
->fbox_diag
[i
];
889 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
890 dx
[i
] += pbc
->fbox_diag
[i
];
895 case epbcdxTRICLINIC
:
898 for(i
=DIM
-1; i
>=0; i
--) {
900 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
902 dx
[j
] -= pbc
->box
[i
][j
];
904 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
906 dx
[j
] += pbc
->box
[i
][j
];
908 d2min
+= dx
[i
]*dx
[i
];
911 if (d2min
> pbc
->max_cutoff2
) {
912 copy_dvec(dx
,dx_start
);
913 /* Now try all possible shifts, when the distance is within max_cutoff
914 * it must be the shortest possible distance.
917 while ((d2min
> pbc
->max_cutoff2
) && (i
< pbc
->ntric_vec
)) {
918 for(j
=0; j
<DIM
; j
++) {
919 trial
[j
] = dx_start
[j
] + pbc
->tric_vec
[i
][j
];
922 for(j
=0; j
<DIM
; j
++) {
924 d2trial
+= trial
[j
]*trial
[j
];
927 if (d2trial
< d2min
) {
935 case epbcdxSCREW_RECT
:
936 /* The shift definition requires x first */
938 while (dx
[XX
] > pbc
->hbox_diag
[XX
]) {
939 dx
[XX
] -= pbc
->fbox_diag
[XX
];
942 while (dx
[XX
] <= pbc
->mhbox_diag
[XX
]) {
943 dx
[XX
] += pbc
->fbox_diag
[YY
];
947 /* Rotate around the x-axis in the middle of the box */
948 dx
[YY
] = pbc
->box
[YY
][YY
] - x1
[YY
] - x2
[YY
];
949 dx
[ZZ
] = pbc
->box
[ZZ
][ZZ
] - x1
[ZZ
] - x2
[ZZ
];
951 /* Normal pbc for y and z */
952 for(i
=YY
; i
<=ZZ
; i
++) {
953 while (dx
[i
] > pbc
->hbox_diag
[i
]) {
954 dx
[i
] -= pbc
->fbox_diag
[i
];
956 while (dx
[i
] <= pbc
->mhbox_diag
[i
]) {
957 dx
[i
] += pbc
->fbox_diag
[i
];
962 case epbcdxUNSUPPORTED
:
965 gmx_fatal(FARGS
,"Internal error in pbc_dx, set_pbc has not been called");
970 bool image_rect(ivec xi
,ivec xj
,ivec box_size
,real rlong2
,int *shift
,real
*r2
)
978 for(m
=0; (m
<DIM
); m
++) {
1002 bool image_cylindric(ivec xi
,ivec xj
,ivec box_size
,real rlong2
,
1003 int *shift
,real
*r2
)
1011 for(m
=0; (m
<DIM
); m
++) {
1038 void calc_shifts(matrix box
,rvec shift_vec
[])
1043 for(m
= -D_BOX_Z
; m
<= D_BOX_Z
; m
++)
1044 for(l
= -D_BOX_Y
; l
<= D_BOX_Y
; l
++)
1045 for(k
= -D_BOX_X
; k
<= D_BOX_X
; k
++,n
++) {
1046 test
= XYZ2IS(k
,l
,m
);
1048 gmx_incons("inconsistent shift numbering");
1049 for(d
= 0; d
< DIM
; d
++)
1050 shift_vec
[n
][d
] = k
*box
[XX
][d
] + l
*box
[YY
][d
] + m
*box
[ZZ
][d
];
1054 void calc_box_center(int ecenter
,matrix box
,rvec box_center
)
1058 clear_rvec(box_center
);
1061 for(m
=0; (m
<DIM
); m
++)
1062 for(d
=0; d
<DIM
; d
++)
1063 box_center
[d
] += 0.5*box
[m
][d
];
1066 for(d
=0; d
<DIM
; d
++)
1067 box_center
[d
] = 0.5*box
[d
][d
];
1072 gmx_fatal(FARGS
,"Unsupported value %d for ecenter",ecenter
);
1076 void calc_triclinic_images(matrix box
,rvec img
[])
1080 /* Calculate 3 adjacent images in the xy-plane */
1081 copy_rvec(box
[0],img
[0]);
1082 copy_rvec(box
[1],img
[1]);
1084 svmul(-1,img
[1],img
[1]);
1085 rvec_sub(img
[1],img
[0],img
[2]);
1087 /* Get the next 3 in the xy-plane as mirror images */
1089 svmul(-1,img
[i
],img
[3+i
]);
1091 /* Calculate the first 4 out of xy-plane images */
1092 copy_rvec(box
[2],img
[6]);
1094 svmul(-1,img
[6],img
[6]);
1096 rvec_add(img
[6],img
[i
+1],img
[7+i
]);
1098 /* Mirror the last 4 from the previous in opposite rotation */
1100 svmul(-1,img
[6 + (2+i
) % 4],img
[10+i
]);
1103 void calc_compact_unitcell_vertices(int ecenter
,matrix box
,rvec vert
[])
1105 rvec img
[NTRICIMG
],box_center
;
1108 calc_triclinic_images(box
,img
);
1111 for(i
=2; i
<=5; i
+=3) {
1119 for(j
=0; j
<4; j
++) {
1120 for(d
=0; d
<DIM
; d
++)
1121 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1125 for(i
=7; i
<=13; i
+=6) {
1133 for(j
=0; j
<4; j
++) {
1134 for(d
=0; d
<DIM
; d
++)
1135 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1139 for(i
=9; i
<=11; i
+=2) {
1150 for(j
=0; j
<4; j
++) {
1151 for(d
=0; d
<DIM
; d
++)
1152 vert
[n
][d
] = img
[i
][d
]+img
[tmp
[j
]][d
]+img
[tmp
[(j
+1)%4]][d
];
1157 calc_box_center(ecenter
,box
,box_center
);
1158 for(i
=0; i
<NCUCVERT
; i
++)
1159 for(d
=0; d
<DIM
; d
++)
1160 vert
[i
][d
] = vert
[i
][d
]*0.25+box_center
[d
];
1163 int *compact_unitcell_edges()
1165 /* this is an index in vert[] (see calc_box_vertices) */
1166 static int edge
[NCUCEDGE
*2];
1167 static int hexcon
[24] = { 0,9, 1,19, 2,15, 3,21,
1168 4,17, 5,11, 6,23, 7,13,
1169 8,20, 10,18, 12,16, 14,22 };
1176 for(j
=0; j
<4; j
++) {
1177 edge
[e
++] = 4*i
+ j
;
1178 edge
[e
++] = 4*i
+ (j
+1) % 4;
1180 for(i
=0; i
<12*2; i
++)
1181 edge
[e
++] = hexcon
[i
];
1189 void put_atom_in_box(matrix box
,rvec x
)
1193 for(m
=DIM
-1; m
>=0; m
--) {
1197 while (x
[m
] >= box
[m
][m
])
1203 void put_atoms_in_box(matrix box
,int natoms
,rvec x
[])
1207 for(i
=0; (i
<natoms
); i
++)
1208 put_atom_in_box(box
,x
[i
]);
1211 void put_atoms_in_triclinic_unitcell(int ecenter
,matrix box
,
1212 int natoms
,rvec x
[])
1214 rvec box_center
,shift_center
;
1215 real shm01
,shm02
,shm12
,shift
;
1218 calc_box_center(ecenter
,box
,box_center
);
1220 /* The product of matrix shm with a coordinate gives the shift vector
1221 which is required determine the periodic cell position */
1222 shm01
= box
[1][0]/box
[1][1];
1223 shm02
= (box
[1][1]*box
[2][0] - box
[2][1]*box
[1][0])/(box
[1][1]*box
[2][2]);
1224 shm12
= box
[2][1]/box
[2][2];
1226 clear_rvec(shift_center
);
1227 for(d
=0; d
<DIM
; d
++)
1228 rvec_inc(shift_center
,box
[d
]);
1229 svmul(0.5,shift_center
,shift_center
);
1230 rvec_sub(box_center
,shift_center
,shift_center
);
1232 shift_center
[0] = shm01
*shift_center
[1] + shm02
*shift_center
[2];
1233 shift_center
[1] = shm12
*shift_center
[2];
1234 shift_center
[2] = 0;
1236 for(i
=0; (i
<natoms
); i
++)
1237 for(m
=DIM
-1; m
>=0; m
--) {
1238 shift
= shift_center
[m
];
1240 shift
+= shm01
*x
[i
][1] + shm02
*x
[i
][2];
1241 } else if (m
== 1) {
1242 shift
+= shm12
*x
[i
][2];
1244 while (x
[i
][m
]-shift
< 0)
1246 x
[i
][d
] += box
[m
][d
];
1247 while (x
[i
][m
]-shift
>= box
[m
][m
])
1249 x
[i
][d
] -= box
[m
][d
];
1254 put_atoms_in_compact_unitcell(int ePBC
,int ecenter
,matrix box
,
1255 int natoms
,rvec x
[])
1261 set_pbc(&pbc
,ePBC
,box
);
1262 calc_box_center(ecenter
,box
,box_center
);
1263 for(i
=0; i
<natoms
; i
++) {
1264 pbc_dx(&pbc
,x
[i
],box_center
,dx
);
1265 rvec_add(box_center
,dx
,x
[i
]);
1268 return pbc
.bLimitDistance
?
1269 "WARNING: Could not put all atoms in the compact unitcell\n"