(git:ed6f26b)
Loading...
Searching...
No Matches
rpa_gw_ic.F
Go to the documentation of this file.
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 image charge corrections
10!> \par History
11!> 06.2019 Moved from rpa_ri_gpw [Frederick Stein]
12! **************************************************************************************************
13MODULE rpa_gw_ic
14 USE cp_dbcsr_api, ONLY: dbcsr_type
15 USE dbt_api, ONLY: &
16 dbt_contract, dbt_copy, dbt_copy_matrix_to_tensor, dbt_create, dbt_destroy, dbt_get_block, &
17 dbt_get_info, dbt_iterator_blocks_left, dbt_iterator_next_block, dbt_iterator_start, &
18 dbt_iterator_stop, dbt_iterator_type, dbt_nblks_total, dbt_pgrid_create, &
19 dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
20 USE kinds, ONLY: dp
21 USE message_passing, ONLY: mp_dims_create,&
22 mp_para_env_type
23 USE physcon, ONLY: evolt
24 USE qs_tensors_types, ONLY: create_2c_tensor
25#include "./base/base_uses.f90"
26
27 IMPLICIT NONE
28
29 PRIVATE
30
31 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_ic'
32
33 PUBLIC :: calculate_ic_correction, apply_ic_corr
34
35CONTAINS
36
37! **************************************************************************************************
38!> \brief ...
39!> \param Eigenval ...
40!> \param mat_SinvVSinv ...
41!> \param t_3c_overl_nnP_ic ...
42!> \param t_3c_overl_nnP_ic_reflected ...
43!> \param gw_corr_lev_tot ...
44!> \param gw_corr_lev_occ ...
45!> \param gw_corr_lev_virt ...
46!> \param homo ...
47!> \param unit_nr ...
48!> \param print_ic_values ...
49!> \param para_env ...
50!> \param do_alpha ...
51!> \param do_beta ...
52! **************************************************************************************************
53 SUBROUTINE calculate_ic_correction(Eigenval, mat_SinvVSinv, &
54 t_3c_overl_nnP_ic, t_3c_overl_nnP_ic_reflected, gw_corr_lev_tot, &
55 gw_corr_lev_occ, gw_corr_lev_virt, homo, unit_nr, &
56 print_ic_values, para_env, &
57 do_alpha, do_beta)
58
59 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: eigenval
60 TYPE(dbcsr_type), INTENT(IN), TARGET :: mat_sinvvsinv
61 TYPE(dbt_type) :: t_3c_overl_nnp_ic, &
62 t_3c_overl_nnp_ic_reflected
63 INTEGER, INTENT(IN) :: gw_corr_lev_tot, gw_corr_lev_occ, &
64 gw_corr_lev_virt, homo, unit_nr
65 LOGICAL, INTENT(IN) :: print_ic_values
66 TYPE(mp_para_env_type), INTENT(IN) :: para_env
67 LOGICAL, INTENT(IN), OPTIONAL :: do_alpha, do_beta
68
69 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_ic_correction'
70
71 CHARACTER(4) :: occ_virt
72 INTEGER :: handle, mo_end, mo_start, n_level_gw, &
73 n_level_gw_ref
74 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_1, dist_2, sizes_ri_split
75 INTEGER, DIMENSION(2) :: pdims
76 LOGICAL :: do_closed_shell, my_do_alpha, my_do_beta
77 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: delta_sigma_neaton
78 TYPE(dbt_pgrid_type) :: pgrid_2d
79 TYPE(dbt_type) :: t_3c_overl_nnp_ic_reflected_ctr, t_sinvvsinv, t_sinvvsinv_tmp
80
81 CALL timeset(routinen, handle)
82
83 IF (PRESENT(do_alpha)) THEN
84 my_do_alpha = do_alpha
85 ELSE
86 my_do_alpha = .false.
87 END IF
88
89 IF (PRESENT(do_beta)) THEN
90 my_do_beta = do_beta
91 ELSE
92 my_do_beta = .false.
93 END IF
94
95 do_closed_shell = .NOT. (my_do_alpha .OR. my_do_beta)
96
97 ALLOCATE (delta_sigma_neaton(gw_corr_lev_tot))
98 delta_sigma_neaton = 0.0_dp
99
100 mo_start = homo - gw_corr_lev_occ + 1
101 mo_end = homo + gw_corr_lev_virt
102 cpassert(mo_end - mo_start + 1 == gw_corr_lev_tot)
103
104 ALLOCATE (sizes_ri_split(dbt_nblks_total(t_3c_overl_nnp_ic_reflected, 1)))
105 CALL dbt_get_info(t_3c_overl_nnp_ic_reflected, blk_size_1=sizes_ri_split)
106
107 CALL dbt_create(mat_sinvvsinv, t_sinvvsinv_tmp)
108 CALL dbt_copy_matrix_to_tensor(mat_sinvvsinv, t_sinvvsinv_tmp)
109 pdims = 0
110 CALL mp_dims_create(para_env%num_pe, pdims)
111 CALL dbt_pgrid_create(para_env, pdims, pgrid_2d)
112 CALL create_2c_tensor(t_sinvvsinv, dist_1, dist_2, pgrid_2d, sizes_ri_split, sizes_ri_split, &
113 name="(RI|RI)")
114 DEALLOCATE (dist_1, dist_2)
115 CALL dbt_pgrid_destroy(pgrid_2d)
116
117 CALL dbt_copy(t_sinvvsinv_tmp, t_sinvvsinv)
118 CALL dbt_destroy(t_sinvvsinv_tmp)
119 CALL dbt_create(t_3c_overl_nnp_ic_reflected, t_3c_overl_nnp_ic_reflected_ctr)
120 CALL dbt_contract(0.5_dp, t_sinvvsinv, t_3c_overl_nnp_ic_reflected, &
121 0.0_dp, t_3c_overl_nnp_ic_reflected_ctr, &
122 contract_1=[2], notcontract_1=[1], &
123 contract_2=[1], notcontract_2=[2, 3], &
124 map_1=[1], map_2=[2, 3])
125
126 CALL trace_ic_gw(t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected_ctr, delta_sigma_neaton, [mo_start, mo_end], para_env)
127
128 delta_sigma_neaton(gw_corr_lev_occ + 1:) = -delta_sigma_neaton(gw_corr_lev_occ + 1:)
129
130 CALL dbt_destroy(t_sinvvsinv)
131 CALL dbt_destroy(t_3c_overl_nnp_ic_reflected_ctr)
132
133 IF (unit_nr > 0) THEN
134
135 WRITE (unit_nr, *) ' '
136
137 IF (do_closed_shell) THEN
138 WRITE (unit_nr, '(T3,A)') 'Single-electron energies with image charge (ic) correction'
139 WRITE (unit_nr, '(T3,A)') '----------------------------------------------------------'
140 ELSE IF (my_do_alpha) THEN
141 WRITE (unit_nr, '(T3,A)') 'Single-electron energies of alpha spins with image charge (ic) correction'
142 WRITE (unit_nr, '(T3,A)') '-------------------------------------------------------------------------'
143 ELSE IF (my_do_beta) THEN
144 WRITE (unit_nr, '(T3,A)') 'Single-electron energies of beta spins with image charge (ic) correction'
145 WRITE (unit_nr, '(T3,A)') '------------------------------------------------------------------------'
146 END IF
147
148 WRITE (unit_nr, *) ' '
149 WRITE (unit_nr, '(T3,A)') 'Reference for the ic: Neaton et al., PRL 97, 216405 (2006)'
150 WRITE (unit_nr, *) ' '
151
152 WRITE (unit_nr, '(T3,A)') ' '
153 WRITE (unit_nr, '(T14,2A)') 'MO E_n before ic corr Delta E_ic', &
154 ' E_n after ic corr'
155
156 DO n_level_gw = 1, gw_corr_lev_tot
157 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
158 IF (n_level_gw <= gw_corr_lev_occ) THEN
159 occ_virt = 'occ'
160 ELSE
161 occ_virt = 'vir'
162 END IF
163
164 WRITE (unit_nr, '(T4,I4,3A,3F21.3)') &
165 n_level_gw_ref, ' ( ', occ_virt, ') ', &
166 eigenval(n_level_gw_ref)*evolt, &
167 delta_sigma_neaton(n_level_gw)*evolt, &
168 (eigenval(n_level_gw_ref) + delta_sigma_neaton(n_level_gw))*evolt
169
170 END DO
171
172 IF (do_closed_shell) THEN
173 WRITE (unit_nr, '(T3,A)') ' '
174 WRITE (unit_nr, '(T3,A,F57.2)') 'IC HOMO-LUMO gap (eV)', (eigenval(homo + 1) + &
175 delta_sigma_neaton(gw_corr_lev_occ + 1) - &
176 eigenval(homo) - &
177 delta_sigma_neaton(gw_corr_lev_occ))*evolt
178 ELSE IF (my_do_alpha) THEN
179 WRITE (unit_nr, '(T3,A)') ' '
180 WRITE (unit_nr, '(T3,A,F51.2)') 'Alpha IC HOMO-LUMO gap (eV)', (eigenval(homo + 1) + &
181 delta_sigma_neaton(gw_corr_lev_occ + 1) - &
182 eigenval(homo) - &
183 delta_sigma_neaton(gw_corr_lev_occ))*evolt
184 ELSE IF (my_do_beta) THEN
185 WRITE (unit_nr, '(T3,A)') ' '
186 WRITE (unit_nr, '(T3,A,F52.2)') 'Beta IC HOMO-LUMO gap (eV)', (eigenval(homo + 1) + &
187 delta_sigma_neaton(gw_corr_lev_occ + 1) - &
188 eigenval(homo) - &
189 delta_sigma_neaton(gw_corr_lev_occ))*evolt
190 END IF
191
192 IF (print_ic_values) THEN
193
194 WRITE (unit_nr, '(T3,A)') ' '
195 WRITE (unit_nr, '(T3,A)') 'Horizontal list for copying the image charge corrections for use as input:'
196 WRITE (unit_nr, '(*(F7.3))') (delta_sigma_neaton(n_level_gw)*evolt, &
197 n_level_gw=1, gw_corr_lev_tot)
198
199 END IF
200
201 END IF
202
203 eigenval(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt) = eigenval(homo - gw_corr_lev_occ + 1: &
204 homo + gw_corr_lev_virt) &
205 + delta_sigma_neaton(1:gw_corr_lev_tot)
206
207 CALL timestop(handle)
208
209 END SUBROUTINE calculate_ic_correction
210
211! **************************************************************************************************
212!> \brief ...
213!> \param t3c_1 ...
214!> \param t3c_2 ...
215!> \param Delta_Sigma_Neaton ...
216!> \param mo_bounds ...
217!> \param para_env ...
218! **************************************************************************************************
219 SUBROUTINE trace_ic_gw(t3c_1, t3c_2, Delta_Sigma_Neaton, mo_bounds, para_env)
220 TYPE(dbt_type), INTENT(INOUT) :: t3c_1, t3c_2
221 REAL(dp), DIMENSION(:), INTENT(INOUT) :: delta_sigma_neaton
222 INTEGER, DIMENSION(2), INTENT(IN) :: mo_bounds
223 TYPE(mp_para_env_type), INTENT(IN) :: para_env
224
225 CHARACTER(LEN=*), PARAMETER :: routinen = 'trace_ic_gw'
226
227 INTEGER :: handle, n, n_end, n_end_block, n_start, &
228 n_start_block
229 INTEGER, DIMENSION(3) :: boff, bsize, ind
230 LOGICAL :: found
231 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: block_1, block_2
232 REAL(kind=dp), &
233 DIMENSION(mo_bounds(2)-mo_bounds(1)+1) :: delta_sigma_neaton_prv
234 TYPE(dbt_iterator_type) :: iter
235
236 CALL timeset(routinen, handle)
237
238 cpassert(SIZE(delta_sigma_neaton_prv) == SIZE(delta_sigma_neaton))
239 delta_sigma_neaton_prv = 0.0_dp
240
241!$OMP PARALLEL DEFAULT(NONE) REDUCTION(+:Delta_Sigma_Neaton_prv) &
242!$OMP SHARED(t3c_1,t3c_2,mo_bounds) &
243!$OMP PRIVATE(iter,ind,bsize,boff,block_1,block_2,found) &
244!$OMP PRIVATE(n,n_start_block,n_start,n_end_block,n_end)
245 CALL dbt_iterator_start(iter, t3c_1)
246 DO WHILE (dbt_iterator_blocks_left(iter))
247 CALL dbt_iterator_next_block(iter, ind, blk_size=bsize, blk_offset=boff)
248 IF (ind(2) /= ind(3)) cycle
249 CALL dbt_get_block(t3c_1, ind, block_1, found)
250 cpassert(found)
251 CALL dbt_get_block(t3c_2, ind, block_2, found)
252 IF (.NOT. found) cycle
253
254 IF (boff(3) < mo_bounds(1)) THEN
255 n_start_block = mo_bounds(1) - boff(3) + 1
256 n_start = 1
257 ELSE
258 n_start_block = 1
259 n_start = boff(3) - mo_bounds(1) + 1
260 END IF
261
262 IF (boff(3) + bsize(3) - 1 > mo_bounds(2)) THEN
263 n_end_block = mo_bounds(2) - boff(3) + 1
264 n_end = mo_bounds(2) - mo_bounds(1) + 1
265 ELSE
266 n_end_block = bsize(3)
267 n_end = boff(3) + bsize(3) - mo_bounds(1)
268 END IF
269
270 delta_sigma_neaton_prv(n_start:n_end) = &
271 delta_sigma_neaton_prv(n_start:n_end) + &
272 (/(dot_product(block_1(:, n, n), &
273 block_2(:, n, n)), &
274 n=n_start_block, n_end_block)/)
275 DEALLOCATE (block_1, block_2)
276 END DO
277 CALL dbt_iterator_stop(iter)
278!$OMP END PARALLEL
279
280 delta_sigma_neaton = delta_sigma_neaton + delta_sigma_neaton_prv
281 CALL para_env%sum(delta_sigma_neaton)
282
283 CALL timestop(handle)
284
285 END SUBROUTINE
286
287! **************************************************************************************************
288!> \brief ...
289!> \param Eigenval ...
290!> \param Eigenval_scf ...
291!> \param ic_corr_list ...
292!> \param gw_corr_lev_occ ...
293!> \param gw_corr_lev_virt ...
294!> \param gw_corr_lev_tot ...
295!> \param homo ...
296!> \param nmo ...
297!> \param unit_nr ...
298!> \param do_alpha ...
299!> \param do_beta ...
300! **************************************************************************************************
301 SUBROUTINE apply_ic_corr(Eigenval, Eigenval_scf, ic_corr_list, &
302 gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, &
303 homo, nmo, unit_nr, do_alpha, do_beta)
304
305 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: eigenval, eigenval_scf
306 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: ic_corr_list
307 INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, &
308 gw_corr_lev_tot, homo, nmo, unit_nr
309 LOGICAL, INTENT(IN), OPTIONAL :: do_alpha, do_beta
310
311 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_ic_corr'
312
313 CHARACTER(4) :: occ_virt
314 INTEGER :: handle, n_level_gw, n_level_gw_ref
315 LOGICAL :: do_closed_shell, my_do_alpha, my_do_beta
316 REAL(kind=dp) :: eigen_diff
317
318 CALL timeset(routinen, handle)
319
320 IF (PRESENT(do_alpha)) THEN
321 my_do_alpha = do_alpha
322 ELSE
323 my_do_alpha = .false.
324 END IF
325
326 IF (PRESENT(do_beta)) THEN
327 my_do_beta = do_beta
328 ELSE
329 my_do_beta = .false.
330 END IF
331
332 do_closed_shell = .NOT. (my_do_alpha .OR. my_do_beta)
333
334 ! check the number of input image charge corrected levels
335 cpassert(SIZE(ic_corr_list) == gw_corr_lev_tot)
336
337 IF (unit_nr > 0) THEN
338
339 WRITE (unit_nr, *) ' '
340
341 IF (do_closed_shell) THEN
342 WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies with image charge (ic) correction'
343 WRITE (unit_nr, '(T3,A)') '-----------------------------------------------------------'
344 ELSE IF (my_do_alpha) THEN
345 WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies of alpha spins with image charge (ic) correction'
346 WRITE (unit_nr, '(T3,A)') '--------------------------------------------------------------------------'
347 ELSE IF (my_do_beta) THEN
348 WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies of beta spins with image charge (ic) correction'
349 WRITE (unit_nr, '(T3,A)') '-------------------------------------------------------------------------'
350 END IF
351
352 WRITE (unit_nr, *) ' '
353
354 DO n_level_gw = 1, gw_corr_lev_tot
355 n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
356 IF (n_level_gw <= gw_corr_lev_occ) THEN
357 occ_virt = 'occ'
358 ELSE
359 occ_virt = 'vir'
360 END IF
361
362 WRITE (unit_nr, '(T4,I4,3A,3F21.3)') &
363 n_level_gw_ref, ' ( ', occ_virt, ') ', &
364 eigenval(n_level_gw_ref)*evolt, &
365 ic_corr_list(n_level_gw)*evolt, &
366 (eigenval(n_level_gw_ref) + ic_corr_list(n_level_gw))*evolt
367
368 END DO
369
370 WRITE (unit_nr, *) ' '
371
372 END IF
373
374 eigenval(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt) = eigenval(homo - gw_corr_lev_occ + 1: &
375 homo + gw_corr_lev_virt) &
376 + ic_corr_list(1:gw_corr_lev_tot)
377
378 eigenval_scf(homo - gw_corr_lev_occ + 1:homo + gw_corr_lev_virt) = eigenval_scf(homo - gw_corr_lev_occ + 1: &
379 homo + gw_corr_lev_virt) &
380 + ic_corr_list(1:gw_corr_lev_tot)
381
382 IF (unit_nr > 0) THEN
383
384 IF (do_closed_shell) THEN
385 WRITE (unit_nr, '(T3,A,F52.2)') 'G0W0 IC HOMO-LUMO gap (eV)', eigenval(homo + 1) - eigenval(homo)
386 ELSE IF (my_do_alpha) THEN
387 WRITE (unit_nr, '(T3,A,F46.2)') 'G0W0 Alpha IC HOMO-LUMO gap (eV)', eigenval(homo + 1) - eigenval(homo)
388 ELSE IF (my_do_beta) THEN
389 WRITE (unit_nr, '(T3,A,F47.2)') 'G0W0 Beta IC HOMO-LUMO gap (eV)', eigenval(homo + 1) - eigenval(homo)
390 END IF
391
392 WRITE (unit_nr, *) ' '
393
394 END IF
395
396 ! for eigenvalue self-consistent GW, all eigenvalues have to be corrected
397 ! 1) the occupied; check if there are occupied MOs not being corrected by the IC model
398 IF (gw_corr_lev_occ < homo .AND. gw_corr_lev_occ > 0) THEN
399
400 ! calculate average IC contribution for occupied orbitals
401 eigen_diff = 0.0_dp
402
403 DO n_level_gw = 1, gw_corr_lev_occ
404 eigen_diff = eigen_diff + ic_corr_list(n_level_gw)
405 END DO
406 eigen_diff = eigen_diff/gw_corr_lev_occ
407
408 ! correct the eigenvalues of the occupied orbitals which have not been corrected by the IC model
409 DO n_level_gw = 1, homo - gw_corr_lev_occ
410 eigenval(n_level_gw) = eigenval(n_level_gw) + eigen_diff
411 eigenval_scf(n_level_gw) = eigenval_scf(n_level_gw) + eigen_diff
412 END DO
413
414 END IF
415
416 ! 2) the virtual: check if there are virtual orbitals not being corrected by the IC model
417 IF (gw_corr_lev_virt < nmo - homo .AND. gw_corr_lev_virt > 0) THEN
418
419 ! calculate average IC correction for virtual orbitals
420 eigen_diff = 0.0_dp
421 DO n_level_gw = gw_corr_lev_occ + 1, gw_corr_lev_tot
422 eigen_diff = eigen_diff + ic_corr_list(n_level_gw)
423 END DO
424 eigen_diff = eigen_diff/gw_corr_lev_virt
425
426 ! correct the eigenvalues of the virtual orbitals which have not been corrected by the IC model
427 DO n_level_gw = homo + gw_corr_lev_virt + 1, nmo
428 eigenval(n_level_gw) = eigenval(n_level_gw) + eigen_diff
429 eigenval_scf(n_level_gw) = eigenval_scf(n_level_gw) + eigen_diff
430 END DO
431
432 END IF
433
434 CALL timestop(handle)
435
436 END SUBROUTINE apply_ic_corr
437
438END MODULE rpa_gw_ic
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_dims_create
subroutine, public mp_dims_create(nodes, dims)
wrapper to MPI_Dims_create
Definition message_passing.F:2132
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_tensors_types
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors_types.F:12
qs_tensors_types::create_2c_tensor
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
Definition qs_tensors_types.F:382
rpa_gw_ic
Routines to calculate image charge corrections.
Definition rpa_gw_ic.F:13
rpa_gw_ic::calculate_ic_correction
subroutine, public calculate_ic_correction(eigenval, mat_sinvvsinv, t_3c_overl_nnp_ic, t_3c_overl_nnp_ic_reflected, gw_corr_lev_tot, gw_corr_lev_occ, gw_corr_lev_virt, homo, unit_nr, print_ic_values, para_env, do_alpha, do_beta)
...
Definition rpa_gw_ic.F:58
rpa_gw_ic::apply_ic_corr
subroutine, public apply_ic_corr(eigenval, eigenval_scf, ic_corr_list, gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, homo, nmo, unit_nr, do_alpha, do_beta)
...
Definition rpa_gw_ic.F:304
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763