Line data Source code
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 : ! **************************************************************************************************
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 60 : 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 60 : CPASSERT(nnodes >= 1)
58 :
59 60 : rscale = (b - a)/REAL(nnodes - 1, kind=dp)
60 840 : DO i = 1, nnodes
61 840 : xnodes(i) = a + rscale*REAL(i - 1, kind=dp)
62 : END DO
63 60 : END SUBROUTINE equidistant_${nametype1}$nodes_a_b
64 : #:endfor
65 :
66 960 : 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 960 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes_angle
79 :
80 960 : CALL timeset(routineN, handle)
81 :
82 1284 : SELECT CASE (shape_id)
83 : CASE (contour_shape_linear)
84 324 : IF (PRESENT(xnodes)) CALL rescale_nodes_linear(nnodes, tnodes, a, b, xnodes)
85 :
86 3220 : IF (PRESENT(weights)) weights(:) = b - a
87 :
88 : CASE (contour_shape_arc)
89 1908 : ALLOCATE (tnodes_angle(nnodes))
90 :
91 25388 : tnodes_angle(:) = tnodes(:)
92 636 : CALL rescale_nodes_pi_phi(a, b, nnodes, tnodes_angle)
93 :
94 636 : IF (PRESENT(xnodes)) CALL rescale_nodes_arc(nnodes, tnodes_angle, a, b, xnodes)
95 :
96 636 : IF (PRESENT(weights)) THEN
97 318 : rscale = (pi - get_arc_smallest_angle(a, b))*get_arc_radius(a, b)
98 :
99 12694 : DO i = 1, nnodes
100 12694 : weights(i) = rscale*CMPLX(SIN(tnodes_angle(i)), -COS(tnodes_angle(i)), kind=dp)
101 : END DO
102 : END IF
103 :
104 636 : DEALLOCATE (tnodes_angle)
105 : CASE DEFAULT
106 960 : CPABORT("Unimplemented integration shape")
107 : END SELECT
108 :
109 960 : CALL timestop(handle)
110 2400 : END SUBROUTINE rescale_normalised_nodes
111 :
112 : ! **************************************************************************************************
113 : !> \brief Compute arc radius.
114 : !> \param a lower bound
115 : !> \param b upper bound
116 : !> \return radius
117 : !> \par History
118 : !> * 05.2017 created [Sergey Chulkov]
119 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
120 : ! c *
121 : ! r * B-------+------
122 : ! a * / . |
123 : ! * r / . | delta
124 : ! * / phi . |
125 : ! A---------*-----------+------
126 : ! <--- r --><-l->
127 : ! <--- r --->
128 : ! **************************************************************************************************
129 636 : PURE FUNCTION get_arc_radius(a, b) RESULT(radius)
130 : COMPLEX(kind=dp), INTENT(in) :: a, b
131 : REAL(kind=dp) :: radius
132 :
133 : COMPLEX(kind=dp) :: b_minus_a
134 :
135 636 : b_minus_a = b - a
136 :
137 : ! l = REAL(B - A); delta = AIMAG(B - A)
138 : ! radius = (l^2 + delta^2) / (2 * l)
139 636 : radius = 0.5_dp*REAL(b_minus_a*CONJG(b_minus_a), kind=dp)/REAL(b_minus_a, kind=dp)
140 636 : END FUNCTION get_arc_radius
141 :
142 : ! **************************************************************************************************
143 : !> \brief Compute the angle phi.
144 : !> \param a lower bound
145 : !> \param b upper bound
146 : !> \return angle
147 : !> \par History
148 : !> * 05.2017 created [Sergey Chulkov]
149 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
150 : ! c *
151 : ! r * B-------+------
152 : ! a * / . |
153 : ! * r / . | delta
154 : ! * / phi . |
155 : ! A---------*-----------+------
156 : ! <--- r --><-l->
157 : ! <--- r --->
158 : ! **************************************************************************************************
159 954 : PURE FUNCTION get_arc_smallest_angle(a, b) RESULT(phi)
160 : COMPLEX(kind=dp), INTENT(in) :: a, b
161 : REAL(kind=dp) :: phi
162 :
163 : COMPLEX(kind=dp) :: b_minus_a
164 : REAL(kind=dp) :: delta2, l2
165 :
166 954 : b_minus_a = b - a
167 :
168 : ! l = REAL(B - A); delta = AIMAG(B - A)
169 : ! phi = arccos((l - radius)/radius) = arccos((l^2 - delta^2) / (l^2 + delta^2))
170 954 : l2 = REAL(b_minus_a, dp)
171 954 : l2 = l2*l2
172 954 : delta2 = AIMAG(b_minus_a)
173 954 : delta2 = delta2*delta2
174 :
175 954 : phi = ACOS((l2 - delta2)/(l2 + delta2))
176 954 : END FUNCTION get_arc_smallest_angle
177 :
178 0 : PURE FUNCTION get_axis_rotation_angle(a, b) RESULT(phi)
179 : COMPLEX(kind=dp), INTENT(in) :: a, b
180 : REAL(kind=dp) :: phi
181 :
182 : COMPLEX(kind=dp) :: b_minus_a
183 :
184 0 : b_minus_a = b - a
185 0 : phi = ACOS(REAL(b_minus_a, dp)/ABS(b_minus_a))
186 0 : END FUNCTION get_axis_rotation_angle
187 :
188 : ! **************************************************************************************************
189 : !> \brief Rescale nodes [pi, phi] -> arc[a, b] .
190 : !> \param nnodes number of nodes
191 : !> \param tnodes_angle parametrically-defined nodes to rescale
192 : !> \param a lower bound
193 : !> \param b upper bound
194 : !> \param xnodes rescaled nodes (initialised on exit)
195 : !> \par History
196 : !> * 05.2017 created [Sergey Chulkov]
197 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
198 : ! **************************************************************************************************
199 318 : SUBROUTINE rescale_nodes_arc(nnodes, tnodes_angle, a, b, xnodes)
200 : INTEGER, INTENT(in) :: nnodes
201 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: tnodes_angle
202 : COMPLEX(kind=dp), INTENT(in) :: a, b
203 : COMPLEX(kind=dp), DIMENSION(:), INTENT(out) :: xnodes
204 :
205 : COMPLEX(kind=dp) :: origin
206 : INTEGER :: i
207 : REAL(kind=dp) :: radius
208 :
209 318 : radius = get_arc_radius(a, b)
210 318 : origin = a + CMPLX(radius, 0.0_dp, kind=dp)
211 :
212 12694 : DO i = 1, nnodes
213 12694 : xnodes(i) = origin + radius*CMPLX(COS(tnodes_angle(i)), SIN(tnodes_angle(i)), kind=dp)
214 : END DO
215 318 : END SUBROUTINE rescale_nodes_arc
216 :
217 : ! **************************************************************************************************
218 : !> \brief Rescale nodes tnodes(i) = cos(pi/2 * (1-tnodes(i))); tnodes \in [-1 .. 1] .
219 : !> \param tnodes parametrically-defined nodes to rescale / rescaled nodes (modified on exit)
220 : !> \par History
221 : !> * 05.2017 created [Sergey Chulkov]
222 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
223 : ! **************************************************************************************************
224 0 : SUBROUTINE rescale_nodes_cos(nnodes, tnodes)
225 : INTEGER, INTENT(in) :: nnodes
226 : REAL(kind=dp), DIMENSION(nnodes), INTENT(inout) :: tnodes
227 :
228 0 : tnodes(:) = COS(0.5_dp*pi*(1.0_dp - tnodes(:)))
229 0 : END SUBROUTINE rescale_nodes_cos
230 :
231 : ! **************************************************************************************************
232 : !> \brief Rescale nodes [-1, 1] -> [a, b] .
233 : !> \param nnodes number of nodes
234 : !> \param tnodes parametrically-defined nodes to rescale
235 : !> \param a lower bound
236 : !> \param b upper bound
237 : !> \param xnodes rescaled nodes (initialised on exit)
238 : !> \par History
239 : !> * 05.2017 created [Sergey Chulkov]
240 : ! **************************************************************************************************
241 162 : SUBROUTINE rescale_nodes_linear(nnodes, tnodes, a, b, xnodes)
242 : INTEGER, INTENT(in) :: nnodes
243 : REAL(kind=dp), DIMENSION(nnodes), INTENT(in) :: tnodes
244 : COMPLEX(kind=dp), INTENT(in) :: a, b
245 : COMPLEX(kind=dp), DIMENSION(nnodes), INTENT(out) :: xnodes
246 :
247 : COMPLEX(kind=dp) :: half_len, median
248 :
249 162 : median = 0.5_dp*(b + a)
250 162 : half_len = 0.5_dp*(b - a)
251 :
252 3058 : xnodes(:) = median + half_len*tnodes(:)
253 162 : END SUBROUTINE rescale_nodes_linear
254 :
255 : ! **************************************************************************************************
256 : !> \brief Rescale nodes [-1, 1] -> [pi, phi] .
257 : !> \param nnodes number of nodes
258 : !> \param a lower bound
259 : !> \param b upper bound
260 : !> \param tnodes parametrically-defined nodes to rescale / rescaled nodes (modified on exit)
261 : !> \par History
262 : !> * 05.2017 created [Sergey Chulkov]
263 : !> \note Assuming Re(a) < Re(b) and Im(a) < Im(b)
264 : ! **************************************************************************************************
265 636 : SUBROUTINE rescale_nodes_pi_phi(a, b, nnodes, tnodes)
266 : COMPLEX(kind=dp), INTENT(in) :: a, b
267 : INTEGER, INTENT(in) :: nnodes
268 : REAL(kind=dp), DIMENSION(nnodes), INTENT(inout) :: tnodes
269 :
270 : REAL(kind=dp) :: half_pi_minus_phi, phi
271 :
272 636 : phi = get_arc_smallest_angle(a, b)
273 636 : half_pi_minus_phi = 0.5_dp*(pi - phi)
274 :
275 25388 : tnodes(:) = phi + half_pi_minus_phi*(1.0_dp - tnodes(:))
276 636 : END SUBROUTINE rescale_nodes_pi_phi
277 : END MODULE negf_integr_utils
|
