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