Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Helper functions for integration routines.
10 : !> \par History
11 : !> * 06.2017 created [Sergey Chulkov]
12 : ! **************************************************************************************************
13 : MODULE negf_integr_utils
14 : USE kinds, ONLY: dp
15 : USE mathconstants, ONLY: pi
16 : #include "./base/base_uses.f90"
17 : #:include 'negf_integr_utils.fypp'
18 : IMPLICIT NONE
19 : PRIVATE
20 :
21 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_integr_utils'
22 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
23 :
24 : PUBLIC :: equidistant_nodes_a_b, rescale_normalised_nodes
25 : PUBLIC :: get_arc_radius, get_arc_smallest_angle
26 : PUBLIC :: rescale_nodes_arc, rescale_nodes_cos, rescale_nodes_linear, rescale_nodes_pi_phi
27 :
28 : INTEGER, PARAMETER, PUBLIC :: contour_shape_linear = 0, &
29 : contour_shape_arc = 1
30 :
31 : INTERFACE equidistant_nodes_a_b
32 : #:for nametype1, type1 in inst_params
33 : MODULE PROCEDURE equidistant_${nametype1}$nodes_a_b
34 : #:endfor
35 : END INTERFACE
36 :
37 : CONTAINS
38 :
39 : #:for nametype1, type1 in inst_params
40 : ! **************************************************************************************************
41 : !> \brief Compute equidistant nodes on an interval [a, b], where a and b are complex numbers.
42 : !> \param a lower bound
43 : !> \param b upper bound
44 : !> \param nnodes number of nodes
45 : !> \param xnodes array to store the nodes
46 : !> \par History
47 : !> * 05.2017 created [Sergey Chulkov]
48 : ! **************************************************************************************************
49 24 : SUBROUTINE equidistant_${nametype1}$nodes_a_b(a, b, nnodes, xnodes)
50 : ${type1}$, INTENT(in) :: a, b
51 : INTEGER, INTENT(in) :: nnodes
52 : ${type1}$, DIMENSION(nnodes), INTENT(out) :: xnodes
53 :
54 : INTEGER :: i
55 : ${type1}$ :: rscale
56 :
57 24 : CPASSERT(nnodes >= 1)
58 :
59 24 : rscale = (b - a)/REAL(nnodes - 1, kind=dp)
60 336 : DO i = 1, nnodes
61 336 : xnodes(i) = a + rscale*REAL(i - 1, kind=dp)
62 : END DO
63 24 : END SUBROUTINE equidistant_${nametype1}$nodes_a_b
64 : #:endfor
65 :
66 260 : SUBROUTINE rescale_normalised_nodes(nnodes, tnodes, a, b, shape_id, xnodes, weights)
67 : INTEGER, INTENT(in) :: nnodes
68 : REAL(kind=dp), DIMENSION(nnodes), INTENT(in) :: tnodes
69 : COMPLEX(kind=dp), INTENT(in) :: a, b
70 : INTEGER, INTENT(in) :: shape_id
71 : COMPLEX(kind=dp), DIMENSION(nnodes), INTENT(out), &
72 : OPTIONAL :: xnodes, weights
73 :
74 : CHARACTER(len=*), PARAMETER :: routineN = 'rescale_normalised_nodes'
75 :
76 : INTEGER :: handle, i
77 : REAL(kind=dp) :: rscale
78 260 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes_angle
79 :
80 260 : CALL timeset(routineN, handle)
81 :
82 344 : SELECT CASE (shape_id)
83 : CASE (contour_shape_linear)
84 84 : IF (PRESENT(xnodes)) &
85 42 : CALL rescale_nodes_linear(nnodes, tnodes, a, b, xnodes)
86 :
87 84 : IF (PRESENT(weights)) &
88 548 : weights(:) = b - a
89 :
90 : CASE (contour_shape_arc)
91 528 : ALLOCATE (tnodes_angle(nnodes))
92 :
93 4516 : tnodes_angle(:) = tnodes(:)
94 176 : CALL rescale_nodes_pi_phi(a, b, nnodes, tnodes_angle)
95 :
96 176 : IF (PRESENT(xnodes)) &
97 88 : CALL rescale_nodes_arc(nnodes, tnodes_angle, a, b, xnodes)
98 :
99 176 : IF (PRESENT(weights)) THEN
100 88 : rscale = (pi - get_arc_smallest_angle(a, b))*get_arc_radius(a, b)
101 :
102 2258 : DO i = 1, nnodes
103 2258 : weights(i) = rscale*CMPLX(SIN(tnodes_angle(i)), -COS(tnodes_angle(i)), kind=dp)
104 : END DO
105 : END IF
106 :
107 176 : DEALLOCATE (tnodes_angle)
108 : CASE DEFAULT
109 260 : CPABORT("Unimplemented integration shape")
110 : END SELECT
111 :
112 260 : CALL timestop(handle)
113 650 : END SUBROUTINE rescale_normalised_nodes
114 :
115 : ! **************************************************************************************************
116 : !> \brief Compute arc radius.
117 : !> \param a lower bound
118 : !> \param b upper bound
119 : !> \return radius
120 : !> \par History
121 : !> * 05.2017 created [Sergey Chulkov]
122 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
123 : ! c *
124 : ! r * B-------+------
125 : ! a * / . |
126 : ! * r / . | delta
127 : ! * / phi . |
128 : ! A---------*-----------+------
129 : ! <--- r --><-l->
130 : ! <--- r --->
131 : ! **************************************************************************************************
132 176 : PURE FUNCTION get_arc_radius(a, b) RESULT(radius)
133 : COMPLEX(kind=dp), INTENT(in) :: a, b
134 : REAL(kind=dp) :: radius
135 :
136 : COMPLEX(kind=dp) :: b_minus_a
137 :
138 176 : b_minus_a = b - a
139 :
140 : ! l = REAL(B - A); delta = AIMAG(B - A)
141 : ! radius = (l^2 + delta^2) / (2 * l)
142 176 : radius = 0.5_dp*REAL(b_minus_a*CONJG(b_minus_a), kind=dp)/REAL(b_minus_a, kind=dp)
143 176 : END FUNCTION get_arc_radius
144 :
145 : ! **************************************************************************************************
146 : !> \brief Compute the angle phi.
147 : !> \param a lower bound
148 : !> \param b upper bound
149 : !> \return angle
150 : !> \par History
151 : !> * 05.2017 created [Sergey Chulkov]
152 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
153 : ! c *
154 : ! r * B-------+------
155 : ! a * / . |
156 : ! * r / . | delta
157 : ! * / phi . |
158 : ! A---------*-----------+------
159 : ! <--- r --><-l->
160 : ! <--- r --->
161 : ! **************************************************************************************************
162 264 : PURE FUNCTION get_arc_smallest_angle(a, b) RESULT(phi)
163 : COMPLEX(kind=dp), INTENT(in) :: a, b
164 : REAL(kind=dp) :: phi
165 :
166 : COMPLEX(kind=dp) :: b_minus_a
167 : REAL(kind=dp) :: delta2, l2
168 :
169 264 : b_minus_a = b - a
170 :
171 : ! l = REAL(B - A); delta = AIMAG(B - A)
172 : ! phi = arccos((l - radius)/radius) = arccos((l^2 - delta^2) / (l^2 + delta^2))
173 264 : l2 = REAL(b_minus_a, dp)
174 264 : l2 = l2*l2
175 264 : delta2 = AIMAG(b_minus_a)
176 264 : delta2 = delta2*delta2
177 :
178 264 : phi = ACOS((l2 - delta2)/(l2 + delta2))
179 264 : END FUNCTION get_arc_smallest_angle
180 :
181 0 : PURE FUNCTION get_axis_rotation_angle(a, b) RESULT(phi)
182 : COMPLEX(kind=dp), INTENT(in) :: a, b
183 : REAL(kind=dp) :: phi
184 :
185 : COMPLEX(kind=dp) :: b_minus_a
186 :
187 0 : b_minus_a = b - a
188 0 : phi = ACOS(REAL(b_minus_a, dp)/ABS(b_minus_a))
189 0 : END FUNCTION get_axis_rotation_angle
190 :
191 : ! **************************************************************************************************
192 : !> \brief Rescale nodes [pi, phi] -> arc[a, b] .
193 : !> \param nnodes number of nodes
194 : !> \param tnodes_angle parametrically-defined nodes to rescale
195 : !> \param a lower bound
196 : !> \param b upper bound
197 : !> \param xnodes rescaled nodes (initialised on exit)
198 : !> \par History
199 : !> * 05.2017 created [Sergey Chulkov]
200 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
201 : ! **************************************************************************************************
202 88 : SUBROUTINE rescale_nodes_arc(nnodes, tnodes_angle, a, b, xnodes)
203 : INTEGER, INTENT(in) :: nnodes
204 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: tnodes_angle
205 : COMPLEX(kind=dp), INTENT(in) :: a, b
206 : COMPLEX(kind=dp), DIMENSION(:), INTENT(out) :: xnodes
207 :
208 : COMPLEX(kind=dp) :: origin
209 : INTEGER :: i
210 : REAL(kind=dp) :: radius
211 :
212 88 : radius = get_arc_radius(a, b)
213 88 : origin = a + CMPLX(radius, 0.0_dp, kind=dp)
214 :
215 2258 : DO i = 1, nnodes
216 2258 : xnodes(i) = origin + radius*CMPLX(COS(tnodes_angle(i)), SIN(tnodes_angle(i)), kind=dp)
217 : END DO
218 88 : END SUBROUTINE rescale_nodes_arc
219 :
220 : ! **************************************************************************************************
221 : !> \brief Rescale nodes tnodes(i) = cos(pi/2 * (1-tnodes(i))); tnodes \in [-1 .. 1] .
222 : !> \param tnodes parametrically-defined nodes to rescale / rescaled nodes (modified on exit)
223 : !> \par History
224 : !> * 05.2017 created [Sergey Chulkov]
225 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
226 : ! **************************************************************************************************
227 0 : SUBROUTINE rescale_nodes_cos(nnodes, tnodes)
228 : INTEGER, INTENT(in) :: nnodes
229 : REAL(kind=dp), DIMENSION(nnodes), INTENT(inout) :: tnodes
230 :
231 0 : tnodes(:) = COS(0.5_dp*pi*(1.0_dp - tnodes(:)))
232 0 : END SUBROUTINE rescale_nodes_cos
233 :
234 : ! **************************************************************************************************
235 : !> \brief Rescale nodes [-1, 1] -> [a, b] .
236 : !> \param nnodes number of nodes
237 : !> \param tnodes parametrically-defined nodes to rescale
238 : !> \param a lower bound
239 : !> \param b upper bound
240 : !> \param xnodes rescaled nodes (initialised on exit)
241 : !> \par History
242 : !> * 05.2017 created [Sergey Chulkov]
243 : ! **************************************************************************************************
244 42 : SUBROUTINE rescale_nodes_linear(nnodes, tnodes, a, b, xnodes)
245 : INTEGER, INTENT(in) :: nnodes
246 : REAL(kind=dp), DIMENSION(nnodes), INTENT(in) :: tnodes
247 : COMPLEX(kind=dp), INTENT(in) :: a, b
248 : COMPLEX(kind=dp), DIMENSION(nnodes), INTENT(out) :: xnodes
249 :
250 : COMPLEX(kind=dp) :: half_len, median
251 :
252 42 : median = 0.5_dp*(b + a)
253 42 : half_len = 0.5_dp*(b - a)
254 :
255 548 : xnodes(:) = median + half_len*tnodes(:)
256 42 : END SUBROUTINE rescale_nodes_linear
257 :
258 : ! **************************************************************************************************
259 : !> \brief Rescale nodes [-1, 1] -> [pi, phi] .
260 : !> \param nnodes number of nodes
261 : !> \param a lower bound
262 : !> \param b upper bound
263 : !> \param tnodes parametrically-defined nodes to rescale / rescaled nodes (modified on exit)
264 : !> \par History
265 : !> * 05.2017 created [Sergey Chulkov]
266 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
267 : ! **************************************************************************************************
268 176 : SUBROUTINE rescale_nodes_pi_phi(a, b, nnodes, tnodes)
269 : COMPLEX(kind=dp), INTENT(in) :: a, b
270 : INTEGER, INTENT(in) :: nnodes
271 : REAL(kind=dp), DIMENSION(nnodes), INTENT(inout) :: tnodes
272 :
273 : REAL(kind=dp) :: half_pi_minus_phi, phi
274 :
275 176 : phi = get_arc_smallest_angle(a, b)
276 176 : half_pi_minus_phi = 0.5_dp*(pi - phi)
277 :
278 4516 : tnodes(:) = phi + half_pi_minus_phi*(1.0_dp - tnodes(:))
279 176 : END SUBROUTINE rescale_nodes_pi_phi
280 : END MODULE negf_integr_utils
|
