Octopus
xc_sic.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!! Copyright (C) 2022 N. Tancogne-Dejean
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
22module xc_sic_oct_m
23 use debug_oct_m
25 use global_oct_m
26 use grid_oct_m
33 use parser_oct_m
36 use space_oct_m
40 use xc_oct_m
41 use xc_f03_lib_m
42 use xc_oep_oct_m
43
44 implicit none
45
46 private
47 public :: &
48 xc_sic_t, &
54
56 integer, parameter, public :: &
57 SIC_NONE = 1, & !< no self-interaction correction
58 sic_pz_oep = 2, &
59 sic_amaldi = 3, &
60 sic_adsic = 4
61
63 type xc_sic_t
64 private
65 integer, public :: level = sic_none
66 real(real64), public :: amaldi_factor
67 type(xc_oep_t), public :: oep
68 end type xc_sic_t
69
70contains
71
72 ! ---------------------------------------------------------
74 !
75 subroutine xc_sic_init(sic, namespace, gr, st, mc, space)
76 type(xc_sic_t), intent(out) :: sic
77 type(namespace_t), intent(in) :: namespace
78 type(grid_t), intent(inout) :: gr
79 type(states_elec_t), intent(in) :: st
80 type(multicomm_t), intent(in) :: mc
81 class(space_t), intent(in) :: space
82
83
84 push_sub(xc_sic_init)
85
86 !%Variable SICCorrection
87 !%Type integer
88 !%Default sic_none
89 !%Section Hamiltonian::XC
90 !%Description
91 !% This variable controls which form of self-interaction correction to use. Note that
92 !% this correction will be applied to the functional chosen by <tt>XCFunctional</tt>.
93 !%Option sic_none 1
94 !% No self-interaction correction.
95 !%Option sic_pz 2
96 !% Perdew-Zunger SIC, handled by the OEP technique.
97 !% J. P. Perdew and Alex Zunger, Phys. Rev. B 23, 5048 (1981)
98 !% Extension to the spinor case follows Tancogne-Dejean et al., J. Chem. Phys. 159, 224110 (2023)
99 !%Option sic_amaldi 3
100 !% Amaldi correction term. Not implemeneted for spinors.
101 !% E. Fermi and E. Amaldi, Mem. Reale Accad. Italia 6, 119 (1934)
102 !%Option sic_adsic 4
103 !% Average-density SIC.
104 !% C. Legrand <i>et al.</i>, <i>J. Phys. B</i> <b>35</b>, 1115 (2002).
105 !% Extension to the spinor case follows Tancogne-Dejean et al., J. Chem. Phys. 159, 224110 (2023)
106 !%End
107 call parse_variable(namespace, 'SICCorrection', sic_none, sic%level)
108 if (.not. varinfo_valid_option('SICCorrection', sic%level)) call messages_input_error(namespace, 'SICCorrection')
109
110 ! check whether we should introduce the Amaldi SIC correction
111 sic%amaldi_factor = m_one
112 if (sic%level == sic_amaldi) then
113 sic%amaldi_factor = (st%qtot - m_one)/st%qtot
114 if(st%d%ispin == spinors) then
115 call messages_not_implemented("Amaldi SIC with non-collinear spins")
116 end if
117 end if
118
119 if(sic%level == sic_pz_oep) then
120 call xc_oep_init(sic%oep, namespace, gr, st, mc, space, oep_type = oep_type_sic)
121
122 if(st%nik > st%d%spin_channels) then
123 call messages_not_implemented("PZ-SIC with k-points")
124 end if
125 end if
126
127 if (allocated(st%rho_core)) then
128 call messages_not_implemented('SIC with nonlinear core corrections')
129 end if
130
131 if (allocated(st%frozen_rho) .and. (sic%level == sic_pz_oep .or. sic%level == sic_amaldi)) then
132 call messages_not_implemented('PZ-SIC with frozen orbitals')
133 end if
134
135 if (space%is_periodic() .and. sic%level /= sic_none) then
136 call messages_not_implemented("SIC corrections in periodic systems")
137 end if
138
139 pop_sub(xc_sic_init)
140 end subroutine xc_sic_init
141
142 ! ---------------------------------------------------------
144 subroutine xc_sic_end(sic)
145 type(xc_sic_t), intent(inout) :: sic
146
147 if (sic%level == sic_none) return
148
149 push_sub(xc_sic_end)
150
151 if(sic%level == sic_pz_oep) call xc_oep_end(sic%oep)
152
153 pop_sub(xc_sic_end)
154 end subroutine xc_sic_end
155
156
157 ! ---------------------------------------------------------
158 subroutine xc_sic_write_info(sic, iunit, namespace)
159 type(xc_sic_t), intent(in) :: sic
160 integer, optional, intent(in) :: iunit
161 type(namespace_t), optional, intent(in) :: namespace
163 if (sic%level == sic_none) return
164
165 push_sub(xc_sic_write_info)
166
167 call messages_print_var_option('SICCorrection', sic%level, iunit=iunit, namespace=namespace)
168
169 pop_sub(xc_sic_write_info)
170 end subroutine xc_sic_write_info
171
172 ! ---------------------------------------------------------
189 subroutine xc_sic_calc_adsic(sic, namespace, space, gr, st, hm, xc, density, vxc, ex, ec)
190 type(xc_sic_t), intent(in) :: sic
191 type(namespace_t), intent(in) :: namespace
192 class(space_t), intent(in) :: space
193 type(grid_t), intent(in) :: gr
194 type(states_elec_t), intent(in) :: st
195 type(hamiltonian_elec_t), intent(in) :: hm
196 type(xc_t), intent(inout) :: xc
197 real(real64), contiguous, intent(in) :: density(:,:)
198 real(real64), contiguous, intent(inout) :: vxc(:,:)
199 real(real64), optional, intent(inout) :: ex, ec
200
201 integer :: ispin, ist, ik, ip
202 real(real64), allocatable :: vxc_sic(:,:), vh_sic(:), rho(:, :)
203 real(real64) :: ex_sic, ec_sic, qsp(2)
204 real(real64) :: dtot, dpol, vpol
205 real(real64) :: nup
206
207 push_sub(xc_sic_calc_adsic)
208
209 assert(sic%level == sic_adsic)
210
211 if (st%d%ispin == spinors .and. .not. in_family(hm%xc%family, [xc_family_lda, xc_family_gga])) then
212 write(message(1),'(a)') 'ADSIC with non-collinear spin is currently only possible'
213 write(message(2),'(a)') 'with LDA and GGA functionals.'
214 call messages_fatal(2, namespace=namespace)
215 end if
216
217 if (xc_is_not_size_consistent(xc, namespace)) then
218 call messages_not_implemented('ADSIC with size inconsistent functionals', namespace=namespace)
219 end if
220
221 ! We compute here the number of electrons per spin channel
222 qsp = m_zero
223 if( .not. allocated(st%frozen_rho)) then
224 select case (st%d%ispin)
226 do ist = 1, st%nst
227 do ik = 1, st%nik
228 ispin = st%d%get_spin_index(ik)
229 qsp(ispin) = qsp(ispin) + st%occ(ist, ik) * st%kweights(ik)
230 end do
231 end do
232 end select
233 else
234 ! In the case of the frozen density, we can only get the charge from the integral
235 ! of the total density, including valence and frozen density
236 qsp(1:st%d%nspin) = dmf_integrate(gr, st%d%nspin, density)
237 end if
238
239 safe_allocate(vxc_sic(1:gr%np, 1:2))
240 safe_allocate(vh_sic(1:gr%np))
241 safe_allocate(rho(1:gr%np, 1:2))
242 ! We first compute the average xc self-interction error and we substract it
243 select case (st%d%ispin)
245 do ispin = 1, st%d%nspin
246 if (abs(qsp(ispin)) <= m_min_occ) cycle
247
248 rho = m_zero
249 vxc_sic = m_zero
250
251 rho(:, ispin) = density(:, ispin) / qsp(ispin)
252 if(present(ex) .and. present(ec)) then
253 ex_sic = m_zero
254 ec_sic = m_zero
255 ! This needs always to be called for the spin-polarized case
256 call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
257 rho, spin_polarized, hm%ions%latt%rcell_volume, vxc_sic, ex = ex_sic, ec = ec_sic)
258 ex = ex - ex_sic * qsp(ispin)
259 ec = ec - ec_sic * qsp(ispin)
260 else
261 ! This needs always to be called for the spin-polarized case
262 call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
263 rho, spin_polarized, hm%ions%latt%rcell_volume, vxc_sic)
264 end if
265
266 call lalg_axpy(gr%np, -m_one, vxc_sic(:, ispin), vxc(:, ispin))
267
268 ! We now substract the averaged Hartree self-interaction error
269 ! See Eq. 15 in [Pietezak and Vieira, Theoretical Chemistry Accounts (2021) 140:130]
270 vh_sic = m_zero
271 call dpoisson_solve(hm%psolver, namespace, vh_sic, rho(:, ispin), all_nodes=.false.)
272 call lalg_axpy(gr%np, -m_one, vh_sic, vxc(:, ispin))
273
274 ! Compute the corresponding energy contribution
275 if(present(ex)) then
276 ex = ex - m_half*dmf_dotp(gr, rho(:,ispin), vh_sic) * qsp(ispin)
277 end if
278
279 end do
280
281 case (spinors)
282 ! Here we only treat the case of LDA/GGA. We rotate the average density in the local frame
283 ! And we then compute the SIC correction from it
284 ! This cannot excerce any xc torque, by construction
285 assert(in_family(hm%xc%family, [xc_family_lda, xc_family_gga]))
286
287 do ispin = 1, 2
288 rho = m_zero
289 vxc_sic = m_zero
290 ! Averaged density in the local frame
291 do ip = 1, gr%np
292 dtot = density(ip, 1) + density(ip, 2)
293 dpol = sqrt((density(ip, 1) - density(ip, 2))**2 + &
294 m_four*(density(ip, 3)**2 + density(ip, 4)**2))
295 if(ispin == 1) then
296 rho(ip, 1) = max(m_half*(dtot + dpol), m_zero)
297 else
298 rho(ip, 2) = max(m_half*(dtot - dpol), m_zero)
299 end if
300 end do
301 nup = dmf_integrate(gr, rho(:,ispin))
302 if (nup <= 1e-14_real64) cycle
303 call lalg_scal(gr%np, m_one/nup, rho(:,ispin))
304
305 ! This needs always to be called for the spin-polarized case
306 if(present(ex) .and. present(ec)) then
307 ex_sic = m_zero
308 ec_sic = m_zero
309 call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
310 rho, spin_polarized, hm%ions%latt%rcell_volume, vxc_sic, ex = ex_sic, ec = ec_sic)
311 ex = ex - ex_sic * nup
312 ec = ec - ec_sic * nup
313 else
314 call xc_get_vxc(gr, xc, st, hm%kpoints, hm%psolver, namespace, space, &
315 rho, spin_polarized, hm%ions%latt%rcell_volume, vxc_sic)
316 end if
317
318 ! Select only the potential correspond to this spin channel
319 if(ispin == 2) then
320 vxc_sic(:, 1) = m_zero
321 else
322 vxc_sic(:, 2) = m_zero
323 end if
324
325 vh_sic = m_zero
326 call dpoisson_solve(hm%psolver, namespace, vh_sic, rho(:, ispin), all_nodes=.false.)
327 call lalg_axpy(gr%np, m_one, vh_sic, vxc_sic(:, ispin))
328 ! Compute the corresponding energy contribution
329 if(present(ex)) then
330 ex = ex - m_half*dmf_dotp(gr, rho(:,ispin), vh_sic) * nup
331 end if
332
333 do ip = 1, gr%np
334 dpol = sqrt((density(ip, 1) - density(ip, 2))**2 + &
335 m_four*(density(ip, 3)**2 + density(ip, 4)**2))
336 vpol = (vxc_sic(ip, 1) - vxc_sic(ip, 2))*(density(ip, 1) - density(ip, 2))/(safe_tol(dpol, xc_tiny))
337
338 vxc(ip, 1) = vxc(ip, 1) - m_half*(vxc_sic(ip, 1) + vxc_sic(ip, 2) + vpol)
339 vxc(ip, 2) = vxc(ip, 2) - m_half*(vxc_sic(ip, 1) + vxc_sic(ip, 2) - vpol)
340 vxc(ip, 3) = vxc(ip, 3) - (vxc_sic(ip, 1) - vxc_sic(ip, 2))*density(ip, 3)/(safe_tol(dpol, xc_tiny))
341 vxc(ip, 4) = vxc(ip, 4) - (vxc_sic(ip, 1) - vxc_sic(ip, 2))*density(ip, 4)/(safe_tol(dpol, xc_tiny))
342 end do
343 end do
344
345
346 end select
347
348 safe_deallocate_a(vxc_sic)
349 safe_deallocate_a(vh_sic)
350 safe_deallocate_a(rho)
351
352 pop_sub(xc_sic_calc_adsic)
353 end subroutine xc_sic_calc_adsic
354
355 ! ---------------------------------------------------------
359 subroutine xc_sic_add_fxc_adsic(namespace, xc, st, gr, rho, fxc)
360 type(namespace_t), intent(in) :: namespace
361 type(xc_t), intent(in) :: xc
362 type(states_elec_t), intent(in) :: st
363 type(grid_t), intent(in) :: gr
364 real(real64), intent(in) :: rho(:,:)
365 real(real64), contiguous, intent(inout) :: fxc(:,:,:)
366
367 real(real64), allocatable :: rho_averaged(:, :)
368 real(real64), allocatable :: fxc_sic(:,:,:)
369 real(real64) :: qtot
370 integer :: ispin
371
372 push_sub(xc_sic_add_fxc_adsic)
373
374 !Check spin and triplets
375 assert(st%d%ispin /= spinors)
376 assert(.not. allocated(st%frozen_rho))
377
378 if (bitand(xc%kernel_family, xc_family_lda) == 0) then
379 message(1) = "fxc calculation with ADSIC not implemented beyond LDAs."
380 call messages_fatal(1, namespace=namespace)
381 end if
382
383 if (xc_is_not_size_consistent(xc, namespace)) then
384 call messages_not_implemented('ADSIC with size inconsistent functionals', namespace=namespace)
385 end if
386
387 if (st%d%ispin == spinors) then
388 call messages_not_implemented('ADSIC fxc with non-collinear spin')
389 end if
390
391 ! This needs always to be called for the spin-polarized case
392 safe_allocate(rho_averaged(1:gr%np, 1:2))
393 safe_allocate(fxc_sic(1:gr%np, 1:2, 1:2))
394
395 do ispin = 1, st%d%nspin
396 rho_averaged = m_zero
397 qtot = dmf_integrate(gr, rho(:, ispin))
398 if (abs(qtot) <= m_min_occ) cycle
399
400 call lalg_copy(gr%np, rho(:, ispin), rho_averaged(:, ispin))
401 call lalg_scal(gr%np, m_one/qtot, rho_averaged(:, ispin))
402
403 fxc_sic = m_zero
404 call xc_get_fxc(xc, gr, namespace, rho_averaged, spin_polarized, fxc_sic)
405
406 call lalg_axpy(gr%np, -m_one/qtot, fxc_sic(:, ispin, ispin), fxc(:, ispin, ispin))
407 end do
408
409 safe_deallocate_a(rho_averaged)
410 safe_deallocate_a(fxc_sic)
411
412 pop_sub(xc_sic_add_fxc_adsic)
413 end subroutine xc_sic_add_fxc_adsic
414
415end module xc_sic_oct_m
416
417!! Local Variables:
418!! mode: f90
419!! coding: utf-8
420!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:173
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
scales a vector by a constant
Definition: lalg_basic.F90:159
double sqrt(double __x) __attribute__((__nothrow__
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_four
Definition: global.F90:195
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public m_one
Definition: global.F90:192
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
Definition: global.F90:214
This module implements the underlying real-space grid.
Definition: grid.F90:119
This module defines various routines, operating on mesh functions.
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:691
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
subroutine, public dpoisson_solve(this, namespace, pot, rho, all_nodes, kernel, reset)
Calculates the Poisson equation. Given the density returns the corresponding potential.
Definition: poisson.F90:871
This module handles spin dimensions of the states and the k-point distribution.
Definition: xc.F90:116
real(real64), parameter, public xc_tiny
Arbitrary definition of tiny, for use in XC context.
Definition: xc.F90:216
subroutine, public xc_get_vxc(gr, xcs, st, kpoints, psolver, namespace, space, rho, ispin, rcell_volume, vxc, ex, ec, deltaxc, vtau, ex_density, ec_density, stress_xc, force_orbitalfree)
Definition: xc.F90:772
logical function, public xc_is_not_size_consistent(xcs, namespace)
Is one of the x or c functional is not size consistent.
Definition: xc.F90:736
subroutine, public xc_get_fxc(xcs, gr, namespace, rho, ispin, fxc, fxc_grad, fxc_grad_spin)
Returns the exchange-correlation kernel.
Definition: xc.F90:2002
pure logical function, public in_family(family, xc_families)
Definition: xc.F90:621
subroutine, public xc_oep_end(oep)
Definition: xc_oep.F90:358
subroutine, public xc_oep_init(oep, namespace, gr, st, mc, space, oep_type)
Definition: xc_oep.F90:219
integer, parameter, public oep_type_sic
Definition: xc_oep.F90:186
subroutine, public xc_sic_write_info(sic, iunit, namespace)
Definition: xc_sic.F90:254
integer, parameter, public sic_adsic
Averaged density SIC.
Definition: xc_sic.F90:151
subroutine, public xc_sic_init(sic, namespace, gr, st, mc, space)
initialize the SIC object
Definition: xc_sic.F90:171
subroutine, public xc_sic_end(sic)
finalize the SIC and, if needed, the included OEP
Definition: xc_sic.F90:240
integer, parameter, public sic_pz_oep
Perdew-Zunger SIC (OEP way)
Definition: xc_sic.F90:151
integer, parameter, public sic_amaldi
Amaldi correction term.
Definition: xc_sic.F90:151
subroutine, public xc_sic_calc_adsic(sic, namespace, space, gr, st, hm, xc, density, vxc, ex, ec)
Computes the ADSIC potential and energy.
Definition: xc_sic.F90:285
subroutine, public xc_sic_add_fxc_adsic(namespace, xc, st, gr, rho, fxc)
Adds to fxc the ADSIC contribution.
Definition: xc_sic.F90:455
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
The states_elec_t class contains all electronic wave functions.
This class contains information about the self-interaction correction.
Definition: xc_sic.F90:158