Line data Source code
1 : MODULE m_tetraWeight
2 :
3 : IMPLICIT NONE
4 :
5 : PRIVATE
6 : PUBLIC :: tetraWeight
7 :
8 : CONTAINS
9 :
10 37029326 : PURE REAL FUNCTION tetraWeight(efermi,etetra,ind,vol,film)
11 :
12 : REAL, INTENT(IN) :: efermi
13 : REAL, INTENT(IN) :: etetra(:)
14 : INTEGER, INTENT(IN) :: ind
15 : REAL, INTENT(IN) :: vol
16 : LOGICAL, INTENT(IN) :: film
17 :
18 37029326 : IF(film) THEN
19 295578 : tetraWeight = tetraWeightFilm(efermi,etetra,ind) * vol/3.0
20 : ELSE
21 36733748 : tetraWeight = tetraWeightBulk(efermi,etetra,ind) * vol/4.0
22 : ENDIF
23 :
24 37029326 : END FUNCTION tetraWeight
25 :
26 36733748 : PURE REAL FUNCTION tetraWeightBulk(efermi,etetra,ind)
27 :
28 : !-------------------------------------------------------------
29 : ! Returns the integration weights for a single tetrahedron
30 : ! at one corner for charge density integration with the
31 : ! given fermi efermi without the volume factor in Bulk
32 : !
33 : ! The energies at the corners etetra are ordered so that
34 : ! etetra(1) <= etetra(2) <= etetra(3) <= etetra(4)
35 : ! ind is the corner index of the kpoint we are interested in
36 : !-------------------------------------------------------------
37 :
38 : REAL, INTENT(IN) :: efermi
39 : REAL, INTENT(IN) :: etetra(:)
40 : INTEGER, INTENT(IN) :: ind
41 :
42 : REAL C1, C2, C3
43 : INTEGER i,j
44 : REAL f(4)
45 : REAL e(4,4)
46 :
47 36733748 : tetraWeightBulk = 0.0
48 :
49 : !Prepare differences between eigenvalues at the corners e(i,j) = ei-ej
50 : !This is to improve readability in the formulas
51 183668740 : DO i = 1, 4
52 771408708 : DO j = 1, 4
53 734674960 : e(i,j) = etetra(i) - etetra(j)
54 : ENDDO
55 : ENDDO
56 :
57 36733748 : IF( efermi>=etetra(4) ) THEN
58 :
59 : tetraWeightBulk = 1.0
60 :
61 22445944 : ELSE IF( efermi>=etetra(3) ) THEN
62 :
63 19343264 : f(1:3) = ( etetra(4)-efermi )/e(4,1:3)
64 :
65 6044770 : SELECT CASE(ind)
66 : CASE(4)
67 1208954 : tetraWeightBulk = 1.0 - f(1) * f(2) * f(3) * ( 4.0 - f(1) - f(2) - f(3))
68 : CASE DEFAULT
69 4835816 : tetraWeightBulk = 1.0 - f(1) * f(2) * f(3) * f(ind)
70 : END SELECT
71 :
72 17610128 : ELSE IF( efermi>=etetra(2) ) THEN
73 :
74 3543196 : C1 = ( efermi-etetra(1) )**2/( e(4,1)*e(3,1) )
75 :
76 3543196 : C2 = ( efermi-etetra(1) )*( efermi-etetra(2) )*( etetra(3)-efermi )/( e(3,1)*e(3,2)*e(4,1) )
77 :
78 3543196 : C3 = ( efermi-etetra(2) )**2*( etetra(4)-efermi )/( e(3,2)*e(4,1)*e(4,2) )
79 :
80 4428995 : SELECT CASE(ind)
81 : CASE(1)
82 : tetraWeightBulk = C1 + (C1 + C2) * (etetra(3)-efermi)/e(3,1) + &
83 885799 : (C1 + C2 + C3) * (etetra(4)-efermi)/e(4,1)
84 : CASE(2)
85 : tetraWeightBulk = C1 + C2 + C3 + (C2 + C3) * (etetra(3)-efermi)/e(3,2) + &
86 885799 : C3 * (etetra(4)-efermi)/e(4,2)
87 : CASE(3)
88 : tetraWeightBulk = (C1 + C2) * (efermi-etetra(1))/e(3,1)+&
89 885799 : (C2 + C3) * (efermi-etetra(2))/e(3,2)
90 : CASE(4)
91 : tetraWeightBulk = (C1 + C2 + C3) * (efermi-etetra(1))/e(4,1) +&
92 3543196 : C3 * (efermi-etetra(2))/e(4,2)
93 : CASE DEFAULT
94 : END SELECT
95 :
96 14066932 : ELSE IF( efermi>=etetra(1) ) THEN
97 :
98 19871280 : f(2:4) = ( efermi-etetra(1) )/e(2:4,1)
99 :
100 6209775 : SELECT CASE(ind)
101 : CASE(1)
102 1241955 : tetraWeightBulk = f(2) * f(3) * f(4) * (4.0 - f(2) - f(3) - f(4))
103 : CASE DEFAULT
104 4967820 : tetraWeightBulk = f(2) * f(3) * f(4) * f(ind)
105 : END SELECT
106 :
107 : END IF
108 :
109 36733748 : END FUNCTION tetraWeightBulk
110 :
111 295578 : PURE REAL FUNCTION tetraWeightFilm(efermi,etetra,ind)
112 :
113 : !-------------------------------------------------------------
114 : ! Returns the integration weights for a single tetrahedron
115 : ! at one corner for charge density integration with the
116 : ! given fermi efermi without the volume factor in Film
117 : !
118 : ! The energies at the corners etetra are ordered so that
119 : ! etetra(1) <= etetra(2) <= etetra(3)
120 : ! ind is the corner index of the kpoint we are interested in
121 : !-------------------------------------------------------------
122 :
123 : REAL, INTENT(IN) :: efermi
124 : REAL, INTENT(IN) :: etetra(:)
125 : INTEGER, INTENT(IN) :: ind
126 :
127 : REAL f(3)
128 :
129 295578 : tetraWeightFilm = 0.0
130 :
131 295578 : IF( efermi>etetra(3) ) THEN
132 :
133 : tetraWeightFilm = 1.0
134 :
135 240882 : ELSE IF( efermi>=etetra(2) ) THEN
136 :
137 352188 : f(1:2) = ( etetra(3)-efermi )/( etetra(3)-etetra(1:2) )
138 :
139 156528 : SELECT CASE(ind)
140 : CASE(3)
141 39132 : tetraWeightFilm = 1.0 - f(1) * f(2) * ( 3.0 - f(1) - f(2) )
142 : CASE DEFAULT
143 117396 : tetraWeightFilm = 1.0 - f(1) * f(2) * f(ind)
144 : END SELECT
145 :
146 123486 : ELSE IF( efermi>=etetra(1) ) THEN
147 :
148 354258 : f(2:3) = ( efermi-etetra(1) )/( etetra(2:3)-etetra(1) )
149 :
150 157448 : SELECT CASE(ind)
151 : CASE(1)
152 39362 : tetraWeightFilm = f(2) * f(3) * ( 3.0 - f(2) - f(3) )
153 : CASE DEFAULT
154 118086 : tetraWeightFilm = f(2) * f(3) * f(ind)
155 : END SELECT
156 :
157 : END IF
158 :
159 295578 : END FUNCTION tetraWeightFilm
160 :
161 : END MODULE m_tetraWeight
|
