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