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