Octopus
hamiltonian_elec.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2020 M. Marques, A. Castro, A. Rubio, G. Bertsch,
2!! N. Tancogne-Dejean, M. Lueders
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
24 use accel_oct_m
26 use batch_oct_m
29 use comm_oct_m
30 use debug_oct_m
33 use energy_oct_m
38 use epot_oct_m
41 use global_oct_m
42 use grid_oct_m
46 use io_oct_m
47 use ions_oct_m
48 use kick_oct_m
49 use, intrinsic :: iso_fortran_env
53 use lasers_oct_m
55 use lda_u_oct_m
58 use math_oct_m
59 use mesh_oct_m
62 use mpi_oct_m
66 use nlcc_oct_m
70 use parser_oct_m
75 use pcm_oct_m
76 use phase_oct_m
79 use space_oct_m
87 use types_oct_m
88 use unit_oct_m
92 use xc_oct_m
93 use xc_cam_oct_m
94 use xc_f03_lib_m
98 use zora_oct_m
99
100 implicit none
101
102 private
103 public :: &
113 magnus, &
114 dvmask, &
115 zvmask, &
130
131
132 type, extends(hamiltonian_abst_t) :: hamiltonian_elec_t
133 ! Components are public by default
134
137 type(space_t), private :: space
138 type(states_elec_dim_t) :: d
139 type(hamiltonian_elec_base_t) :: hm_base
140 type(phase_t) :: phase
141 type(energy_t), allocatable :: energy
142 type(absorbing_boundaries_t) :: abs_boundaries
143 type(ks_potential_t) :: ks_pot
144 real(real64), allocatable :: vberry(:,:)
145
146 type(derivatives_t), pointer, private :: der
147
148 type(nonlocal_pseudopotential_t) :: vnl
149
150 type(ions_t), pointer :: ions
151 real(real64) :: exx_coef
152
153 type(poisson_t) :: psolver
154
156 logical :: self_induced_magnetic
157 real(real64), allocatable :: a_ind(:, :)
158 real(real64), allocatable :: b_ind(:, :)
159
160 integer :: theory_level
161 type(xc_t), pointer :: xc
162 type(xc_photons_t), pointer :: xc_photons
163
164 type(epot_t) :: ep
165 type(pcm_t) :: pcm
166
168 logical, private :: adjoint
169
171 real(real64), private :: mass
172
174 logical, private :: inh_term
175 type(states_elec_t) :: inh_st
176
179 type(oct_exchange_t) :: oct_exchange
180
181 type(scissor_t) :: scissor
182
183 real(real64) :: current_time
184 logical, private :: is_applied_packed
185
187 type(lda_u_t) :: lda_u
188 integer :: lda_u_level
189
190 logical, public :: time_zero
191
192 type(exchange_operator_t), public :: exxop
193
194 type(kpoints_t), pointer, public :: kpoints => null()
195
196 type(partner_list_t) :: external_potentials
197 real(real64), allocatable, public :: v_ext_pot(:)
198 real(real64), allocatable, public :: v_static(:)
199
200 type(ion_electron_local_potential_t) :: v_ie_loc
201 type(nlcc_t) :: nlcc
202
203 type(magnetic_constrain_t) :: magnetic_constrain
204
206 type(kick_t) :: kick
207
209 type(mxll_coupling_t) :: mxll
210 type(zora_t), pointer :: zora => null()
211
212 contains
213 procedure :: update => hamiltonian_elec_update
214 procedure :: apply_packed => hamiltonian_elec_apply_packed
215 procedure :: update_span => hamiltonian_elec_span
216 procedure :: dapply => dhamiltonian_elec_apply
217 procedure :: zapply => zhamiltonian_elec_apply
218 procedure :: dmagnus_apply => dhamiltonian_elec_magnus_apply
219 procedure :: zmagnus_apply => zhamiltonian_elec_magnus_apply
220 procedure :: is_hermitian => hamiltonian_elec_hermitian
221 procedure :: needs_mgga_term => hamiltonian_elec_needs_mgga_term
222 procedure :: set_mass => hamiltonian_elec_set_mass
223 end type hamiltonian_elec_t
224
225 integer, public, parameter :: &
226 LENGTH = 1, &
228
229
230contains
231
232 ! ---------------------------------------------------------
233 subroutine hamiltonian_elec_init(hm, namespace, space, gr, ions, ext_partners, st, theory_level, xc, &
234 mc, kpoints, need_exchange, xc_photons)
235 type(hamiltonian_elec_t), target, intent(inout) :: hm
236 type(namespace_t), intent(in) :: namespace
237 class(space_t), intent(in) :: space
238 type(grid_t), target, intent(inout) :: gr
239 type(ions_t), target, intent(inout) :: ions
240 type(partner_list_t), intent(inout) :: ext_partners
241 type(states_elec_t), target, intent(inout) :: st
242 integer, intent(in) :: theory_level
243 type(xc_t), target, intent(in) :: xc
244 type(multicomm_t), intent(in) :: mc
245 type(kpoints_t), target, intent(in) :: kpoints
246 logical, optional, intent(in) :: need_exchange
247 type(xc_photons_t), optional, target, intent(in) :: xc_photons
249
250 logical :: need_exchange_
251 real(real64) :: rashba_coupling
254 call profiling_in('HAMILTONIAN_ELEC_INIT')
256 ! make a couple of local copies
257 hm%space = space
258 hm%theory_level = theory_level
259 call states_elec_dim_copy(hm%d, st%d)
261 hm%kpoints => kpoints
262
263 !%Variable ParticleMass
264 !%Type float
265 !%Default 1.0
266 !%Section Hamiltonian
267 !%Description
268 !% It is possible to make calculations for a particle with a mass
269 !% different from one (atomic unit of mass, or mass of the electron).
270 !% This is useful to describe non-electronic systems, or for
271 !% esoteric purposes.
272 !%End
273 call parse_variable(namespace, 'ParticleMass', m_one, hm%mass)
275 !%Variable RashbaSpinOrbitCoupling
276 !%Type float
277 !%Default 0.0
278 !%Section Hamiltonian
279 !%Description
280 !% (Experimental.) For systems described in 2D (electrons confined to 2D in semiconductor structures), one
281 !% may add the Bychkov-Rashba spin-orbit coupling term [Bychkov and Rashba, <i>J. Phys. C: Solid
282 !% State Phys.</i> <b>17</b>, 6031 (1984)]. This variable determines the strength
283 !% of this perturbation, and has dimensions of energy times length.
284 !%End
285 call parse_variable(namespace, 'RashbaSpinOrbitCoupling', m_zero, rashba_coupling, units_inp%energy*units_inp%length)
286 if (parse_is_defined(namespace, 'RashbaSpinOrbitCoupling')) then
287 if (space%dim /= 2) then
288 write(message(1),'(a)') 'Rashba spin-orbit coupling can only be used for two-dimensional systems.'
289 call messages_fatal(1, namespace=namespace)
290 end if
291 call messages_experimental('RashbaSpinOrbitCoupling', namespace=namespace)
292 end if
294 call hm%hm_base%init(hm%d%nspin, hm%mass, rashba_coupling)
295 call hm%vnl%init()
297 assert(associated(gr%der%lapl))
298 hm%hm_base%kinetic => gr%der%lapl
299
300 safe_allocate(hm%energy)
302 !Keep pointers to derivatives, geometry and xc
303 hm%der => gr%der
304 hm%ions => ions
305 hm%xc => xc
306
307 if(present(xc_photons)) then
308 hm%xc_photons => xc_photons
309 else
310 hm%xc_photons => null()
311 end if
313 ! allocate potentials and density of the cores
314 ! In the case of spinors, vxc_11 = hm%vxc(:, 1), vxc_22 = hm%vxc(:, 2), Re(vxc_12) = hm%vxc(:. 3);
315 ! Im(vxc_12) = hm%vxc(:, 4)
316 call hm%ks_pot%init(gr%der, gr%np, gr%np_part, hm%d%nspin, hm%theory_level, family_is_mgga_with_exc(hm%xc))
318 !Initialize Poisson solvers
319 call poisson_init(hm%psolver, namespace, space, gr%der, mc, gr%stencil, st%qtot)
321 ! Initialize external potential
322 call epot_init(hm%ep, namespace, gr, hm%ions, hm%psolver, hm%d%ispin, hm%xc%family, hm%kpoints)
323 call kick_init(hm%kick, namespace, space, hm%kpoints, hm%d%ispin)
324
325 hm%zora => zora_t(namespace, hm%der, hm%d, hm%ep, hm%mass)
326
327 !Temporary construction of the ion-electron interactions
328 call hm%v_ie_loc%init(gr, hm%psolver, hm%ions, namespace)
329 if (hm%ep%nlcc) then
330 call hm%nlcc%init(gr, hm%ions)
331 safe_allocate(st%rho_core(1:gr%np))
332 st%rho_core(:) = m_zero
333 end if
334
335 !Static magnetic field or rashba spin-orbit interaction requires complex wavefunctions
336 if (parse_is_defined(namespace, 'StaticMagneticField') .or. list_has_gauge_field(ext_partners) .or. &
337 parse_is_defined(namespace, 'RashbaSpinOrbitCoupling')) then
338 call states_set_complex(st)
339 end if
340
341 !%Variable CalculateSelfInducedMagneticField
342 !%Type logical
343 !%Default no
344 !%Section Hamiltonian
345 !%Description
346 !% The existence of an electronic current implies the creation of a self-induced magnetic
347 !% field, which may in turn back-react on the system. Of course, a fully consistent treatment
348 !% of this kind of effect should be done in QED theory, but we will attempt a first
349 !% approximation to the problem by considering the lowest-order relativistic terms
350 !% plugged into the normal Hamiltonian equations (spin-other-orbit coupling terms, etc.).
351 !% For the moment being, none of this is done, but a first step is taken by calculating
352 !% the induced magnetic field of a system that has a current, by considering the magnetostatic
353 !% approximation and Biot-Savart law:
354 !%
355 !% <math> \nabla^2 \vec{A} + 4\pi\alpha \vec{J} = 0</math>
356 !%
357 !% <math> \vec{B} = \vec{\nabla} \times \vec{A}</math>
358 !%
359 !% If <tt>CalculateSelfInducedMagneticField</tt> is set to yes, this <i>B</i> field is
360 !% calculated at the end of a <tt>gs</tt> calculation (nothing is done -- yet -- in the <tt>td</tt>case)
361 !% and printed out, if the <tt>Output</tt> variable contains the <tt>potential</tt> keyword (the prefix
362 !% of the output files is <tt>Bind</tt>).
363 !%End
364 call parse_variable(namespace, 'CalculateSelfInducedMagneticField', .false., hm%self_induced_magnetic)
365 if (hm%self_induced_magnetic) then
366 safe_allocate(hm%a_ind(1:gr%np_part, 1:space%dim))
367 safe_allocate(hm%b_ind(1:gr%np_part, 1:space%dim))
368
369 !(for dim = we could save some memory, but it is better to keep it simple)
370 end if
371
372 ! Absorbing boundaries
373 call absorbing_boundaries_init(hm%abs_boundaries, namespace, space, gr)
374
375 hm%inh_term = .false.
376 call oct_exchange_remove(hm%oct_exchange)
377
378 hm%adjoint = .false.
379
380 call hm%phase%init(gr, hm%d%kpt, hm%kpoints, st%d, space)
381
382 !%Variable DFTULevel
383 !%Type integer
384 !%Default no
385 !%Section Hamiltonian::XC
386 !%Description
387 !% This variable selects which DFT+U expression is added to the Hamiltonian.
388 !%Option dft_u_none 0
389 !% No +U term is not applied.
390 !%Option dft_u_empirical 1
391 !% An empiricial Hubbard U is added on the orbitals specified in the block species
392 !% with hubbard_l and hubbard_u
393 !%Option dft_u_acbn0 2
394 !% Octopus determines the effective U term using the
395 !% ACBN0 functional as defined in PRX 5, 011006 (2015)
396 !%End
397 call parse_variable(namespace, 'DFTULevel', dft_u_none, hm%lda_u_level)
398 call messages_print_var_option('DFTULevel', hm%lda_u_level, namespace=namespace)
399 if (hm%lda_u_level /= dft_u_none) then
400 call lda_u_init(hm%lda_u, namespace, space, hm%lda_u_level, gr, ions, st, mc, &
401 hm%kpoints, hm%phase%is_allocated())
402
403 !In the present implementation of DFT+U, in case of spinors, we have off-diagonal terms
404 !in spin space which break the assumption of the generalized Bloch theorem
405 if (kick_get_type(hm%kick) == kick_magnon_mode .and. gr%der%boundaries%spiral) then
406 call messages_not_implemented("DFT+U with generalized Bloch theorem and magnon kick", namespace=namespace)
407 end if
408
409 ! We rebuild the phase for the orbital projection, similarly to the one of the pseudopotentials
410 if(hm%lda_u_level /= dft_u_none .and. hm%phase%is_allocated()) then
411 call lda_u_build_phase_correction(hm%lda_u, space, hm%d, gr%der%boundaries, namespace, hm%kpoints)
412 end if
413 end if
414
415 !%Variable HamiltonianApplyPacked
416 !%Type logical
417 !%Default yes
418 !%Section Execution::Optimization
419 !%Description
420 !% If set to yes (the default), Octopus will 'pack' the
421 !% wave-functions when operating with them. This might involve some
422 !% additional copying but makes operations more efficient.
423 !% See also the related <tt>StatesPack</tt> variable.
424 !%End
425 call parse_variable(namespace, 'HamiltonianApplyPacked', .true., hm%is_applied_packed)
426
427 if (hm%theory_level == hartree_fock .and. st%parallel_in_states) then
428 call messages_experimental('Hartree-Fock parallel in states', namespace=namespace)
429 end if
430
431 if (hm%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(hm%xc) &
432 .and. st%parallel_in_states) then
433 call messages_experimental('Hybrid functionals parallel in states', namespace=namespace)
434 end if
435
436 !%Variable TimeZero
437 !%Type logical
438 !%Default no
439 !%Section Hamiltonian
440 !%Description
441 !% (Experimental) If set to yes, the ground state and other time
442 !% dependent calculation will assume that they are done at time
443 !% zero, so that all time depedent field at that time will be
444 !% included.
445 !%End
446 call parse_variable(namespace, 'TimeZero', .false., hm%time_zero)
447 if (hm%time_zero) call messages_experimental('TimeZero', namespace=namespace)
448
449 !Cam parameters are irrelevant here and are updated later
450 need_exchange_ = optional_default(need_exchange, .false.)
451 if (hm%xc%compute_exchange(hm%theory_level) .or. need_exchange_) then
452 !We test Slater before OEP, as Slater is treated as OEP for the moment....
453 if (hm%xc%functional(func_x,1)%id == xc_oep_x_slater) then
454 call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
455 hm%kpoints, cam_exact_exchange)
456 else if (bitand(hm%xc%family, xc_family_oep) /= 0 .or. hm%theory_level == rdmft) then
457 call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
458 hm%kpoints, hm%xc%cam)
459 if (hm%theory_level == rdmft) hm%exxop%useACE = .false.
460 else
461 call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
462 hm%kpoints, cam_exact_exchange)
463 end if
464 end if
465
466 if (hm%is_applied_packed .and. accel_is_enabled()) then
467 ! Check if we can actually apply the hamiltonian packed
468 if (gr%use_curvilinear) then
469 if (accel_allow_cpu_only()) then
470 hm%is_applied_packed = .false.
471 call messages_write('Cannot use CUDA or OpenCL as curvilinear coordinates are used.')
472 call messages_warning(namespace=namespace)
473 else
474 call messages_write('Cannot use CUDA or OpenCL as curvilinear coordinates are used.', new_line = .true.)
475 call messages_write('Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
476 call messages_fatal(namespace=namespace)
477 end if
478 end if
479 end if
480
481 !We are building the list of external potentials
482 !This is done here at the moment, because we pass directly the mesh
483 !TODO: Once the abstract Hamiltonian knows about an abstract basis, we might move this to the
484 ! abstract Hamiltonian
485 call load_external_potentials(hm%external_potentials, namespace)
486
487 !Some checks which are electron specific, like k-points
489
490 !At the moment we do only have static external potential, so we never update them
492
493 !Build the resulting interactions
494 !TODO: This will be moved to the actual interactions
495 call build_interactions()
496
497 ! Constrained DFT for noncollinear magnetism
498 if (hm%theory_level /= independent_particles) then
499 call magnetic_constrain_init(hm%magnetic_constrain, namespace, gr, st%d, ions%natoms, ions%min_distance())
500 end if
501
502 ! init maxwell-electrons coupling
503 call mxll_coupling_init(hm%mxll, st%d, gr, namespace, hm%mass)
504
505 if (associated(hm%xc_photons)) then
506 if (hm%xc_photons%wants_to_renormalize_mass()) then
507 ! remornalize the electron mass due to light-matter interaction; here we only deal with it in free space
508 call hm%set_mass(namespace, hm%xc_photons%get_renormalized_mass())
509 end if
510 end if
511
512 if (hm%xc%compute_exchange(hm%theory_level) .or. need_exchange_) call hm%exxop%write_info(namespace)
513
514 call profiling_out('HAMILTONIAN_ELEC_INIT')
515 pop_sub(hamiltonian_elec_init)
516
517 contains
518
519 ! ---------------------------------------------------------
520 subroutine build_external_potentials()
521 type(list_iterator_t) :: iter
522 class(*), pointer :: potential
523 integer :: iop
524
526
527 safe_allocate(hm%v_ext_pot(1:gr%np))
528 hm%v_ext_pot(1:gr%np) = m_zero
529
530 call iter%start(hm%external_potentials)
531 do while (iter%has_next())
532 potential => iter%get_next()
533 select type (potential)
534 class is (external_potential_t)
535
536 call potential%allocate_memory(gr)
537 call potential%calculate(namespace, gr, hm%psolver)
538 !To preserve the old behavior, we are adding the various potentials
539 !to the corresponding arrays
540 select case (potential%type)
542 call lalg_axpy(gr%np, m_one, potential%pot, hm%v_ext_pot)
543
545 if (states_are_real(st)) then
546 message(1) = "Cannot use static magnetic field with real wavefunctions"
547 call messages_fatal(1, namespace=namespace)
548 end if
549
550 if (.not. allocated(hm%ep%b_field)) then
551 safe_allocate(hm%ep%b_field(1:3)) !Cannot be space%dim
552 hm%ep%b_field(1:3) = m_zero
553 end if
554 hm%ep%b_field(1:3) = hm%ep%b_field(1:3) + potential%b_field(1:3)
555
556 if (.not. allocated(hm%ep%a_static)) then
557 safe_allocate(hm%ep%a_static(1:gr%np, 1:space%dim))
558 hm%ep%a_static(1:gr%np, 1:space%dim) = m_zero
559 end if
560 call lalg_axpy(gr%np, space%dim, m_one, potential%a_static, hm%ep%a_static)
561
563 if (.not. allocated(hm%ep%e_field)) then
564 safe_allocate(hm%ep%e_field(1:space%dim))
565 hm%ep%e_field(1:space%dim) = m_zero
566 end if
567 hm%ep%e_field(1:space%dim) = hm%ep%e_field(1:space%dim) + potential%e_field(1:space%dim)
568
569 !In the fully periodic case, we use Berry phases
570 if (space%periodic_dim < space%dim) then
571 if (.not. allocated(hm%v_static)) then
572 safe_allocate(hm%v_static(1:gr%np))
573 hm%v_static(1:gr%np) = m_zero
574 end if
575 if (.not. allocated(hm%ep%v_ext)) then
576 safe_allocate(hm%ep%v_ext(1:gr%np_part))
577 hm%ep%v_ext(1:gr%np_part) = m_zero
578 end if
579 call lalg_axpy(gr%np, m_one, potential%pot, hm%v_static)
580 call lalg_axpy(gr%np, m_one, potential%v_ext, hm%ep%v_ext)
581 end if
582
583 if (hm%kpoints%use_symmetries) then
584 do iop = 1, symmetries_number(hm%kpoints%symm)
585 if (iop == symmetries_identity_index(hm%kpoints%symm)) cycle
586 if (.not. symm_op_invariant_cart(hm%kpoints%symm%ops(iop), hm%ep%e_field, 1e-5_real64)) then
587 message(1) = "The StaticElectricField breaks (at least) one of the symmetries used to reduce the k-points."
588 message(2) = "Set SymmetryBreakDir equal to StaticElectricField."
589 call messages_fatal(2, namespace=namespace)
590 end if
591 end do
592 end if
593
594 end select
595 call potential%deallocate_memory()
596
597 class default
598 assert(.false.)
599 end select
600 end do
601
603 end subroutine build_external_potentials
604
605 ! ---------------------------------------------------------
606 subroutine external_potentials_checks()
607 type(list_iterator_t) :: iter
608 class(*), pointer :: potential
609
611
612 call iter%start(hm%external_potentials)
613 do while (iter%has_next())
614 potential => iter%get_next()
615 select type (potential)
616 class is (external_potential_t)
617
618 if (potential%type == external_pot_static_efield .and. hm%kpoints%reduced%npoints > 1) then
619 message(1) = "Applying StaticElectricField in a periodic direction is only accurate for large supercells."
620 message(2) = "Single-point Berry phase is not appropriate when k-point sampling is needed."
621 call messages_warning(2, namespace=namespace)
622 end if
623
624 class default
625 assert(.false.)
626 end select
627 end do
628
630 end subroutine external_potentials_checks
631
632
633 !The code in this routines needs to know about the external potentials.
634 !This will be treated in the future by the interactions directly.
635 subroutine build_interactions()
636 logical :: external_potentials_present
637 logical :: kick_present
638
640
641 if (allocated(hm%ep%e_field) .and. space%is_periodic() .and. .not. list_has_gauge_field(ext_partners)) then
642 ! only need vberry if there is a field in a periodic direction
643 ! and we are not setting a gauge field
644 if (any(abs(hm%ep%e_field(1:space%periodic_dim)) > m_epsilon)) then
645 safe_allocate(hm%vberry(1:gr%np, 1:hm%d%nspin))
646 hm%vberry = m_zero
647 end if
648 end if
649
650 external_potentials_present = epot_have_external_potentials(hm%ep) .or. &
651 list_has_lasers(ext_partners) .or. allocated(hm%v_static)
652
653
654 kick_present = hamiltonian_elec_has_kick(hm)
655
656 call pcm_init(hm%pcm, namespace, space, ions, gr, st%qtot, st%val_charge, external_potentials_present, kick_present)
657 if (hm%pcm%run_pcm) then
658 if (hm%theory_level /= kohn_sham_dft) call messages_not_implemented("PCM for TheoryLevel /= kohn_sham", namespace=namespace)
659 end if
660
662
663 end subroutine build_interactions
664
665
666 end subroutine hamiltonian_elec_init
667
668
669
670
671 ! ---------------------------------------------------------
672 subroutine hamiltonian_elec_end(hm)
673 type(hamiltonian_elec_t), target, intent(inout) :: hm
674
675 type(partner_iterator_t) :: iter
676 class(interaction_partner_t), pointer :: potential
677
678 push_sub(hamiltonian_elec_end)
679
680 call hm%hm_base%end()
681 call hm%vnl%end()
682
683 call hm%phase%end()
684
685 call hm%ks_pot%end()
686 safe_deallocate_a(hm%vberry)
687 safe_deallocate_a(hm%a_ind)
688 safe_deallocate_a(hm%b_ind)
689 safe_deallocate_a(hm%v_ext_pot)
690
691 safe_deallocate_p(hm%zora)
692
693 call poisson_end(hm%psolver)
694
695 nullify(hm%xc)
696
697 call kick_end(hm%kick)
698 call epot_end(hm%ep)
699 nullify(hm%ions)
700
701 call absorbing_boundaries_end(hm%abs_boundaries)
702
703 call states_elec_dim_end(hm%d)
704
705 if (hm%scissor%apply) call scissor_end(hm%scissor)
706
707 call exchange_operator_end(hm%exxop)
708 call lda_u_end(hm%lda_u)
709
710 safe_deallocate_a(hm%energy)
711
712 if (hm%pcm%run_pcm) call pcm_end(hm%pcm)
713
714 call hm%v_ie_loc%end()
715 call hm%nlcc%end()
716
717 call iter%start(hm%external_potentials)
718 do while (iter%has_next())
719 potential => iter%get_next()
720 safe_deallocate_p(potential)
721 end do
722 call hm%external_potentials%empty()
723 safe_deallocate_a(hm%v_static)
724
725 call magnetic_constrain_end(hm%magnetic_constrain)
726
727 call mxll_coupling_end(hm%mxll)
728
729 pop_sub(hamiltonian_elec_end)
730 end subroutine hamiltonian_elec_end
731
732
733 ! ---------------------------------------------------------
734 ! True if the Hamiltonian is Hermitian, false otherwise
735 logical function hamiltonian_elec_hermitian(hm)
736 class(hamiltonian_elec_t), intent(in) :: hm
737
739 hamiltonian_elec_hermitian = .not.((hm%abs_boundaries%abtype == imaginary_absorbing) .or. &
740 oct_exchange_enabled(hm%oct_exchange))
741
743 end function hamiltonian_elec_hermitian
744
745 ! ---------------------------------------------------------
746 ! True if the Hamiltonian needs to apply the mGGA term
747 pure logical function hamiltonian_elec_needs_mgga_term(hm, terms)
748 class(hamiltonian_elec_t), intent(in) :: hm
749 integer, intent(in) :: terms
750
752
753 if (bitand(term_mgga, terms) /= 0 .and. family_is_mgga_with_exc(hm%xc) &
754 .and. hm%theory_level == generalized_kohn_sham_dft) then
756 end if
757
758 !For OEP
759 if(terms == term_mgga .and. family_is_mgga_with_exc(hm%xc) .and. hm%theory_level == kohn_sham_dft) then
761 end if
762
764
765
766
767 ! ---------------------------------------------------------
768 subroutine hamiltonian_elec_span(hm, delta, emin, namespace)
769 class(hamiltonian_elec_t), intent(inout) :: hm
770 real(real64), intent(in) :: delta(:)
771 real(real64), intent(in) :: emin
772 type(namespace_t), intent(in) :: namespace
773
774 real(real64) :: emax
775
776 push_sub(hamiltonian_elec_span)
777
778 ! estimate maximum energy of discrete kinetic operator
779 ! this neglects possible contributions from the non-local part of the pseudopotentials
781
782 hm%spectral_middle_point = (emax + emin) / m_two
783 hm%spectral_half_span = (emax - emin) / m_two
784
785 pop_sub(hamiltonian_elec_span)
786 end subroutine hamiltonian_elec_span
787
788
789 ! ---------------------------------------------------------
790 pure logical function hamiltonian_elec_inh_term(hm) result(inh)
791 type(hamiltonian_elec_t), intent(in) :: hm
792
793 inh = hm%inh_term
794 end function hamiltonian_elec_inh_term
795
796
797 ! ---------------------------------------------------------
798 subroutine hamiltonian_elec_set_inh(hm, st)
799 type(hamiltonian_elec_t), intent(inout) :: hm
800 type(states_elec_t), intent(in) :: st
801
803
804 if (hm%inh_term) call states_elec_end(hm%inh_st)
805 call states_elec_copy(hm%inh_st, st)
806 hm%inh_term = .true.
807
809 end subroutine hamiltonian_elec_set_inh
810
811
812 ! ---------------------------------------------------------
813 subroutine hamiltonian_elec_remove_inh(hm)
814 type(hamiltonian_elec_t), intent(inout) :: hm
815
817
818 if (hm%inh_term) then
819 call states_elec_end(hm%inh_st)
820 hm%inh_term = .false.
821 end if
822
824 end subroutine hamiltonian_elec_remove_inh
825
826 ! ---------------------------------------------------------
827 subroutine hamiltonian_elec_adjoint(hm)
828 type(hamiltonian_elec_t), intent(inout) :: hm
829
831
832 if (.not. hm%adjoint) then
833 hm%adjoint = .true.
834 if (hm%abs_boundaries%abtype == imaginary_absorbing) then
835 hm%abs_boundaries%mf = -hm%abs_boundaries%mf
836 end if
837 end if
838
840 end subroutine hamiltonian_elec_adjoint
841
843 ! ---------------------------------------------------------
844 subroutine hamiltonian_elec_not_adjoint(hm)
845 type(hamiltonian_elec_t), intent(inout) :: hm
846
848
849 if (hm%adjoint) then
850 hm%adjoint = .false.
851 if (hm%abs_boundaries%abtype == imaginary_absorbing) then
852 hm%abs_boundaries%mf = -hm%abs_boundaries%mf
853 end if
854 end if
855
857 end subroutine hamiltonian_elec_not_adjoint
858
859
860 ! ---------------------------------------------------------
862 subroutine hamiltonian_elec_update(this, mesh, namespace, space, ext_partners, time)
863 class(hamiltonian_elec_t), intent(inout) :: this
864 class(mesh_t), intent(in) :: mesh
865 type(namespace_t), intent(in) :: namespace
866 class(space_t), intent(in) :: space
867 type(partner_list_t), intent(in) :: ext_partners
868 real(real64), optional, intent(in) :: time
869
870 integer :: ispin, ip, idir, iatom, ilaser
871 real(real64) :: aa(space%dim), time_
872 real(real64), allocatable :: vp(:,:)
873 type(lasers_t), pointer :: lasers
874 type(gauge_field_t), pointer :: gfield
875 real(real64) :: am(space%dim)
876
878 call profiling_in("HAMILTONIAN_ELEC_UPDATE")
879
880 this%current_time = m_zero
881 if (present(time)) this%current_time = time
882
883 time_ = optional_default(time, 0.0_real64)
884
885 ! set everything to zero
886 call this%hm_base%clear(mesh%np)
887
888 ! alllocate the scalar potential for the xc, hartree and external potentials
889 call this%hm_base%allocate_field(mesh, field_potential, &
890 complex_potential = this%abs_boundaries%abtype == imaginary_absorbing)
891
892 ! the lasers
893 if (present(time) .or. this%time_zero) then
894
895 lasers => list_get_lasers(ext_partners)
896 if(associated(lasers)) then
897 do ilaser = 1, lasers%no_lasers
898 select case (laser_kind(lasers%lasers(ilaser)))
900 do ispin = 1, this%d%spin_channels
901 call laser_potential(lasers%lasers(ilaser), mesh, &
902 this%hm_base%potential(:, ispin), time_)
903 end do
904 case (e_field_magnetic)
905 call this%hm_base%allocate_field(mesh, field_vector_potential + field_uniform_magnetic_field, &
906 .false.)
907 ! get the vector potential
908 safe_allocate(vp(1:mesh%np, 1:space%dim))
909 vp(1:mesh%np, 1:space%dim) = m_zero
910 call laser_vector_potential(lasers%lasers(ilaser), mesh, vp, time_)
911 !$omp parallel do private(idir) schedule(static)
912 do ip = 1, mesh%np
913 do idir = 1, space%dim
914 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + vp(ip, idir)/p_c
915 end do
916 end do
917 ! and the magnetic field
918 call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
919 safe_deallocate_a(vp)
921 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
922 ! get the uniform vector potential associated with a magnetic field
923 aa = m_zero
924 call laser_field(lasers%lasers(ilaser), aa, time_)
925 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
926 end select
927 end do
928
929 if (lasers_with_nondipole_field(lasers)) then
930 assert( allocated(this%hm_base%uniform_vector_potential))
931 call lasers_nondipole_laser_field_step(lasers, am, time_)
932 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - am/p_c
933 end if
934 end if
935
936 ! the gauge field
937 gfield => list_get_gauge_field(ext_partners)
938 if (associated(gfield)) then
939 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
940 call gauge_field_get_vec_pot(gfield, aa)
941 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
942 end if
943
944 ! the electric field for a periodic system through the gauge field
945 if (allocated(this%ep%e_field) .and. associated(gfield)) then
946 this%hm_base%uniform_vector_potential(1:space%periodic_dim) = &
947 this%hm_base%uniform_vector_potential(1:space%periodic_dim) - time_*this%ep%e_field(1:space%periodic_dim)
948 end if
949
950 ! add the photon-free mean-field vector potential
951 if (associated(this%xc_photons)) then
952 if(this%xc_photons%lpfmf .and. allocated(this%xc_photons%mf_vector_potential)) then
953 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
954 ! here we put a minus sign in front of the mean field term to get the right answer (need to check the formula)
955 this%hm_base%uniform_vector_potential(1:space%dim) = &
956 this%hm_base%uniform_vector_potential(1:space%dim) - this%xc_photons%mf_vector_potential(1:space%dim)/p_c
957 end if
958 end if
959
960 end if
961
962 ! the vector potential of a static magnetic field
963 if (allocated(this%ep%a_static)) then
964 call this%hm_base%allocate_field(mesh, field_vector_potential, .false.)
965 !ep%a_static contains 1/c A(r)
966 !$omp parallel do private(idir) schedule(static)
967 do ip = 1, mesh%np
968 do idir = 1, space%dim
969 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + this%ep%a_static(ip, idir)
970 end do
971 end do
972 end if
973
974 ! add Maxwell coupling to Hamiltonian and sets the magnetic field for the Zeeman term added below
975 call mxll_coupling_calc(this%mxll, this%hm_base, mesh, this%d, space)
976
977 !The electric field was added to the KS potential
978 call this%hm_base%accel_copy_pot(mesh)
979
980 ! and the static magnetic field
981 if (allocated(this%ep%b_field)) then
982 call this%hm_base%allocate_field(mesh, field_uniform_magnetic_field, .false.)
983 do idir = 1, 3
984 this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
985 end do
986 end if
987
988 ! Combine the uniform and non-uniform fields and compute the Zeeman term
989 call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
990
991 ! This needs to be called at the end as the zeeman term enters the potential
992 call hamiltonian_elec_update_pot(this, mesh, accumulate = .true.)
993
994 if (this%mxll%test_equad) then
995 call set_electric_quadrupole_pot(this%mxll, mesh)
996 end if
997
998 call build_phase()
999
1000 call profiling_out("HAMILTONIAN_ELEC_UPDATE")
1002
1003 contains
1004
1005 subroutine build_phase()
1006 integer :: ik, imat, nmat, max_npoints, offset
1007 integer :: ip
1008 integer :: iphase, nphase
1009
1011
1012 if ((.not. this%kpoints%gamma_only()) .or. allocated(this%hm_base%uniform_vector_potential)) then
1013
1014 call profiling_in('UPDATE_PHASES')
1015 ! now regenerate the phases for the pseudopotentials
1016 do iatom = 1, this%ep%natoms
1017 call projector_init_phases(this%ep%proj(iatom), space%dim, this%d, this%der%boundaries, this%kpoints, &
1018 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1019 end do
1020
1021 call profiling_out('UPDATE_PHASES')
1022 end if
1023
1024 if (allocated(this%hm_base%uniform_vector_potential)) then
1025
1026 call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
1027
1028 ! We rebuild the phase for the orbital projection, similarly to the one of the pseudopotentials
1029 if (this%lda_u_level /= dft_u_none) then
1030 call lda_u_build_phase_correction(this%lda_u, space, this%d, this%der%boundaries, namespace, this%kpoints, &
1031 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1032 end if
1033 end if
1034
1035 max_npoints = this%vnl%max_npoints
1036 nmat = this%vnl%nprojector_matrices
1037
1038
1039 if (this%phase%is_allocated() .and. allocated(this%vnl%projector_matrices)) then
1040
1041 nphase = 1
1042 if (this%der%boundaries%spiralBC) nphase = 3
1043
1044 if (.not. allocated(this%vnl%projector_phases)) then
1045 safe_allocate(this%vnl%projector_phases(1:max_npoints, 1:nphase, nmat, this%d%kpt%start:this%d%kpt%end))
1046 if (accel_is_enabled()) then
1047 call accel_create_buffer(this%vnl%buff_projector_phases, accel_mem_read_only, &
1048 type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
1049 ! We need to save nphase, with which the array has been build,
1050 ! as the number might change throughout the run
1051 this%vnl%nphase = nphase
1052 end if
1053 end if
1054
1055 offset = 0
1056 do ik = this%d%kpt%start, this%d%kpt%end
1057 do imat = 1, this%vnl%nprojector_matrices
1058 iatom = this%vnl%projector_to_atom(imat)
1059 do iphase = 1, nphase
1060 !$omp parallel do simd schedule(static)
1061 do ip = 1, this%vnl%projector_matrices(imat)%npoints
1062 this%vnl%projector_phases(ip, iphase, imat, ik) = this%ep%proj(iatom)%phase(ip, iphase, ik)
1063 end do
1064
1065 if (accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0) then
1066 call accel_write_buffer(this%vnl%buff_projector_phases, &
1067 this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
1068 offset = offset, async=.true.)
1069 end if
1070 offset = offset + this%vnl%projector_matrices(imat)%npoints
1071 end do
1072 end do
1073 end do
1074
1075 end if
1076
1077 call accel_finish()
1078
1080 end subroutine build_phase
1081
1082 end subroutine hamiltonian_elec_update
1083
1084
1085 !----------------------------------------------------------------
1088 ! TODO: See Issue #1064
1089 subroutine hamiltonian_elec_update_pot(this, mesh, accumulate)
1090 type(hamiltonian_elec_t), intent(inout) :: this
1091 class(mesh_t), intent(in) :: mesh
1092 logical, optional, intent(in) :: accumulate
1093
1094 integer :: ispin, ip
1095
1097
1098 ! By default we nullify first the result
1099 if (.not. optional_default(accumulate, .false.)) then
1100 !$omp parallel private(ip, ispin)
1101 do ispin = 1, this%d%nspin
1102 !$omp do simd schedule(static)
1103 do ip = 1, mesh%np
1104 this%hm_base%potential(ip, ispin) = m_zero
1105 end do
1106 end do
1107 !$omp end parallel
1108 end if
1109
1110 !$omp parallel private(ip, ispin)
1111 do ispin = 1, this%d%nspin
1112 if (ispin <= 2) then
1113 !$omp do simd schedule(static)
1114 ! this%vhxc(ip, ispin) is added after the calculation of ZORA potential
1115 do ip = 1, mesh%np
1116 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%ep%vpsl(ip) + this%v_ext_pot(ip)
1117 end do
1118
1120 if (this%pcm%run_pcm) then
1121 if (this%pcm%solute) then
1122 !$omp do simd schedule(static)
1123 do ip = 1, mesh%np
1124 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
1125 this%pcm%v_e_rs(ip) + this%pcm%v_n_rs(ip)
1126 end do
1127 !$omp end do simd nowait
1128 end if
1129 if (this%pcm%localf) then
1130 !$omp do simd schedule(static)
1131 do ip = 1, mesh%np
1132 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
1133 this%pcm%v_ext_rs(ip)
1134 end do
1135 !$omp end do simd nowait
1136 end if
1137 end if
1138
1140 if (this%abs_boundaries%abtype == imaginary_absorbing) then
1141 !$omp do simd schedule(static)
1142 do ip = 1, mesh%np
1143 this%hm_base%Impotential(ip, ispin) = this%hm_base%Impotential(ip, ispin) + this%abs_boundaries%mf(ip)
1144 end do
1145 !$omp end do simd nowait
1146 end if
1147 end if
1148 end do
1149 !$omp end parallel
1150
1151 ! scalar relativistic ZORA contribution
1152 ! \boldsymbol{p} \frac{c^2}{2c^2 - V} \boldsymbol{p} \Phi^\mathrm{ZORA}
1153 if (this%ep%reltype == scalar_relativistic_zora .or. this%ep%reltype == fully_relativistic_zora) then
1154 call this%zora%update(this%der, this%hm_base%potential)
1155 end if
1156
1157 !$omp parallel private(ip, ispin)
1158 do ispin = 1, this%d%nspin
1159 ! Adding Zeeman potential to hm_base%potential
1160 if (allocated(this%hm_base%zeeman_pot)) then
1161 !$omp do simd schedule(static)
1162 do ip = 1, mesh%np
1163 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%hm_base%zeeman_pot(ip, ispin)
1164 end do
1165 !$omp end do simd nowait
1166 end if
1167
1168 ! Adding Quadrupole potential from static E-field (test)
1169 if (this%mxll%test_equad) then
1170 !$omp do simd schedule(static)
1171 do ip = 1, mesh%np
1172 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%mxll%e_quadrupole_pot(ip)
1173 end do
1174 !$omp end do simd
1175 end if
1176 end do
1177 !$omp end parallel
1178
1179
1180 ! Add the Hartree and KS potential
1181 call this%ks_pot%add_vhxc(this%hm_base%potential)
1182
1183 call this%hm_base%accel_copy_pot(mesh)
1186 end subroutine hamiltonian_elec_update_pot
1187
1188 ! ---------------------------------------------------------
1189 subroutine hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
1190 type(hamiltonian_elec_t), intent(inout) :: this
1191 type(namespace_t), intent(in) :: namespace
1192 class(electron_space_t), intent(in) :: space
1193 type(grid_t), intent(in) :: gr
1194 type(ions_t), target, intent(inout) :: ions
1195 type(partner_list_t), intent(in) :: ext_partners
1196 type(states_elec_t), intent(inout) :: st
1197 real(real64), optional, intent(in) :: time
1198
1200
1201 this%ions => ions
1202 call epot_generate(this%ep, namespace, gr, this%ions, this%d)
1203
1204 ! Interation terms are treated below
1205
1206 ! First we add the static electric field
1207 if (allocated(this%ep%e_field) .and. space%periodic_dim < space%dim) then
1208 call lalg_axpy(gr%np, m_one, this%v_static, this%ep%vpsl)
1209 end if
1210
1211 ! Here we need to pass this again, else test are failing.
1212 ! This is not a real problem, as the multisystem framework will indeed to this anyway
1213 this%v_ie_loc%atoms_dist => ions%atoms_dist
1214 this%v_ie_loc%atom => ions%atom
1215 call this%v_ie_loc%calculate()
1216
1217 ! At the moment we need to add this to ep%vpsl, to keep the behavior of the code
1218 call lalg_axpy(gr%np, m_one, this%v_ie_loc%potential(:,1), this%ep%vpsl)
1219
1220 ! Here we need to pass this again, else test are failing.
1221 ! This is not a real problem, as the multisystem framework will indeed to this anyway
1222 if (this%ep%nlcc) then
1223 this%nlcc%atoms_dist => ions%atoms_dist
1224 this%nlcc%atom => ions%atom
1225 call this%nlcc%calculate()
1226 call lalg_copy(gr%np, this%nlcc%density(:,1), st%rho_core)
1227 end if
1228
1229 call this%vnl%build(space, gr, this%ep)
1230 call hamiltonian_elec_update(this, gr, namespace, space, ext_partners, time)
1231
1232 ! Check if projectors are still compatible with apply_packed on GPU
1233 if (this%is_applied_packed .and. accel_is_enabled()) then
1234 if (this%ep%non_local .and. .not. this%vnl%apply_projector_matrices) then
1235 if (accel_allow_cpu_only()) then
1236 call messages_write('Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.')
1237 call messages_warning(namespace=namespace)
1238 else
1239 call messages_write('Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.',&
1240 new_line = .true.)
1241 call messages_write('Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
1242 call messages_fatal(namespace=namespace)
1243 end if
1244 end if
1245
1246 end if
1247
1248 if (this%pcm%run_pcm) then
1251 if (this%pcm%solute) then
1252 call pcm_calc_pot_rs(this%pcm, gr, this%psolver, ions = ions)
1253 end if
1254
1257 ! Interpolation is needed, hence gr%np_part -> 1:gr%np
1258 if (this%pcm%localf .and. allocated(this%v_static)) then
1259 call pcm_calc_pot_rs(this%pcm, gr, this%psolver, v_ext = this%ep%v_ext(1:gr%np_part))
1260 end if
1261
1262 end if
1263
1264 call lda_u_update_basis(this%lda_u, space, gr, ions, st, this%psolver, namespace, this%kpoints, &
1265 this%phase%is_allocated())
1266
1268 end subroutine hamiltonian_elec_epot_generate
1269
1270 ! -----------------------------------------------------------------
1271
1272 real(real64) function hamiltonian_elec_get_time(this) result(time)
1273 type(hamiltonian_elec_t), intent(inout) :: this
1274
1275 time = this%current_time
1276 end function hamiltonian_elec_get_time
1277
1278 ! -----------------------------------------------------------------
1279
1280 pure logical function hamiltonian_elec_apply_packed(this) result(apply)
1281 class(hamiltonian_elec_t), intent(in) :: this
1282
1283 apply = this%is_applied_packed
1286
1287
1288 ! -----------------------------------------------------------------
1289 subroutine zhamiltonian_elec_apply_atom (hm, namespace, space, latt, species, pos, ia, mesh, psi, vpsi)
1290 type(hamiltonian_elec_t), intent(in) :: hm
1291 type(namespace_t), intent(in) :: namespace
1292 class(space_t), intent(in) :: space
1293 type(lattice_vectors_t), intent(in) :: latt
1294 class(species_t), intent(in) :: species
1295 real(real64), intent(in) :: pos(1:space%dim)
1296 integer, intent(in) :: ia
1297 class(mesh_t), intent(in) :: mesh
1298 complex(real64), intent(in) :: psi(:,:)
1299 complex(real64), intent(out) :: vpsi(:,:)
1300
1301 integer :: idim
1302 real(real64), allocatable :: vlocal(:)
1304
1305 safe_allocate(vlocal(1:mesh%np_part))
1306 vlocal = m_zero
1307 call epot_local_potential(hm%ep, namespace, space, latt, mesh, species, pos, ia, vlocal)
1308
1309 do idim = 1, hm%d%dim
1310 vpsi(1:mesh%np, idim) = vlocal(1:mesh%np) * psi(1:mesh%np, idim)
1311 end do
1312
1313 safe_deallocate_a(vlocal)
1315 end subroutine zhamiltonian_elec_apply_atom
1316
1317 ! ---------------------------------------------------------
1322 subroutine hamiltonian_elec_update_with_ext_pot(this, mesh, space, ext_partners, time, mu)
1323 type(hamiltonian_elec_t), intent(inout) :: this
1324 class(space_t), intent(in) :: space
1325 class(mesh_t), intent(in) :: mesh
1326 type(partner_list_t), intent(in) :: ext_partners
1327 real(real64), intent(in) :: time(1:2)
1328 real(real64), intent(in) :: mu(1:2)
1329
1330 integer :: ispin, ip, idir, iatom, ilaser, itime
1331 real(real64) :: aa(space%dim), time_
1332 real(real64), allocatable :: vp(:,:)
1333 real(real64), allocatable :: velectric(:)
1334 type(lasers_t), pointer :: lasers
1335 type(gauge_field_t), pointer :: gfield
1336
1338 call profiling_in("HAMILTONIAN_ELEC_UPDATE_EXT_POT")
1339
1340 this%current_time = m_zero
1341 this%current_time = time(1)
1342
1343 ! set everything to zero
1344 call this%hm_base%clear(mesh%np)
1345
1346 ! the xc, hartree and external potentials
1347 call this%hm_base%allocate_field(mesh, field_potential, &
1348 complex_potential = this%abs_boundaries%abtype == imaginary_absorbing)
1349
1350 do itime = 1, 2
1351 time_ = time(itime)
1352
1353 lasers => list_get_lasers(ext_partners)
1354 if(associated(lasers)) then
1355 do ilaser = 1, lasers%no_lasers
1356 select case (laser_kind(lasers%lasers(ilaser)))
1357 case (e_field_scalar_potential, e_field_electric)
1358 safe_allocate(velectric(1:mesh%np))
1359 do ispin = 1, this%d%spin_channels
1360 velectric = m_zero
1361 call laser_potential(lasers%lasers(ilaser), mesh, velectric, time_)
1362 !$omp parallel do simd schedule(static)
1363 do ip = 1, mesh%np
1364 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + mu(itime) * velectric(ip)
1365 end do
1366 end do
1367 safe_deallocate_a(velectric)
1368 case (e_field_magnetic)
1369 call this%hm_base%allocate_field(mesh, field_vector_potential + field_uniform_magnetic_field, .false.)
1370 ! get the vector potential
1371 safe_allocate(vp(1:mesh%np, 1:space%dim))
1372 vp(1:mesh%np, 1:space%dim) = m_zero
1373 call laser_vector_potential(lasers%lasers(ilaser), mesh, vp, time_)
1374 do idir = 1, space%dim
1375 !$omp parallel do schedule(static)
1376 do ip = 1, mesh%np
1377 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) &
1378 - mu(itime) * vp(ip, idir)/p_c
1379 end do
1380 end do
1381 ! and the magnetic field
1382 call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
1383 safe_deallocate_a(vp)
1384 case (e_field_vector_potential)
1385 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
1386 ! get the uniform vector potential associated with a magnetic field
1387 aa = m_zero
1388 call laser_field(lasers%lasers(ilaser), aa, time_)
1389 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) &
1390 - mu(itime) * aa/p_c
1391 end select
1392 end do
1393 end if
1394
1395 ! the gauge field
1396 gfield => list_get_gauge_field(ext_partners)
1397 if (associated(gfield)) then
1398 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
1399 call gauge_field_get_vec_pot(gfield, aa)
1400 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
1401 end if
1402
1403 ! the electric field for a periodic system through the gauge field
1404 ! TODO: The condition is wrong here: the e_field should be in non-periodic dims as E field
1405 ! and as a gauge field in the periodic dim, unless we use a Bery phase, in which, we do not use it
1406 ! this way. But this is unrelated to the gauge field
1407 if (allocated(this%ep%e_field) .and. associated(gfield)) then
1408 this%hm_base%uniform_vector_potential(1:space%periodic_dim) = &
1409 this%hm_base%uniform_vector_potential(1:space%periodic_dim) - time_*this%ep%e_field(1:space%periodic_dim)
1410 end if
1411
1412 end do
1413
1414 ! the vector potential of a static magnetic field
1415 if (allocated(this%ep%a_static)) then
1416 call this%hm_base%allocate_field(mesh, field_vector_potential, .false.)
1417 !ep%a_static contains 1/c A(r)
1418 !$omp parallel do schedule(static) private(idir)
1419 do ip = 1, mesh%np
1420 do idir = 1, space%dim
1421 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + this%ep%a_static(ip, idir)
1422 end do
1423 end do
1424 end if
1425
1426 !The electric field is added to the KS potential
1427 call this%hm_base%accel_copy_pot(mesh)
1428
1429 ! and the static magnetic field
1430 if (allocated(this%ep%b_field)) then
1431 call this%hm_base%allocate_field(mesh, field_uniform_magnetic_field, .false.)
1432 do idir = 1, 3
1433 this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
1434 end do
1435 end if
1436
1437 call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
1438
1439 call hamiltonian_elec_update_pot(this, mesh)
1440
1441 call build_phase()
1442
1443 call profiling_out("HAMILTONIAN_ELEC_UPDATE_EXT_POT")
1445
1446 contains
1447
1448 subroutine build_phase()
1449 integer :: ik, imat, nmat, max_npoints, offset, iphase, nphase
1450
1452
1453 if ((.not. this%kpoints%gamma_only()) .or. allocated(this%hm_base%uniform_vector_potential)) then
1454
1455 call profiling_in('UPDATE_PHASES')
1456 ! now regenerate the phases for the pseudopotentials
1457 do iatom = 1, this%ep%natoms
1458 call projector_init_phases(this%ep%proj(iatom), space%dim, this%d, this%der%boundaries, this%kpoints, &
1459 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1460 end do
1461
1462 call profiling_out('UPDATE_PHASES')
1463 end if
1464
1465 if (allocated(this%hm_base%uniform_vector_potential)) then
1466 call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
1467 end if
1468
1469 max_npoints = this%vnl%max_npoints
1470 nmat = this%vnl%nprojector_matrices
1471
1472
1473 if (this%phase%is_allocated() .and. allocated(this%vnl%projector_matrices)) then
1474
1475 nphase = 1
1476 if (this%der%boundaries%spiralBC) nphase = 3
1477
1478 if (.not. allocated(this%vnl%projector_phases)) then
1479 safe_allocate(this%vnl%projector_phases(1:max_npoints, nphase, nmat, this%d%kpt%start:this%d%kpt%end))
1480 if (accel_is_enabled()) then
1481 call accel_create_buffer(this%vnl%buff_projector_phases, accel_mem_read_only, &
1482 type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
1483 end if
1484 end if
1485
1486 offset = 0
1487 do ik = this%d%kpt%start, this%d%kpt%end
1488 do imat = 1, this%vnl%nprojector_matrices
1489 iatom = this%vnl%projector_to_atom(imat)
1490 do iphase = 1, nphase
1491 !$omp parallel do schedule(static)
1492 do ip = 1, this%vnl%projector_matrices(imat)%npoints
1493 this%vnl%projector_phases(ip, imat, iphase, ik) = this%ep%proj(iatom)%phase(ip, iphase, ik)
1494 end do
1495
1496 if (accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0) then
1497 call accel_write_buffer(this%vnl%buff_projector_phases, &
1498 this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
1499 offset = offset)
1500 end if
1501 offset = offset + this%vnl%projector_matrices(imat)%npoints
1502 end do
1503 end do
1504 end do
1505
1506 end if
1507
1509 end subroutine build_phase
1510
1512
1513 logical function hamiltonian_elec_needs_current(hm, states_are_real)
1514 type(hamiltonian_elec_t), intent(in) :: hm
1515 logical, intent(in) :: states_are_real
1516
1518
1519 if (hm%self_induced_magnetic) then
1520 if (.not. states_are_real) then
1522 else
1523 message(1) = 'No current density for real states since it is identically zero.'
1524 call messages_warning(1)
1525 end if
1526 end if
1527
1529
1530 ! ---------------------------------------------------------
1531 subroutine zhamiltonian_elec_apply_all(hm, namespace, gr, st, hst)
1532 type(hamiltonian_elec_t), intent(inout) :: hm
1533 type(namespace_t), intent(in) :: namespace
1534 type(grid_t), intent(in) :: gr
1535 type(states_elec_t), intent(inout) :: st
1536 type(states_elec_t), intent(inout) :: hst
1537
1538 integer :: ik, ib, ist
1539 complex(real64), allocatable :: psi(:, :)
1540 complex(real64), allocatable :: psiall(:, :, :, :)
1541
1543
1544 do ik = st%d%kpt%start, st%d%kpt%end
1545 do ib = st%group%block_start, st%group%block_end
1546 call zhamiltonian_elec_apply_batch(hm, namespace, gr, st%group%psib(ib, ik), hst%group%psib(ib, ik))
1547 end do
1548 end do
1549
1550 if (oct_exchange_enabled(hm%oct_exchange)) then
1551
1552 safe_allocate(psiall(gr%np_part, 1:hst%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
1553
1554 call states_elec_get_state(st, gr, psiall)
1555
1556 call oct_exchange_prepare(hm%oct_exchange, gr, psiall, hm%xc, hm%psolver, namespace)
1557
1558 safe_deallocate_a(psiall)
1559
1560 safe_allocate(psi(gr%np_part, 1:hst%d%dim))
1561
1562 do ik = 1, st%nik
1563 do ist = 1, st%nst
1564 call states_elec_get_state(hst, gr, ist, ik, psi)
1565 call oct_exchange_operator(hm%oct_exchange, namespace, gr, psi, ist, ik)
1566 call states_elec_set_state(hst, gr, ist, ik, psi)
1567 end do
1568 end do
1569
1570 safe_deallocate_a(psi)
1571
1572 end if
1573
1575 end subroutine zhamiltonian_elec_apply_all
1576
1577
1578 ! ---------------------------------------------------------
1579
1580 subroutine magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase)
1581 type(hamiltonian_elec_t), intent(in) :: hm
1582 type(namespace_t), intent(in) :: namespace
1583 class(mesh_t), intent(in) :: mesh
1584 complex(real64), contiguous, intent(inout) :: psi(:,:)
1585 complex(real64), contiguous, intent(out) :: hpsi(:,:)
1586 integer, intent(in) :: ik
1587 real(real64), intent(in) :: vmagnus(:, :, :)
1588 logical, optional, intent(in) :: set_phase
1589
1590 complex(real64), allocatable :: auxpsi(:, :), aux2psi(:, :)
1591 integer :: idim, ispin
1592
1593 push_sub(magnus)
1594
1595 ! We will assume, for the moment, no spinors.
1596 if (hm%d%dim /= 1) then
1597 call messages_not_implemented("Magnus with spinors", namespace=namespace)
1598 end if
1599
1600 safe_allocate( auxpsi(1:mesh%np_part, 1:hm%d%dim))
1601 safe_allocate(aux2psi(1:mesh%np, 1:hm%d%dim))
1602
1603 ispin = hm%d%get_spin_index(ik)
1604
1605 ! Compute (T + Vnl)|psi> and store it
1606 call zhamiltonian_elec_apply_single(hm, namespace, mesh, psi, auxpsi, 1, ik, &
1607 terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
1609 ! H|psi> = (T + Vnl)|psi> + Vpsl|psi> + Vmagnus(t2)|psi> + Vborders
1610 do idim = 1, hm%d%dim
1611 call lalg_copy(mesh%np, auxpsi(:, idim), hpsi(:, idim))
1612 hpsi(1:mesh%np, idim) = hpsi(1:mesh%np, idim) + hm%ep%Vpsl(1:mesh%np)*psi(1:mesh%np,idim)
1613 call vborders(mesh, hm, psi(:, idim), hpsi(:, idim))
1614 end do
1615 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + vmagnus(1:mesh%np, ispin, 2)*psi(1:mesh%np, 1)
1616
1617 ! Add first term of the commutator: - i Vmagnus(t1) (T + Vnl) |psi>
1618 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) - m_zi*vmagnus(1:mesh%np, ispin, 1)*auxpsi(1:mesh%np, 1)
1619
1620 ! Add second term of commutator: i (T + Vnl) Vmagnus(t1) |psi>
1621 auxpsi(1:mesh%np, 1) = vmagnus(1:mesh%np, ispin, 1)*psi(1:mesh%np, 1)
1622 call zhamiltonian_elec_apply_single(hm, namespace, mesh, auxpsi, aux2psi, 1, ik, &
1623 terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
1624 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + m_zi*aux2psi(1:mesh%np, 1)
1625
1626 safe_deallocate_a(auxpsi)
1627 safe_deallocate_a(aux2psi)
1628 pop_sub(magnus)
1629 end subroutine magnus
1630
1631 ! ---------------------------------------------------------
1632 subroutine vborders (mesh, hm, psi, hpsi)
1633 class(mesh_t), intent(in) :: mesh
1634 type(hamiltonian_elec_t), intent(in) :: hm
1635 complex(real64), intent(in) :: psi(:)
1636 complex(real64), intent(inout) :: hpsi(:)
1637
1638 integer :: ip
1639
1640 push_sub(vborders)
1641
1642 if (hm%abs_boundaries%abtype == imaginary_absorbing) then
1643 do ip = 1, mesh%np
1644 hpsi(ip) = hpsi(ip) + m_zi*hm%abs_boundaries%mf(ip)*psi(ip)
1645 end do
1646 end if
1647
1648 pop_sub(vborders)
1649 end subroutine vborders
1650
1651 ! ---------------------------------------------------------
1652 logical function hamiltonian_elec_has_kick(hm)
1653 type(hamiltonian_elec_t), intent(in) :: hm
1654
1656
1657 hamiltonian_elec_has_kick = (abs(hm%kick%delta_strength) > m_epsilon)
1658
1660 end function hamiltonian_elec_has_kick
1661
1663 !
1664 subroutine hamiltonian_elec_set_mass(this, namespace, mass)
1665 class(hamiltonian_elec_t) , intent(inout) :: this
1666 type(namespace_t), intent(in) :: namespace
1667 real(real64), intent(in) :: mass
1668
1670
1671 if (parse_is_defined(namespace, 'ParticleMass')) then
1672 message(1) = 'Attempting to redefine a non-unit electron mass'
1673 call messages_fatal(1)
1674 else
1675 this%mass = mass
1676 end if
1677
1679 end subroutine hamiltonian_elec_set_mass
1680
1681 !----------------------------------------------------------
1688 subroutine hamiltonian_elec_copy_and_set_phase(hm, gr, kpt, psib, psib_with_phase)
1689 type(hamiltonian_elec_t), intent(in) :: hm
1690 type(grid_t), intent(in) :: gr
1691 type(distributed_t), intent(in) :: kpt
1692 type(wfs_elec_t), intent(in) :: psib
1693 type(wfs_elec_t), intent(out) :: psib_with_phase
1694
1695 integer :: k_offset, n_boundary_points
1696
1698
1699 call psib%copy_to(psib_with_phase)
1700 if (hm%phase%is_allocated()) then
1701 call hm%phase%apply_to(gr, gr%np, conjugate = .false., psib = psib_with_phase, src = psib, async=.true.)
1702 ! apply phase correction while setting boundary -> memory needs to be
1703 ! accessed only once
1704 k_offset = psib%ik - kpt%start
1705 n_boundary_points = int(gr%np_part - gr%np)
1706 call boundaries_set(gr%der%boundaries, gr, psib_with_phase, phase_correction = hm%phase%phase_corr(:, psib%ik), &
1707 buff_phase_corr = hm%phase%buff_phase_corr, offset=k_offset*n_boundary_points, async=.true.)
1708 else
1709 call psib%copy_data_to(gr%np, psib_with_phase)
1710 call boundaries_set(gr%der%boundaries, gr, psib_with_phase)
1711 end if
1712
1713 call psib_with_phase%do_pack(copy = .true.)
1714
1717
1718 ! ---------------------------------------------------------
1719 subroutine hamiltonian_elec_diagonal (hm, mesh, diag, ik)
1720 type(hamiltonian_elec_t), intent(in) :: hm
1721 class(mesh_t), intent(in) :: mesh
1722 real(real64), contiguous, intent(out) :: diag(:,:)
1723 integer, intent(in) :: ik
1724
1725 integer :: idim, ip, ispin
1726
1727 real(real64), allocatable :: ldiag(:)
1728
1730
1731 safe_allocate(ldiag(1:mesh%np))
1732
1733 diag = m_zero
1734
1735 call derivatives_lapl_diag(hm%der, ldiag)
1736
1737 do idim = 1, hm%d%dim
1738 diag(1:mesh%np, idim) = -m_half/hm%mass*ldiag(1:mesh%np)
1739 end do
1740
1741 select case (hm%d%ispin)
1742
1743 case (unpolarized, spin_polarized)
1744 ispin = hm%d%get_spin_index(ik)
1745 diag(1:mesh%np, 1) = diag(1:mesh%np, 1) + hm%ep%vpsl(1:mesh%np)
1746
1747 case (spinors)
1748 do ip = 1, mesh%np
1749 diag(ip, 1) = diag(ip, 1) + hm%ep%vpsl(ip)
1750 diag(ip, 2) = diag(ip, 2) + hm%ep%vpsl(ip)
1751 end do
1752
1753 end select
1754
1755 call hm%ks_pot%add_vhxc(diag)
1756
1758 end subroutine hamiltonian_elec_diagonal
1760
1761
1762#include "undef.F90"
1763#include "real.F90"
1764#include "hamiltonian_elec_inc.F90"
1765
1766#include "undef.F90"
1767#include "complex.F90"
1768#include "hamiltonian_elec_inc.F90"
1769
1770end module hamiltonian_elec_oct_m
1771
1772!! Local Variables:
1773!! mode: f90
1774!! coding: utf-8
1775!! End:
subroutine build_external_potentials()
subroutine build_phase()
subroutine external_potentials_checks()
subroutine build_interactions()
constant times a vector plus a vector
Definition: lalg_basic.F90:173
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
integer, parameter, public imaginary_absorbing
subroutine, public absorbing_boundaries_end(this)
subroutine, public absorbing_boundaries_init(this, namespace, space, gr)
pure logical function, public accel_allow_cpu_only()
Definition: accel.F90:428
subroutine, public accel_finish()
Definition: accel.F90:984
pure logical function, public accel_is_enabled()
Definition: accel.F90:418
integer, parameter, public accel_mem_read_only
Definition: accel.F90:195
This module implements batches of mesh functions.
Definition: batch.F90:135
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:118
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:124
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
real(real64) function, public derivatives_lapl_get_max_eigenvalue(this)
Get maximum eigenvalue of discrete Laplacian. For the star and star_general stencils,...
logical function, public epot_have_external_potentials(ep)
Definition: epot.F90:604
integer, parameter, public scalar_relativistic_zora
Definition: epot.F90:168
subroutine, public epot_end(ep)
Definition: epot.F90:383
integer, parameter, public fully_relativistic_zora
Definition: epot.F90:168
subroutine, public epot_init(ep, namespace, gr, ions, psolver, ispin, xc_family, kpoints)
Definition: epot.F90:220
subroutine, public epot_generate(ep, namespace, mesh, ions, st_d)
Definition: epot.F90:419
subroutine, public exchange_operator_init(this, namespace, space, st, der, mc, stencil, kpoints, cam)
subroutine, public exchange_operator_end(this)
logical function, public list_has_gauge_field(partners)
type(gauge_field_t) function, pointer, public list_get_gauge_field(partners)
logical function, public list_has_lasers(partners)
type(lasers_t) function, pointer, public list_get_lasers(partners)
integer, parameter, public external_pot_from_file
potential, defined in a file
subroutine, public load_external_potentials(external_potentials, namespace)
integer, parameter, public external_pot_charge_density
user-defined function for charge density
integer, parameter, public external_pot_usdef
user-defined function for local potential
integer, parameter, public external_pot_static_efield
Static electric field.
integer, parameter, public external_pot_static_bfield
Static magnetic field.
subroutine, public gauge_field_get_vec_pot(this, vec_pot)
real(real64), parameter, public m_two
Definition: global.F90:192
real(real64), parameter, public m_zero
Definition: global.F90:190
integer, parameter, public rdmft
Definition: global.F90:236
integer, parameter, public hartree_fock
Definition: global.F90:236
integer, parameter, public independent_particles
Theory level.
Definition: global.F90:236
integer, parameter, public generalized_kohn_sham_dft
Definition: global.F90:236
integer, parameter, public kohn_sham_dft
Definition: global.F90:236
real(real64), parameter, public m_epsilon
Definition: global.F90:206
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
Definition: global.F90:228
real(real64), parameter, public m_one
Definition: global.F90:191
This module implements the underlying real-space grid.
Definition: grid.F90:119
This module defines an abstract class for Hamiltonians.
integer, parameter, public field_uniform_magnetic_field
integer, parameter, public field_uniform_vector_potential
integer, parameter, public field_vector_potential
integer, parameter, public term_mgga
integer, parameter, public field_potential
subroutine vborders(mesh, hm, psi, hpsi)
pure logical function, public hamiltonian_elec_apply_packed(this)
subroutine, public hamiltonian_elec_set_inh(hm, st)
subroutine, public zvmask(mesh, hm, st)
subroutine, public zhamiltonian_elec_apply_batch(hm, namespace, mesh, psib, hpsib, terms, set_bc)
subroutine, public hamiltonian_elec_adjoint(hm)
subroutine, public hamiltonian_elec_end(hm)
subroutine, public hamiltonian_elec_copy_and_set_phase(hm, gr, kpt, psib, psib_with_phase)
Copy a batch to another batch and apply the Bloch phase to it.
subroutine, public zhamiltonian_elec_apply_single(hm, namespace, mesh, psi, hpsi, ist, ik, terms, set_bc, set_phase)
pure logical function hamiltonian_elec_needs_mgga_term(hm, terms)
logical function, public hamiltonian_elec_has_kick(hm)
subroutine, public magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase)
logical function hamiltonian_elec_hermitian(hm)
subroutine, public dhamiltonian_elec_apply_single(hm, namespace, mesh, psi, hpsi, ist, ik, terms, set_bc, set_phase)
subroutine, public hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
real(real64) function, public hamiltonian_elec_get_time(this)
subroutine, public dvmask(mesh, hm, st)
logical function, public hamiltonian_elec_needs_current(hm, states_are_real)
subroutine, public hamiltonian_elec_remove_inh(hm)
subroutine zhamiltonian_elec_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
subroutine, public zhamiltonian_elec_apply_atom(hm, namespace, space, latt, species, pos, ia, mesh, psi, vpsi)
subroutine, public zhamiltonian_elec_apply_all(hm, namespace, gr, st, hst)
subroutine, public hamiltonian_elec_diagonal(hm, mesh, diag, ik)
subroutine, public hamiltonian_elec_update_pot(this, mesh, accumulate)
Update the KS potential of the electronic Hamiltonian.
integer, parameter, public velocity
subroutine dhamiltonian_elec_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
subroutine hamiltonian_elec_update(this, mesh, namespace, space, ext_partners, time)
(re-)build the Hamiltonian for the next application:
subroutine hamiltonian_elec_span(hm, delta, emin, namespace)
subroutine, public hamiltonian_elec_init(hm, namespace, space, gr, ions, ext_partners, st, theory_level, xc, mc, kpoints, need_exchange, xc_photons)
subroutine dhamiltonian_elec_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
pure logical function, public hamiltonian_elec_inh_term(hm)
subroutine hamiltonian_elec_set_mass(this, namespace, mass)
set the effective electron mass, checking whether it was previously redefined.
subroutine, public dhamiltonian_elec_apply_batch(hm, namespace, mesh, psib, hpsib, terms, set_bc)
subroutine, public hamiltonian_elec_update_with_ext_pot(this, mesh, space, ext_partners, time, mu)
This is an extension of "hamiltonian_elec_update_pot" to be used by the CFM4 propagator....
subroutine, public hamiltonian_elec_not_adjoint(hm)
subroutine zhamiltonian_elec_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
This module defines classes and functions for interaction partners.
Definition: io.F90:116
integer, parameter, public kick_magnon_mode
Definition: kick.F90:165
subroutine, public kick_end(kick)
Definition: kick.F90:796
subroutine, public kick_init(kick, namespace, space, kpoints, nspin)
Definition: kick.F90:225
pure integer function, public kick_get_type(kick)
Definition: kick.F90:1365
A module to handle KS potential, without the external potential.
subroutine, public laser_vector_potential(laser, mesh, aa, time)
Definition: lasers.F90:1080
subroutine, public lasers_nondipole_laser_field_step(this, field, time)
Retrieves the NDSFA vector_potential correction. The nondipole field is obtained for consecutive time...
Definition: lasers.F90:1152
logical function, public lasers_with_nondipole_field(lasers)
Check if a nondipole SFA correction should be computed for the given laser.
Definition: lasers.F90:741
integer, parameter, public e_field_electric
Definition: lasers.F90:179
integer, parameter, public e_field_vector_potential
Definition: lasers.F90:179
subroutine, public laser_potential(laser, mesh, pot, time)
Definition: lasers.F90:1045
integer, parameter, public e_field_scalar_potential
Definition: lasers.F90:179
integer pure elemental function, public laser_kind(laser)
Definition: lasers.F90:717
subroutine, public laser_field(laser, field, time)
Retrieves the value of either the electric or the magnetic field. If the laser is given by a scalar p...
Definition: lasers.F90:1117
integer, parameter, public e_field_magnetic
Definition: lasers.F90:179
integer, parameter, public dft_u_none
Definition: lda_u.F90:203
subroutine, public lda_u_update_basis(this, space, gr, ions, st, psolver, namespace, kpoints, has_phase)
Definition: lda_u.F90:699
subroutine, public lda_u_build_phase_correction(this, space, std, boundaries, namespace, kpoints, vec_pot, vec_pot_var)
Build the phase correction to the global phase for all orbitals.
Definition: lda_u.F90:827
subroutine, public lda_u_end(this)
Definition: lda_u.F90:657
subroutine, public lda_u_init(this, namespace, space, level, gr, ions, st, mc, kpoints, has_phase)
Definition: lda_u.F90:284
This module implements fully polymorphic linked lists, and some specializations thereof.
This modules implements the routines for doing constrain DFT for noncollinear magnetism.
subroutine, public magnetic_constrain_end(this)
Releases memory of the magnetic constrain.
subroutine, public magnetic_constrain_init(this, namespace, mesh, std, natoms, min_dist)
Initilializes the magnetic_constrain_t object.
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1097
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:531
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:416
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1069
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
subroutine, public mxll_coupling_init(this, d, gr, namespace, mass)
Parse variables and initialize Maxwell coupling.
subroutine, public set_electric_quadrupole_pot(this, mesh)
Computes the electric quadrupole potential where .
subroutine, public mxll_coupling_end(this)
Finalize and deallocate Maxwell coupling arrays.
subroutine, public mxll_coupling_calc(this, hm_base, mesh, d, space)
Add the Maxwell coupling to the electronic Hamiltonian.
logical function, public oct_exchange_enabled(this)
subroutine, public oct_exchange_remove(this)
this module contains the low-level part of the output system
Definition: output_low.F90:117
Some general things and nomenclature:
Definition: par_vec.F90:173
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:505
subroutine, public pcm_calc_pot_rs(pcm, mesh, psolver, ions, v_h, v_ext, kick, time_present, kick_time)
Definition: pcm.F90:1218
subroutine, public pcm_end(pcm)
Definition: pcm.F90:3073
real(real64), dimension(:,:), allocatable delta
D_E matrix in JCP 139, 024105 (2013).
Definition: pcm.F90:271
subroutine, public pcm_init(pcm, namespace, space, ions, grid, qtot, val_charge, external_potentials_present, kick_present)
Initializes the PCM calculation: reads the VdW molecular cavity and generates the PCM response matrix...
Definition: pcm.F90:297
subroutine, public poisson_init(this, namespace, space, der, mc, stencil, qtot, label, solver, verbose, force_serial, force_cmplx)
Definition: poisson.F90:236
subroutine, public poisson_end(this)
Definition: poisson.F90:688
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
subroutine, public projector_init_phases(this, dim, std, bnd, kpoints, vec_pot, vec_pot_var)
Definition: projector.F90:264
subroutine, public scissor_end(this)
Definition: scissor.F90:241
subroutine, public states_set_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_dim_copy(dout, din)
subroutine, public states_elec_dim_end(dim)
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
This module provides routines for communicating states when using states parallelization.
integer pure function, public symmetries_identity_index(this)
Definition: symmetries.F90:591
integer pure function, public symmetries_number(this)
Definition: symmetries.F90:545
type(type_t), public type_cmplx
Definition: types.F90:136
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
This module defines the unit system, used for input and output.
type(unit_system_t), public units_inp
the units systems for reading and writing
type(xc_cam_t), parameter, public cam_exact_exchange
Use only Hartree Fock exact exchange.
Definition: xc_cam.F90:155
integer, parameter, public xc_oep_x_slater
Slater approximation to the exact exchange.
integer, parameter, public func_x
Definition: xc.F90:116
logical pure function, public family_is_mgga_with_exc(xcs)
Is the xc function part of the mGGA family with an energy functional.
Definition: xc.F90:620
logical pure function, public family_is_hybrid(xcs)
Returns true if the functional is an hybrid functional.
Definition: xc.F90:635
This module implements the "photon-free" electron-photon exchange-correlation functional.
Definition: xc_photons.F90:123
This module implements the ZORA terms for the Hamoiltonian.
Definition: zora.F90:118
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
The abstract Hamiltonian class defines a skeleton for specific implementations.
abstract class for general interaction partners
Describes mesh distribution to nodes.
Definition: mesh.F90:187
Stores all communicators and groups.
Definition: multicomm.F90:208
The states_elec_t class contains all electronic wave functions.
This class described the 'photon-exchange' electron-photon xc functionals, based on QEDFT.
Definition: xc_photons.F90:159
This class is responsible for calculating and applying the ZORA.
Definition: zora.F90:147
int true(void)