(git:ccc2433)
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-2024 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 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  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 
438 END MODULE rpa_gw_ic
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