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 Routines to calculate RI-RPA energy and Sigma correction to the RPA energies
10 : !> using the cubic spline based on eigen values of Q(w).
11 : !> \par History
12 : ! **************************************************************************************************
13 : MODULE rpa_sigma_functional
14 : USE cp_fm_diag, ONLY: choose_eigv_solver
15 : USE cp_fm_struct, ONLY: cp_fm_struct_type
16 : USE cp_fm_types, ONLY: cp_fm_create,&
17 : cp_fm_get_info,&
18 : cp_fm_release,&
19 : cp_fm_to_fm,&
20 : cp_fm_type
21 : USE input_constants, ONLY: sigma_PBE0_S1,&
22 : sigma_PBE0_S2,&
23 : sigma_PBE_S1,&
24 : sigma_PBE_S2,&
25 : sigma_none
26 : USE kinds, ONLY: dp
27 : USE machine, ONLY: m_flush
28 : USE mathconstants, ONLY: pi
29 : USE message_passing, ONLY: mp_comm_type,&
30 : mp_para_env_type
31 : #include "./base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_sigma_functional'
38 :
39 : PUBLIC :: rpa_sigma_matrix_spectral, rpa_sigma_create, rpa_sigma_type, finalize_rpa_sigma
40 :
41 : TYPE rpa_sigma_type
42 : PRIVATE
43 : REAL(KIND=dp) :: e_sigma_corr = 0.0_dp
44 : REAL(KIND=dp) :: e_rpa_by_eig_val = 0.0_dp
45 : INTEGER :: sigma_param = 0
46 : TYPE(cp_fm_type) :: mat_Q_diagonal = cp_fm_type()
47 : TYPE(cp_fm_type) :: fm_evec = cp_fm_type()
48 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: sigma_eigenvalue
49 : INTEGER :: dimen_RI_red = 0
50 : END TYPE rpa_sigma_type
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief ... Collect the Q(w) (fm_mat_Q) matrix to create rpa_sigma a derived type variable.
56 : !> and write out the choosen parametrization for the cubic spline.
57 : !> \param rpa_sigma ...
58 : !> \param sigma_param ...
59 : !> \param fm_mat_Q ...
60 : !> \param unit_nr ...
61 : !> \param para_env ...
62 : ! **************************************************************************************************
63 10 : SUBROUTINE rpa_sigma_create(rpa_sigma, sigma_param, fm_mat_Q, unit_nr, para_env)
64 :
65 : TYPE(rpa_sigma_type), INTENT(OUT) :: rpa_sigma
66 : INTEGER :: sigma_param
67 : TYPE(cp_fm_type) :: fm_mat_Q
68 : INTEGER :: unit_nr
69 :
70 : CLASS(mp_comm_type), INTENT(IN) :: para_env
71 :
72 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
73 :
74 : ! Getting information about the Q matrix and creating initializing two matrices to pass it to the diagonalising driver.
75 10 : CALL cp_fm_get_info(fm_mat_Q, matrix_struct=matrix_struct, nrow_global=rpa_sigma%dimen_RI_red)
76 :
77 30 : ALLOCATE (rpa_sigma%sigma_eigenvalue(rpa_sigma%dimen_RI_red))
78 :
79 10 : CALL cp_fm_create(rpa_sigma%fm_evec, matrix_struct)
80 10 : CALL cp_fm_create(rpa_sigma%mat_Q_diagonal, matrix_struct)
81 :
82 10 : rpa_sigma%sigma_param = sigma_param
83 :
84 0 : SELECT CASE (rpa_sigma%sigma_param)
85 :
86 : CASE (sigma_none)
87 : ! There is nothing to do
88 : CASE DEFAULT
89 0 : CPABORT("Unknown parameterization")
90 :
91 : CASE (sigma_PBE0_S1)
92 0 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A)") &
93 0 : "SIGMA_INFO| Sigma eigenvalues parameterized with PBE0_S1 reference"
94 :
95 : CASE (sigma_PBE0_S2)
96 10 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A)") &
97 5 : "SIGMA_INFO| Sigma eigenvalues parameterized with PBE0_S2 reference"
98 :
99 : CASE (sigma_PBE_S1)
100 0 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A)") &
101 0 : "SIGMA_INFO| Sigma eigenvalues parameterized with PBE_S1 reference"
102 :
103 : CASE (sigma_PBE_S2)
104 0 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3, A)") &
105 10 : "SIGMA_INFO| Sigma eigenvalues parameterized with PBE_S2 reference"
106 : END SELECT
107 10 : IF (unit_nr > 0) CALL m_flush(unit_nr)
108 10 : CALL para_env%sync()
109 :
110 10 : END SUBROUTINE rpa_sigma_create
111 :
112 : ! **************************************************************************************************
113 : !> \brief ... Memory cleanup routine
114 : !> \param rpa_sigma ...
115 : ! **************************************************************************************************
116 10 : SUBROUTINE rpa_sigma_cleanup(rpa_sigma)
117 :
118 : TYPE(rpa_sigma_type), INTENT(INOUT) :: rpa_sigma
119 :
120 10 : CALL cp_fm_release(rpa_sigma%mat_Q_diagonal)
121 10 : CALL cp_fm_release(rpa_sigma%fm_evec)
122 10 : DEALLOCATE (rpa_sigma%sigma_eigenvalue)
123 :
124 10 : END SUBROUTINE rpa_sigma_cleanup
125 :
126 : ! **************************************************************************************************
127 : !> \brief ... Diagonalize and store the eigenvalues of fm_mat_Q in rpa_sigma%sigma_eigenvalue.
128 : !> \param rpa_sigma ...
129 : !> \param fm_mat_Q ...
130 : !> \param wj ...
131 : !> \param para_env_RPA ...
132 : ! **************************************************************************************************
133 30 : SUBROUTINE rpa_sigma_matrix_spectral(rpa_sigma, fm_mat_Q, wj, para_env_RPA)
134 :
135 : TYPE(rpa_sigma_type) :: rpa_sigma
136 : TYPE(cp_fm_type) :: fm_mat_Q
137 : REAL(KIND=dp) :: wj
138 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA
139 :
140 : ! copy the Q matrix into the dummy matrix to avoid changing it.
141 30 : CALL cp_fm_to_fm(fm_mat_Q, rpa_sigma%mat_Q_diagonal)
142 :
143 : !diagonalising driver
144 30 : CALL choose_eigv_solver(rpa_sigma%mat_Q_diagonal, rpa_sigma%fm_evec, rpa_sigma%sigma_eigenvalue)
145 :
146 : ! Computing the integration to calculate the sigma correction.
147 30 : CALL compute_e_sigma_corr_by_freq_int(rpa_sigma, wj, para_env_RPA)
148 :
149 30 : END SUBROUTINE rpa_sigma_matrix_spectral
150 : ! **************************************************************************************************
151 :
152 : ! **************************************************************************************************
153 : !> \brief ... To compute the e_sigma_corr and e_rpa_by_eig_val by freq integration over the eigenvalues of Q(w)
154 : !> e_sigma_corr = - H(sigma) & e_rpa_by_eig_val = log(1+sigma)-sigma
155 : !> \param rpa_sigma ...
156 : !> \param wj ...
157 : !> \param para_env_RPA ...
158 : ! **************************************************************************************************
159 30 : SUBROUTINE compute_e_sigma_corr_by_freq_int(rpa_sigma, wj, para_env_RPA)
160 : TYPE(rpa_sigma_type), INTENT(INOUT) :: rpa_sigma
161 : REAL(KIND=dp), INTENT(IN) :: wj
162 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_RPA
163 :
164 : INTEGER :: iaux
165 : REAL(KIND=dp) :: dedw_rpa, dedw_sigma
166 :
167 30 : dedw_sigma = 0.0_dp
168 30 : dedw_rpa = 0.0_dp
169 :
170 : ! Loop which take each eigenvalue to: i) get E_RPA & ii) integrates the spline to get E_c correction
171 1692 : DO iaux = 1, rpa_sigma%dimen_RI_red
172 1692 : IF (rpa_sigma%sigma_eigenvalue(iaux) > 0.0_dp) THEN
173 1501 : IF (MODULO(iaux, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE
174 1501 : dedw_rpa = dedw_rpa + LOG(1.0_dp + rpa_sigma%sigma_eigenvalue(iaux)) - rpa_sigma%sigma_eigenvalue(iaux)
175 : IF (MODULO(iaux, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE
176 1501 : dedw_sigma = dedw_sigma - cubic_spline_integr(rpa_sigma%sigma_eigenvalue(iaux), rpa_sigma%sigma_param)
177 : END IF
178 : END DO
179 :
180 : ! (use 2.0_dp its better for compilers)
181 30 : rpa_sigma%e_sigma_corr = rpa_sigma%e_sigma_corr + (wj*dedw_sigma/(2.0_dp*pi*2.0_dp))
182 30 : rpa_sigma%e_rpa_by_eig_val = rpa_sigma%e_rpa_by_eig_val + (wj*dedw_rpa/(2.0_dp*pi*2.0_dp))
183 :
184 30 : END SUBROUTINE compute_e_sigma_corr_by_freq_int
185 :
186 : ! **************************************************************************************************
187 : !> \brief ... Save the calculated value of E_c correction to the global variable and memory clean.
188 : !> \param rpa_sigma ...
189 : !> \param unit_nr ...
190 : !> \param e_sigma_corr ...
191 : !> \param para_env ...
192 : !> \param do_minimax_quad ...
193 : ! **************************************************************************************************
194 10 : SUBROUTINE finalize_rpa_sigma(rpa_sigma, unit_nr, e_sigma_corr, para_env, do_minimax_quad)
195 : TYPE(rpa_sigma_type), INTENT(INOUT) :: rpa_sigma
196 : INTEGER :: unit_nr
197 : REAL(KIND=dp), INTENT(OUT) :: e_sigma_corr
198 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
199 : LOGICAL, INTENT(IN) :: do_minimax_quad
200 :
201 10 : IF (do_minimax_quad) rpa_sigma%e_rpa_by_eig_val = rpa_sigma%e_rpa_by_eig_val/2.0_dp
202 10 : CALL para_env%sum(rpa_sigma%e_rpa_by_eig_val)
203 10 : e_sigma_corr = rpa_sigma%e_sigma_corr
204 10 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') &
205 5 : 'RI-RPA energy from eigenvalues of Q(w) = ', &
206 10 : rpa_sigma%e_rpa_by_eig_val
207 :
208 10 : CALL rpa_sigma_cleanup(rpa_sigma)
209 :
210 10 : END SUBROUTINE finalize_rpa_sigma
211 :
212 : ! **************************************************************************************************
213 : !> \brief ... integrates cubic spline to get eigenvalue sigma based E_c correction.
214 : !> \param sigma ...
215 : !> \param sigma_param ...
216 : !> \return ... Integration value H(sigma) wrt to coupling constant eq 35 in Trushin et al JCP 2021
217 : ! **************************************************************************************************
218 1501 : FUNCTION cubic_spline_integr(sigma, sigma_param) RESULT(integral)
219 : REAL(KIND=dp), INTENT(in) :: sigma
220 : INTEGER :: sigma_param
221 : REAL(KIND=dp) :: integral
222 :
223 : INTEGER :: i, m, n
224 : REAL(KIND=dp) :: h
225 1501 : REAL(KIND=dp), ALLOCATABLE :: coeff(:, :), x_coord(:)
226 :
227 1501 : SELECT CASE (sigma_param)
228 : CASE (sigma_PBE0_S1)
229 0 : n = 21
230 0 : ALLOCATE (x_coord(n))
231 0 : ALLOCATE (coeff(4, n))
232 :
233 0 : coeff(1, 1) = 0.000000000000D+00
234 0 : coeff(1, 2) = 0.000000000000D+00
235 0 : coeff(1, 3) = -0.149500660756D-03
236 0 : coeff(1, 4) = -0.353017276233D-02
237 0 : coeff(1, 5) = -0.109810247734D-01
238 0 : coeff(1, 6) = -0.231246943777D-01
239 0 : coeff(1, 7) = -0.268999962858D-01
240 0 : coeff(1, 8) = -0.634751994007D-03
241 0 : coeff(1, 9) = 0.118792892470D-01
242 0 : coeff(1, 10) = -0.473431931326D-01
243 0 : coeff(1, 11) = -0.817589390539D-01
244 0 : coeff(1, 12) = 0.125726011069D-01
245 0 : coeff(1, 13) = 0.108028492092D+00
246 0 : coeff(1, 14) = 0.193548206759D+00
247 0 : coeff(1, 15) = 0.358395561305D-01
248 0 : coeff(1, 16) = -0.497714974829D-01
249 0 : coeff(1, 17) = 0.341059348835D-01
250 0 : coeff(1, 18) = 0.341050720155D-01
251 0 : coeff(1, 19) = 0.785549033229D-01
252 0 : coeff(1, 20) = 0.000000000000D+00
253 0 : coeff(1, 21) = 0.000000000000D+00
254 0 : coeff(2, 1) = 0.000000000000D+00
255 0 : coeff(2, 2) = 0.000000000000D+00
256 0 : coeff(2, 3) = -0.208376539581D+01
257 0 : coeff(2, 4) = -0.469755869285D+01
258 0 : coeff(2, 5) = -0.565503803415D+01
259 0 : coeff(2, 6) = -0.135502867642D+01
260 0 : coeff(2, 7) = 0.000000000000D+00
261 0 : coeff(2, 8) = 0.284340701746D+01
262 0 : coeff(2, 9) = 0.000000000000D+00
263 0 : coeff(2, 10) = -0.342695931351D+01
264 0 : coeff(2, 11) = 0.000000000000D+00
265 0 : coeff(2, 12) = 0.358739081268D+01
266 0 : coeff(2, 13) = 0.203368806130D+01
267 0 : coeff(2, 14) = 0.000000000000D+00
268 0 : coeff(2, 15) = -0.901387663218D+00
269 0 : coeff(2, 16) = 0.000000000000D+00
270 0 : coeff(2, 17) = 0.000000000000D+00
271 0 : coeff(2, 18) = 0.000000000000D+00
272 0 : coeff(2, 19) = 0.000000000000D+00
273 0 : coeff(2, 20) = 0.000000000000D+00
274 0 : coeff(2, 21) = 0.000000000000D+00
275 0 : coeff(3, 1) = -0.000000000000D+00
276 0 : coeff(3, 2) = -0.322176662524D+05
277 0 : coeff(3, 3) = -0.267090835643D+04
278 0 : coeff(3, 4) = -0.373532067350D+04
279 0 : coeff(3, 5) = -0.797121299000D+03
280 0 : coeff(3, 6) = 0.111299540119D+03
281 0 : coeff(3, 7) = 0.299284621116D+04
282 0 : coeff(3, 8) = -0.319333485618D+02
283 0 : coeff(3, 9) = -0.140910103454D+04
284 0 : coeff(3, 10) = -0.848330431187D+01
285 0 : coeff(3, 11) = 0.435025012278D+03
286 0 : coeff(3, 12) = -0.700327539634D+01
287 0 : coeff(3, 13) = 0.545486142353D+01
288 0 : coeff(3, 14) = -0.453346282407D+02
289 0 : coeff(3, 15) = 0.371921027910D+00
290 0 : coeff(3, 16) = 0.464101795796D+01
291 0 : coeff(3, 17) = -0.190069531714D-04
292 0 : coeff(3, 18) = 0.514345336660D-01
293 0 : coeff(3, 19) = -0.431543078188D-02
294 0 : coeff(3, 20) = -0.000000000000D+00
295 0 : coeff(3, 21) = 0.000000000000D+00
296 0 : coeff(4, 1) = 0.000000000000D+00
297 0 : coeff(4, 2) = 0.152897717268D+09
298 0 : coeff(4, 3) = 0.902815532735D+06
299 0 : coeff(4, 4) = 0.191760493084D+07
300 0 : coeff(4, 5) = 0.445372471512D+06
301 0 : coeff(4, 6) = 0.188362654331D+04
302 0 : coeff(4, 7) = -0.383203258784D+06
303 0 : coeff(4, 8) = -0.170027418959D+05
304 0 : coeff(4, 9) = 0.819629330224D+05
305 0 : coeff(4, 10) = 0.560228610945D+04
306 0 : coeff(4, 11) = -0.108203002413D+05
307 0 : coeff(4, 12) = -0.363378668069D+03
308 0 : coeff(4, 13) = -0.260332257619D+03
309 0 : coeff(4, 14) = 0.291068208088D+03
310 0 : coeff(4, 15) = 0.122322834276D+02
311 0 : coeff(4, 16) = -0.132875656470D+02
312 0 : coeff(4, 17) = 0.343356030115D-04
313 0 : coeff(4, 18) = -0.212958640167D-01
314 0 : coeff(4, 19) = 0.389311916174D-03
315 0 : coeff(4, 20) = 0.000000000000D+00
316 0 : coeff(4, 21) = 0.000000000000D+00
317 0 : x_coord(1) = 0.000000000000D+00
318 0 : x_coord(2) = 0.100000000000D-04
319 0 : x_coord(3) = 0.100000000000D-03
320 0 : x_coord(4) = 0.100000000000D-02
321 0 : x_coord(5) = 0.215443469000D-02
322 0 : x_coord(6) = 0.464158883000D-02
323 0 : x_coord(7) = 0.100000000000D-01
324 0 : x_coord(8) = 0.146779926762D-01
325 0 : x_coord(9) = 0.215443469003D-01
326 0 : x_coord(10) = 0.316227766017D-01
327 0 : x_coord(11) = 0.464158883361D-01
328 0 : x_coord(12) = 0.681292069058D-01
329 0 : x_coord(13) = 0.100000000000D+00
330 0 : x_coord(14) = 0.158489319246D+00
331 0 : x_coord(15) = 0.251188643151D+00
332 0 : x_coord(16) = 0.398107170553D+00
333 0 : x_coord(17) = 0.630957344480D+00
334 0 : x_coord(18) = 0.100000000000D+01
335 0 : x_coord(19) = 0.261015721568D+01
336 0 : x_coord(20) = 0.100000000000D+02
337 0 : x_coord(21) = 0.215443469000D+02
338 :
339 : CASE (sigma_PBE0_S2)
340 1501 : n = 21
341 1501 : ALLOCATE (x_coord(n))
342 1501 : ALLOCATE (coeff(4, n))
343 :
344 1501 : coeff(1, 1) = 0.000000000000D+00
345 1501 : coeff(1, 2) = 0.000000000000D+00
346 1501 : coeff(1, 3) = -0.431405252048D-04
347 1501 : coeff(1, 4) = -0.182874853131D-02
348 1501 : coeff(1, 5) = -0.852003132762D-02
349 1501 : coeff(1, 6) = -0.218177403992D-01
350 1501 : coeff(1, 7) = -0.305777654735D-01
351 1501 : coeff(1, 8) = -0.870882903969D-02
352 1501 : coeff(1, 9) = 0.137878988102D-01
353 1501 : coeff(1, 10) = -0.284352007440D-01
354 1501 : coeff(1, 11) = -0.798812002431D-01
355 1501 : coeff(1, 12) = -0.334010771574D-02
356 1501 : coeff(1, 13) = 0.934182748715D-01
357 1501 : coeff(1, 14) = 0.204960802253D+00
358 1501 : coeff(1, 15) = 0.213204380281D-01
359 1501 : coeff(1, 16) = -0.401220283037D-01
360 1501 : coeff(1, 17) = 0.321629738336D-01
361 1501 : coeff(1, 18) = 0.321618301891D-01
362 1501 : coeff(1, 19) = 0.808763912948D-01
363 1501 : coeff(1, 20) = 0.000000000000D+00
364 1501 : coeff(1, 21) = 0.000000000000D+00
365 1501 : coeff(2, 1) = 0.000000000000D+00
366 1501 : coeff(2, 2) = 0.000000000000D+00
367 1501 : coeff(2, 3) = -0.661870777583D+00
368 1501 : coeff(2, 4) = -0.289752912590D+01
369 1501 : coeff(2, 5) = -0.558979946652D+01
370 1501 : coeff(2, 6) = -0.267765704540D+01
371 1501 : coeff(2, 7) = 0.000000000000D+00
372 1501 : coeff(2, 8) = 0.389592612611D+01
373 1501 : coeff(2, 9) = 0.000000000000D+00
374 1501 : coeff(2, 10) = -0.382296397421D+01
375 1501 : coeff(2, 11) = 0.000000000000D+00
376 1501 : coeff(2, 12) = 0.327772498106D+01
377 1501 : coeff(2, 13) = 0.239633724310D+01
378 1501 : coeff(2, 14) = 0.000000000000D+00
379 1501 : coeff(2, 15) = -0.726304793204D+00
380 1501 : coeff(2, 16) = 0.000000000000D+00
381 1501 : coeff(2, 17) = 0.000000000000D+00
382 1501 : coeff(2, 18) = 0.000000000000D+00
383 1501 : coeff(2, 19) = 0.000000000000D+00
384 1501 : coeff(2, 20) = 0.000000000000D+00
385 1501 : coeff(2, 21) = 0.000000000000D+00
386 1501 : coeff(3, 1) = -0.000000000000D+00
387 1501 : coeff(3, 2) = -0.862385254713D+04
388 1501 : coeff(3, 3) = -0.192306222883D+04
389 1501 : coeff(3, 4) = -0.520047462362D+04
390 1501 : coeff(3, 5) = -0.877473657666D+03
391 1501 : coeff(3, 6) = 0.841408344046D+02
392 1501 : coeff(3, 7) = 0.216516760964D+04
393 1501 : coeff(3, 8) = 0.296702212913D+03
394 1501 : coeff(3, 9) = -0.867733655494D+03
395 1501 : coeff(3, 10) = -0.188410055380D+03
396 1501 : coeff(3, 11) = 0.336084151111D+03
397 1501 : coeff(3, 12) = 0.489746728744D+01
398 1501 : coeff(3, 13) = 0.158746877181D+02
399 1501 : coeff(3, 14) = -0.562764882273D+02
400 1501 : coeff(3, 15) = 0.134759277149D+01
401 1501 : coeff(3, 16) = 0.399959778866D+01
402 1501 : coeff(3, 17) = -0.251917983154D-04
403 1501 : coeff(3, 18) = 0.563694092760D-01
404 1501 : coeff(3, 19) = -0.444296223097D-02
405 1501 : coeff(3, 20) = -0.000000000000D+00
406 1501 : coeff(3, 21) = 0.000000000000D+00
407 1501 : coeff(4, 1) = 0.000000000000D+00
408 1501 : coeff(4, 2) = 0.366429086790D+08
409 1501 : coeff(4, 3) = 0.504466528222D+06
410 1501 : coeff(4, 4) = 0.232980923705D+07
411 1501 : coeff(4, 5) = 0.392124287301D+06
412 1501 : coeff(4, 6) = 0.206173887726D+05
413 1501 : coeff(4, 7) = -0.249217659838D+06
414 1501 : coeff(4, 8) = -0.563519876566D+05
415 1501 : coeff(4, 9) = 0.448530826095D+05
416 1501 : coeff(4, 10) = 0.143140667434D+05
417 1501 : coeff(4, 11) = -0.800144415404D+04
418 1501 : coeff(4, 12) = -0.391685311241D+03
419 1501 : coeff(4, 13) = -0.414433988077D+03
420 1501 : coeff(4, 14) = 0.376550449117D+03
421 1501 : coeff(4, 15) = 0.510124747789D+01
422 1501 : coeff(4, 16) = -0.114511339236D+02
423 1501 : coeff(4, 17) = 0.455083767664D-04
424 1501 : coeff(4, 18) = -0.233390912502D-01
425 1501 : coeff(4, 19) = 0.400817027790D-03
426 1501 : coeff(4, 20) = 0.000000000000D+00
427 1501 : coeff(4, 21) = 0.000000000000D+00
428 1501 : x_coord(1) = 0.000000000000D+00
429 1501 : x_coord(2) = 0.100000000000D-04
430 1501 : x_coord(3) = 0.100000000000D-03
431 1501 : x_coord(4) = 0.100000000000D-02
432 1501 : x_coord(5) = 0.215443469000D-02
433 1501 : x_coord(6) = 0.464158883000D-02
434 1501 : x_coord(7) = 0.100000000000D-01
435 1501 : x_coord(8) = 0.146779926762D-01
436 1501 : x_coord(9) = 0.215443469003D-01
437 1501 : x_coord(10) = 0.316227766017D-01
438 1501 : x_coord(11) = 0.464158883361D-01
439 1501 : x_coord(12) = 0.681292069058D-01
440 1501 : x_coord(13) = 0.100000000000D+00
441 1501 : x_coord(14) = 0.158489319246D+00
442 1501 : x_coord(15) = 0.251188643151D+00
443 1501 : x_coord(16) = 0.398107170553D+00
444 1501 : x_coord(17) = 0.630957344480D+00
445 1501 : x_coord(18) = 0.100000000000D+01
446 1501 : x_coord(19) = 0.261015721568D+01
447 1501 : x_coord(20) = 0.100000000000D+02
448 1501 : x_coord(21) = 0.215443469000D+02
449 :
450 : CASE (sigma_PBE_S1)
451 0 : n = 22
452 0 : ALLOCATE (x_coord(n))
453 0 : ALLOCATE (coeff(4, n))
454 :
455 0 : coeff(1, 1) = 0.000000000000D+00
456 0 : coeff(1, 2) = 0.000000000000D+00
457 0 : coeff(1, 3) = -0.493740326815D-04
458 0 : coeff(1, 4) = -0.136110637329D-02
459 0 : coeff(1, 5) = -0.506905111755D-02
460 0 : coeff(1, 6) = -0.127411222930D-01
461 0 : coeff(1, 7) = -0.220144968504D-01
462 0 : coeff(1, 8) = -0.239939034695D-01
463 0 : coeff(1, 9) = -0.436386416290D-01
464 0 : coeff(1, 10) = -0.117890214262D+00
465 0 : coeff(1, 11) = -0.141123921668D+00
466 0 : coeff(1, 12) = 0.865524876740D-01
467 0 : coeff(1, 13) = 0.179390274565D+00
468 0 : coeff(1, 14) = 0.269368658116D+00
469 0 : coeff(1, 15) = 0.785040456996D-01
470 0 : coeff(1, 16) = 0.490248637276D-01
471 0 : coeff(1, 17) = -0.111571911794D+00
472 0 : coeff(1, 18) = -0.197712184164D-01
473 0 : coeff(1, 19) = -0.197716870218D-01
474 0 : coeff(1, 20) = -0.372253617253D-01
475 0 : coeff(1, 21) = 0.000000000000D+00
476 0 : coeff(1, 22) = 0.000000000000D+00
477 0 : coeff(2, 1) = 0.000000000000D+00
478 0 : coeff(2, 2) = 0.000000000000D+00
479 0 : coeff(2, 3) = -0.709484897949D+00
480 0 : coeff(2, 4) = -0.197447407686D+01
481 0 : coeff(2, 5) = -0.315478745349D+01
482 0 : coeff(2, 6) = -0.229603163128D+01
483 0 : coeff(2, 7) = -0.670801534786D+00
484 0 : coeff(2, 8) = -0.704199644986D+00
485 0 : coeff(2, 9) = -0.400987325224D+01
486 0 : coeff(2, 10) = -0.269982990241D+01
487 0 : coeff(2, 11) = 0.000000000000D+00
488 0 : coeff(2, 12) = 0.472814414167D+01
489 0 : coeff(2, 13) = 0.207638470052D+01
490 0 : coeff(2, 14) = 0.000000000000D+00
491 0 : coeff(2, 15) = -0.389846972557D+00
492 0 : coeff(2, 16) = -0.298496119087D+00
493 0 : coeff(2, 17) = 0.000000000000D+00
494 0 : coeff(2, 18) = 0.000000000000D+00
495 0 : coeff(2, 19) = -0.601781536636D-06
496 0 : coeff(2, 20) = 0.000000000000D+00
497 0 : coeff(2, 21) = 0.000000000000D+00
498 0 : coeff(2, 22) = 0.000000000000D+00
499 0 : coeff(3, 1) = -0.000000000000D+00
500 0 : coeff(3, 2) = -0.104035132381D+05
501 0 : coeff(3, 3) = -0.108777473624D+04
502 0 : coeff(3, 4) = -0.219328637518D+04
503 0 : coeff(3, 5) = -0.260711341283D+03
504 0 : coeff(3, 6) = 0.132509852177D+02
505 0 : coeff(3, 7) = 0.165970301474D+03
506 0 : coeff(3, 8) = -0.460909893146D+03
507 0 : coeff(3, 9) = -0.112939707971D+04
508 0 : coeff(3, 10) = 0.465035067500D+02
509 0 : coeff(3, 11) = 0.123097490767D+04
510 0 : coeff(3, 12) = -0.876616265219D+02
511 0 : coeff(3, 13) = 0.790484996078D+01
512 0 : coeff(3, 14) = -0.624281400584D+02
513 0 : coeff(3, 15) = 0.324152775194D+01
514 0 : coeff(3, 16) = -0.632212496608D+01
515 0 : coeff(3, 17) = 0.202215332970D+01
516 0 : coeff(3, 18) = -0.308693235932D-06
517 0 : coeff(3, 19) = -0.495067060383D-02
518 0 : coeff(3, 20) = 0.116855980641D-02
519 0 : coeff(3, 21) = -0.000000000000D+00
520 0 : coeff(3, 22) = 0.000000000000D+00
521 0 : coeff(4, 1) = 0.000000000000D+00
522 0 : coeff(4, 2) = 0.478661516427D+08
523 0 : coeff(4, 3) = 0.285187385316D+06
524 0 : coeff(4, 4) = 0.971371823345D+06
525 0 : coeff(4, 5) = 0.116156741398D+06
526 0 : coeff(4, 6) = 0.172191903906D+05
527 0 : coeff(4, 7) = -0.241613612898D+05
528 0 : coeff(4, 8) = 0.213790845631D+05
529 0 : coeff(4, 9) = 0.790063233314D+05
530 0 : coeff(4, 10) = 0.201667888760D+04
531 0 : coeff(4, 11) = -0.344519214370D+05
532 0 : coeff(4, 12) = 0.963471669433D+03
533 0 : coeff(4, 13) = -0.292417702205D+03
534 0 : coeff(4, 14) = 0.433842720035D+03
535 0 : coeff(4, 15) = -0.132982468090D+02
536 0 : coeff(4, 16) = 0.199358142858D+02
537 0 : coeff(4, 17) = -0.365297127483D+01
538 0 : coeff(4, 18) = 0.434041376596D-07
539 0 : coeff(4, 19) = 0.101490424907D-02
540 0 : coeff(4, 20) = -0.796902275213D-04
541 0 : coeff(4, 21) = 0.000000000000D+00
542 0 : coeff(4, 22) = 0.000000000000D+00
543 0 : x_coord(1) = 0.000000000000D+00
544 0 : x_coord(2) = 0.100000000000D-04
545 0 : x_coord(3) = 0.100000000000D-03
546 0 : x_coord(4) = 0.100000000000D-02
547 0 : x_coord(5) = 0.215443469000D-02
548 0 : x_coord(6) = 0.464158883000D-02
549 0 : x_coord(7) = 0.100000000000D-01
550 0 : x_coord(8) = 0.146779926762D-01
551 0 : x_coord(9) = 0.215443469003D-01
552 0 : x_coord(10) = 0.316227766017D-01
553 0 : x_coord(11) = 0.464158883361D-01
554 0 : x_coord(12) = 0.681292069058D-01
555 0 : x_coord(13) = 0.100000000000D+00
556 0 : x_coord(14) = 0.158489319246D+00
557 0 : x_coord(15) = 0.251188643151D+00
558 0 : x_coord(16) = 0.398107170553D+00
559 0 : x_coord(17) = 0.630957344480D+00
560 0 : x_coord(18) = 0.100000000000D+01
561 0 : x_coord(19) = 0.237137370566D+01
562 0 : x_coord(20) = 0.562341325000D+01
563 0 : x_coord(21) = 0.153992652606D+02
564 0 : x_coord(22) = 0.316227766000D+02
565 :
566 : CASE (sigma_PBE_S2)
567 0 : n = 22
568 0 : ALLOCATE (x_coord(n))
569 0 : ALLOCATE (coeff(4, n))
570 :
571 0 : coeff(1, 1) = 0.000000000000D+00
572 0 : coeff(1, 2) = 0.000000000000D+00
573 0 : coeff(1, 3) = -0.156157535801D-03
574 0 : coeff(1, 4) = -0.365199003270D-02
575 0 : coeff(1, 5) = -0.108302033233D-01
576 0 : coeff(1, 6) = -0.203436953346D-01
577 0 : coeff(1, 7) = -0.214330355346D-01
578 0 : coeff(1, 8) = 0.109617244934D-03
579 0 : coeff(1, 9) = 0.813969827075D-02
580 0 : coeff(1, 10) = -0.701367130014D-01
581 0 : coeff(1, 11) = -0.162002361715D+00
582 0 : coeff(1, 12) = 0.337288711362D-01
583 0 : coeff(1, 13) = 0.140348429629D+00
584 0 : coeff(1, 14) = 0.271234417677D+00
585 0 : coeff(1, 15) = 0.780732751240D-01
586 0 : coeff(1, 16) = 0.436066976238D-01
587 0 : coeff(1, 17) = -0.106097689688D+00
588 0 : coeff(1, 18) = -0.133141637069D-01
589 0 : coeff(1, 19) = -0.133143525246D-01
590 0 : coeff(1, 20) = -0.430994711278D-01
591 0 : coeff(1, 21) = 0.000000000000D+00
592 0 : coeff(1, 22) = 0.000000000000D+00
593 0 : coeff(2, 1) = 0.000000000000D+00
594 0 : coeff(2, 2) = 0.000000000000D+00
595 0 : coeff(2, 3) = -0.217211651544D+01
596 0 : coeff(2, 4) = -0.473638379726D+01
597 0 : coeff(2, 5) = -0.487821808504D+01
598 0 : coeff(2, 6) = -0.433631413905D+00
599 0 : coeff(2, 7) = 0.000000000000D+00
600 0 : coeff(2, 8) = 0.193813387881D+01
601 0 : coeff(2, 9) = 0.000000000000D+00
602 0 : coeff(2, 10) = -0.695060290528D+01
603 0 : coeff(2, 11) = 0.000000000000D+00
604 0 : coeff(2, 12) = 0.502541925806D+01
605 0 : coeff(2, 13) = 0.273498669354D+01
606 0 : coeff(2, 14) = 0.000000000000D+00
607 0 : coeff(2, 15) = -0.448708826169D+00
608 0 : coeff(2, 16) = -0.332102918195D+00
609 0 : coeff(2, 17) = 0.000000000000D+00
610 0 : coeff(2, 18) = 0.000000000000D+00
611 0 : coeff(2, 19) = -0.242488141082D-06
612 0 : coeff(2, 20) = 0.000000000000D+00
613 0 : coeff(2, 21) = 0.000000000000D+00
614 0 : coeff(2, 22) = 0.000000000000D+00
615 0 : coeff(3, 1) = -0.000000000000D+00
616 0 : coeff(3, 2) = -0.337014964214D+05
617 0 : coeff(3, 3) = -0.285795351280D+04
618 0 : coeff(3, 4) = -0.372723918347D+04
619 0 : coeff(3, 5) = -0.516689374427D+03
620 0 : coeff(3, 6) = 0.480322803175D+02
621 0 : coeff(3, 7) = 0.253894893657D+04
622 0 : coeff(3, 8) = -0.535684993409D+02
623 0 : coeff(3, 9) = -0.162223464755D+04
624 0 : coeff(3, 10) = -0.319667723139D+03
625 0 : coeff(3, 11) = 0.101401359817D+04
626 0 : coeff(3, 12) = -0.862770702569D+02
627 0 : coeff(3, 13) = 0.212578002151D+02
628 0 : coeff(3, 14) = -0.625949163782D+02
629 0 : coeff(3, 15) = 0.357838438707D+01
630 0 : coeff(3, 16) = -0.543078279308D+01
631 0 : coeff(3, 17) = 0.204380282001D+01
632 0 : coeff(3, 18) = -0.124376927880D-06
633 0 : coeff(3, 19) = -0.844892173480D-02
634 0 : coeff(3, 20) = 0.135295689023D-02
635 0 : coeff(3, 21) = -0.000000000000D+00
636 0 : coeff(3, 22) = 0.000000000000D+00
637 0 : coeff(4, 1) = 0.000000000000D+00
638 0 : coeff(4, 2) = 0.160253203309D+09
639 0 : coeff(4, 3) = 0.106174857663D+07
640 0 : coeff(4, 4) = 0.211694319531D+07
641 0 : coeff(4, 5) = 0.377995031873D+06
642 0 : coeff(4, 6) = -0.941771032564D+03
643 0 : coeff(4, 7) = -0.332306990184D+06
644 0 : coeff(4, 8) = -0.850176312292D+04
645 0 : coeff(4, 9) = 0.844978829514D+05
646 0 : coeff(4, 10) = 0.249933770676D+05
647 0 : coeff(4, 11) = -0.275803550484D+05
648 0 : coeff(4, 12) = 0.105308484492D+04
649 0 : coeff(4, 13) = -0.508788318128D+03
650 0 : coeff(4, 14) = 0.432758845019D+03
651 0 : coeff(4, 15) = -0.144367801061D+02
652 0 : coeff(4, 16) = 0.175904487062D+02
653 0 : coeff(4, 17) = -0.369208055752D+01
654 0 : coeff(4, 18) = 0.174842962107D-07
655 0 : coeff(4, 19) = 0.173203285755D-02
656 0 : coeff(4, 20) = -0.922652326542D-04
657 0 : coeff(4, 21) = 0.000000000000D+00
658 0 : coeff(4, 22) = 0.000000000000D+00
659 0 : x_coord(1) = 0.000000000000D+00
660 0 : x_coord(2) = 0.100000000000D-04
661 0 : x_coord(3) = 0.100000000000D-03
662 0 : x_coord(4) = 0.100000000000D-02
663 0 : x_coord(5) = 0.215443469000D-02
664 0 : x_coord(6) = 0.464158883000D-02
665 0 : x_coord(7) = 0.100000000000D-01
666 0 : x_coord(8) = 0.146779926762D-01
667 0 : x_coord(9) = 0.215443469003D-01
668 0 : x_coord(10) = 0.316227766017D-01
669 0 : x_coord(11) = 0.464158883361D-01
670 0 : x_coord(12) = 0.681292069058D-01
671 0 : x_coord(13) = 0.100000000000D+00
672 0 : x_coord(14) = 0.158489319246D+00
673 0 : x_coord(15) = 0.251188643151D+00
674 0 : x_coord(16) = 0.398107170553D+00
675 0 : x_coord(17) = 0.630957344480D+00
676 0 : x_coord(18) = 0.100000000000D+01
677 0 : x_coord(19) = 0.237137370566D+01
678 0 : x_coord(20) = 0.562341325000D+01
679 0 : x_coord(21) = 0.153992652606D+02
680 1501 : x_coord(22) = 0.316227766000D+02
681 :
682 : END SELECT
683 :
684 : ! determine to which interval sigma eigenvalue belongs
685 1501 : m = intervalnum(x_coord, n, sigma)
686 :
687 : ! Numerically evaluate integral
688 1501 : integral = 0.0_dp
689 1501 : IF (m == 1) THEN
690 378 : integral = 0.5_dp*coeff(2, 1)*sigma
691 : END IF
692 1501 : IF ((m > 1) .AND. (m < n)) THEN
693 1123 : h = sigma - x_coord(m)
694 : integral = 0.5_dp*coeff(2, 1)*x_coord(2)**2/sigma &
695 : + (coeff(1, m)*h + coeff(2, m)/2.0_dp*h**2 + &
696 1123 : coeff(3, m)/3.0_dp*h**3 + coeff(4, m)/4.0_dp*h**4)/sigma
697 8042 : DO i = 2, m - 1
698 6919 : h = x_coord(i + 1) - x_coord(i)
699 : integral = integral &
700 : + (coeff(1, i)*h + coeff(2, i)/2.0_dp*h**2 + &
701 8042 : coeff(3, i)/3.0_dp*h**3 + coeff(4, i)/4.0_dp*h**4)/sigma
702 : END DO
703 : END IF
704 1501 : IF (m == n) THEN
705 0 : integral = 0.5_dp*coeff(2, 1)*x_coord(2)**2/sigma
706 0 : DO i = 2, m - 1
707 0 : h = x_coord(i + 1) - x_coord(i)
708 : integral = integral &
709 : + (coeff(1, i)*h + coeff(2, i)/2.0_dp*h**2 &
710 0 : + coeff(3, i)/3.0_dp*h**3 + coeff(4, i)/4.0_dp*h**4)/sigma
711 : END DO
712 0 : integral = integral + coeff(1, n)*(1.0_dp - x_coord(n)/sigma)
713 : END IF
714 1501 : integral = integral*sigma
715 :
716 1501 : DEALLOCATE (x_coord, coeff)
717 :
718 1501 : END FUNCTION cubic_spline_integr
719 :
720 : ! **************************************************************************************************
721 : !> \brief ... Determine the interval which contains the sigma eigenvalue.
722 : !> \param x_coord ...
723 : !> \param n ...
724 : !> \param sigma ...
725 : !> \return ... an integer m
726 : ! **************************************************************************************************
727 1501 : INTEGER FUNCTION intervalnum(x_coord, n, sigma) RESULT(inum)
728 :
729 : INTEGER :: n
730 : REAL(KIND=dp), INTENT(in) :: x_coord(n), sigma
731 :
732 : INTEGER :: i
733 :
734 1501 : IF (sigma <= 0.0_dp) CPABORT('intervalnum: sigma should be positive')
735 :
736 1501 : inum = -1
737 1501 : IF (sigma > x_coord(n)) inum = n
738 31521 : DO i = 1, n - 1
739 31521 : IF ((sigma > x_coord(i)) .AND. (sigma <= x_coord(i + 1))) inum = i
740 : END DO
741 :
742 1501 : IF (inum == -1) CPABORT('interval: something was wrong')
743 :
744 1501 : END FUNCTION intervalnum
745 :
746 0 : END MODULE rpa_sigma_functional
|
