(git:374b731)
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-2024 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)) &
105 CALL rescale_nodes_linear(nnodes, tnodes, a, b, xnodes)
106
107 IF (PRESENT(weights)) &
108 weights(:) = b - a
109
110 CASE (contour_shape_arc)
111 ALLOCATE (tnodes_angle(nnodes))
112
113 tnodes_angle(:) = tnodes(:)
114 CALL rescale_nodes_pi_phi(a, b, nnodes, tnodes_angle)
115
116 IF (PRESENT(xnodes)) &
117 CALL rescale_nodes_arc(nnodes, tnodes_angle, a, b, xnodes)
118
119 IF (PRESENT(weights)) THEN
120 rscale = (pi - get_arc_smallest_angle(a, b))*get_arc_radius(a, b)
121
122 DO i = 1, nnodes
123 weights(i) = rscale*cmplx(sin(tnodes_angle(i)), -cos(tnodes_angle(i)), kind=dp)
124 END DO
125 END IF
126
127 DEALLOCATE (tnodes_angle)
128 CASE DEFAULT
129 cpabort("Unimplemented integration shape")
130 END SELECT
131
132 CALL timestop(handle)
133 END SUBROUTINE rescale_normalised_nodes
134
135! **************************************************************************************************
136!> \brief Compute arc radius.
137!> \param a lower bound
138!> \param b upper bound
139!> \return radius
140!> \par History
141!> * 05.2017 created [Sergey Chulkov]
142!> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
143! c *
144! r * B-------+------
145! a * / . |
146! * r / . | delta
147! * / phi . |
148! A---------*-----------+------
149! <--- r --><-l->
150! <--- r --->
151! **************************************************************************************************
152 PURE FUNCTION get_arc_radius(a, b) RESULT(radius)
153 COMPLEX(kind=dp), INTENT(in) :: a, b
154 REAL(kind=dp) :: radius
155
156 COMPLEX(kind=dp) :: b_minus_a
157
158 b_minus_a = b - a
159
160 ! l = REAL(B - A); delta = AIMAG(B - A)
161 ! radius = (l^2 + delta^2) / (2 * l)
162 radius = 0.5_dp*real(b_minus_a*conjg(b_minus_a), kind=dp)/real(b_minus_a, kind=dp)
163 END FUNCTION get_arc_radius
164
165! **************************************************************************************************
166!> \brief Compute the angle phi.
167!> \param a lower bound
168!> \param b upper bound
169!> \return angle
170!> \par History
171!> * 05.2017 created [Sergey Chulkov]
172!> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
173! c *
174! r * B-------+------
175! a * / . |
176! * r / . | delta
177! * / phi . |
178! A---------*-----------+------
179! <--- r --><-l->
180! <--- r --->
181! **************************************************************************************************
182 PURE FUNCTION get_arc_smallest_angle(a, b) RESULT(phi)
183 COMPLEX(kind=dp), INTENT(in) :: a, b
184 REAL(kind=dp) :: phi
185
186 COMPLEX(kind=dp) :: b_minus_a
187 REAL(kind=dp) :: delta2, l2
188
189 b_minus_a = b - a
190
191 ! l = REAL(B - A); delta = AIMAG(B - A)
192 ! phi = arccos((l - radius)/radius) = arccos((l^2 - delta^2) / (l^2 + delta^2))
193 l2 = real(b_minus_a, dp)
194 l2 = l2*l2
195 delta2 = aimag(b_minus_a)
196 delta2 = delta2*delta2
197
198 phi = acos((l2 - delta2)/(l2 + delta2))
199 END FUNCTION get_arc_smallest_angle
200
201 PURE FUNCTION get_axis_rotation_angle(a, b) RESULT(phi)
202 COMPLEX(kind=dp), INTENT(in) :: a, b
203 REAL(kind=dp) :: phi
204
205 COMPLEX(kind=dp) :: b_minus_a
206
207 b_minus_a = b - a
208 phi = acos(real(b_minus_a, dp)/abs(b_minus_a))
209 END FUNCTION get_axis_rotation_angle
210
211! **************************************************************************************************
212!> \brief Rescale nodes [pi, phi] -> arc[a, b] .
213!> \param nnodes number of nodes
214!> \param tnodes_angle parametrically-defined nodes to rescale
215!> \param a lower bound
216!> \param b upper bound
217!> \param xnodes rescaled nodes (initialised on exit)
218!> \par History
219!> * 05.2017 created [Sergey Chulkov]
220!> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
221! **************************************************************************************************
222 SUBROUTINE rescale_nodes_arc(nnodes, tnodes_angle, a, b, xnodes)
223 INTEGER, INTENT(in) :: nnodes
224 REAL(kind=dp), DIMENSION(:), INTENT(in) :: tnodes_angle
225 COMPLEX(kind=dp), INTENT(in) :: a, b
226 COMPLEX(kind=dp), DIMENSION(:), INTENT(out) :: xnodes
227
228 COMPLEX(kind=dp) :: origin
229 INTEGER :: i
230 REAL(kind=dp) :: radius
231
232 radius = get_arc_radius(a, b)
233 origin = a + cmplx(radius, 0.0_dp, kind=dp)
234
235 DO i = 1, nnodes
236 xnodes(i) = origin + radius*cmplx(cos(tnodes_angle(i)), sin(tnodes_angle(i)), kind=dp)
237 END DO
238 END SUBROUTINE rescale_nodes_arc
239
240! **************************************************************************************************
241!> \brief Rescale nodes tnodes(i) = cos(pi/2 * (1-tnodes(i))); tnodes \in [-1 .. 1] .
242!> \param tnodes parametrically-defined nodes to rescale / rescaled nodes (modified on exit)
243!> \par History
244!> * 05.2017 created [Sergey Chulkov]
245!> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
246! **************************************************************************************************
247 SUBROUTINE rescale_nodes_cos(nnodes, tnodes)
248 INTEGER, INTENT(in) :: nnodes
249 REAL(kind=dp), DIMENSION(nnodes), INTENT(inout) :: tnodes
250
251 tnodes(:) = cos(0.5_dp*pi*(1.0_dp - tnodes(:)))
252 END SUBROUTINE rescale_nodes_cos
253
254! **************************************************************************************************
255!> \brief Rescale nodes [-1, 1] -> [a, b] .
256!> \param nnodes number of nodes
257!> \param tnodes parametrically-defined nodes to rescale
258!> \param a lower bound
259!> \param b upper bound
260!> \param xnodes rescaled nodes (initialised on exit)
261!> \par History
262!> * 05.2017 created [Sergey Chulkov]
263! **************************************************************************************************
264 SUBROUTINE rescale_nodes_linear(nnodes, tnodes, a, b, xnodes)
265 INTEGER, INTENT(in) :: nnodes
266 REAL(kind=dp), DIMENSION(nnodes), INTENT(in) :: tnodes
267 COMPLEX(kind=dp), INTENT(in) :: a, b
268 COMPLEX(kind=dp), DIMENSION(nnodes), INTENT(out) :: xnodes
269
270 COMPLEX(kind=dp) :: half_len, median
271
272 median = 0.5_dp*(b + a)
273 half_len = 0.5_dp*(b - a)
274
275 xnodes(:) = median + half_len*tnodes(:)
276 END SUBROUTINE rescale_nodes_linear
277
278! **************************************************************************************************
279!> \brief Rescale nodes [-1, 1] -> [pi, phi] .
280!> \param nnodes number of nodes
281!> \param a lower bound
282!> \param b upper bound
283!> \param tnodes parametrically-defined nodes to rescale / rescaled nodes (modified on exit)
284!> \par History
285!> * 05.2017 created [Sergey Chulkov]
286!> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
287! **************************************************************************************************
288 SUBROUTINE rescale_nodes_pi_phi(a, b, nnodes, tnodes)
289 COMPLEX(kind=dp), INTENT(in) :: a, b
290 INTEGER, INTENT(in) :: nnodes
291 REAL(kind=dp), DIMENSION(nnodes), INTENT(inout) :: tnodes
292
293 REAL(kind=dp) :: half_pi_minus_phi, phi
294
295 phi = get_arc_smallest_angle(a, b)
296 half_pi_minus_phi = 0.5_dp*(pi - phi)
297
298 tnodes(:) = phi + half_pi_minus_phi*(1.0_dp - tnodes(:))
299 END SUBROUTINE rescale_nodes_pi_phi
300END 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:223
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:183
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:248
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:153
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:289
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:265