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, &
53
55 integer, parameter, public :: &
56 SIC_NONE = 1, & !< no self-interaction correction
57 sic_pz_oep = 2, &
58 sic_amaldi = 3, &
59 sic_adsic = 4
60
62 type xc_sic_t
63 private
64 integer, public :: level = sic_none
65 real(real64), public :: amaldi_factor
66 type(xc_oep_t), public :: oep
67 end type xc_sic_t
68
69contains
70
71 ! ---------------------------------------------------------
73 !
74 subroutine xc_sic_init(sic, namespace, gr, st, mc, space)
75 type(xc_sic_t), intent(out) :: sic
76 type(namespace_t), intent(in) :: namespace
77 type(grid_t), intent(inout) :: gr
78 type(states_elec_t), intent(in) :: st
79 type(multicomm_t), intent(in) :: mc
80 class(space_t), intent(in) :: space
81
82
83 push_sub(xc_sic_init)
84
85 !%Variable SICCorrection
86 !%Type integer
87 !%Default sic_none
88 !%Section Hamiltonian::XC
89 !%Description
90 !% This variable controls which form of self-interaction correction to use. Note that
91 !% this correction will be applied to the functional chosen by <tt>XCFunctional</tt>.
92 !%Option sic_none 1
93 !% No self-interaction correction.
94 !%Option sic_pz 2
95 !% Perdew-Zunger SIC, handled by the OEP technique.
96 !% J. P. Perdew and Alex Zunger, Phys. Rev. B 23, 5048 (1981)
97 !% Extension to the spinor case follows Tancogne-Dejean et al., J. Chem. Phys. 159, 224110 (2023)
98 !%Option sic_amaldi 3
99 !% Amaldi correction term. Not implemeneted for spinors.
100 !% E. Fermi and E. Amaldi, Mem. Reale Accad. Italia 6, 119 (1934)
101 !%Option sic_adsic 4
102 !% Average-density SIC.
103 !% C. Legrand <i>et al.</i>, <i>J. Phys. B</i> <b>35</b>, 1115 (2002).
104 !% Extension to the spinor case follows Tancogne-Dejean et al., J. Chem. Phys. 159, 224110 (2023)
105 !%End
106 call parse_variable(namespace, 'SICCorrection', sic_none, sic%level)
107 if (.not. varinfo_valid_option('SICCorrection', sic%level)) call messages_input_error(namespace, 'SICCorrection')
108
109 ! check whether we should introduce the Amaldi SIC correction
110 sic%amaldi_factor = m_one
111 if (sic%level == sic_amaldi) then
112 sic%amaldi_factor = (st%qtot - m_one)/st%qtot
113 if(st%d%ispin == spinors) then
114 call messages_not_implemented("Amaldi SIC with non-collinear spins")
115 end if
116 end if
118 if(sic%level == sic_pz_oep) then
119 call xc_oep_init(sic%oep, namespace, gr, st, mc, space, oep_type = oep_type_sic)
120
121 if(st%nik > st%d%spin_channels) then
122 call messages_not_implemented("PZ-SIC with k-points")
123 end if
124 end if
125
126 if (allocated(st%rho_core)) then
127 call messages_not_implemented('SIC with nonlinear core corrections')
128 end if
129
130 if (allocated(st%frozen_rho) .and. (sic%level == sic_pz_oep .or. sic%level == sic_amaldi)) then
131 call messages_not_implemented('PZ-SIC with frozen orbitals')
132 end if
133
134 if (space%is_periodic() .and. sic%level /= sic_none) then
135 call messages_not_implemented("SIC corrections in periodic systems")
136 end if
137
138 pop_sub(xc_sic_init)
139 end subroutine xc_sic_init
140
141 ! ---------------------------------------------------------
143 subroutine xc_sic_end(sic)
144 type(xc_sic_t), intent(inout) :: sic
145
146 if (sic%level == sic_none) return
147
148 push_sub(xc_sic_end)
149
150 if(sic%level == sic_pz_oep) call xc_oep_end(sic%oep)
151
152 pop_sub(xc_sic_end)
153 end subroutine xc_sic_end
154
155
156 ! ---------------------------------------------------------
157 subroutine xc_sic_write_info(sic, iunit, namespace)
158 type(xc_sic_t), intent(in) :: sic
159 integer, optional, intent(in) :: iunit
160 type(namespace_t), optional, intent(in) :: namespace
162 if (sic%level == sic_none) return
163
164 push_sub(xc_sic_write_info)
165
166 call messages_print_var_option('SICCorrection', sic%level, iunit=iunit, namespace=namespace)
167
168 pop_sub(xc_sic_write_info)
169 end subroutine xc_sic_write_info
170
171 ! ---------------------------------------------------------
188 subroutine xc_sic_calc_adsic(sic, namespace, space, gr, st, hm, xc, density, vxc, ex, ec)
189 type(xc_sic_t), intent(in) :: sic
190 type(namespace_t), intent(in) :: namespace
191 class(space_t), intent(in) :: space
192 type(grid_t), intent(in) :: gr
193 type(states_elec_t), intent(in) :: st
194 type(hamiltonian_elec_t), intent(in) :: hm
195 type(xc_t), intent(inout) :: xc
196 real(real64), contiguous, intent(in) :: density(:,:)
197 real(real64), contiguous, intent(inout) :: vxc(:,:)
198 real(real64), optional, intent(inout) :: ex, ec
199
200 integer :: ispin, ist, ik, ip
201 real(real64), allocatable :: vxc_sic(:,:), vh_sic(:), rho(:, :), qsp(:)
202 real(real64) :: ex_sic, ec_sic
203 real(real64) :: dtot, dpol, vpol
204 real(real64) :: nup
205
206 push_sub(xc_sic_calc_adsic)
207
208 assert(sic%level == sic_adsic)
209
210 if (st%d%ispin == spinors .and. .not. in_family(hm%xc%family, [xc_family_lda, xc_family_gga])) then
211 write(message(1),'(a)') 'ADSIC with non-collinear spin is currently only possible'
212 write(message(2),'(a)') 'with LDA and GGA functionals.'
213 call messages_fatal(2, namespace=namespace)
214 end if
215
216 if (xc_is_not_size_consistent(xc, namespace)) then
217 call messages_not_implemented('ADSIC with size inconsistent functionals', namespace=namespace)
218 end if
219
220 ! We compute here the number of electrons per spin channel
221 safe_allocate(qsp(1:2))
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 = dmf_integrate(gr, st%d%nspin, density)
237 end if
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(qsp)
349 safe_deallocate_a(vxc_sic)
350 safe_deallocate_a(vh_sic)
351 safe_deallocate_a(rho)
352
353 pop_sub(xc_sic_calc_adsic)
354 end subroutine xc_sic_calc_adsic
355
356end module xc_sic_oct_m
357
358!! Local Variables:
359!! mode: f90
360!! coding: utf-8
361!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:173
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:190
real(real64), parameter, public m_four
Definition: global.F90:194
real(real64), parameter, public m_half
Definition: global.F90:196
real(real64), parameter, public m_one
Definition: global.F90:191
real(real64), parameter, public m_min_occ
Minimal occupation that is considered to be non-zero.
Definition: global.F90:213
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:1097
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:416
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:697
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:799
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:763
pure logical function, public in_family(family, xc_families)
Definition: xc.F90:648
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:253
integer, parameter, public sic_adsic
Averaged density SIC.
Definition: xc_sic.F90:150
subroutine, public xc_sic_init(sic, namespace, gr, st, mc, space)
initialize the SIC object
Definition: xc_sic.F90:170
subroutine, public xc_sic_end(sic)
finalize the SIC and, if needed, the included OEP
Definition: xc_sic.F90:239
integer, parameter, public sic_pz_oep
Perdew-Zunger SIC (OEP way)
Definition: xc_sic.F90:150
integer, parameter, public sic_amaldi
Amaldi correction term.
Definition: xc_sic.F90:150
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:284
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:157