(git:e7e05ae)
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 ! **************************************************************************************************
13 MODULE 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 
35 CONTAINS
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
300 END MODULE negf_integr_utils
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