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, ks%gr, st, space, ks%calc%energy%correlation, vxc = ks%calc%vxc)
965
966 end if
967
968 end if
969
970 if (bitand(ks%xc_family, xc_family_ks_inversion) /= 0) then
971 ! Also treat KS inversion separately (not part of libxc)
972 call xc_ks_inversion_calc(ks%ks_inversion, namespace, space, ks%gr, hm, ext_partners, st, vxc = ks%calc%vxc, &
973 time = ks%calc%time)
974 end if
975
976 ! compute the photon-free photon exchange potential and energy
977 if (ks%xc_photon /= 0) then
978
979 call ks%xc_photons%v_ks(namespace, ks%calc%total_density, ks%gr, space, hm%psolver, st)
980
981 ! add the photon-free px potential into the xc potential
982 do ispin = 1, hm%d%spin_channels
983 call lalg_axpy(ks%gr%np, m_one, ks%xc_photons%vpx(1:ks%gr%np), ks%calc%vxc(1:ks%gr%np, ispin) )
984 end do
985
986 ! photon-exchange energy
987 ks%calc%energy%photon_exchange = ks%xc_photons%ex
988 end if
989
990 end if
991
992 if (ks%calc%calc_energy) then
993 ! MGGA vtau contribution is done after copying vtau to hm%vtau
994
995 call v_ks_update_dftu_energy(ks, namespace, hm, st, ks%calc%energy%int_dft_u)
996 end if
997
998 call profiling_out("XC")
999 pop_sub(v_ks_calc_start.v_a_xc)
1000 end subroutine v_a_xc
1001
1002 end subroutine v_ks_calc_start
1003 ! ---------------------------------------------------------
1004
1005 subroutine v_ks_calc_finish(ks, hm, namespace, space, latt, st, ext_partners, force_semilocal)
1006 type(v_ks_t), target, intent(inout) :: ks
1007 type(hamiltonian_elec_t), intent(inout) :: hm
1008 type(namespace_t), intent(in) :: namespace
1009 class(space_t), intent(in) :: space
1010 type(lattice_vectors_t), intent(in) :: latt
1011 type(states_elec_t), intent(inout) :: st
1012 type(partner_list_t), intent(in) :: ext_partners
1013 logical, optional, intent(in) :: force_semilocal
1014
1015 integer :: ip, ispin
1016 type(states_elec_t) :: xst
1018 real(real64) :: exx_energy
1019 real(real64) :: factor
1020
1021 push_sub(v_ks_calc_finish)
1022
1023 assert(ks%calc%calculating)
1024 ks%calc%calculating = .false.
1025
1026 if (ks%frozen_hxc) then
1027 pop_sub(v_ks_calc_finish)
1028 return
1029 end if
1030
1031 !change the pointer to the energy object
1032 safe_deallocate_a(hm%energy)
1033 call move_alloc(ks%calc%energy, hm%energy)
1034
1035 if (hm%self_induced_magnetic) then
1036 hm%a_ind(1:ks%gr%np, 1:space%dim) = ks%calc%a_ind(1:ks%gr%np, 1:space%dim)
1037 hm%b_ind(1:ks%gr%np, 1:space%dim) = ks%calc%b_ind(1:ks%gr%np, 1:space%dim)
1038
1039 safe_deallocate_a(ks%calc%a_ind)
1040 safe_deallocate_a(ks%calc%b_ind)
1041 end if
1042
1043 if (allocated(hm%v_static)) then
1044 hm%energy%intnvstatic = dmf_dotp(ks%gr, ks%calc%total_density, hm%v_static)
1045 else
1046 hm%energy%intnvstatic = m_zero
1047 end if
1048
1049 if (ks%theory_level == independent_particles .or. abs(ks%sic%amaldi_factor) <= m_epsilon) then
1050
1051 hm%ks_pot%vhxc = m_zero
1052 hm%energy%intnvxc = m_zero
1053 hm%energy%hartree = m_zero
1054 hm%energy%exchange = m_zero
1055 hm%energy%exchange_hf = m_zero
1056 hm%energy%correlation = m_zero
1057 else
1058
1059 hm%energy%hartree = m_zero
1060 call v_ks_hartree(namespace, ks, space, hm, ext_partners)
1061
1062 if (.not. optional_default(force_semilocal, .false.)) then
1063 !PZ-SIC
1064 if(ks%sic%level == sic_pz_oep) then
1065 if (states_are_real(st)) then
1066 call dxc_oep_calc(ks%sic%oep, namespace, ks%xc, ks%gr, hm, st, space, &
1067 latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1068 else
1069 call zxc_oep_calc(ks%sic%oep, namespace, ks%xc, ks%gr, hm, st, space, &
1070 latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1071 end if
1072 end if
1073
1074 ! OEP for exchange ad MGGAs (within Kohn-Sham DFT)
1075 if (ks%theory_level == kohn_sham_dft .and. ks%oep%level /= oep_level_none) then
1076 ! The OEP family has to be handled specially
1077 if (ks%xc%functional(func_x,1)%id == xc_oep_x .or. family_is_mgga_with_exc(ks%xc)) then
1078 if (states_are_real(st)) then
1079 call dxc_oep_calc(ks%oep, namespace, ks%xc, ks%gr, hm, st, space, &
1080 latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1081 else
1082 call zxc_oep_calc(ks%oep, namespace, ks%xc, ks%gr, hm, st, space, &
1083 latt%rcell_volume, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1084 end if
1085 end if
1086 end if
1087 end if
1088
1089 if (ks%theory_level == kohn_sham_dft .and. ks%oep_photon%level /= oep_level_none) then
1090 if (states_are_real(st)) then
1091 call dxc_oep_photon_calc(ks%oep_photon, namespace, ks%xc, ks%gr, &
1092 hm, st, space, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1093 else
1094 call zxc_oep_photon_calc(ks%oep_photon, namespace, ks%xc, ks%gr, &
1095 hm, st, space, hm%energy%exchange, hm%energy%correlation, vxc = ks%calc%vxc)
1096 end if
1097 hm%energy%photon_exchange = ks%oep_photon%pt%ex
1098 end if
1099
1101 if (ks%calc%calc_energy) then
1102 ! Now we calculate Int[n vxc] = energy%intnvxc
1103 hm%energy%intnvxc = m_zero
1104
1105 if (ks%theory_level /= independent_particles .and. ks%theory_level /= hartree .and. ks%theory_level /= rdmft) then
1106 do ispin = 1, hm%d%nspin
1107 if (ispin <= 2) then
1108 factor = m_one
1109 else
1110 factor = m_two
1111 end if
1112 hm%energy%intnvxc = hm%energy%intnvxc + &
1113 factor*dmf_dotp(ks%gr, st%rho(:, ispin), ks%calc%vxc(:, ispin), reduce = .false.)
1114 end do
1115 call ks%gr%allreduce(hm%energy%intnvxc)
1116 end if
1117 end if
1118
1119
1120 if (ks%theory_level /= hartree .and. ks%theory_level /= rdmft) then
1121 ! move allocation of vxc from ks%calc to hm
1122 safe_deallocate_a(hm%ks_pot%vxc)
1123 call move_alloc(ks%calc%vxc, hm%ks_pot%vxc)
1124
1125 if (family_is_mgga_with_exc(hm%xc)) then
1126 call hm%ks_pot%set_vtau(ks%calc%vtau)
1127 safe_deallocate_a(ks%calc%vtau)
1128
1129 ! We need to evaluate the energy after copying vtau to hm%vtau
1130 if (ks%theory_level == generalized_kohn_sham_dft .and. ks%calc%calc_energy) then
1131 ! MGGA vtau contribution
1132 if (states_are_real(st)) then
1133 hm%energy%intnvxc = hm%energy%intnvxc &
1134 + denergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_mgga)
1135 else
1136 hm%energy%intnvxc = hm%energy%intnvxc &
1137 + zenergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_mgga)
1138 end if
1139 end if
1140 end if
1141
1142 else
1143 hm%ks_pot%vxc = m_zero
1144 end if
1145
1146 if (.not. ks%xc_photon_include_hartree) then
1147 hm%energy%hartree = m_zero
1148 hm%ks_pot%vhartree = m_zero
1149 end if
1150
1151 ! Build Hartree + XC potential
1152
1153 do ip = 1, ks%gr%np
1154 hm%ks_pot%vhxc(ip, 1) = hm%ks_pot%vxc(ip, 1) + hm%ks_pot%vhartree(ip)
1155 end do
1156 if (allocated(hm%vberry)) then
1157 do ip = 1, ks%gr%np
1158 hm%ks_pot%vhxc(ip, 1) = hm%ks_pot%vhxc(ip, 1) + hm%vberry(ip, 1)
1159 end do
1160 end if
1161
1162 if (hm%d%ispin > unpolarized) then
1163 do ip = 1, ks%gr%np
1164 hm%ks_pot%vhxc(ip, 2) = hm%ks_pot%vxc(ip, 2) + hm%ks_pot%vhartree(ip)
1165 end do
1166 if (allocated(hm%vberry)) then
1167 do ip = 1, ks%gr%np
1168 hm%ks_pot%vhxc(ip, 2) = hm%ks_pot%vhxc(ip, 2) + hm%vberry(ip, 2)
1169 end do
1170 end if
1171 end if
1172
1173 if (hm%d%ispin == spinors) then
1174 do ispin=3, 4
1175 do ip = 1, ks%gr%np
1176 hm%ks_pot%vhxc(ip, ispin) = hm%ks_pot%vxc(ip, ispin)
1177 end do
1178 end do
1179 end if
1180
1181 ! Note: this includes hybrids calculated with the Fock operator instead of OEP
1182 hm%energy%exchange_hf = m_zero
1183 if (ks%theory_level == hartree .or. ks%theory_level == hartree_fock &
1184 .or. ks%theory_level == rdmft &
1185 .or. (ks%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(ks%xc))) then
1186
1187 ! swap the states object
1188 if (.not. hm%exxop%useACE) then
1189 ! We also close the MPI remote memory access to the old object
1190 if (associated(hm%exxop%st)) then
1192 call states_elec_end(hm%exxop%st)
1193 safe_deallocate_p(hm%exxop%st)
1194 end if
1195 ! We activate the MPI remote memory access for ks%calc%hf_st
1196 ! This allows to have all calls to exchange_operator_apply_standard to access
1197 ! the states over MPI
1199 end if
1200
1201 ! The exchange operator will use ks%calc%hf_st
1202 ! For the ACE case, this is the same as st
1203 if (.not. optional_default(force_semilocal, .false.)) then
1204 select case (ks%theory_level)
1206 if (family_is_hybrid(ks%xc)) then
1207 call exchange_operator_reinit(hm%exxop, ks%xc%cam, ks%calc%hf_st)
1208 end if
1209 case (hartree_fock)
1210 call exchange_operator_reinit(hm%exxop, ks%xc%cam, ks%calc%hf_st)
1211 case (hartree, rdmft)
1212 call exchange_operator_reinit(hm%exxop, cam_exact_exchange, ks%calc%hf_st)
1213 end select
1214
1215 !This should be changed and the CAM parameters should also be obtained from the restart information
1216 !Maybe the parameters should be mixed too.
1217 exx_energy = m_zero
1218 if (hm%exxop%useACE) then
1219 call xst%nullify()
1220 if (states_are_real(ks%calc%hf_st)) then
1221 ! TODO(Alex) Clean up nested if statements
1222 if (hm%exxop%with_isdf) then
1223 ! Find interpolation points from density, or read from file
1224 ! TODO(Alex) Issue #1195 Extend ISDF to spin-polarised systems
1225 call hm%exxop%isdf%get_interpolation_points(namespace, space, ks%gr, st%rho(1:ks%gr%np, 1))
1226 call isdf_ace_compute_potentials(hm%exxop, namespace, space, ks%gr, &
1227 ks%calc%hf_st, xst, hm%kpoints)
1228 else
1229 call dexchange_operator_compute_potentials(hm%exxop, namespace, space, ks%gr, &
1230 ks%calc%hf_st, xst, hm%kpoints)
1231 endif
1232 exx_energy = dexchange_operator_compute_ex(ks%gr, ks%calc%hf_st, xst)
1233 call dexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst)
1234 else
1235 call zexchange_operator_compute_potentials(hm%exxop, namespace, space, ks%gr, &
1236 ks%calc%hf_st, xst, hm%kpoints)
1237 exx_energy = zexchange_operator_compute_ex(ks%gr, ks%calc%hf_st, xst)
1238 if (hm%phase%is_allocated()) then
1239 call zexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst, hm%phase)
1240 else
1241 call zexchange_operator_ace(hm%exxop, namespace, ks%gr, ks%calc%hf_st, xst)
1242 end if
1243 end if
1244 call states_elec_end(xst)
1245 exx_energy = exx_energy + hm%exxop%singul%energy
1246 end if
1247
1248 ! Add the energy only the ACE case. In the non-ACE case, the singularity energy is added in energy_calc.F90
1249 select case (ks%theory_level)
1251 if (family_is_hybrid(ks%xc)) then
1252 hm%energy%exchange_hf = hm%energy%exchange_hf + exx_energy
1253 end if
1254 case (hartree_fock)
1255 hm%energy%exchange_hf = hm%energy%exchange_hf + exx_energy
1256 end select
1257 else
1258 ! If we ask for semilocal, we deactivate the exchange operator entirely
1259 call exchange_operator_reinit(hm%exxop, cam_null, ks%calc%hf_st)
1260 end if
1261 end if
1262
1263 end if
1264
1265 ! Because of the intent(in) in v_ks_calc_start, we need to update the parameters of hybrids for OEP
1266 ! here
1267 if (ks%theory_level == kohn_sham_dft .and. bitand(ks%xc_family, xc_family_oep) /= 0) then
1268 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
1269 call exchange_operator_reinit(hm%exxop, ks%xc%cam)
1270 end if
1271 end if
1272
1273 if (ks%has_photons .and. (ks%xc_photon == 0)) then
1274 if (associated(ks%pt_mx%vmf)) then
1275 forall(ip = 1:ks%gr%np) hm%ks_pot%vhxc(ip, 1) = hm%ks_pot%vhxc(ip, 1) + ks%pt_mx%vmf(ip)
1276 if (hm%d%ispin > unpolarized) then
1277 forall(ip = 1:ks%gr%np) hm%ks_pot%vhxc(ip, 2) = hm%ks_pot%vhxc(ip, 2) + ks%pt_mx%vmf(ip)
1278 end if
1279 end if
1280 hm%ep%photon_forces(1:space%dim) = ks%pt_mx%fmf(1:space%dim)
1281 end if
1282
1283 if (ks%vdw%vdw_correction /= option__vdwcorrection__none) then
1284 assert(allocated(ks%vdw%forces))
1285 hm%ep%vdw_forces(:, :) = ks%vdw%forces(:, :)
1286 hm%ep%vdw_stress = ks%vdw%stress
1287 safe_deallocate_a(ks%vdw%forces)
1288 else
1289 hm%ep%vdw_forces = 0.0_real64
1290 end if
1291
1292 if (ks%calc%time_present .or. hm%time_zero) then
1293 call hm%update(ks%gr, namespace, space, ext_partners, time = ks%calc%time)
1294 else
1295 call hamiltonian_elec_update_pot(hm, ks%gr)
1296 end if
1297
1298
1299 safe_deallocate_a(ks%calc%density)
1300 if (ks%calc%total_density_alloc) then
1301 safe_deallocate_p(ks%calc%total_density)
1302 end if
1303 nullify(ks%calc%total_density)
1304
1305 pop_sub(v_ks_calc_finish)
1306 end subroutine v_ks_calc_finish
1307
1308
1311 !
1312 !! TODO(Alex) Once the implementation is finalised and benchmarked
1313 !! remove the serial version, and get rid of this routine.
1314 subroutine isdf_ace_compute_potentials(exxop, namespace, space, gr, hf_st, xst, kpoints)
1315 type(exchange_operator_t), intent(in ) :: exxop
1316 type(namespace_t), intent(in ) :: namespace
1317 class(space_t), intent(in ) :: space
1318 class(mesh_t), intent(in ) :: gr
1319 type(states_elec_t), intent(in ) :: hf_st
1320 type(kpoints_t), intent(in ) :: kpoints
1321
1322 type(states_elec_t), intent(inout) :: xst
1323
1324 if (exxop%isdf%use_serial) then
1325 call isdf_serial_ace_compute_potentials(exxop, namespace, space, gr, &
1326 hf_st, xst, kpoints)
1327 else
1328 call isdf_parallel_ace_compute_potentials(exxop, namespace, space, gr, &
1329 hf_st, xst, kpoints)
1330 endif
1331
1332 end subroutine isdf_ace_compute_potentials
1333
1334 ! ---------------------------------------------------------
1335 !
1339 !
1340 subroutine v_ks_hartree(namespace, ks, space, hm, ext_partners)
1341 type(namespace_t), intent(in) :: namespace
1342 type(v_ks_t), intent(inout) :: ks
1343 class(space_t), intent(in) :: space
1344 type(hamiltonian_elec_t), intent(inout) :: hm
1345 type(partner_list_t), intent(in) :: ext_partners
1346
1347 push_sub(v_ks_hartree)
1348
1349 if (.not. poisson_is_async(hm%psolver)) then
1350 ! solve the Poisson equation
1351 call dpoisson_solve(hm%psolver, namespace, hm%ks_pot%vhartree, ks%calc%total_density, reset=.false.)
1352 else
1353 ! The calculation was started by v_ks_calc_start.
1354 call dpoisson_solve_finish(hm%psolver, hm%ks_pot%vhartree)
1355 end if
1356
1357 if (ks%calc%calc_energy) then
1358 ! Get the Hartree energy
1359 hm%energy%hartree = m_half*dmf_dotp(ks%gr, ks%calc%total_density, hm%ks_pot%vhartree)
1360 end if
1361
1363 if(ks%calc%time_present) then
1364 if(hamiltonian_elec_has_kick(hm)) then
1365 call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
1366 ks%calc%total_density, hm%energy%pcm_corr, kick=hm%kick, time=ks%calc%time)
1367 else
1368 call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
1369 ks%calc%total_density, hm%energy%pcm_corr, time=ks%calc%time)
1370 end if
1371 else
1372 if(hamiltonian_elec_has_kick(hm)) then
1373 call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
1374 ks%calc%total_density, hm%energy%pcm_corr, kick=hm%kick)
1375 else
1376 call pcm_hartree_potential(hm%pcm, space, ks%gr, hm%psolver, ext_partners, hm%ks_pot%vhartree, &
1377 ks%calc%total_density, hm%energy%pcm_corr)
1378 end if
1379 end if
1380
1381 pop_sub(v_ks_hartree)
1382 end subroutine v_ks_hartree
1383 ! ---------------------------------------------------------
1384
1385
1386 ! ---------------------------------------------------------
1387 subroutine v_ks_freeze_hxc(ks)
1388 type(v_ks_t), intent(inout) :: ks
1389
1390 push_sub(v_ks_freeze_hxc)
1391
1392 ks%frozen_hxc = .true.
1393
1394 pop_sub(v_ks_freeze_hxc)
1395 end subroutine v_ks_freeze_hxc
1396 ! ---------------------------------------------------------
1397
1398 subroutine v_ks_calculate_current(this, calc_cur)
1399 type(v_ks_t), intent(inout) :: this
1400 logical, intent(in) :: calc_cur
1401
1402 push_sub(v_ks_calculate_current)
1403
1404 this%calculate_current = calc_cur
1405
1406 pop_sub(v_ks_calculate_current)
1407 end subroutine v_ks_calculate_current
1408
1410 subroutine v_ks_update_dftu_energy(ks, namespace, hm, st, int_dft_u)
1411 type(v_ks_t), intent(inout) :: ks
1412 type(hamiltonian_elec_t), intent(in) :: hm
1413 type(namespace_t), intent(in) :: namespace
1414 type(states_elec_t), intent(inout) :: st
1415 real(real64), intent(out) :: int_dft_u
1416
1417 if (hm%lda_u_level == dft_u_none) return
1418
1419 push_sub(v_ks_update_dftu_energy)
1420
1421 if (states_are_real(st)) then
1422 int_dft_u = denergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_dft_u)
1423 else
1424 int_dft_u = zenergy_calc_electronic(namespace, hm, ks%gr%der, st, terms = term_dft_u)
1425 end if
1426
1428 end subroutine v_ks_update_dftu_energy
1429end module v_ks_oct_m
1430
1431!! Local Variables:
1432!! mode: f90
1433!! coding: utf-8
1434!! 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:419
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:193
real(real64), parameter, public m_zero
Definition: global.F90:191
integer, parameter, public rdmft
Definition: global.F90:237
integer, parameter, public hartree_fock
Definition: global.F90:237
integer, parameter, public independent_particles
Theory level.
Definition: global.F90:237
integer, parameter, public generalized_kohn_sham_dft
Definition: global.F90:237
integer, parameter, public kohn_sham_dft
Definition: global.F90:237
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_half
Definition: global.F90:197
real(real64), parameter, public m_one
Definition: global.F90:192
integer, parameter, public hartree
Definition: global.F90:237
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:898
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1023
subroutine, public messages_new_line()
Definition: messages.F90:1112
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
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1063
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
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:455
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:280
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:1436
subroutine, public v_ks_calc_finish(ks, hm, namespace, space, latt, st, ext_partners, force_semilocal)
Definition: v_ks.F90:1101
subroutine, public v_ks_freeze_hxc(ks)
Definition: v_ks.F90:1483
subroutine, public v_ks_end(ks)
Definition: v_ks.F90:628
subroutine, public v_ks_calculate_current(this, calc_cur)
Definition: v_ks.F90:1494
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:1506
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, gr, st, space, ec, vxc)
Sturm-Liouville version of the FBE local-density correlation functional.
Definition: xc_fbe.F90:456
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:276
pure logical function, public family_is_mgga(family, only_collinear)
Is the xc function part of the mGGA family.
Definition: xc.F90:580
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:593
subroutine, public xc_end(xcs)
Definition: xc.F90:518
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 pure function, public family_is_hybrid(xcs)
Returns true if the functional is an hybrid functional.
Definition: xc.F90:608
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:2530
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:151
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
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