Octopus
v_ks.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module v_ks_oct_m
22 use accel_oct_m
23 use comm_oct_m
25 use debug_oct_m
28 use energy_oct_m
32 use global_oct_m
33 use grid_oct_m
37 use ions_oct_m
38 use, intrinsic :: iso_fortran_env
39 use isdf_oct_m, only: isdf_parallel_ace_compute_potentials => isdf_ace_compute_potentials
44 use lda_u_oct_m
48 use mesh_oct_m
51 use mpi_oct_m
54 use parser_oct_m
59 use pseudo_oct_m
62 use sort_oct_m
63 use space_oct_m
71 use xc_cam_oct_m
72 use xc_oct_m
73 use xc_f03_lib_m
74 use xc_fbe_oct_m
78 use xc_oep_oct_m
79 use xc_sic_oct_m
81 use xc_vdw_oct_m
83
84 ! from the dftd3 library
85 use dftd3_api
86
87 implicit none
88
89 private
90 public :: &
91 v_ks_t, &
92 v_ks_init, &
93 v_ks_end, &
96 v_ks_calc, &
103
104 type v_ks_calc_t
105 private
106 logical :: calculating
107 logical :: time_present
108 real(real64) :: time
109 real(real64), allocatable :: density(:, :)
110 logical :: total_density_alloc
111 real(real64), pointer, contiguous :: total_density(:)
112 type(energy_t), allocatable :: energy
113
114 type(states_elec_t), pointer :: hf_st
118
119 real(real64), allocatable :: vxc(:, :)
120 real(real64), allocatable :: vtau(:, :)
121 real(real64), allocatable :: axc(:, :, :)
122 real(real64), allocatable :: a_ind(:, :)
123 real(real64), allocatable :: b_ind(:, :)
124 logical :: calc_energy
125 end type v_ks_calc_t
126
127 type v_ks_t
128 private
129 integer, public :: theory_level = -1
130 logical, public :: frozen_hxc = .false.
131
132 integer, public :: xc_family = 0
133 integer, public :: xc_flags = 0
134 integer, public :: xc_photon = 0
135 type(xc_t), public :: xc
136 type(xc_photons_t), public :: xc_photons
137 type(xc_oep_t), public :: oep
138 type(xc_oep_photon_t), public :: oep_photon
139 type(xc_ks_inversion_t), public :: ks_inversion
140 type(xc_sic_t), public :: sic
141 type(xc_vdw_t), public :: vdw
142 type(grid_t), pointer, public :: gr
143 type(v_ks_calc_t) :: calc
144 logical :: calculate_current = .false.
145 type(current_t) :: current_calculator
146 logical :: include_td_field = .false.
147 logical, public :: has_photons = .false.
148 logical :: xc_photon_include_hartree = .true.
149
150 real(real64), public :: stress_xc_gga(3, 3)
151 type(photon_mode_t), pointer, public :: pt => null()
152 type(mf_t), public :: pt_mx
153 end type v_ks_t
154
155contains
156
157 ! ---------------------------------------------------------
158 subroutine v_ks_init(ks, namespace, gr, st, ions, mc, space, kpoints)
159 type(v_ks_t), intent(inout) :: ks
160 type(namespace_t), intent(in) :: namespace
161 type(grid_t), target, intent(inout) :: gr
162 type(states_elec_t), intent(in) :: st
163 type(ions_t), intent(inout) :: ions
164 type(multicomm_t), intent(in) :: mc
165 class(space_t), intent(in) :: space
166 type(kpoints_t), intent(in) :: kpoints
167
168 integer :: x_id, c_id, xk_id, ck_id, default, val
169 logical :: parsed_theory_level, using_hartree_fock
170 integer :: pseudo_x_functional, pseudo_c_functional
171 integer :: oep_type
172
173 push_sub(v_ks_init)
174
175 ! We need to parse TheoryLevel and XCFunctional, this is
176 ! complicated because they are interdependent.
177
178 !%Variable TheoryLevel
179 !%Type integer
180 !%Section Hamiltonian
181 !%Description
182 !% The calculations can be run with different "theory levels" that
183 !% control how electrons are simulated. The default is
184 !% <tt>dft</tt>. When hybrid functionals are requested, through
185 !% the <tt>XCFunctional</tt> variable, the default is
186 !% <tt>hartree_fock</tt>.
187 !%Option independent_particles 2
188 !% Particles will be considered as independent, <i>i.e.</i> as non-interacting.
189 !% This mode is mainly used for testing purposes, as the code is usually
190 !% much faster with <tt>independent_particles</tt>.
191 !%Option hartree 1
192 !% Calculation within the Hartree method (experimental). Note that, contrary to popular
193 !% belief, the Hartree potential is self-interaction-free. Therefore, this run
194 !% mode will not yield the same result as <tt>kohn-sham</tt> without exchange-correlation.
195 !%Option hartree_fock 3
196 !% This is the traditional Hartree-Fock scheme. Like the Hartree scheme, it is fully
197 !% self-interaction-free.
198 !%Option kohn_sham 4
199 !% This is the default density-functional theory scheme. Note that you can also use
200 !% hybrid functionals in this scheme, but they will be handled the "DFT" way, <i>i.e.</i>,
201 !% solving the OEP equation.
202 !%Option generalized_kohn_sham 5
203 !% This is similar to the <tt>kohn-sham</tt> scheme, except that this allows for nonlocal operators.
204 !% This is the default mode to run hybrid functionals, meta-GGA functionals, or DFT+U.
205 !% It can be more convenient to use <tt>kohn-sham</tt> DFT within the OEP scheme to get similar (but not the same) results.
206 !% Note that within this scheme you can use a correlation functional, or a hybrid
207 !% functional (see <tt>XCFunctional</tt>). In the latter case, you will be following the
208 !% quantum-chemistry recipe to use hybrids.
209 !%Option rdmft 7
210 !% (Experimental) Reduced Density Matrix functional theory.
211 !%End
212
213 ks%xc_family = xc_family_none
214 ks%sic%level = sic_none
215 ks%oep%level = oep_level_none
216 ks%oep_photon%level = oep_level_none
218 ks%theory_level = kohn_sham_dft
219 parsed_theory_level = .false.
220
221 ! the user knows what he wants, give her that
222 if (parse_is_defined(namespace, 'TheoryLevel')) then
223 call parse_variable(namespace, 'TheoryLevel', kohn_sham_dft, ks%theory_level)
224 if (.not. varinfo_valid_option('TheoryLevel', ks%theory_level)) call messages_input_error(namespace, 'TheoryLevel')
226 parsed_theory_level = .true.
227 end if
229 ! parse the XC functional
231 call get_functional_from_pseudos(pseudo_x_functional, pseudo_c_functional)
233 default = 0
234 if (ks%theory_level == kohn_sham_dft .or. ks%theory_level == generalized_kohn_sham_dft) then
235 default = xc_get_default_functional(space%dim, pseudo_x_functional, pseudo_c_functional)
236 end if
238 if (.not. parse_is_defined(namespace, 'XCFunctional') &
239 .and. (pseudo_x_functional /= pseudo_exchange_any .or. pseudo_c_functional /= pseudo_correlation_any)) then
240 call messages_write('Info: the XCFunctional has been selected to match the pseudopotentials', new_line = .true.)
241 call messages_write(' used in the calculation.')
242 call messages_info(namespace=namespace)
243 end if
244
245 ! The description of this variable can be found in file src/xc/functionals_list.F90
246 call parse_variable(namespace, 'XCFunctional', default, val)
248 ! the first 3 digits of the number indicate the X functional and
249 ! the next 3 the C functional.
250 c_id = val / libxc_c_index
251 x_id = val - c_id * libxc_c_index
252
253 if ((x_id /= pseudo_x_functional .and. pseudo_x_functional /= pseudo_exchange_any) .or. &
254 (c_id /= pseudo_c_functional .and. pseudo_c_functional /= pseudo_exchange_any)) then
255 call messages_write('The XCFunctional that you selected does not match the one used', new_line = .true.)
256 call messages_write('to generate the pseudopotentials.')
257 call messages_warning(namespace=namespace)
258 end if
259
260 ! FIXME: we rarely need this. We should only parse when necessary.
261
262 !%Variable XCKernel
263 !%Type integer
264 !%Default -1
265 !%Section Hamiltonian::XC
266 !%Description
267 !% Defines the exchange-correlation kernel. Only LDA kernels are available currently.
268 !% The options are the same as <tt>XCFunctional</tt>.
269 !% Note: the kernel is only needed for Casida, Sternheimer, or optimal-control calculations.
270 !%Option xc_functional -1
271 !% The same functional defined by <tt>XCFunctional</tt>. By default, this is the case.
272 !%End
273 call parse_variable(namespace, 'XCKernel', -1, val)
274 if (-1 == val) then
275 ck_id = c_id
276 xk_id = x_id
277 else
278 ck_id = val / libxc_c_index
279 xk_id = val - ck_id * libxc_c_index
280 end if
281
282 call messages_obsolete_variable(namespace, 'XFunctional', 'XCFunctional')
283 call messages_obsolete_variable(namespace, 'CFunctional', 'XCFunctional')
284
285 !%Variable XCPhotonFunctional
286 !%Type integer
287 !%Default 0
288 !%Section Hamiltonian::XC
289 !%Description
290 !% Defines the exchange and correlation functionals to be used for the QEDFT
291 !% description of the electron-photon system.
292 !%Option none 0
293 !% No functional is used
294 !%Option photon_x_lda 10
295 !% Exchange-only local density approcimation
296 !%Option photon_xc_lda 11
297 !% Exchange-correlation local density approcimation
298 !%Option photon_x_wfn 20
299 !% Exchange-only based on wave functions
300 !%Option photon_xc_wfn 21
301 !% Exchange-correlation based on wave functions
302 !%End
303
304 call parse_variable(namespace, 'XCPhotonFunctional', option__xcphotonfunctional__none, ks%xc_photon)
305
306 !%Variable XCPhotonIncludeHartree
307 !%Type logical
308 !%Default yes
309 !%Section Hamiltonian::XC
310 !%Description
311 !% Use the Hartree potential and energy in calculations
312 !%End
313
314 call parse_variable(namespace, 'XCPhotonIncludeHartree', .true., ks%xc_photon_include_hartree)
315
316 if (.not. ks%xc_photon_include_hartree) then
317 call messages_write('turn off hartree potential and energy')
318 call messages_warning(namespace=namespace)
319 end if
320
321 ! initialize XC modules
322
323 ! This is a bit ugly, theory_level might not be generalized KS or HF now
324 ! but it might become generalized KS or HF later. This is safe because it
325 ! becomes generalized KS in the cases where the functional is hybrid
326 ! and the ifs inside check for both conditions.
327 using_hartree_fock = (ks%theory_level == hartree_fock) &
328 .or. (ks%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(ks%xc))
329 call xc_init(ks%xc, namespace, space%dim, space%periodic_dim, st%qtot, &
330 x_id, c_id, xk_id, ck_id, hartree_fock = using_hartree_fock, ispin=st%d%ispin)
331
332 ks%xc_family = ks%xc%family
333 ks%xc_flags = ks%xc%flags
334
335 if (.not. parsed_theory_level) then
336 default = kohn_sham_dft
337
338 ! the functional is a hybrid, use Hartree-Fock as theory level by default
339 if (family_is_hybrid(ks%xc) .or. family_is_mgga_with_exc(ks%xc)) then
341 end if
342
343 ! In principle we do not need to parse. However we do it for consistency
344 call parse_variable(namespace, 'TheoryLevel', default, ks%theory_level)
345 if (.not. varinfo_valid_option('TheoryLevel', ks%theory_level)) call messages_input_error(namespace, 'TheoryLevel')
346
347 end if
348
349 ! In case we need OEP, we need to find which type of OEP it is
350 oep_type = -1
351 if (family_is_mgga_with_exc(ks%xc)) then
352 call messages_experimental('MGGA energy functionals')
353
354 if (ks%theory_level == kohn_sham_dft) then
355 call messages_experimental("MGGA within the Kohn-Sham scheme")
356 ks%xc_family = ior(ks%xc_family, xc_family_oep)
357 oep_type = oep_type_mgga
358 end if
359 end if
360
361 call messages_obsolete_variable(namespace, 'NonInteractingElectrons', 'TheoryLevel')
362 call messages_obsolete_variable(namespace, 'HartreeFock', 'TheoryLevel')
363
364 ! Due to how the code is made, we need to set this to have theory level other than DFT
365 ! correct...
366 ks%sic%amaldi_factor = m_one
367
368 select case (ks%theory_level)
370
371 case (hartree)
372 call messages_experimental("Hartree theory level")
373 if (space%periodic_dim == space%dim) then
374 call messages_experimental("Hartree in fully periodic system")
375 end if
376 if (kpoints%full%npoints > 1) then
377 call messages_not_implemented("Hartree with k-points", namespace=namespace)
378 end if
379
380 case (hartree_fock)
381 if (kpoints%full%npoints > 1) then
382 call messages_experimental("Hartree-Fock with k-points")
383 end if
384
386 if (kpoints%full%npoints > 1 .and. family_is_hybrid(ks%xc)) then
387 call messages_experimental("Hybrid functionals with k-points")
388 end if
389
390 case (rdmft)
391 call messages_experimental('RDMFT theory level')
392
393 case (kohn_sham_dft)
394
395 ! check for SIC
396 if (bitand(ks%xc_family, xc_family_lda + xc_family_gga) /= 0) then
397 call xc_sic_init(ks%sic, namespace, gr, st, mc, space)
398 end if
399
400 if (bitand(ks%xc_family, xc_family_oep) /= 0) then
401 select case (ks%xc%functional(func_x,1)%id)
402 case (xc_oep_x_slater)
403 if (kpoints%reduced%npoints > 1) then
404 call messages_not_implemented("Slater with k-points", namespace=namespace)
405 end if
406 ks%oep%level = oep_level_none
407 case (xc_oep_x_fbe)
408 if (kpoints%reduced%npoints > 1) then
409 call messages_not_implemented("FBE functional with k-points", namespace=namespace)
410 end if
411 ks%oep%level = oep_level_none
412 case default
413 if((.not. ks%has_photons) .or. (ks%xc_photon /= 0)) then
414 if(oep_type == -1) then ! Else we have a MGGA
415 oep_type = oep_type_exx
416 end if
417 call xc_oep_init(ks%oep, namespace, gr, st, mc, space, oep_type)
418 end if
419 end select
420 else
421 ks%oep%level = oep_level_none
422 end if
423
424 if (bitand(ks%xc_family, xc_family_ks_inversion) /= 0) then
425 call xc_ks_inversion_init(ks%ks_inversion, namespace, gr, ions, st, ks%xc, mc, space, kpoints)
426 end if
427
428 end select
429
430 if (ks%theory_level /= kohn_sham_dft .and. parse_is_defined(namespace, "SICCorrection")) then
431 message(1) = "SICCorrection can only be used with Kohn-Sham DFT"
432 call messages_fatal(1, namespace=namespace)
433 end if
434
435 if (st%d%ispin == spinors) then
436 if (bitand(ks%xc_family, xc_family_mgga + xc_family_hyb_mgga) /= 0) then
437 call messages_not_implemented("MGGA with spinors", namespace=namespace)
438 end if
439 end if
440
441 ks%frozen_hxc = .false.
442
443 call v_ks_write_info(ks, namespace=namespace)
444
445 ks%gr => gr
446 ks%calc%calculating = .false.
447
448 !The value of ks%calculate_current is set to false or true by Output
449 call current_init(ks%current_calculator, namespace)
450
451 call ks%vdw%init(namespace, space, gr, ks%xc, ions, x_id, c_id)
452 if (ks%vdw%vdw_correction /= option__vdwcorrection__none .and. ks%theory_level == rdmft) then
453 message(1) = "VDWCorrection and RDMFT are not compatible"
454 call messages_fatal(1, namespace=namespace)
455 end if
456 if (ks%vdw%vdw_correction /= option__vdwcorrection__none .and. ks%theory_level == independent_particles) then
457 message(1) = "VDWCorrection and independent particles are not compatible"
458 call messages_fatal(1, namespace=namespace)
459 end if
460
461 if (ks%xc_photon /= 0) then
462 ! initilize the photon free variables
463 call ks%xc_photons%init(namespace, ks%xc_photon , space, gr, st)
464 ! remornalize the electron mass due to light-matter interaction; here we only deal with it in free space
465 ks%oep_photon%level = oep_level_none
466 end if
467
468
469 pop_sub(v_ks_init)
470
471 contains
472
474 subroutine get_functional_from_pseudos(x_functional, c_functional)
475 integer, intent(out) :: x_functional
476 integer, intent(out) :: c_functional
477
478 integer :: xf, cf, ispecies
479 logical :: warned_inconsistent
480
481 x_functional = pseudo_exchange_any
482 c_functional = pseudo_correlation_any
483
484 warned_inconsistent = .false.
485 do ispecies = 1, ions%nspecies
486 select type(spec=>ions%species(ispecies)%s)
487 class is(pseudopotential_t)
488 xf = spec%x_functional()
489 cf = spec%c_functional()
490
491 if (xf == pseudo_exchange_unknown .or. cf == pseudo_correlation_unknown) then
492 call messages_write("Unknown XC functional for species '"//trim(ions%species(ispecies)%s%get_label())//"'")
493 call messages_warning(namespace=namespace)
494 cycle
495 end if
496
497 if (x_functional == pseudo_exchange_any) then
498 x_functional = xf
499 else
500 if (xf /= x_functional .and. .not. warned_inconsistent) then
501 call messages_write('Inconsistent XC functional detected between species')
502 call messages_warning(namespace=namespace)
503 warned_inconsistent = .true.
504 end if
505 end if
506
507 if (c_functional == pseudo_correlation_any) then
508 c_functional = cf
509 else
510 if (cf /= c_functional .and. .not. warned_inconsistent) then
511 call messages_write('Inconsistent XC functional detected between species')
512 call messages_warning(namespace=namespace)
513 warned_inconsistent = .true.
514 end if
515 end if
516
517 class default
520 end select
521
522 end do
523
524 assert(x_functional /= pseudo_exchange_unknown)
525 assert(c_functional /= pseudo_correlation_unknown)
526
527 end subroutine get_functional_from_pseudos
528 end subroutine v_ks_init
529 ! ---------------------------------------------------------
530
531 ! ---------------------------------------------------------
532 subroutine v_ks_end(ks)
533 type(v_ks_t), intent(inout) :: ks
534
535 push_sub(v_ks_end)
536
537 call ks%vdw%end()
538
539 select case (ks%theory_level)
540 case (kohn_sham_dft)
541 if (bitand(ks%xc_family, xc_family_ks_inversion) /= 0) then
542 call xc_ks_inversion_end(ks%ks_inversion)
543 end if
544 if (bitand(ks%xc_family, xc_family_oep) /= 0) then
545 call xc_oep_end(ks%oep)
546 end if
547 call xc_end(ks%xc)
549 call xc_end(ks%xc)
550 end select
551
552 call xc_sic_end(ks%sic)
553
554 if (ks%xc_photon /= 0) then
555 call ks%xc_photons%end()
556 end if
557
558 pop_sub(v_ks_end)
559 end subroutine v_ks_end
560 ! ---------------------------------------------------------
561
562
563 ! ---------------------------------------------------------
564 subroutine v_ks_write_info(ks, iunit, namespace)
565 type(v_ks_t), intent(in) :: ks
566 integer, optional, intent(in) :: iunit
567 type(namespace_t), optional, intent(in) :: namespace
568
570
571 call messages_print_with_emphasis(msg="Theory Level", iunit=iunit, namespace=namespace)
572 call messages_print_var_option("TheoryLevel", ks%theory_level, iunit=iunit, namespace=namespace)
573
574 select case (ks%theory_level)
576 call messages_info(iunit=iunit, namespace=namespace)
577 call xc_write_info(ks%xc, iunit, namespace)
578
579 case (kohn_sham_dft)
580 call messages_info(iunit=iunit, namespace=namespace)
581 call xc_write_info(ks%xc, iunit, namespace)
582
583 call messages_info(iunit=iunit, namespace=namespace)
584
585 call xc_sic_write_info(ks%sic, iunit, namespace)
586 call xc_oep_write_info(ks%oep, iunit, namespace)
587 call xc_ks_inversion_write_info(ks%ks_inversion, iunit, namespace)
588
589 end select
590
591 call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
592
593 pop_sub(v_ks_write_info)
594 end subroutine v_ks_write_info
595 ! ---------------------------------------------------------
596
597
598 !----------------------------------------------------------
599 subroutine v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, calc_eigenval, calc_current)
600 type(namespace_t), intent(in) :: namespace
601 type(electron_space_t), intent(in) :: space
602 type(grid_t), intent(in) :: gr
603 type(ions_t), intent(in) :: ions
604 type(partner_list_t), intent(in) :: ext_partners
605 type(states_elec_t), intent(inout) :: st
606 type(v_ks_t), intent(inout) :: ks
607 type(hamiltonian_elec_t), intent(inout) :: hm
608 logical, optional, intent(in) :: calc_eigenval
609 logical, optional, intent(in) :: calc_current
610
611 integer, allocatable :: ind(:)
612 integer :: ist, ik
613 real(real64), allocatable :: copy_occ(:)
614 logical :: calc_eigenval_
615 logical :: calc_current_
616
617 push_sub(v_ks_h_setup)
618
619 calc_eigenval_ = optional_default(calc_eigenval, .true.)
620 calc_current_ = optional_default(calc_current, .true.)
621 call states_elec_fermi(st, namespace, gr)
622 call density_calc(st, gr, st%rho)
623 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
624 calc_eigenval = calc_eigenval_, calc_current = calc_current_) ! get potentials
625
626 if (st%restart_reorder_occs .and. .not. st%fromScratch) then
627 message(1) = "Reordering occupations for restart."
628 call messages_info(1, namespace=namespace)
629
630 safe_allocate(ind(1:st%nst))
631 safe_allocate(copy_occ(1:st%nst))
632
633 do ik = 1, st%nik
634 call sort(st%eigenval(:, ik), ind)
635 copy_occ(1:st%nst) = st%occ(1:st%nst, ik)
636 do ist = 1, st%nst
637 st%occ(ist, ik) = copy_occ(ind(ist))
638 end do
639 end do
640
641 safe_deallocate_a(ind)
642 safe_deallocate_a(copy_occ)
643 end if
644
645 if (calc_eigenval_) call states_elec_fermi(st, namespace, gr) ! occupations
646 call energy_calc_total(namespace, space, hm, gr, st, ext_partners)
647
648 pop_sub(v_ks_h_setup)
649 end subroutine v_ks_h_setup
650
651 ! ---------------------------------------------------------
652 subroutine v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
653 calc_eigenval, time, calc_energy, calc_current, force_semilocal)
654 type(v_ks_t), intent(inout) :: ks
655 type(namespace_t), intent(in) :: namespace
656 type(electron_space_t), intent(in) :: space
657 type(hamiltonian_elec_t), intent(inout) :: hm
658 type(states_elec_t), intent(inout) :: st
659 type(ions_t), intent(in) :: ions
660 type(partner_list_t), intent(in) :: ext_partners
661 logical, optional, intent(in) :: calc_eigenval
662 real(real64), optional, intent(in) :: time
663 logical, optional, intent(in) :: calc_energy
664 logical, optional, intent(in) :: calc_current
665 logical, optional, intent(in) :: force_semilocal
666
667 logical :: calc_current_
668
669 push_sub(v_ks_calc)
670
671 calc_current_ = optional_default(calc_current, .true.) &
672 .and. ks%calculate_current &
673 .and. states_are_complex(st) &
675
676 if (calc_current_) then
677 call states_elec_allocate_current(st, space, ks%gr)
678 call current_calculate(ks%current_calculator, namespace, ks%gr, hm, space, st)
679 end if
680
681 call v_ks_calc_start(ks, namespace, space, hm, st, ions, hm%kpoints%latt, ext_partners, time, &
682 calc_energy, force_semilocal=force_semilocal)
683 call v_ks_calc_finish(ks, hm, namespace, space, hm%kpoints%latt, st, &
684 ext_partners, force_semilocal=force_semilocal)
685
686 if (optional_default(calc_eigenval, .false.)) then
687 call energy_calc_eigenvalues(namespace, hm, ks%gr%der, st)
688 end if
689
690 ! Update the magnetic constrain
691 call magnetic_constrain_update(hm%magnetic_constrain, ks%gr, st%d, space, hm%kpoints%latt, ions%pos, st%rho)
692 ! We add the potential to vxc, as this way the potential gets mixed together with vxc
693 ! While this is not ideal, this is a simple practical solution
694 if (hm%magnetic_constrain%level /= constrain_none) then
695 call lalg_axpy(ks%gr%np, st%d%nspin, m_one, hm%magnetic_constrain%pot, hm%ks_pot%vhxc)
696 end if
697
698 pop_sub(v_ks_calc)
699 end subroutine v_ks_calc
700
701 ! ---------------------------------------------------------
702
707 subroutine v_ks_calc_start(ks, namespace, space, hm, st, ions, latt, ext_partners, time, &
708 calc_energy, force_semilocal)
709 type(v_ks_t), target, intent(inout) :: ks
710 type(namespace_t), intent(in) :: namespace
711 class(space_t), intent(in) :: space
712 type(hamiltonian_elec_t), target, intent(in) :: hm
713 type(states_elec_t), target, intent(inout) :: st
714 type(ions_t), intent(in) :: ions
715 type(lattice_vectors_t), intent(in) :: latt
716 type(partner_list_t), intent(in) :: ext_partners
717 real(real64), optional, intent(in) :: time
718 logical, optional, intent(in) :: calc_energy
719 logical, optional, intent(in) :: force_semilocal
720
721 push_sub(v_ks_calc_start)
722
723 call profiling_in("KOHN_SHAM_CALC")
724
725 assert(.not. ks%calc%calculating)
726 ks%calc%calculating = .true.
727
728 write(message(1), '(a)') 'Debug: Calculating Kohn-Sham potential.'
729 call messages_info(1, namespace=namespace, debug_only=.true.)
730
731 ks%calc%time_present = present(time)
732 ks%calc%time = optional_default(time, m_zero)
733
734 ks%calc%calc_energy = optional_default(calc_energy, .true.)
735
736 ! If the Hxc term is frozen, there is nothing more to do (WARNING: MISSING ks%calc%energy%intnvxc)
737 if (ks%frozen_hxc) then
738 call profiling_out("KOHN_SHAM_CALC")
739 pop_sub(v_ks_calc_start)
740 return
741 end if
742
743 allocate(ks%calc%energy)
744
745 call energy_copy(hm%energy, ks%calc%energy)
746
747 ks%calc%energy%intnvxc = m_zero
748
749 nullify(ks%calc%total_density)
750
751 if (ks%theory_level /= independent_particles .and. abs(ks%sic%amaldi_factor) > m_epsilon) then
752
753 call calculate_density()
754
755 if (poisson_is_async(hm%psolver)) then
756 call dpoisson_solve_start(hm%psolver, ks%calc%total_density)
757 end if
758
759 if (ks%theory_level /= hartree .and. ks%theory_level /= rdmft) call v_a_xc(hm, force_semilocal)
760 else
761 ks%calc%total_density_alloc = .false.
762 end if
763
764 ! The exchange operator is computed from the states of the previous iteration
765 ! This is done by copying the state object to ks%calc%hf_st
766 ! For ACE, the states are the same in ks%calc%hf_st and st, as we compute the
767 ! ACE potential in v_ks_finish, so the copy is not needed
768 nullify(ks%calc%hf_st)
769 if (ks%theory_level == hartree .or. ks%theory_level == hartree_fock &
770 .or. ks%theory_level == rdmft .or. (ks%theory_level == generalized_kohn_sham_dft &
771 .and. family_is_hybrid(ks%xc))) then
772
773 if (st%parallel_in_states) then
774 if (accel_is_enabled()) then
775 call messages_write('State parallelization of Hartree-Fock exchange is not supported')
776 call messages_new_line()
777 call messages_write('when running with OpenCL/CUDA. Please use domain parallelization')
778 call messages_new_line()
779 call messages_write("or disable acceleration using 'DisableAccel = yes'.")
780 call messages_fatal(namespace=namespace)
781 end if
782 end if
783
784 if (hm%exxop%useACE) then
785 ks%calc%hf_st => st
786 else
787 safe_allocate(ks%calc%hf_st)
788 call states_elec_copy(ks%calc%hf_st, st)
789 end if
790 end if
791
792 ! Calculate the vector potential induced by the electronic current.
793 ! WARNING: calculating the self-induced magnetic field here only makes
794 ! sense if it is going to be used in the Hamiltonian, which does not happen
795 ! now. Otherwise one could just calculate it at the end of the calculation.
796 if (hm%self_induced_magnetic) then
797 safe_allocate(ks%calc%a_ind(1:ks%gr%np_part, 1:space%dim))
798 safe_allocate(ks%calc%b_ind(1:ks%gr%np_part, 1:space%dim))
799 call magnetic_induced(namespace, ks%gr, st, hm%psolver, hm%kpoints, ks%calc%a_ind, ks%calc%b_ind)
800 end if
801
802 if ((ks%has_photons) .and. (ks%calc%time_present) .and. (ks%xc_photon == 0) ) then
803 call mf_calc(ks%pt_mx, ks%gr, st, ions, ks%pt, time)
804 end if
805
806 ! if (ks%has_vibrations) then
807 ! call vibrations_eph_coup(ks%vib, ks%gr, hm, ions, st)
808 ! end if
809
810 call profiling_out("KOHN_SHAM_CALC")
811 pop_sub(v_ks_calc_start)
812
813 contains
814
815 subroutine calculate_density()
816 integer :: ip
817
819
820 ! get density taking into account non-linear core corrections
821 safe_allocate(ks%calc%density(1:ks%gr%np, 1:st%d%nspin))
822 call states_elec_total_density(st, ks%gr, ks%calc%density)
823
824 ! Amaldi correction
825 if (ks%sic%level == sic_amaldi) then
826 call lalg_scal(ks%gr%np, st%d%nspin, ks%sic%amaldi_factor, ks%calc%density)
827 end if
828
829 nullify(ks%calc%total_density)
830 if (allocated(st%rho_core) .or. hm%d%spin_channels > 1) then
831 ks%calc%total_density_alloc = .true.
832
833 safe_allocate(ks%calc%total_density(1:ks%gr%np))
834
835 do ip = 1, ks%gr%np
836 ks%calc%total_density(ip) = sum(ks%calc%density(ip, 1:hm%d%spin_channels))
837 end do
838
839 ! remove non-local core corrections
840 if (allocated(st%rho_core)) then
841 call lalg_axpy(ks%gr%np, -ks%sic%amaldi_factor, st%rho_core, ks%calc%total_density)
842 end if
843 else
844 ks%calc%total_density_alloc = .false.
845 ks%calc%total_density => ks%calc%density(:, 1)
846 end if
847
849 end subroutine calculate_density
850
851 ! ---------------------------------------------------------
852 subroutine v_a_xc(hm, force_semilocal)
853 type(hamiltonian_elec_t), intent(in) :: hm
854 logical, optional, intent(in) :: force_semilocal
855
856 integer :: ispin
857
858 push_sub(v_ks_calc_start.v_a_xc)
859 call profiling_in("XC")
860
861 ks%calc%energy%exchange = m_zero
862 ks%calc%energy%correlation = m_zero
863 ks%calc%energy%xc_j = m_zero
864 ks%calc%energy%vdw = m_zero
865
866 allocate(ks%calc%vxc(1:ks%gr%np, 1:st%d%nspin))
867 ks%calc%vxc = m_zero
868
869 if (family_is_mgga_with_exc(hm%xc)) then
870 safe_allocate(ks%calc%vtau(1:ks%gr%np, 1:st%d%nspin))
871 ks%calc%vtau = m_zero
872 end if
873
874 ! Get the *local* XC term
875 if (ks%calc%calc_energy) then
876 if (family_is_mgga_with_exc(hm%xc)) then
877 call xc_get_vxc(ks%gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, ks%calc%density, st%d%ispin, &
878 latt%rcell_volume, ks%calc%vxc, ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation, &
879 deltaxc = ks%calc%energy%delta_xc, vtau = ks%calc%vtau, force_orbitalfree=force_semilocal)
880 else
881 call xc_get_vxc(ks%gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, ks%calc%density, st%d%ispin, &
882 latt%rcell_volume, ks%calc%vxc, ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation, &
883 deltaxc = ks%calc%energy%delta_xc, stress_xc=ks%stress_xc_gga, force_orbitalfree=force_semilocal)
884 end if
885 else
886 if (family_is_mgga_with_exc(hm%xc)) then
887 call xc_get_vxc(ks%gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, ks%calc%density, &
888 st%d%ispin, latt%rcell_volume, ks%calc%vxc, vtau = ks%calc%vtau, force_orbitalfree=force_semilocal)
889 else
890 call xc_get_vxc(ks%gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, ks%calc%density, &
891 st%d%ispin, latt%rcell_volume, ks%calc%vxc, stress_xc=ks%stress_xc_gga, force_orbitalfree=force_semilocal)
892 end if
893 end if
894
895 !Noncollinear functionals
896 if (bitand(hm%xc%family, xc_family_nc_lda + xc_family_nc_mgga) /= 0) then
897 if (st%d%ispin /= spinors) then
898 message(1) = "Noncollinear functionals can only be used with spinor wavefunctions."
899 call messages_fatal(1)
900 end if
901
902 if (optional_default(force_semilocal, .false.)) then
903 message(1) = "Cannot perform LCAO for noncollinear MGGAs."
904 message(2) = "Please perform a LDA calculation first."
905 call messages_fatal(2)
906 end if
907
908 if (ks%calc%calc_energy) then
909 if (family_is_mgga_with_exc(hm%xc)) then
910 call xc_get_nc_vxc(ks%gr, ks%xc, st, hm%kpoints, space, namespace, ks%calc%density, ks%calc%vxc, &
911 vtau = ks%calc%vtau, ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation)
912 else
913 call xc_get_nc_vxc(ks%gr, ks%xc, st, hm%kpoints, space, namespace, ks%calc%density, ks%calc%vxc, &
914 ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation)
915 end if
916 else
917 if (family_is_mgga_with_exc(hm%xc)) then
918 call xc_get_nc_vxc(ks%gr, ks%xc, st, hm%kpoints, space, namespace, ks%calc%density, &
919 ks%calc%vxc, vtau = ks%calc%vtau)
920 else
921 call xc_get_nc_vxc(ks%gr, ks%xc, st, hm%kpoints, space, namespace, ks%calc%density, ks%calc%vxc)
922 end if
923 end if
924 end if
925
926 call ks%vdw%calc(namespace, space, latt, ions%atom, ions%natoms, ions%pos, &
927 ks%gr, st, ks%calc%energy%vdw, ks%calc%vxc)
928
929 if (optional_default(force_semilocal, .false.)) then
930 call profiling_out("XC")
931 pop_sub(v_ks_calc_start.v_a_xc)
932 return
933 end if
934
935 ! ADSIC correction
936 if (ks%sic%level == sic_adsic) then
937 if (family_is_mgga(hm%xc%family)) then
938 call messages_not_implemented('ADSIC with MGGAs', namespace=namespace)
939 end if
940 if (ks%calc%calc_energy) then
941 call xc_sic_calc_adsic(ks%sic, namespace, space, ks%gr, st, hm, ks%xc, ks%calc%density, &
942 ks%calc%vxc, ex = ks%calc%energy%exchange, ec = ks%calc%energy%correlation)
943 else
944 call xc_sic_calc_adsic(ks%sic, namespace, space, ks%gr, st, hm, ks%xc, ks%calc%density, &
945 ks%calc%vxc)
946 end if
947 end if
948 !PZ SIC is done in the finish routine as OEP full needs to update the Hamiltonian
949
950 if (ks%theory_level == kohn_sham_dft) then
951 ! The OEP family has to be handled specially
952 ! Note that OEP is done in the finish state, as it requires updating the Hamiltonian and needs the new Hartre and vxc term
953 if (bitand(ks%xc_family, xc_family_oep) /= 0 .or. family_is_mgga_with_exc(ks%xc)) then
954
955 if (ks%xc%functional(func_x,1)%id == xc_oep_x_slater) then
956 call x_slater_calc(namespace, ks%gr, space, hm%exxop, st, hm%kpoints, ks%calc%energy%exchange, &
957 vxc = ks%calc%vxc)
958 else if (ks%xc%functional(func_x,1)%id == xc_oep_x_fbe .or. ks%xc%functional(func_x,1)%id == xc_oep_x_fbe_sl) then
959 call x_fbe_calc(ks%xc%functional(func_x,1)%id, namespace, hm%psolver, ks%gr, st, space, &
960 ks%calc%energy%exchange, vxc = ks%calc%vxc)
961
962 else if (ks%xc%functional(func_c,1)%id == xc_lda_c_fbe_sl) then
963
964 call fbe_c_lda_sl(namespace, hm%psolver, ks%gr, st, space, &
965 ks%calc%energy%correlation, vxc = ks%calc%vxc)
966
967 end if
968
969 end if
970
971 if (bitand(ks%xc_family, xc_family_ks_inversion) /= 0) then
972 ! Also treat KS inversion separately (not part of libxc)
973 call xc_ks_inversion_calc(ks%ks_inversion, namespace, space, ks%gr, hm, ext_partners, st, vxc = ks%calc%vxc, &
974 time = ks%calc%time)
975 end if
976
977 ! compute the photon-free photon exchange potential and energy
978 if (ks%xc_photon /= 0) then
979
980 call ks%xc_photons%v_ks(namespace, ks%calc%total_density, ks%gr, space, hm%psolver, hm%ep, st)
981
982 ! add the photon-free px potential into the xc potential
983 do ispin = 1, hm%d%spin_channels
984 call lalg_axpy(ks%gr%np, m_one, ks%xc_photons%vpx(1:ks%gr%np), ks%calc%vxc(1:ks%gr%np, ispin) )
985 end do
986
987 ! photon-exchange energy
988 ks%calc%energy%photon_exchange = ks%xc_photons%ex
989 end if
990
991 end if
992
993 if (ks%calc%calc_energy) then
994 ! MGGA vtau contribution is done after copying vtau to hm%vtau
995
996 call v_ks_update_dftu_energy(ks, namespace, hm, st, ks%calc%energy%int_dft_u)
997 end if
998
999 call profiling_out("XC")
1000 pop_sub(v_ks_calc_start.v_a_xc)
1001 end subroutine v_a_xc
1002
1003 end subroutine v_ks_calc_start
1004 ! ---------------------------------------------------------
1005
1006 subroutine v_ks_calc_finish(ks, hm, namespace, space, latt, st, ext_partners, force_semilocal)
1007 type(v_ks_t), target, intent(inout) :: ks
1008 type(hamiltonian_elec_t), intent(inout) :: hm
1009 type(namespace_t), intent(in) :: namespace
1010 class(space_t), intent(in) :: space
1011 type(lattice_vectors_t), intent(in) :: latt
1012 type(states_elec_t), intent(inout) :: st
1013 type(partner_list_t), intent(in) :: ext_partners
1014 logical, optional, intent(in) :: force_semilocal
1015
1016 integer :: ip, ispin
1017 type(states_elec_t) :: xst
1019 real(real64) :: exx_energy
1020 real(real64) :: factor
1021
1022 push_sub(v_ks_calc_finish)
1023
1024 assert(ks%calc%calculating)
1025 ks%calc%calculating = .false.
1026
1027 if (ks%frozen_hxc) then
1028 pop_sub(v_ks_calc_finish)
1029 return
1030 end if
1031
1032 !change the pointer to the energy object
1033 safe_deallocate_a(hm%energy)
1034 call move_alloc(ks%calc%energy, hm%energy)
1035
1036 if (hm%self_induced_magnetic) then
1037 hm%a_ind(1:ks%gr%np, 1:space%dim) = ks%calc%a_ind(1:ks%gr%np, 1:space%dim)
1038 hm%b_ind(1:ks%gr%np, 1:space%dim) = ks%calc%b_ind(1:ks%gr%np, 1:space%dim)
1039
1040 safe_deallocate_a(ks%calc%a_ind)
1041 safe_deallocate_a(ks%calc%b_ind)
1042 end if
1043
1044 if (allocated(hm%v_static)) then
1045 hm%energy%intnvstatic = dmf_dotp(ks%gr, ks%calc%total_density, hm%v_static)
1046 else
1047 hm%energy%intnvstatic = m_zero
1048 end if
1049
1050 if (ks%theory_level == independent_particles .or. abs(ks%sic%amaldi_factor) <= m_epsilon) then
1051
1052 hm%ks_pot%vhxc = m_zero
1053 hm%energy%intnvxc = m_zero
1054 hm%energy%hartree = m_zero
1055 hm%energy%exchange = m_zero
1056 hm%energy%exchange_hf = m_zero
1057 hm%energy%correlation = m_zero
1058 else
1059
1060 hm%energy%hartree = m_zero
1061 call v_ks_hartree(namespace, ks, space, hm, ext_partners)
1062
1063 if (.not. optional_default(force_semilocal, .false.)) then
1064 !PZ-SIC
1065 if(ks%sic%level == sic_pz_oep) then
1066 if (states_are_real(st)) then
1067 call dxc_oep_calc(ks%sic%oep, namespace, ks%xc, ks%gr, hm, st, space, &
1068 latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1069 else
1070 call zxc_oep_calc(ks%sic%oep, namespace, ks%xc, ks%gr, hm, st, space, &
1071 latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1072 end if
1073 end if
1074
1075 ! OEP for exchange ad MGGAs (within Kohn-Sham DFT)
1076 if (ks%theory_level == kohn_sham_dft .and. ks%oep%level /= oep_level_none) then
1077 ! The OEP family has to be handled specially
1078 if (ks%xc%functional(func_x,1)%id == xc_oep_x .or. family_is_mgga_with_exc(ks%xc)) then
1079 if (states_are_real(st)) then
1080 call dxc_oep_calc(ks%oep, namespace, ks%xc, ks%gr, hm, st, space, &
1081 latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1082 else
1083 call zxc_oep_calc(ks%oep, namespace, ks%xc, ks%gr, hm, st, space, &
1084 latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1085 end if
1086 end if
1087 end if
1088 end if
1089
1090 if (ks%theory_level == kohn_sham_dft .and. ks%oep_photon%level /= oep_level_none) then
1091 if (states_are_real(st)) then
1092 call dxc_oep_photon_calc(ks%oep_photon, namespace, ks%xc, ks%gr, &
1093 hm, st, space, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1094 else
1095 call zxc_oep_photon_calc(ks%oep_photon, namespace, ks%xc, ks%gr, &
1096 hm, st, space, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1097 end if
1098 hm%energy%photon_exchange = ks%oep_photon%pt%ex
1099 end if
1100
1102 if (ks%calc%calc_energy) then
1103 ! Now we calculate Int[n vxc] = energy%intnvxc
1104 hm%energy%intnvxc = m_zero
1105
1106 if (ks%theory_level /= independent_particles .and. ks%theory_level /= hartree .and. ks%theory_level /= rdmft) then
1107 do ispin = 1, hm%d%nspin
1108 if (ispin <= 2) then
1109 factor = m_one
1110 else
1111 factor = m_two
1112 end if
1113 hm%energy%intnvxc = hm%energy%intnvxc + &
1114 factor*dmf_dotp(ks%gr, st%rho(:, ispin), ks%calc%vxc(:, ispin), reduce = .false.)
1115 end do
1116 call ks%gr%allreduce(hm%energy%intnvxc)
1117 end if
1118 end if
1119
1120
1121 if (ks%theory_level /= hartree .and. ks%theory_level /= rdmft) then
1122 ! move allocation of vxc from ks%calc to hm
1123 safe_deallocate_a(hm%ks_pot%vxc)
1124 call move_alloc(ks%calc%vxc, hm%ks_pot%vxc)
1125
1126 if (family_is_mgga_with_exc(hm%xc)) then
1127 call hm%ks_pot%set_vtau(ks%calc%vtau)
1128 safe_deallocate_a(ks%calc%vtau)
1129
1130 ! We need to evaluate the energy after copying vtau to hm%vtau
1131 if (ks%theory_level == generalized_kohn_sham_dft .and. ks%calc%calc_energy) then
1132 ! MGGA vtau contribution
1133 if (states_are_real(st)) then
1134 hm%energy%intnvxc = hm%energy%intnvxc &
1135 + denergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_mgga)
1136 else
1137 hm%energy%intnvxc = hm%energy%intnvxc &
1138 + zenergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_mgga)
1139 end if
1140 end if
1141 end if
1142
1143 else
1144 hm%ks_pot%vxc = m_zero
1145 end if
1146
1147 if (.not. ks%xc_photon_include_hartree) then
1148 hm%energy%hartree = m_zero
1149 hm%ks_pot%vhartree = m_zero
1150 end if
1151
1152 ! Build Hartree + XC potential
1153
1154 do ip = 1, ks%gr%np
1155 hm%ks_pot%vhxc(ip, 1) = hm%ks_pot%vxc(ip, 1) + hm%ks_pot%vhartree(ip)
1156 end do
1157 if (allocated(hm%vberry)) then
1158 do ip = 1, ks%gr%np
1159 hm%ks_pot%vhxc(ip, 1) = hm%ks_pot%vhxc(ip, 1) + hm%vberry(ip, 1)
1160 end do
1161 end if
1162
1163 if (hm%d%ispin > unpolarized) then
1164 do ip = 1, ks%gr%np
1165 hm%ks_pot%vhxc(ip, 2) = hm%ks_pot%vxc(ip, 2) + hm%ks_pot%vhartree(ip)
1166 end do
1167 if (allocated(hm%vberry)) then
1168 do ip = 1, ks%gr%np
1169 hm%ks_pot%vhxc(ip, 2) = hm%ks_pot%vhxc(ip, 2) + hm%vberry(ip, 2)
1170 end do
1171 end if
1172 end if
1173
1174 if (hm%d%ispin == spinors) then
1175 do ispin=3, 4
1176 do ip = 1, ks%gr%np
1177 hm%ks_pot%vhxc(ip, ispin) = hm%ks_pot%vxc(ip, ispin)
1178 end do
1179 end do
1180 end if
1181
1182 ! Note: this includes hybrids calculated with the Fock operator instead of OEP
1183 hm%energy%exchange_hf = m_zero
1184 if (ks%theory_level == hartree .or. ks%theory_level == hartree_fock &
1185 .or. ks%theory_level == rdmft &
1186 .or. (ks%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(ks%xc))) then
1187
1188 ! swap the states object
1189 if (.not. hm%exxop%useACE) then
1190 ! We also close the MPI remote memory access to the old object
1191 if (associated(hm%exxop%st)) then
1193 call states_elec_end(hm%exxop%st)
1194 safe_deallocate_p(hm%exxop%st)
1195 end if
1196 ! We activate the MPI remote memory access for ks%calc%hf_st
1197 ! This allows to have all calls to exchange_operator_apply_standard to access
1198 ! the states over MPI
1200 end if
1201
1202 ! The exchange operator will use ks%calc%hf_st
1203 ! For the ACE case, this is the same as st
1204 if (.not. optional_default(force_semilocal, .false.)) then
1205 select case (ks%theory_level)
1207 if (family_is_hybrid(ks%xc)) then
1208 call exchange_operator_reinit(hm%exxop, ks%xc%cam, ks%calc%hf_st)
1209 end if
1210 case (hartree_fock)
1211 call exchange_operator_reinit(hm%exxop, ks%xc%cam, ks%calc%hf_st)
1212 case (hartree, rdmft)
1213 call exchange_operator_reinit(hm%exxop, cam_exact_exchange, ks%calc%hf_st)
1214 end select
1215
1216 !This should be changed and the CAM parameters should also be obtained from the restart information
1217 !Maybe the parameters should be mixed too.
1218 exx_energy = m_zero
1219 if (hm%exxop%useACE) then
1220 call xst%nullify()
1221 if (states_are_real(ks%calc%hf_st)) then
1222 ! TODO(Alex) Clean up nested if statements
1223 if (hm%exxop%with_isdf) then
1224 ! Find interpolation points from density, or read from file
1225 ! TODO(Alex) Issue #1195 Extend ISDF to spin-polarised systems
1226 call hm%exxop%isdf%get_interpolation_points(namespace, space, ks%gr, st%rho(1:ks%gr%np, 1))
1227 call isdf_ace_compute_potentials(hm%exxop, namespace, space, ks%gr, &
1228 ks%calc%hf_st, xst, hm%kpoints)
1229 else
1230 call dexchange_operator_compute_potentials(hm%exxop, namespace, space, ks%gr, &
1231 ks%calc%hf_st, xst, hm%kpoints)
1232 endif
1233 exx_energy = dexchange_operator_compute_ex(ks%gr, ks%calc%hf_st, xst)
1234 call dexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst)
1235 else
1236 call zexchange_operator_compute_potentials(hm%exxop, namespace, space, ks%gr, &
1237 ks%calc%hf_st, xst, hm%kpoints)
1238 exx_energy = zexchange_operator_compute_ex(ks%gr, ks%calc%hf_st, xst)
1239 if (hm%phase%is_allocated()) then
1240 call zexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst, hm%phase)
1241 else
1242 call zexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst)
1243 end if
1244 end if
1245 call states_elec_end(xst)
1246 exx_energy = exx_energy + hm%exxop%singul%energy
1247 end if
1248
1249 ! Add the energy only the ACE case. In the non-ACE case, the singularity energy is added in energy_calc.F90
1250 select case (ks%theory_level)
1252 if (family_is_hybrid(ks%xc)) then
1253 hm%energy%exchange_hf = hm%energy%exchange_hf + exx_energy
1254 end if
1255 case (hartree_fock)
1256 hm%energy%exchange_hf = hm%energy%exchange_hf + exx_energy
1257 end select
1258 else
1259 ! If we ask for semilocal, we deactivate the exchange operator entirely
1260 call exchange_operator_reinit(hm%exxop, cam_null, ks%calc%hf_st)
1261 end if
1262 end if
1263
1264 end if
1265
1266 ! Because of the intent(in) in v_ks_calc_start, we need to update the parameters of hybrids for OEP
1267 ! here
1268 if (ks%theory_level == kohn_sham_dft .and. bitand(ks%xc_family, xc_family_oep) /= 0) then
1269 if (ks%xc%functional(func_x,1)%id /= xc_oep_x_slater .and. ks%xc%functional(func_x,1)%id /= xc_oep_x_fbe) then
1270 call exchange_operator_reinit(hm%exxop, ks%xc%cam)
1271 end if
1272 end if
1273
1274 if (ks%has_photons .and. (ks%xc_photon == 0)) then
1275 if (associated(ks%pt_mx%vmf)) then
1276 forall(ip = 1:ks%gr%np) hm%ks_pot%vhxc(ip, 1) = hm%ks_pot%vhxc(ip, 1) + ks%pt_mx%vmf(ip)
1277 if (hm%d%ispin > unpolarized) then
1278 forall(ip = 1:ks%gr%np) hm%ks_pot%vhxc(ip, 2) = hm%ks_pot%vhxc(ip, 2) + ks%pt_mx%vmf(ip)
1279 end if
1280 end if
1281 hm%ep%photon_forces(1:space%dim) = ks%pt_mx%fmf(1:space%dim)
1282 end if
1283
1284 if (ks%vdw%vdw_correction /= option__vdwcorrection__none) then
1285 assert(allocated(ks%vdw%forces))
1286 hm%ep%vdw_forces(:, :) = ks%vdw%forces(:, :)
1287 hm%ep%vdw_stress = ks%vdw%stress
1288 safe_deallocate_a(ks%vdw%forces)
1289 else
1290 hm%ep%vdw_forces = 0.0_real64
1291 end if
1292
1293 if (ks%calc%time_present .or. hm%time_zero) then
1294 call hm%update(ks%gr, namespace, space, ext_partners, time = ks%calc%time)
1295 else
1296 call hamiltonian_elec_update_pot(hm, ks%gr)
1297 end if
1298
1299
1300 safe_deallocate_a(ks%calc%density)
1301 if (ks%calc%total_density_alloc) then
1302 safe_deallocate_p(ks%calc%total_density)
1303 end if
1304 nullify(ks%calc%total_density)
1305
1306 pop_sub(v_ks_calc_finish)
1307 end subroutine v_ks_calc_finish
1308
1309
1312 !
1313 !! TODO(Alex) Once the implementation is finalised and benchmarked
1314 !! remove the serial version, and get rid of this routine.
1315 subroutine isdf_ace_compute_potentials(exxop, namespace, space, gr, hf_st, xst, kpoints)
1316 type(exchange_operator_t), intent(in ) :: exxop
1317 type(namespace_t), intent(in ) :: namespace
1318 class(space_t), intent(in ) :: space
1319 class(mesh_t), intent(in ) :: gr
1320 type(states_elec_t), intent(in ) :: hf_st
1321 type(kpoints_t), intent(in ) :: kpoints
1322
1323 type(states_elec_t), intent(inout) :: xst
1324
1325 if (exxop%isdf%use_serial) then
1326 call isdf_serial_ace_compute_potentials(exxop, namespace, space, gr, &
1327 hf_st, xst, kpoints)
1328 else
1329 call isdf_parallel_ace_compute_potentials(exxop, namespace, space, gr, &
1330 hf_st, xst, kpoints)
1331 endif
1332
1333 end subroutine isdf_ace_compute_potentials
1334
1335 ! ---------------------------------------------------------
1336 !
1340 !
1341 subroutine v_ks_hartree(namespace, ks, space, hm, ext_partners)
1342 type(namespace_t), intent(in) :: namespace
1343 type(v_ks_t), intent(inout) :: ks
1344 class(space_t), intent(in) :: space
1345 type(hamiltonian_elec_t), intent(inout) :: hm
1346 type(partner_list_t), intent(in) :: ext_partners
1347
1348 push_sub(v_ks_hartree)
1349
1350 if (.not. poisson_is_async(hm%psolver)) then
1351 ! solve the Poisson equation
1352 call dpoisson_solve(hm%psolver, namespace, hm%ks_pot%vhartree, ks%calc%total_density, reset=.false.)
1353 else
1354 ! The calculation was started by v_ks_calc_start.
1355 call dpoisson_solve_finish(hm%psolver, hm%ks_pot%vhartree)
1356 end if
1357
1358 if (ks%calc%calc_energy) then
1359 ! Get the Hartree energy
1360 hm%energy%hartree = m_half*dmf_dotp(ks%gr, ks%calc%total_density, hm%ks_pot%vhartree)
1361 end if
1362
1364 if(ks%calc%time_present) then
1365 if(hamiltonian_elec_has_kick(hm)) then
1366 call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
1367 ks%calc%total_density, hm%energy%pcm_corr, kick=hm%kick, time=ks%calc%time)
1368 else
1369 call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
1370 ks%calc%total_density, hm%energy%pcm_corr, time=ks%calc%time)
1371 end if
1372 else
1373 if(hamiltonian_elec_has_kick(hm)) then
1374 call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
1375 ks%calc%total_density, hm%energy%pcm_corr, kick=hm%kick)
1376 else
1377 call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
1378 ks%calc%total_density, hm%energy%pcm_corr)
1379 end if
1380 end if
1381
1382 pop_sub(v_ks_hartree)
1383 end subroutine v_ks_hartree
1384 ! ---------------------------------------------------------
1385
1386
1387 ! ---------------------------------------------------------
1388 subroutine v_ks_freeze_hxc(ks)
1389 type(v_ks_t), intent(inout) :: ks
1390
1391 push_sub(v_ks_freeze_hxc)
1392
1393 ks%frozen_hxc = .true.
1394
1395 pop_sub(v_ks_freeze_hxc)
1396 end subroutine v_ks_freeze_hxc
1397 ! ---------------------------------------------------------
1398
1399 subroutine v_ks_calculate_current(this, calc_cur)
1400 type(v_ks_t), intent(inout) :: this
1401 logical, intent(in) :: calc_cur
1402
1403 push_sub(v_ks_calculate_current)
1404
1405 this%calculate_current = calc_cur
1406
1407 pop_sub(v_ks_calculate_current)
1408 end subroutine v_ks_calculate_current
1409
1411 subroutine v_ks_update_dftu_energy(ks, namespace, hm, st, int_dft_u)
1412 type(v_ks_t), intent(inout) :: ks
1413 type(hamiltonian_elec_t), intent(in) :: hm
1414 type(namespace_t), intent(in) :: namespace
1415 type(states_elec_t), intent(inout) :: st
1416 real(real64), intent(out) :: int_dft_u
1417
1418 if (hm%lda_u_level == dft_u_none) return
1419
1420 push_sub(v_ks_update_dftu_energy)
1421
1422 if (states_are_real(st)) then
1423 int_dft_u = denergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_dft_u)
1424 else
1425 int_dft_u = zenergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_dft_u)
1426 end if
1427
1429 end subroutine v_ks_update_dftu_energy
1430end module v_ks_oct_m
1431
1432!! Local Variables:
1433!! mode: f90
1434!! coding: utf-8
1435!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:173
scales a vector by a constant
Definition: lalg_basic.F90:159
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:151
pure logical function, public accel_is_enabled()
Definition: accel.F90:418
subroutine, public current_calculate(this, namespace, gr, hm, space, st)
Compute total electronic current density.
Definition: current.F90:371
subroutine, public current_init(this, namespace)
Definition: current.F90:181
This module implements a calculator for the density and defines related functions.
Definition: density.F90:122
subroutine, public states_elec_total_density(st, mesh, total_rho)
This routine calculates the total electronic density.
Definition: density.F90:853
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
Definition: density.F90:612
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
subroutine, public energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full)
This subroutine calculates the total energy of the system. Basically, it adds up the KS eigenvalues,...
real(real64) function, public zenergy_calc_electronic(namespace, hm, der, st, terms)
real(real64) function, public denergy_calc_electronic(namespace, hm, der, st, terms)
subroutine, public energy_calc_eigenvalues(namespace, hm, der, st)
subroutine, public energy_copy(ein, eout)
Definition: energy.F90:170
subroutine, public dexchange_operator_ace(this, namespace, mesh, st, xst, phase)
subroutine, public zexchange_operator_compute_potentials(this, namespace, space, gr, st, xst, kpoints, F_out)
subroutine, public exchange_operator_reinit(this, cam, st)
subroutine, public dexchange_operator_compute_potentials(this, namespace, space, gr, st, xst, kpoints, F_out)
subroutine, public zexchange_operator_ace(this, namespace, mesh, st, xst, phase)
real(real64) function, public dexchange_operator_compute_ex(mesh, st, xst)
Compute the exact exchange energy.
real(real64) function, public zexchange_operator_compute_ex(mesh, st, xst)
Compute the exact exchange energy.
real(real64), parameter, public m_two
Definition: global.F90:192
real(real64), parameter, public m_zero
Definition: global.F90:190
integer, parameter, public rdmft
Definition: global.F90:236
integer, parameter, public hartree_fock
Definition: global.F90:236
integer, parameter, public independent_particles
Theory level.
Definition: global.F90:236
integer, parameter, public generalized_kohn_sham_dft
Definition: global.F90:236
integer, parameter, public kohn_sham_dft
Definition: global.F90:236
real(real64), parameter, public m_epsilon
Definition: global.F90:206
real(real64), parameter, public m_half
Definition: global.F90:196
real(real64), parameter, public m_one
Definition: global.F90:191
integer, parameter, public hartree
Definition: global.F90:236
This module implements the underlying real-space grid.
Definition: grid.F90:119
integer, parameter, public term_mgga
integer, parameter, public term_dft_u
logical function, public hamiltonian_elec_has_kick(hm)
logical function, public hamiltonian_elec_needs_current(hm, states_are_real)
subroutine, public hamiltonian_elec_update_pot(this, mesh, accumulate)
Update the KS potential of the electronic Hamiltonian.
This module defines classes and functions for interaction partners.
Interoperable Separable Density Fitting (ISDF) molecular implementation.
Definition: isdf.F90:116
subroutine, public isdf_ace_compute_potentials(exxop, namespace, space, mesh, st, Vx_on_st, kpoints)
ISDF wrapper computing interpolation points and vectors, which are used to build the potential used ...
Definition: isdf.F90:162
Serial prototype for benchmarking and validating ISDF implementation.
subroutine, public isdf_serial_ace_compute_potentials(exxop, namespace, space, mesh, st, Vx_on_st, kpoints)
ISDF wrapper computing interpolation points and vectors, which are used to build the potential used ...
A module to handle KS potential, without the external potential.
integer, parameter, public dft_u_none
Definition: lda_u.F90:203
This modules implements the routines for doing constrain DFT for noncollinear magnetism.
integer, parameter, public constrain_none
subroutine, public magnetic_constrain_update(this, mesh, std, space, latt, pos, rho)
Recomputes the magnetic contraining potential.
subroutine, public magnetic_induced(namespace, gr, st, psolver, kpoints, a_ind, b_ind)
This subroutine receives as input a current, and produces as an output the vector potential that it i...
Definition: magnetic.F90:371
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:904
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1097
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:531
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1029
subroutine, public messages_new_line()
Definition: messages.F90:1118
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
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1069
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:600
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:505
subroutine, public pcm_hartree_potential(pcm, space, mesh, psolver, ext_partners, vhartree, density, pcm_corr, kick, time)
PCM reaction field due to the electronic density.
subroutine, public mf_calc(this, gr, st, ions, pt_mode, time)
subroutine, public dpoisson_solve_start(this, rho)
Definition: poisson.F90:2025
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
subroutine, public dpoisson_solve_finish(this, pot)
Definition: poisson.F90:2033
logical pure function, public poisson_is_async(this)
Definition: poisson.F90:1158
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
integer, parameter, public pseudo_exchange_unknown
Definition: pseudo.F90:190
integer, parameter, public pseudo_correlation_unknown
Definition: pseudo.F90:194
integer, parameter, public pseudo_correlation_any
Definition: pseudo.F90:194
integer, parameter, public pseudo_exchange_any
Definition: pseudo.F90:190
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:119
integer, parameter, private libxc_c_index
Definition: species.F90:278
pure logical function, public states_are_complex(st)
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_fermi(st, namespace, mesh, compute_spin)
calculate the Fermi level for the states in this object
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
subroutine, public states_elec_allocate_current(st, space, mesh)
This module provides routines for communicating states when using states parallelization.
subroutine, public states_elec_parallel_remote_access_stop(this)
stop remote memory access for states on other processors
subroutine, public states_elec_parallel_remote_access_start(this)
start remote memory access for states on other processors
subroutine v_ks_hartree(namespace, ks, space, hm, ext_partners)
Hartree contribution to the KS potential. This function is designed to be used by v_ks_calc_finish an...
Definition: v_ks.F90:1437
subroutine, public v_ks_calc_finish(ks, hm, namespace, space, latt, st, ext_partners, force_semilocal)
Definition: v_ks.F90:1102
subroutine, public v_ks_freeze_hxc(ks)
Definition: v_ks.F90:1484
subroutine, public v_ks_end(ks)
Definition: v_ks.F90:628
subroutine, public v_ks_calculate_current(this, calc_cur)
Definition: v_ks.F90:1495
subroutine, public v_ks_write_info(ks, iunit, namespace)
Definition: v_ks.F90:660
subroutine, public v_ks_update_dftu_energy(ks, namespace, hm, st, int_dft_u)
Update the value of <\psi | V_U | \psi>, where V_U is the DFT+U potential.
Definition: v_ks.F90:1507
subroutine, public v_ks_calc_start(ks, namespace, space, hm, st, ions, latt, ext_partners, time, calc_energy, force_semilocal)
This routine starts the calculation of the Kohn-Sham potential. The routine v_ks_calc_finish must be ...
Definition: v_ks.F90:804
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Definition: v_ks.F90:749
subroutine, public v_ks_h_setup(namespace, space, gr, ions, ext_partners, st, ks, hm, calc_eigenval, calc_current)
Definition: v_ks.F90:695
subroutine, public v_ks_init(ks, namespace, gr, st, ions, mc, space, kpoints)
Definition: v_ks.F90:254
subroutine, public x_slater_calc(namespace, gr, space, exxop, st, kpoints, ex, vxc)
Interface to X(slater_calc)
Definition: x_slater.F90:147
type(xc_cam_t), parameter, public cam_null
All CAM parameters set to zero.
Definition: xc_cam.F90:152
type(xc_cam_t), parameter, public cam_exact_exchange
Use only Hartree Fock exact exchange.
Definition: xc_cam.F90:155
subroutine, public x_fbe_calc(id, namespace, psolver, gr, st, space, ex, vxc)
Interface to X(x_fbe_calc) Two possible run modes possible: adiabatic and Sturm-Liouville....
Definition: xc_fbe.F90:172
subroutine, public fbe_c_lda_sl(namespace, psolver, gr, st, space, ec, vxc)
Sturm-Liouville version of the FBE local-density correlation functional.
Definition: xc_fbe.F90:455
integer, parameter, public xc_family_ks_inversion
declaring 'family' constants for 'functionals' not handled by libxc careful not to use a value define...
integer function, public xc_get_default_functional(dim, pseudo_x_functional, pseudo_c_functional)
Returns the default functional given the one parsed from the pseudopotentials and the space dimension...
integer, parameter, public xc_family_nc_mgga
integer, parameter, public xc_oep_x
Exact exchange.
integer, parameter, public xc_lda_c_fbe_sl
LDA correlation based ib the force-balance equation - Sturm-Liouville version.
integer, parameter, public xc_family_nc_lda
integer, parameter, public xc_oep_x_fbe_sl
Exchange approximation based on the force balance equation - Sturn-Liouville version.
integer, parameter, public xc_oep_x_fbe
Exchange approximation based on the force balance equation.
integer, parameter, public xc_oep_x_slater
Slater approximation to the exact exchange.
integer, parameter, public func_c
integer, parameter, public func_x
subroutine, public xc_ks_inversion_end(ks_inv)
subroutine, public xc_ks_inversion_write_info(ks_inversion, iunit, namespace)
subroutine, public xc_ks_inversion_init(ks_inv, namespace, gr, ions, st, xc, mc, space, kpoints)
subroutine, public xc_ks_inversion_calc(ks_inversion, namespace, space, gr, hm, ext_partners, st, vxc, time)
Definition: xc.F90:116
subroutine, public xc_write_info(xcs, iunit, namespace)
Definition: xc.F90:226
subroutine, public xc_init(xcs, namespace, ndim, periodic_dim, nel, x_id, c_id, xk_id, ck_id, hartree_fock, ispin)
Definition: xc.F90:256
pure logical function, public family_is_mgga(family, only_collinear)
Is the xc function part of the mGGA family.
Definition: xc.F90:607
logical pure function, public family_is_mgga_with_exc(xcs)
Is the xc function part of the mGGA family with an energy functional.
Definition: xc.F90:620
subroutine, public xc_end(xcs)
Definition: xc.F90:545
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 pure function, public family_is_hybrid(xcs)
Returns true if the functional is an hybrid functional.
Definition: xc.F90:635
subroutine, public xc_get_nc_vxc(gr, xcs, st, kpoints, space, namespace, rho, vxc, ex, ec, vtau, ex_density, ec_density)
This routines is similar to xc_get_vxc but for noncollinear functionals, which are not implemented in...
Definition: xc.F90:2557
integer, parameter, public oep_type_mgga
Definition: xc_oep.F90:186
integer, parameter, public oep_level_none
the OEP levels
Definition: xc_oep.F90:174
subroutine, public xc_oep_end(oep)
Definition: xc_oep.F90:358
subroutine, public zxc_oep_calc(oep, namespace, xcs, gr, hm, st, space, rcell_volume, ex, ec, vxc)
This file handles the evaluation of the OEP potential, in the KLI or full OEP as described in S....
Definition: xc_oep.F90:2306
subroutine, public dxc_oep_calc(oep, namespace, xcs, gr, hm, st, space, rcell_volume, ex, ec, vxc)
This file handles the evaluation of the OEP potential, in the KLI or full OEP as described in S....
Definition: xc_oep.F90:1447
subroutine, public xc_oep_write_info(oep, iunit, namespace)
Definition: xc_oep.F90:380
integer, parameter, public oep_type_exx
The different types of OEP that we can work with.
Definition: xc_oep.F90:186
subroutine, public xc_oep_init(oep, namespace, gr, st, mc, space, oep_type)
Definition: xc_oep.F90:219
subroutine, public zxc_oep_photon_calc(oep, namespace, xcs, gr, hm, st, space, ex, ec, vxc)
This file handles the evaluation of the OEP potential, in the KLI or full OEP as described in S....
subroutine, public dxc_oep_photon_calc(oep, namespace, xcs, gr, hm, st, space, ex, ec, vxc)
This file handles the evaluation of the OEP potential, in the KLI or full OEP as described in S....
This module implements the "photon-free" electron-photon exchange-correlation functional.
Definition: xc_photons.F90:123
integer, parameter, public sic_none
no self-interaction correction
Definition: xc_sic.F90:150
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
A module that takes care of xc contribution from vdW interactions.
Definition: xc_vdw.F90:118
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
Describes mesh distribution to nodes.
Definition: mesh.F90:187
The states_elec_t class contains all electronic wave functions.
int true(void)
subroutine get_functional_from_pseudos(x_functional, c_functional)
Tries to find out the functional from the pseudopotential.
Definition: v_ks.F90:570
subroutine v_a_xc(hm, force_semilocal)
Definition: v_ks.F90:948
subroutine calculate_density()
Definition: v_ks.F90:911