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 image charge corrections
10 : !> \par History
11 : !> 06.2019 Moved from rpa_ri_gpw [Frederick Stein]
12 : ! **************************************************************************************************
13 : MODULE 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 :
35 : CONTAINS
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 2 : 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 2 : 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 6 : TYPE(dbt_pgrid_type) :: pgrid_2d
79 38 : TYPE(dbt_type) :: t_3c_overl_nnP_ic_reflected_ctr, t_SinvVSinv, t_SinvVSinv_tmp
80 :
81 2 : CALL timeset(routineN, handle)
82 :
83 2 : IF (PRESENT(do_alpha)) THEN
84 0 : my_do_alpha = do_alpha
85 : ELSE
86 : my_do_alpha = .FALSE.
87 : END IF
88 :
89 2 : IF (PRESENT(do_beta)) THEN
90 0 : my_do_beta = do_beta
91 : ELSE
92 : my_do_beta = .FALSE.
93 : END IF
94 :
95 2 : do_closed_shell = .NOT. (my_do_alpha .OR. my_do_beta)
96 :
97 6 : ALLOCATE (Delta_Sigma_Neaton(gw_corr_lev_tot))
98 14 : Delta_Sigma_Neaton = 0.0_dp
99 :
100 2 : mo_start = homo - gw_corr_lev_occ + 1
101 2 : mo_end = homo + gw_corr_lev_virt
102 2 : CPASSERT(mo_end - mo_start + 1 == gw_corr_lev_tot)
103 :
104 6 : ALLOCATE (sizes_RI_split(dbt_nblks_total(t_3c_overl_nnP_ic_reflected, 1)))
105 2 : CALL dbt_get_info(t_3c_overl_nnP_ic_reflected, blk_size_1=sizes_RI_split)
106 :
107 2 : CALL dbt_create(mat_SinvVSinv, t_SinvVSinv_tmp)
108 2 : CALL dbt_copy_matrix_to_tensor(mat_SinvVSinv, t_SinvVSinv_tmp)
109 2 : pdims = 0
110 2 : CALL mp_dims_create(para_env%num_pe, pdims)
111 2 : 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 2 : DEALLOCATE (dist_1, dist_2)
115 2 : CALL dbt_pgrid_destroy(pgrid_2d)
116 :
117 2 : CALL dbt_copy(t_SinvVSinv_tmp, t_SinvVSinv)
118 2 : CALL dbt_destroy(t_SinvVSinv_tmp)
119 2 : 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 2 : map_1=[1], map_2=[2, 3])
125 :
126 6 : 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 8 : Delta_Sigma_Neaton(gw_corr_lev_occ + 1:) = -Delta_Sigma_Neaton(gw_corr_lev_occ + 1:)
129 :
130 2 : CALL dbt_destroy(t_SinvVSinv)
131 2 : CALL dbt_destroy(t_3c_overl_nnP_ic_reflected_ctr)
132 :
133 2 : IF (unit_nr > 0) THEN
134 :
135 1 : WRITE (unit_nr, *) ' '
136 :
137 1 : IF (do_closed_shell) THEN
138 1 : WRITE (unit_nr, '(T3,A)') 'Single-electron energies with image charge (ic) correction'
139 1 : WRITE (unit_nr, '(T3,A)') '----------------------------------------------------------'
140 0 : ELSE IF (my_do_alpha) THEN
141 0 : WRITE (unit_nr, '(T3,A)') 'Single-electron energies of alpha spins with image charge (ic) correction'
142 0 : WRITE (unit_nr, '(T3,A)') '-------------------------------------------------------------------------'
143 0 : ELSE IF (my_do_beta) THEN
144 0 : WRITE (unit_nr, '(T3,A)') 'Single-electron energies of beta spins with image charge (ic) correction'
145 0 : WRITE (unit_nr, '(T3,A)') '------------------------------------------------------------------------'
146 : END IF
147 :
148 1 : WRITE (unit_nr, *) ' '
149 1 : WRITE (unit_nr, '(T3,A)') 'Reference for the ic: Neaton et al., PRL 97, 216405 (2006)'
150 1 : WRITE (unit_nr, *) ' '
151 :
152 1 : WRITE (unit_nr, '(T3,A)') ' '
153 1 : WRITE (unit_nr, '(T14,2A)') 'MO E_n before ic corr Delta E_ic', &
154 2 : ' E_n after ic corr'
155 :
156 7 : DO n_level_gw = 1, gw_corr_lev_tot
157 6 : n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
158 6 : IF (n_level_gw <= gw_corr_lev_occ) THEN
159 3 : occ_virt = 'occ'
160 : ELSE
161 3 : occ_virt = 'vir'
162 : END IF
163 :
164 : WRITE (unit_nr, '(T4,I4,3A,3F21.3)') &
165 6 : n_level_gw_ref, ' ( ', occ_virt, ') ', &
166 6 : Eigenval(n_level_gw_ref)*evolt, &
167 6 : Delta_Sigma_Neaton(n_level_gw)*evolt, &
168 13 : (Eigenval(n_level_gw_ref) + Delta_Sigma_Neaton(n_level_gw))*evolt
169 :
170 : END DO
171 :
172 1 : IF (do_closed_shell) THEN
173 1 : WRITE (unit_nr, '(T3,A)') ' '
174 1 : 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 2 : Delta_Sigma_Neaton(gw_corr_lev_occ))*evolt
178 0 : ELSE IF (my_do_alpha) THEN
179 0 : WRITE (unit_nr, '(T3,A)') ' '
180 0 : 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 0 : Delta_Sigma_Neaton(gw_corr_lev_occ))*evolt
184 0 : ELSE IF (my_do_beta) THEN
185 0 : WRITE (unit_nr, '(T3,A)') ' '
186 0 : 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 0 : Delta_Sigma_Neaton(gw_corr_lev_occ))*evolt
190 : END IF
191 :
192 1 : IF (print_ic_values) THEN
193 :
194 0 : WRITE (unit_nr, '(T3,A)') ' '
195 0 : WRITE (unit_nr, '(T3,A)') 'Horizontal list for copying the image charge corrections for use as input:'
196 0 : WRITE (unit_nr, '(*(F7.3))') (Delta_Sigma_Neaton(n_level_gw)*evolt, &
197 0 : 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 14 : + Delta_Sigma_Neaton(1:gw_corr_lev_tot)
206 :
207 2 : CALL timestop(handle)
208 :
209 6 : 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 2 : 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 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: block_1, block_2
232 : REAL(KIND=dp), &
233 4 : DIMENSION(mo_bounds(2)-mo_bounds(1)+1) :: Delta_Sigma_Neaton_prv
234 : TYPE(dbt_iterator_type) :: iter
235 :
236 2 : CALL timeset(routineN, handle)
237 :
238 2 : CPASSERT(SIZE(Delta_Sigma_Neaton_prv) == SIZE(Delta_Sigma_Neaton))
239 14 : 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 2 : !$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 14 : Delta_Sigma_Neaton = Delta_Sigma_Neaton + Delta_Sigma_Neaton_prv
281 26 : CALL para_env%sum(Delta_Sigma_Neaton)
282 :
283 2 : CALL timestop(handle)
284 :
285 4 : 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 0 : 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 0 : CALL timeset(routineN, handle)
319 :
320 0 : IF (PRESENT(do_alpha)) THEN
321 0 : my_do_alpha = do_alpha
322 : ELSE
323 : my_do_alpha = .FALSE.
324 : END IF
325 :
326 0 : IF (PRESENT(do_beta)) THEN
327 0 : my_do_beta = do_beta
328 : ELSE
329 : my_do_beta = .FALSE.
330 : END IF
331 :
332 0 : do_closed_shell = .NOT. (my_do_alpha .OR. my_do_beta)
333 :
334 : ! check the number of input image charge corrected levels
335 0 : CPASSERT(SIZE(ic_corr_list) == gw_corr_lev_tot)
336 :
337 0 : IF (unit_nr > 0) THEN
338 :
339 0 : WRITE (unit_nr, *) ' '
340 :
341 0 : IF (do_closed_shell) THEN
342 0 : WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies with image charge (ic) correction'
343 0 : WRITE (unit_nr, '(T3,A)') '-----------------------------------------------------------'
344 0 : ELSE IF (my_do_alpha) THEN
345 0 : WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies of alpha spins with image charge (ic) correction'
346 0 : WRITE (unit_nr, '(T3,A)') '--------------------------------------------------------------------------'
347 0 : ELSE IF (my_do_beta) THEN
348 0 : WRITE (unit_nr, '(T3,A)') 'GW quasiparticle energies of beta spins with image charge (ic) correction'
349 0 : WRITE (unit_nr, '(T3,A)') '-------------------------------------------------------------------------'
350 : END IF
351 :
352 0 : WRITE (unit_nr, *) ' '
353 :
354 0 : DO n_level_gw = 1, gw_corr_lev_tot
355 0 : n_level_gw_ref = n_level_gw + homo - gw_corr_lev_occ
356 0 : IF (n_level_gw <= gw_corr_lev_occ) THEN
357 0 : occ_virt = 'occ'
358 : ELSE
359 0 : occ_virt = 'vir'
360 : END IF
361 :
362 : WRITE (unit_nr, '(T4,I4,3A,3F21.3)') &
363 0 : n_level_gw_ref, ' ( ', occ_virt, ') ', &
364 0 : Eigenval(n_level_gw_ref)*evolt, &
365 0 : ic_corr_list(n_level_gw)*evolt, &
366 0 : (Eigenval(n_level_gw_ref) + ic_corr_list(n_level_gw))*evolt
367 :
368 : END DO
369 :
370 0 : 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 0 : + 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 0 : + ic_corr_list(1:gw_corr_lev_tot)
381 :
382 0 : IF (unit_nr > 0) THEN
383 :
384 0 : IF (do_closed_shell) THEN
385 0 : WRITE (unit_nr, '(T3,A,F52.2)') 'G0W0 IC HOMO-LUMO gap (eV)', Eigenval(homo + 1) - Eigenval(homo)
386 0 : ELSE IF (my_do_alpha) THEN
387 0 : WRITE (unit_nr, '(T3,A,F46.2)') 'G0W0 Alpha IC HOMO-LUMO gap (eV)', Eigenval(homo + 1) - Eigenval(homo)
388 0 : ELSE IF (my_do_beta) THEN
389 0 : WRITE (unit_nr, '(T3,A,F47.2)') 'G0W0 Beta IC HOMO-LUMO gap (eV)', Eigenval(homo + 1) - Eigenval(homo)
390 : END IF
391 :
392 0 : 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 0 : 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 0 : DO n_level_gw = 1, gw_corr_lev_occ
404 0 : eigen_diff = eigen_diff + ic_corr_list(n_level_gw)
405 : END DO
406 0 : 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 0 : DO n_level_gw = 1, homo - gw_corr_lev_occ
410 0 : Eigenval(n_level_gw) = Eigenval(n_level_gw) + eigen_diff
411 0 : 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 0 : IF (gw_corr_lev_virt < nmo - homo .AND. gw_corr_lev_virt > 0) THEN
418 :
419 : ! calculate average IC correction for virtual orbitals
420 0 : eigen_diff = 0.0_dp
421 0 : DO n_level_gw = gw_corr_lev_occ + 1, gw_corr_lev_tot
422 0 : eigen_diff = eigen_diff + ic_corr_list(n_level_gw)
423 : END DO
424 0 : 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 0 : DO n_level_gw = homo + gw_corr_lev_virt + 1, nmo
428 0 : Eigenval(n_level_gw) = Eigenval(n_level_gw) + eigen_diff
429 0 : Eigenval_scf(n_level_gw) = Eigenval_scf(n_level_gw) + eigen_diff
430 : END DO
431 :
432 : END IF
433 :
434 0 : CALL timestop(handle)
435 :
436 0 : END SUBROUTINE apply_ic_corr
437 :
438 : END MODULE rpa_gw_ic
|
