Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_boxdim
8 : CONTAINS
9 51695 : SUBROUTINE boxdim( &
10 : bmat, &
11 : arltv1, arltv2, arltv3)
12 : !*********************************************************************
13 : ! This subroutine determines the maximum number L, M and N
14 : ! nrgvl, nrgvm, nrgvn, of reciprocal lattice vectors G(L,M,N)
15 : ! along the directions G(1), G(2), G(3), respectively, for which
16 :
17 : ! | G(L,M,N) | <= GMAX.
18 :
19 : ! This equation defines a "sphere" and the nrgvl,m,n define
20 : ! the DIMension of the BOX in which the sphere is placed.
21 :
22 : ! In reality the G(i)'s, do not form a carteasian coordinate
23 : ! system. Therefore the "sphere" is not a sphere, but an
24 : ! ellipsoid. For this ellipsoid the largest L, M and N component
25 : ! is determined as arltv1, arltv2, arltv3 to construct the boxes and
26 :
27 : ! nrgvl = int( gmax/arltv1 ) + 1
28 : ! nrgvm = int( gmax/arltv2 ) + 1
29 : ! nrgvn = int( gmax/arltv3 ) + 1
30 :
31 : ! G(i,xyz) is stored in bmat(i,xyz)
32 :
33 : ! routine by s.bluegel from carpar-program
34 :
35 : ! S. Bl"ugel, IFF, 13. Nov. 97
36 : ! tested by S. Heinze , IFF,
37 : !*********************************************************************
38 : USE m_juDFT
39 : ! SE m_constants
40 :
41 : ! .. Parameters ..
42 : IMPLICIT NONE
43 :
44 : ! .. Scalar Arguments ..
45 : REAL, INTENT(OUT) :: arltv1, arltv2, arltv3
46 : ! ..
47 : ! .. Array Arguments ..
48 : REAL, INTENT(IN) :: bmat(3, 3)
49 : ! ..
50 : ! .. Local Scalars ..
51 : INTEGER :: ixyz, j, k
52 : REAL :: denom, eps, one, zero
53 : ! ..
54 : ! .. Local Arrays ..
55 : REAL :: det(3, 3), rr(3, 3)
56 : ! ..
57 : DATA one, eps, zero/1.0, 1e-10, 0.0/
58 :
59 : !---> build up quadratic form for ellipsoid
60 :
61 206780 : DO j = 1, 3
62 672035 : DO k = 1, 3
63 465255 : rr(k, j) = zero
64 2016105 : DO ixyz = 1, 3
65 1861020 : rr(k, j) = rr(k, j) + bmat(k, ixyz)*bmat(j, ixyz)
66 : ENDDO
67 : ENDDO
68 : ENDDO
69 :
70 : !---> build determinants for Cramer's rule
71 :
72 51695 : det(1, 1) = rr(2, 2)*rr(3, 3) - rr(3, 2)*rr(2, 3)
73 51695 : det(1, 2) = rr(2, 1)*rr(3, 3) - rr(3, 1)*rr(2, 3)
74 51695 : det(1, 3) = rr(2, 1)*rr(3, 2) - rr(3, 1)*rr(2, 2)
75 51695 : det(2, 1) = rr(1, 2)*rr(3, 3) - rr(3, 2)*rr(1, 3)
76 51695 : det(2, 2) = rr(1, 1)*rr(3, 3) - rr(3, 1)*rr(1, 3)
77 51695 : det(2, 3) = rr(1, 1)*rr(3, 2) - rr(3, 1)*rr(1, 2)
78 51695 : det(3, 1) = rr(1, 2)*rr(2, 3) - rr(2, 2)*rr(1, 3)
79 51695 : det(3, 2) = rr(1, 1)*rr(2, 3) - rr(2, 1)*rr(1, 3)
80 51695 : det(3, 3) = rr(1, 1)*rr(2, 2) - rr(2, 1)*rr(1, 2)
81 :
82 : !---> check on the zeros of some determinants
83 :
84 206780 : DO j = 1, 3
85 206780 : IF (det(j, j) < eps) THEN
86 : CALL juDFT_error("boxdim: determinant. Problem with det(" &
87 0 : //int2str(j) // ","//int2str(j) // ")", calledby="boxdim")
88 : END IF
89 : ENDDO
90 :
91 : !---> scale determinants
92 :
93 206780 : DO k = 1, 3
94 155085 : denom = one/det(k, k)
95 620340 : DO j = 1, 3
96 620340 : det(k, j) = denom*det(k, j)
97 : ENDDO
98 206780 : det(k, k) = one/denom
99 : ENDDO
100 :
101 : !---> calculate the maximum l, m, n components of the ellipsoid
102 :
103 : arltv1 = sqrt(rr(1, 1) + rr(2, 2)*det(1, 2)*det(1, 2) &
104 : + rr(3, 3)*det(1, 3)*det(1, 3) &
105 : - (rr(1, 2) + rr(1, 2))*det(1, 2) &
106 : + (rr(1, 3) + rr(1, 3))*det(1, 3) &
107 51695 : - (rr(2, 3) + rr(2, 3))*det(1, 2)*det(1, 3))
108 : arltv2 = sqrt(rr(2, 2) + rr(1, 1)*det(2, 1)*det(2, 1) &
109 : + rr(3, 3)*det(2, 3)*det(2, 3) &
110 : - (rr(1, 2) + rr(1, 2))*det(2, 1) &
111 : + (rr(1, 3) + rr(1, 3))*det(2, 1)*det(2, 3) &
112 51695 : - (rr(2, 3) + rr(2, 3))*det(2, 3))
113 : arltv3 = sqrt(rr(3, 3) + rr(1, 1)*det(3, 1)*det(3, 1) &
114 : + rr(2, 2)*det(3, 2)*det(3, 2) &
115 : - (rr(1, 2) + rr(1, 2))*det(3, 1)*det(3, 2) &
116 : + (rr(1, 3) + rr(1, 3))*det(3, 1) &
117 51695 : - (rr(2, 3) + rr(2, 3))*det(3, 2))
118 :
119 51695 : RETURN
120 : end SUBROUTINE boxdim
121 : end module m_boxdim
|
