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, hm%kpoints)
401
402 !In the present implementation of DFT+U, in case of spinors, we have off-diagonal terms
403 !in spin space which break the assumption of the generalized Bloch theorem
404 if (kick_get_type(hm%kick) == kick_magnon_mode .and. gr%der%boundaries%spiral) then
405 call messages_not_implemented("DFT+U with generalized Bloch theorem and magnon kick", namespace=namespace)
406 end if
407
408 ! We rebuild the phase for the orbital projection, similarly to the one of the pseudopotentials
409 if(hm%lda_u_level /= dft_u_none .and. hm%phase%is_allocated()) then
410 call lda_u_build_phase_correction(hm%lda_u, space, hm%d, gr%der%boundaries, namespace, hm%kpoints)
411 end if
412 end if
413
414 !%Variable HamiltonianApplyPacked
415 !%Type logical
416 !%Default yes
417 !%Section Execution::Optimization
418 !%Description
419 !% If set to yes (the default), Octopus will 'pack' the
420 !% wave-functions when operating with them. This might involve some
421 !% additional copying but makes operations more efficient.
422 !% See also the related <tt>StatesPack</tt> variable.
423 !%End
424 call parse_variable(namespace, 'HamiltonianApplyPacked', .true., hm%is_applied_packed)
425
426 if (hm%theory_level == hartree_fock .and. st%parallel_in_states) then
427 call messages_experimental('Hartree-Fock parallel in states', namespace=namespace)
428 end if
429
430 if (hm%theory_level == generalized_kohn_sham_dft .and. family_is_hybrid(hm%xc) &
431 .and. st%parallel_in_states) then
432 call messages_experimental('Hybrid functionals parallel in states', namespace=namespace)
433 end if
434
435 !%Variable TimeZero
436 !%Type logical
437 !%Default no
438 !%Section Hamiltonian
439 !%Description
440 !% (Experimental) If set to yes, the ground state and other time
441 !% dependent calculation will assume that they are done at time
442 !% zero, so that all time depedent field at that time will be
443 !% included.
444 !%End
445 call parse_variable(namespace, 'TimeZero', .false., hm%time_zero)
446 if (hm%time_zero) call messages_experimental('TimeZero', namespace=namespace)
447
448 !Cam parameters are irrelevant here and are updated later
449 need_exchange_ = optional_default(need_exchange, .false.)
450 if (hm%xc%compute_exchange(hm%theory_level) .or. need_exchange_) then
451 !We test Slater before OEP, as Slater is treated as OEP for the moment....
452 if (hm%xc%functional(func_x,1)%id == xc_oep_x_slater) then
453 call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
454 hm%kpoints, cam_exact_exchange)
455 else if (bitand(hm%xc%family, xc_family_oep) /= 0 .or. hm%theory_level == rdmft) then
456 call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
457 hm%kpoints, hm%xc%cam)
458 if (hm%theory_level == rdmft) hm%exxop%useACE = .false.
459 else
460 call exchange_operator_init(hm%exxop, namespace, space, st, gr%der, mc, gr%stencil, &
461 hm%kpoints, cam_exact_exchange)
462 end if
463 end if
464
465 if (hm%is_applied_packed .and. accel_is_enabled()) then
466 ! Check if we can actually apply the hamiltonian packed
467 if (gr%use_curvilinear) then
468 if (accel_allow_cpu_only()) then
469 hm%is_applied_packed = .false.
470 call messages_write('Cannot use CUDA or OpenCL as curvilinear coordinates are used.')
471 call messages_warning(namespace=namespace)
472 else
473 call messages_write('Cannot use CUDA or OpenCL as curvilinear coordinates are used.', new_line = .true.)
474 call messages_write('Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
475 call messages_fatal(namespace=namespace)
476 end if
477 end if
478 end if
479
480 !We are building the list of external potentials
481 !This is done here at the moment, because we pass directly the mesh
482 !TODO: Once the abstract Hamiltonian knows about an abstract basis, we might move this to the
483 ! abstract Hamiltonian
484 call load_external_potentials(hm%external_potentials, namespace)
485
486 !Some checks which are electron specific, like k-points
488
489 !At the moment we do only have static external potential, so we never update them
491
492 !Build the resulting interactions
493 !TODO: This will be moved to the actual interactions
494 call build_interactions()
495
496 ! Constrained DFT for noncollinear magnetism
497 if (hm%theory_level /= independent_particles) then
498 call magnetic_constrain_init(hm%magnetic_constrain, namespace, gr, st%d, ions%natoms, ions%min_distance())
499 end if
500
501 ! init maxwell-electrons coupling
502 call mxll_coupling_init(hm%mxll, st%d, gr, namespace, hm%mass)
503
504 if (associated(hm%xc_photons)) then
505 if (hm%xc_photons%wants_to_renormalize_mass()) then
506 ! remornalize the electron mass due to light-matter interaction; here we only deal with it in free space
507 call hm%set_mass(namespace, hm%xc_photons%get_renormalized_mass())
508 end if
509 end if
510
511 if (hm%xc%compute_exchange(hm%theory_level) .or. need_exchange_) call hm%exxop%write_info(namespace)
512
513 call profiling_out('HAMILTONIAN_ELEC_INIT')
514 pop_sub(hamiltonian_elec_init)
515
516 contains
517
518 ! ---------------------------------------------------------
519 subroutine build_external_potentials()
520 type(list_iterator_t) :: iter
521 class(*), pointer :: potential
522 integer :: iop
523
525
526 safe_allocate(hm%v_ext_pot(1:gr%np))
527 hm%v_ext_pot(1:gr%np) = m_zero
528
529 call iter%start(hm%external_potentials)
530 do while (iter%has_next())
531 potential => iter%get_next()
532 select type (potential)
533 class is (external_potential_t)
534
535 call potential%allocate_memory(gr)
536 call potential%calculate(namespace, gr, hm%psolver)
537 !To preserve the old behavior, we are adding the various potentials
538 !to the corresponding arrays
539 select case (potential%type)
541 call lalg_axpy(gr%np, m_one, potential%pot, hm%v_ext_pot)
542
544 if (states_are_real(st)) then
545 message(1) = "Cannot use static magnetic field with real wavefunctions"
546 call messages_fatal(1, namespace=namespace)
547 end if
548
549 if (.not. allocated(hm%ep%b_field)) then
550 safe_allocate(hm%ep%b_field(1:3)) !Cannot be space%dim
551 hm%ep%b_field(1:3) = m_zero
552 end if
553 hm%ep%b_field(1:3) = hm%ep%b_field(1:3) + potential%b_field(1:3)
554
555 if (.not. allocated(hm%ep%a_static)) then
556 safe_allocate(hm%ep%a_static(1:gr%np, 1:space%dim))
557 hm%ep%a_static(1:gr%np, 1:space%dim) = m_zero
558 end if
559 call lalg_axpy(gr%np, space%dim, m_one, potential%a_static, hm%ep%a_static)
560
562 if (.not. allocated(hm%ep%e_field)) then
563 safe_allocate(hm%ep%e_field(1:space%dim))
564 hm%ep%e_field(1:space%dim) = m_zero
565 end if
566 hm%ep%e_field(1:space%dim) = hm%ep%e_field(1:space%dim) + potential%e_field(1:space%dim)
567
568 !In the fully periodic case, we use Berry phases
569 if (space%periodic_dim < space%dim) then
570 if (.not. allocated(hm%v_static)) then
571 safe_allocate(hm%v_static(1:gr%np))
572 hm%v_static(1:gr%np) = m_zero
573 end if
574 if (.not. allocated(hm%ep%v_ext)) then
575 safe_allocate(hm%ep%v_ext(1:gr%np_part))
576 hm%ep%v_ext(1:gr%np_part) = m_zero
577 end if
578 call lalg_axpy(gr%np, m_one, potential%pot, hm%v_static)
579 call lalg_axpy(gr%np, m_one, potential%v_ext, hm%ep%v_ext)
580 end if
581
582 if (hm%kpoints%use_symmetries) then
583 do iop = 1, symmetries_number(hm%kpoints%symm)
584 if (iop == symmetries_identity_index(hm%kpoints%symm)) cycle
585 if (.not. symm_op_invariant_cart(hm%kpoints%symm%ops(iop), hm%ep%e_field, 1e-5_real64)) then
586 message(1) = "The StaticElectricField breaks (at least) one of the symmetries used to reduce the k-points."
587 message(2) = "Set SymmetryBreakDir equal to StaticElectricField."
588 call messages_fatal(2, namespace=namespace)
589 end if
590 end do
591 end if
592
593 end select
594 call potential%deallocate_memory()
595
596 class default
597 assert(.false.)
598 end select
599 end do
600
602 end subroutine build_external_potentials
603
604 ! ---------------------------------------------------------
605 subroutine external_potentials_checks()
606 type(list_iterator_t) :: iter
607 class(*), pointer :: potential
608
610
611 call iter%start(hm%external_potentials)
612 do while (iter%has_next())
613 potential => iter%get_next()
614 select type (potential)
615 class is (external_potential_t)
616
617 if (potential%type == external_pot_static_efield .and. hm%kpoints%reduced%npoints > 1) then
618 message(1) = "Applying StaticElectricField in a periodic direction is only accurate for large supercells."
619 message(2) = "Single-point Berry phase is not appropriate when k-point sampling is needed."
620 call messages_warning(2, namespace=namespace)
621 end if
622
623 class default
624 assert(.false.)
625 end select
626 end do
627
629 end subroutine external_potentials_checks
630
631
632 !The code in this routines needs to know about the external potentials.
633 !This will be treated in the future by the interactions directly.
634 subroutine build_interactions()
635 logical :: external_potentials_present
636 logical :: kick_present
637
639
640 if (allocated(hm%ep%e_field) .and. space%is_periodic() .and. .not. list_has_gauge_field(ext_partners)) then
641 ! only need vberry if there is a field in a periodic direction
642 ! and we are not setting a gauge field
643 if (any(abs(hm%ep%e_field(1:space%periodic_dim)) > m_epsilon)) then
644 safe_allocate(hm%vberry(1:gr%np, 1:hm%d%nspin))
645 hm%vberry = m_zero
646 end if
647 end if
648
649 external_potentials_present = epot_have_external_potentials(hm%ep) .or. &
650 list_has_lasers(ext_partners) .or. allocated(hm%v_static)
651
652
653 kick_present = hamiltonian_elec_has_kick(hm)
654
655 call pcm_init(hm%pcm, namespace, space, ions, gr, st%qtot, st%val_charge, external_potentials_present, kick_present)
656 if (hm%pcm%run_pcm) then
657 if (hm%theory_level /= kohn_sham_dft) call messages_not_implemented("PCM for TheoryLevel /= kohn_sham", namespace=namespace)
658 end if
659
661
662 end subroutine build_interactions
663
664
665 end subroutine hamiltonian_elec_init
666
667
668
669
670 ! ---------------------------------------------------------
671 subroutine hamiltonian_elec_end(hm)
672 type(hamiltonian_elec_t), target, intent(inout) :: hm
673
674 type(partner_iterator_t) :: iter
675 class(interaction_partner_t), pointer :: potential
676
677 push_sub(hamiltonian_elec_end)
678
679 call hm%hm_base%end()
680 call hm%vnl%end()
681
682 call hm%phase%end()
683
684 call hm%ks_pot%end()
685 safe_deallocate_a(hm%vberry)
686 safe_deallocate_a(hm%a_ind)
687 safe_deallocate_a(hm%b_ind)
688 safe_deallocate_a(hm%v_ext_pot)
689
690 safe_deallocate_p(hm%zora)
691
692 call poisson_end(hm%psolver)
693
694 nullify(hm%xc)
695
696 call kick_end(hm%kick)
697 call epot_end(hm%ep)
698 nullify(hm%ions)
699
700 call absorbing_boundaries_end(hm%abs_boundaries)
701
702 call states_elec_dim_end(hm%d)
703
704 if (hm%scissor%apply) call scissor_end(hm%scissor)
705
706 call exchange_operator_end(hm%exxop)
707 call lda_u_end(hm%lda_u)
708
709 safe_deallocate_a(hm%energy)
710
711 if (hm%pcm%run_pcm) call pcm_end(hm%pcm)
712
713 call hm%v_ie_loc%end()
714 call hm%nlcc%end()
715
716 call iter%start(hm%external_potentials)
717 do while (iter%has_next())
718 potential => iter%get_next()
719 safe_deallocate_p(potential)
720 end do
721 call hm%external_potentials%empty()
722 safe_deallocate_a(hm%v_static)
723
724 call magnetic_constrain_end(hm%magnetic_constrain)
725
726 call mxll_coupling_end(hm%mxll)
727
728 pop_sub(hamiltonian_elec_end)
729 end subroutine hamiltonian_elec_end
730
731
732 ! ---------------------------------------------------------
733 ! True if the Hamiltonian is Hermitian, false otherwise
734 logical function hamiltonian_elec_hermitian(hm)
735 class(hamiltonian_elec_t), intent(in) :: hm
736
738 hamiltonian_elec_hermitian = .not.((hm%abs_boundaries%abtype == imaginary_absorbing) .or. &
739 oct_exchange_enabled(hm%oct_exchange))
740
742 end function hamiltonian_elec_hermitian
743
744 ! ---------------------------------------------------------
745 ! True if the Hamiltonian needs to apply the mGGA term
746 pure logical function hamiltonian_elec_needs_mgga_term(hm, terms)
747 class(hamiltonian_elec_t), intent(in) :: hm
748 integer, intent(in) :: terms
749
751
752 if (bitand(term_mgga, terms) /= 0 .and. family_is_mgga_with_exc(hm%xc) &
753 .and. hm%theory_level == generalized_kohn_sham_dft) then
755 end if
756
757 !For OEP
758 if(terms == term_mgga .and. family_is_mgga_with_exc(hm%xc) .and. hm%theory_level == kohn_sham_dft) then
760 end if
761
763
764
765
766 ! ---------------------------------------------------------
767 subroutine hamiltonian_elec_span(hm, delta, emin, namespace)
768 class(hamiltonian_elec_t), intent(inout) :: hm
769 real(real64), intent(in) :: delta(:)
770 real(real64), intent(in) :: emin
771 type(namespace_t), intent(in) :: namespace
772
773 real(real64) :: emax
774
775 push_sub(hamiltonian_elec_span)
776
777 ! estimate maximum energy of discrete kinetic operator
778 ! this neglects possible contributions from the non-local part of the pseudopotentials
780
781 hm%spectral_middle_point = (emax + emin) / m_two
782 hm%spectral_half_span = (emax - emin) / m_two
783
784 pop_sub(hamiltonian_elec_span)
785 end subroutine hamiltonian_elec_span
786
787
788 ! ---------------------------------------------------------
789 pure logical function hamiltonian_elec_inh_term(hm) result(inh)
790 type(hamiltonian_elec_t), intent(in) :: hm
791
792 inh = hm%inh_term
793 end function hamiltonian_elec_inh_term
794
795
796 ! ---------------------------------------------------------
797 subroutine hamiltonian_elec_set_inh(hm, st)
798 type(hamiltonian_elec_t), intent(inout) :: hm
799 type(states_elec_t), intent(in) :: st
800
802
803 if (hm%inh_term) call states_elec_end(hm%inh_st)
804 call states_elec_copy(hm%inh_st, st)
805 hm%inh_term = .true.
806
808 end subroutine hamiltonian_elec_set_inh
809
810
811 ! ---------------------------------------------------------
812 subroutine hamiltonian_elec_remove_inh(hm)
813 type(hamiltonian_elec_t), intent(inout) :: hm
814
816
817 if (hm%inh_term) then
818 call states_elec_end(hm%inh_st)
819 hm%inh_term = .false.
820 end if
821
823 end subroutine hamiltonian_elec_remove_inh
824
825 ! ---------------------------------------------------------
826 subroutine hamiltonian_elec_adjoint(hm)
827 type(hamiltonian_elec_t), intent(inout) :: hm
828
830
831 if (.not. hm%adjoint) then
832 hm%adjoint = .true.
833 if (hm%abs_boundaries%abtype == imaginary_absorbing) then
834 hm%abs_boundaries%mf = -hm%abs_boundaries%mf
835 end if
836 end if
837
839 end subroutine hamiltonian_elec_adjoint
840
842 ! ---------------------------------------------------------
843 subroutine hamiltonian_elec_not_adjoint(hm)
844 type(hamiltonian_elec_t), intent(inout) :: hm
845
847
848 if (hm%adjoint) then
849 hm%adjoint = .false.
850 if (hm%abs_boundaries%abtype == imaginary_absorbing) then
851 hm%abs_boundaries%mf = -hm%abs_boundaries%mf
852 end if
853 end if
854
856 end subroutine hamiltonian_elec_not_adjoint
857
858
859 ! ---------------------------------------------------------
861 subroutine hamiltonian_elec_update(this, mesh, namespace, space, ext_partners, time)
862 class(hamiltonian_elec_t), intent(inout) :: this
863 class(mesh_t), intent(in) :: mesh
864 type(namespace_t), intent(in) :: namespace
865 class(space_t), intent(in) :: space
866 type(partner_list_t), intent(in) :: ext_partners
867 real(real64), optional, intent(in) :: time
868
869 integer :: ispin, ip, idir, iatom, ilaser
870 real(real64) :: aa(space%dim), time_
871 real(real64), allocatable :: vp(:,:)
872 type(lasers_t), pointer :: lasers
873 type(gauge_field_t), pointer :: gfield
874 real(real64) :: am(space%dim)
875
877 call profiling_in("HAMILTONIAN_ELEC_UPDATE")
878
879 this%current_time = m_zero
880 if (present(time)) this%current_time = time
881
882 time_ = optional_default(time, 0.0_real64)
883
884 ! set everything to zero
885 call this%hm_base%clear(mesh%np)
886
887 ! alllocate the scalar potential for the xc, hartree and external potentials
888 call this%hm_base%allocate_field(mesh, field_potential, &
889 complex_potential = this%abs_boundaries%abtype == imaginary_absorbing)
890
891 ! the lasers
892 if (present(time) .or. this%time_zero) then
893
894 lasers => list_get_lasers(ext_partners)
895 if(associated(lasers)) then
896 do ilaser = 1, lasers%no_lasers
897 select case (laser_kind(lasers%lasers(ilaser)))
899 do ispin = 1, this%d%spin_channels
900 call laser_potential(lasers%lasers(ilaser), mesh, &
901 this%hm_base%potential(:, ispin), time_)
902 end do
903 case (e_field_magnetic)
904 call this%hm_base%allocate_field(mesh, field_vector_potential + field_uniform_magnetic_field, &
905 .false.)
906 ! get the vector potential
907 safe_allocate(vp(1:mesh%np, 1:space%dim))
908 vp(1:mesh%np, 1:space%dim) = m_zero
909 call laser_vector_potential(lasers%lasers(ilaser), mesh, vp, time_)
910 !$omp parallel do private(idir) schedule(static)
911 do ip = 1, mesh%np
912 do idir = 1, space%dim
913 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + vp(ip, idir)/p_c
914 end do
915 end do
916 ! and the magnetic field
917 call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
918 safe_deallocate_a(vp)
920 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
921 ! get the uniform vector potential associated with a magnetic field
922 aa = m_zero
923 call laser_field(lasers%lasers(ilaser), aa, time_)
924 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
925 end select
926 end do
927
928 if (lasers_with_nondipole_field(lasers)) then
929 assert( allocated(this%hm_base%uniform_vector_potential))
930 call lasers_nondipole_laser_field_step(lasers, am, time_)
931 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - am/p_c
932 end if
933 end if
934
935 ! the gauge field
936 gfield => list_get_gauge_field(ext_partners)
937 if (associated(gfield)) then
938 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
939 call gauge_field_get_vec_pot(gfield, aa)
940 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
941 end if
942
943 ! the electric field for a periodic system through the gauge field
944 if (allocated(this%ep%e_field) .and. associated(gfield)) then
945 this%hm_base%uniform_vector_potential(1:space%periodic_dim) = &
946 this%hm_base%uniform_vector_potential(1:space%periodic_dim) - time_*this%ep%e_field(1:space%periodic_dim)
947 end if
948
949 ! add the photon-free mean-field vector potential
950 if (associated(this%xc_photons)) then
951 if(this%xc_photons%lpfmf .and. allocated(this%xc_photons%mf_vector_potential)) then
952 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
953 ! here we put a minus sign in front of the mean field term to get the right answer (need to check the formula)
954 this%hm_base%uniform_vector_potential(1:space%dim) = &
955 this%hm_base%uniform_vector_potential(1:space%dim) - this%xc_photons%mf_vector_potential(1:space%dim)/p_c
956 end if
957 end if
958
959 end if
960
961 ! the vector potential of a static magnetic field
962 if (allocated(this%ep%a_static)) then
963 call this%hm_base%allocate_field(mesh, field_vector_potential, .false.)
964 !ep%a_static contains 1/c A(r)
965 !$omp parallel do private(idir) schedule(static)
966 do ip = 1, mesh%np
967 do idir = 1, space%dim
968 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + this%ep%a_static(ip, idir)
969 end do
970 end do
971 end if
972
973 ! add Maxwell coupling to Hamiltonian and sets the magnetic field for the Zeeman term added below
974 call mxll_coupling_calc(this%mxll, this%hm_base, mesh, this%d, space)
975
976 !The electric field was added to the KS potential
977 call this%hm_base%accel_copy_pot(mesh)
978
979 ! and the static magnetic field
980 if (allocated(this%ep%b_field)) then
981 call this%hm_base%allocate_field(mesh, field_uniform_magnetic_field, .false.)
982 do idir = 1, 3
983 this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
984 end do
985 end if
986
987 ! Combine the uniform and non-uniform fields and compute the Zeeman term
988 call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
989
990 ! This needs to be called at the end as the zeeman term enters the potential
991 call hamiltonian_elec_update_pot(this, mesh, accumulate = .true.)
992
993 if (this%mxll%test_equad) then
994 call set_electric_quadrupole_pot(this%mxll, mesh)
995 end if
996
997 call build_phase()
998
999 call profiling_out("HAMILTONIAN_ELEC_UPDATE")
1001
1002 contains
1003
1004 subroutine build_phase()
1005 integer :: ik, imat, nmat, max_npoints, offset
1006 integer :: ip
1007 integer :: iphase, nphase
1008
1010
1011 if ((.not. this%kpoints%gamma_only()) .or. allocated(this%hm_base%uniform_vector_potential)) then
1012
1013 call profiling_in('UPDATE_PHASES')
1014 ! now regenerate the phases for the pseudopotentials
1015 do iatom = 1, this%ep%natoms
1016 call projector_init_phases(this%ep%proj(iatom), space%dim, this%d, this%der%boundaries, this%kpoints, &
1017 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1018 end do
1019
1020 call profiling_out('UPDATE_PHASES')
1021 end if
1022
1023 if (allocated(this%hm_base%uniform_vector_potential)) then
1024
1025 call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
1026
1027 ! We rebuild the phase for the orbital projection, similarly to the one of the pseudopotentials
1028 if (this%lda_u_level /= dft_u_none) then
1029 call lda_u_build_phase_correction(this%lda_u, space, this%d, this%der%boundaries, namespace, this%kpoints, &
1030 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1031 end if
1032 end if
1033
1034 max_npoints = this%vnl%max_npoints
1035 nmat = this%vnl%nprojector_matrices
1036
1037
1038 if (this%phase%is_allocated() .and. allocated(this%vnl%projector_matrices)) then
1039
1040 nphase = 1
1041 if (this%der%boundaries%spiralBC) nphase = 3
1042
1043 if (.not. allocated(this%vnl%projector_phases)) then
1044 safe_allocate(this%vnl%projector_phases(1:max_npoints, 1:nphase, nmat, this%d%kpt%start:this%d%kpt%end))
1045 if (accel_is_enabled()) then
1046 call accel_create_buffer(this%vnl%buff_projector_phases, accel_mem_read_only, &
1047 type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
1048 ! We need to save nphase, with which the array has been build,
1049 ! as the number might change throughout the run
1050 this%vnl%nphase = nphase
1051 end if
1052 end if
1053
1054 offset = 0
1055 do ik = this%d%kpt%start, this%d%kpt%end
1056 do imat = 1, this%vnl%nprojector_matrices
1057 iatom = this%vnl%projector_to_atom(imat)
1058 do iphase = 1, nphase
1059 !$omp parallel do simd schedule(static)
1060 do ip = 1, this%vnl%projector_matrices(imat)%npoints
1061 this%vnl%projector_phases(ip, iphase, imat, ik) = this%ep%proj(iatom)%phase(ip, iphase, ik)
1062 end do
1063
1064 if (accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0) then
1065 call accel_write_buffer(this%vnl%buff_projector_phases, &
1066 this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
1067 offset = offset, async=.true.)
1068 end if
1069 offset = offset + this%vnl%projector_matrices(imat)%npoints
1070 end do
1071 end do
1072 end do
1073
1074 end if
1075
1076 call accel_finish()
1077
1079 end subroutine build_phase
1080
1081 end subroutine hamiltonian_elec_update
1082
1083
1084 !----------------------------------------------------------------
1087 ! TODO: See Issue #1064
1088 subroutine hamiltonian_elec_update_pot(this, mesh, accumulate)
1089 type(hamiltonian_elec_t), intent(inout) :: this
1090 class(mesh_t), intent(in) :: mesh
1091 logical, optional, intent(in) :: accumulate
1092
1093 integer :: ispin, ip
1094
1096
1097 ! By default we nullify first the result
1098 if (.not. optional_default(accumulate, .false.)) then
1099 !$omp parallel private(ip, ispin)
1100 do ispin = 1, this%d%nspin
1101 !$omp do simd schedule(static)
1102 do ip = 1, mesh%np
1103 this%hm_base%potential(ip, ispin) = m_zero
1104 end do
1105 end do
1106 !$omp end parallel
1107 end if
1108
1109 !$omp parallel private(ip, ispin)
1110 do ispin = 1, this%d%nspin
1111 if (ispin <= 2) then
1112 !$omp do simd schedule(static)
1113 ! this%vhxc(ip, ispin) is added after the calculation of ZORA potential
1114 do ip = 1, mesh%np
1115 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%ep%vpsl(ip) + this%v_ext_pot(ip)
1116 end do
1117
1119 if (this%pcm%run_pcm) then
1120 if (this%pcm%solute) then
1121 !$omp do simd schedule(static)
1122 do ip = 1, mesh%np
1123 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
1124 this%pcm%v_e_rs(ip) + this%pcm%v_n_rs(ip)
1125 end do
1126 !$omp end do simd nowait
1127 end if
1128 if (this%pcm%localf) then
1129 !$omp do simd schedule(static)
1130 do ip = 1, mesh%np
1131 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
1132 this%pcm%v_ext_rs(ip)
1133 end do
1134 !$omp end do simd nowait
1135 end if
1136 end if
1137
1139 if (this%abs_boundaries%abtype == imaginary_absorbing) then
1140 !$omp do simd schedule(static)
1141 do ip = 1, mesh%np
1142 this%hm_base%Impotential(ip, ispin) = this%hm_base%Impotential(ip, ispin) + this%abs_boundaries%mf(ip)
1143 end do
1144 !$omp end do simd nowait
1145 end if
1146 end if
1147 end do
1148 !$omp end parallel
1149
1150 ! scalar relativistic ZORA contribution
1151 ! \boldsymbol{p} \frac{c^2}{2c^2 - V} \boldsymbol{p} \Phi^\mathrm{ZORA}
1152 if (this%ep%reltype == scalar_relativistic_zora .or. this%ep%reltype == fully_relativistic_zora) then
1153 call this%zora%update(this%der, this%hm_base%potential)
1154 end if
1155
1156 !$omp parallel private(ip, ispin)
1157 do ispin = 1, this%d%nspin
1158 ! Adding Zeeman potential to hm_base%potential
1159 if (allocated(this%hm_base%zeeman_pot)) then
1160 !$omp do simd schedule(static)
1161 do ip = 1, mesh%np
1162 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%hm_base%zeeman_pot(ip, ispin)
1163 end do
1164 !$omp end do simd nowait
1165 end if
1166
1167 ! Adding Quadrupole potential from static E-field (test)
1168 if (this%mxll%test_equad) then
1169 !$omp do simd schedule(static)
1170 do ip = 1, mesh%np
1171 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%mxll%e_quadrupole_pot(ip)
1172 end do
1173 !$omp end do simd
1174 end if
1175 end do
1176 !$omp end parallel
1177
1178
1179 ! Add the Hartree and KS potential
1180 call this%ks_pot%add_vhxc(this%hm_base%potential)
1181
1182 call this%hm_base%accel_copy_pot(mesh)
1185 end subroutine hamiltonian_elec_update_pot
1186
1187 ! ---------------------------------------------------------
1188 subroutine hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
1189 type(hamiltonian_elec_t), intent(inout) :: this
1190 type(namespace_t), intent(in) :: namespace
1191 class(electron_space_t), intent(in) :: space
1192 type(grid_t), intent(in) :: gr
1193 type(ions_t), target, intent(inout) :: ions
1194 type(partner_list_t), intent(in) :: ext_partners
1195 type(states_elec_t), intent(inout) :: st
1196 real(real64), optional, intent(in) :: time
1197
1199
1200 this%ions => ions
1201 call epot_generate(this%ep, namespace, gr, this%ions, this%d)
1202
1203 ! Interation terms are treated below
1204
1205 ! First we add the static electric field
1206 if (allocated(this%ep%e_field) .and. space%periodic_dim < space%dim) then
1207 call lalg_axpy(gr%np, m_one, this%v_static, this%ep%vpsl)
1208 end if
1209
1210 ! Here we need to pass this again, else test are failing.
1211 ! This is not a real problem, as the multisystem framework will indeed to this anyway
1212 this%v_ie_loc%atoms_dist => ions%atoms_dist
1213 this%v_ie_loc%atom => ions%atom
1214 call this%v_ie_loc%calculate()
1215
1216 ! At the moment we need to add this to ep%vpsl, to keep the behavior of the code
1217 call lalg_axpy(gr%np, m_one, this%v_ie_loc%potential(:,1), this%ep%vpsl)
1218
1219 ! Here we need to reinit the NLCC object
1220 ! This is not a real problem, as the multisystem framework will indeed to this anyway
1221 if (this%ep%nlcc) then
1222 call this%nlcc%end()
1223 call this%nlcc%init(gr, ions)
1224 call this%nlcc%calculate()
1225 call lalg_copy(gr%np, this%nlcc%density(:,1), st%rho_core)
1226 end if
1227
1228 call this%vnl%build(space, gr, this%ep)
1229 call hamiltonian_elec_update(this, gr, namespace, space, ext_partners, time)
1230
1231 ! Check if projectors are still compatible with apply_packed on GPU
1232 if (this%is_applied_packed .and. accel_is_enabled()) then
1233 if (this%ep%non_local .and. .not. this%vnl%apply_projector_matrices) then
1234 if (accel_allow_cpu_only()) then
1235 call messages_write('Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.')
1236 call messages_warning(namespace=namespace)
1237 else
1238 call messages_write('Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.',&
1239 new_line = .true.)
1240 call messages_write('Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
1241 call messages_fatal(namespace=namespace)
1242 end if
1243 end if
1244
1245 end if
1246
1247 if (this%pcm%run_pcm) then
1250 if (this%pcm%solute) then
1251 call pcm_calc_pot_rs(this%pcm, gr, this%psolver, ions = ions)
1252 end if
1253
1256 ! Interpolation is needed, hence gr%np_part -> 1:gr%np
1257 if (this%pcm%localf .and. allocated(this%v_static)) then
1258 call pcm_calc_pot_rs(this%pcm, gr, this%psolver, v_ext = this%ep%v_ext(1:gr%np_part))
1259 end if
1260
1261 end if
1262
1263 call lda_u_update_basis(this%lda_u, space, gr, ions, st, this%psolver, namespace, this%kpoints, &
1264 this%phase%is_allocated())
1265
1267 end subroutine hamiltonian_elec_epot_generate
1268
1269 ! -----------------------------------------------------------------
1270
1271 real(real64) function hamiltonian_elec_get_time(this) result(time)
1272 type(hamiltonian_elec_t), intent(inout) :: this
1273
1274 time = this%current_time
1275 end function hamiltonian_elec_get_time
1276
1277 ! -----------------------------------------------------------------
1278
1279 pure logical function hamiltonian_elec_apply_packed(this) result(apply)
1280 class(hamiltonian_elec_t), intent(in) :: this
1281
1282 apply = this%is_applied_packed
1285
1286
1287 ! -----------------------------------------------------------------
1288 subroutine zhamiltonian_elec_apply_atom (hm, namespace, space, latt, species, pos, ia, mesh, psi, vpsi)
1289 type(hamiltonian_elec_t), intent(in) :: hm
1290 type(namespace_t), intent(in) :: namespace
1291 class(space_t), intent(in) :: space
1292 type(lattice_vectors_t), intent(in) :: latt
1293 class(species_t), intent(in) :: species
1294 real(real64), intent(in) :: pos(1:space%dim)
1295 integer, intent(in) :: ia
1296 class(mesh_t), intent(in) :: mesh
1297 complex(real64), intent(in) :: psi(:,:)
1298 complex(real64), intent(out) :: vpsi(:,:)
1299
1300 integer :: idim
1301 real(real64), allocatable :: vlocal(:)
1303
1304 safe_allocate(vlocal(1:mesh%np_part))
1305 vlocal = m_zero
1306 call epot_local_potential(hm%ep, namespace, space, latt, mesh, species, pos, ia, vlocal)
1307
1308 do idim = 1, hm%d%dim
1309 vpsi(1:mesh%np, idim) = vlocal(1:mesh%np) * psi(1:mesh%np, idim)
1310 end do
1311
1312 safe_deallocate_a(vlocal)
1314 end subroutine zhamiltonian_elec_apply_atom
1315
1316 ! ---------------------------------------------------------
1321 subroutine hamiltonian_elec_update_with_ext_pot(this, mesh, space, ext_partners, time, mu)
1322 type(hamiltonian_elec_t), intent(inout) :: this
1323 class(space_t), intent(in) :: space
1324 class(mesh_t), intent(in) :: mesh
1325 type(partner_list_t), intent(in) :: ext_partners
1326 real(real64), intent(in) :: time(1:2)
1327 real(real64), intent(in) :: mu(1:2)
1328
1329 integer :: ispin, ip, idir, iatom, ilaser, itime
1330 real(real64) :: aa(space%dim), time_
1331 real(real64), allocatable :: vp(:,:)
1332 real(real64), allocatable :: velectric(:)
1333 type(lasers_t), pointer :: lasers
1334 type(gauge_field_t), pointer :: gfield
1335
1337 call profiling_in("HAMILTONIAN_ELEC_UPDATE_EXT_POT")
1338
1339 this%current_time = m_zero
1340 this%current_time = time(1)
1341
1342 ! set everything to zero
1343 call this%hm_base%clear(mesh%np)
1344
1345 ! the xc, hartree and external potentials
1346 call this%hm_base%allocate_field(mesh, field_potential, &
1347 complex_potential = this%abs_boundaries%abtype == imaginary_absorbing)
1348
1349 do itime = 1, 2
1350 time_ = time(itime)
1351
1352 lasers => list_get_lasers(ext_partners)
1353 if(associated(lasers)) then
1354 do ilaser = 1, lasers%no_lasers
1355 select case (laser_kind(lasers%lasers(ilaser)))
1356 case (e_field_scalar_potential, e_field_electric)
1357 safe_allocate(velectric(1:mesh%np))
1358 do ispin = 1, this%d%spin_channels
1359 velectric = m_zero
1360 call laser_potential(lasers%lasers(ilaser), mesh, velectric, time_)
1361 !$omp parallel do simd schedule(static)
1362 do ip = 1, mesh%np
1363 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + mu(itime) * velectric(ip)
1364 end do
1365 end do
1366 safe_deallocate_a(velectric)
1367 case (e_field_magnetic)
1368 call this%hm_base%allocate_field(mesh, field_vector_potential + field_uniform_magnetic_field, .false.)
1369 ! get the vector potential
1370 safe_allocate(vp(1:mesh%np, 1:space%dim))
1371 vp(1:mesh%np, 1:space%dim) = m_zero
1372 call laser_vector_potential(lasers%lasers(ilaser), mesh, vp, time_)
1373 do idir = 1, space%dim
1374 !$omp parallel do schedule(static)
1375 do ip = 1, mesh%np
1376 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) &
1377 - mu(itime) * vp(ip, idir)/p_c
1378 end do
1379 end do
1380 ! and the magnetic field
1381 call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
1382 safe_deallocate_a(vp)
1383 case (e_field_vector_potential)
1384 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
1385 ! get the uniform vector potential associated with a magnetic field
1386 aa = m_zero
1387 call laser_field(lasers%lasers(ilaser), aa, time_)
1388 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) &
1389 - mu(itime) * aa/p_c
1390 end select
1391 end do
1392 end if
1393
1394 ! the gauge field
1395 gfield => list_get_gauge_field(ext_partners)
1396 if (associated(gfield)) then
1397 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
1398 call gauge_field_get_vec_pot(gfield, aa)
1399 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
1400 end if
1401
1402 ! the electric field for a periodic system through the gauge field
1403 ! TODO: The condition is wrong here: the e_field should be in non-periodic dims as E field
1404 ! and as a gauge field in the periodic dim, unless we use a Bery phase, in which, we do not use it
1405 ! this way. But this is unrelated to the gauge field
1406 if (allocated(this%ep%e_field) .and. associated(gfield)) then
1407 this%hm_base%uniform_vector_potential(1:space%periodic_dim) = &
1408 this%hm_base%uniform_vector_potential(1:space%periodic_dim) - time_*this%ep%e_field(1:space%periodic_dim)
1409 end if
1410
1411 end do
1412
1413 ! the vector potential of a static magnetic field
1414 if (allocated(this%ep%a_static)) then
1415 call this%hm_base%allocate_field(mesh, field_vector_potential, .false.)
1416 !ep%a_static contains 1/c A(r)
1417 !$omp parallel do schedule(static) private(idir)
1418 do ip = 1, mesh%np
1419 do idir = 1, space%dim
1420 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + this%ep%a_static(ip, idir)
1421 end do
1422 end do
1423 end if
1424
1425 !The electric field is added to the KS potential
1426 call this%hm_base%accel_copy_pot(mesh)
1427
1428 ! and the static magnetic field
1429 if (allocated(this%ep%b_field)) then
1430 call this%hm_base%allocate_field(mesh, field_uniform_magnetic_field, .false.)
1431 do idir = 1, 3
1432 this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
1433 end do
1434 end if
1435
1436 call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
1437
1438 call hamiltonian_elec_update_pot(this, mesh)
1439
1440 call build_phase()
1441
1442 call profiling_out("HAMILTONIAN_ELEC_UPDATE_EXT_POT")
1444
1445 contains
1446
1447 subroutine build_phase()
1448 integer :: ik, imat, nmat, max_npoints, offset, iphase, nphase
1449
1451
1452 if ((.not. this%kpoints%gamma_only()) .or. allocated(this%hm_base%uniform_vector_potential)) then
1453
1454 call profiling_in('UPDATE_PHASES')
1455 ! now regenerate the phases for the pseudopotentials
1456 do iatom = 1, this%ep%natoms
1457 call projector_init_phases(this%ep%proj(iatom), space%dim, this%d, this%der%boundaries, this%kpoints, &
1458 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1459 end do
1460
1461 call profiling_out('UPDATE_PHASES')
1462 end if
1463
1464 if (allocated(this%hm_base%uniform_vector_potential)) then
1465 call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
1466 end if
1467
1468 max_npoints = this%vnl%max_npoints
1469 nmat = this%vnl%nprojector_matrices
1470
1471
1472 if (this%phase%is_allocated() .and. allocated(this%vnl%projector_matrices)) then
1473
1474 nphase = 1
1475 if (this%der%boundaries%spiralBC) nphase = 3
1476
1477 if (.not. allocated(this%vnl%projector_phases)) then
1478 safe_allocate(this%vnl%projector_phases(1:max_npoints, nphase, nmat, this%d%kpt%start:this%d%kpt%end))
1479 if (accel_is_enabled()) then
1480 call accel_create_buffer(this%vnl%buff_projector_phases, accel_mem_read_only, &
1481 type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
1482 end if
1483 end if
1484
1485 offset = 0
1486 do ik = this%d%kpt%start, this%d%kpt%end
1487 do imat = 1, this%vnl%nprojector_matrices
1488 iatom = this%vnl%projector_to_atom(imat)
1489 do iphase = 1, nphase
1490 !$omp parallel do schedule(static)
1491 do ip = 1, this%vnl%projector_matrices(imat)%npoints
1492 this%vnl%projector_phases(ip, imat, iphase, ik) = this%ep%proj(iatom)%phase(ip, iphase, ik)
1493 end do
1494
1495 if (accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0) then
1496 call accel_write_buffer(this%vnl%buff_projector_phases, &
1497 this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
1498 offset = offset)
1499 end if
1500 offset = offset + this%vnl%projector_matrices(imat)%npoints
1501 end do
1502 end do
1503 end do
1504
1505 end if
1506
1508 end subroutine build_phase
1509
1511
1512 logical function hamiltonian_elec_needs_current(hm, states_are_real)
1513 type(hamiltonian_elec_t), intent(in) :: hm
1514 logical, intent(in) :: states_are_real
1515
1517
1518 if (hm%self_induced_magnetic) then
1519 if (.not. states_are_real) then
1521 else
1522 message(1) = 'No current density for real states since it is identically zero.'
1523 call messages_warning(1)
1524 end if
1525 end if
1526
1528
1529 ! ---------------------------------------------------------
1530 subroutine zhamiltonian_elec_apply_all(hm, namespace, gr, st, hst)
1531 type(hamiltonian_elec_t), intent(inout) :: hm
1532 type(namespace_t), intent(in) :: namespace
1533 type(grid_t), intent(in) :: gr
1534 type(states_elec_t), intent(inout) :: st
1535 type(states_elec_t), intent(inout) :: hst
1536
1537 integer :: ik, ib, ist
1538 complex(real64), allocatable :: psi(:, :)
1539 complex(real64), allocatable :: psiall(:, :, :, :)
1540
1542
1543 do ik = st%d%kpt%start, st%d%kpt%end
1544 do ib = st%group%block_start, st%group%block_end
1545 call zhamiltonian_elec_apply_batch(hm, namespace, gr, st%group%psib(ib, ik), hst%group%psib(ib, ik))
1546 end do
1547 end do
1548
1549 if (oct_exchange_enabled(hm%oct_exchange)) then
1550
1551 safe_allocate(psiall(gr%np_part, 1:hst%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
1552
1553 call states_elec_get_state(st, gr, psiall)
1554
1555 call oct_exchange_prepare(hm%oct_exchange, gr, psiall, hm%xc, hm%psolver, namespace)
1556
1557 safe_deallocate_a(psiall)
1558
1559 safe_allocate(psi(gr%np_part, 1:hst%d%dim))
1560
1561 do ik = 1, st%nik
1562 do ist = 1, st%nst
1563 call states_elec_get_state(hst, gr, ist, ik, psi)
1564 call oct_exchange_operator(hm%oct_exchange, namespace, gr, psi, ist, ik)
1565 call states_elec_set_state(hst, gr, ist, ik, psi)
1566 end do
1567 end do
1568
1569 safe_deallocate_a(psi)
1570
1571 end if
1572
1574 end subroutine zhamiltonian_elec_apply_all
1575
1576
1577 ! ---------------------------------------------------------
1578
1579 subroutine magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase)
1580 type(hamiltonian_elec_t), intent(in) :: hm
1581 type(namespace_t), intent(in) :: namespace
1582 class(mesh_t), intent(in) :: mesh
1583 complex(real64), contiguous, intent(inout) :: psi(:,:)
1584 complex(real64), contiguous, intent(out) :: hpsi(:,:)
1585 integer, intent(in) :: ik
1586 real(real64), intent(in) :: vmagnus(:, :, :)
1587 logical, optional, intent(in) :: set_phase
1588
1589 complex(real64), allocatable :: auxpsi(:, :), aux2psi(:, :)
1590 integer :: idim, ispin
1591
1592 push_sub(magnus)
1593
1594 ! We will assume, for the moment, no spinors.
1595 if (hm%d%dim /= 1) then
1596 call messages_not_implemented("Magnus with spinors", namespace=namespace)
1597 end if
1598
1599 safe_allocate( auxpsi(1:mesh%np_part, 1:hm%d%dim))
1600 safe_allocate(aux2psi(1:mesh%np, 1:hm%d%dim))
1601
1602 ispin = hm%d%get_spin_index(ik)
1603
1604 ! Compute (T + Vnl)|psi> and store it
1605 call zhamiltonian_elec_apply_single(hm, namespace, mesh, psi, auxpsi, 1, ik, &
1606 terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
1608 ! H|psi> = (T + Vnl)|psi> + Vpsl|psi> + Vmagnus(t2)|psi> + Vborders
1609 do idim = 1, hm%d%dim
1610 call lalg_copy(mesh%np, auxpsi(:, idim), hpsi(:, idim))
1611 hpsi(1:mesh%np, idim) = hpsi(1:mesh%np, idim) + hm%ep%Vpsl(1:mesh%np)*psi(1:mesh%np,idim)
1612 call vborders(mesh, hm, psi(:, idim), hpsi(:, idim))
1613 end do
1614 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + vmagnus(1:mesh%np, ispin, 2)*psi(1:mesh%np, 1)
1615
1616 ! Add first term of the commutator: - i Vmagnus(t1) (T + Vnl) |psi>
1617 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) - m_zi*vmagnus(1:mesh%np, ispin, 1)*auxpsi(1:mesh%np, 1)
1618
1619 ! Add second term of commutator: i (T + Vnl) Vmagnus(t1) |psi>
1620 auxpsi(1:mesh%np, 1) = vmagnus(1:mesh%np, ispin, 1)*psi(1:mesh%np, 1)
1621 call zhamiltonian_elec_apply_single(hm, namespace, mesh, auxpsi, aux2psi, 1, ik, &
1622 terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
1623 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + m_zi*aux2psi(1:mesh%np, 1)
1624
1625 safe_deallocate_a(auxpsi)
1626 safe_deallocate_a(aux2psi)
1627 pop_sub(magnus)
1628 end subroutine magnus
1629
1630 ! ---------------------------------------------------------
1631 subroutine vborders (mesh, hm, psi, hpsi)
1632 class(mesh_t), intent(in) :: mesh
1633 type(hamiltonian_elec_t), intent(in) :: hm
1634 complex(real64), intent(in) :: psi(:)
1635 complex(real64), intent(inout) :: hpsi(:)
1636
1637 integer :: ip
1638
1639 push_sub(vborders)
1640
1641 if (hm%abs_boundaries%abtype == imaginary_absorbing) then
1642 do ip = 1, mesh%np
1643 hpsi(ip) = hpsi(ip) + m_zi*hm%abs_boundaries%mf(ip)*psi(ip)
1644 end do
1645 end if
1646
1647 pop_sub(vborders)
1648 end subroutine vborders
1649
1650 ! ---------------------------------------------------------
1651 logical function hamiltonian_elec_has_kick(hm)
1652 type(hamiltonian_elec_t), intent(in) :: hm
1653
1655
1656 hamiltonian_elec_has_kick = (abs(hm%kick%delta_strength) > m_epsilon)
1657
1659 end function hamiltonian_elec_has_kick
1660
1662 !
1663 subroutine hamiltonian_elec_set_mass(this, namespace, mass)
1664 class(hamiltonian_elec_t) , intent(inout) :: this
1665 type(namespace_t), intent(in) :: namespace
1666 real(real64), intent(in) :: mass
1667
1669
1670 if (parse_is_defined(namespace, 'ParticleMass')) then
1671 message(1) = 'Attempting to redefine a non-unit electron mass'
1672 call messages_fatal(1)
1673 else
1674 this%mass = mass
1675 end if
1676
1678 end subroutine hamiltonian_elec_set_mass
1679
1680 !----------------------------------------------------------
1687 subroutine hamiltonian_elec_copy_and_set_phase(hm, gr, kpt, psib, psib_with_phase)
1688 type(hamiltonian_elec_t), intent(in) :: hm
1689 type(grid_t), intent(in) :: gr
1690 type(distributed_t), intent(in) :: kpt
1691 type(wfs_elec_t), intent(in) :: psib
1692 type(wfs_elec_t), intent(out) :: psib_with_phase
1693
1694 integer :: k_offset, n_boundary_points
1695
1697
1698 call psib%copy_to(psib_with_phase)
1699 if (hm%phase%is_allocated()) then
1700 call hm%phase%apply_to(gr, gr%np, conjugate = .false., psib = psib_with_phase, src = psib, async=.true.)
1701 ! apply phase correction while setting boundary -> memory needs to be
1702 ! accessed only once
1703 k_offset = psib%ik - kpt%start
1704 n_boundary_points = int(gr%np_part - gr%np)
1705 call boundaries_set(gr%der%boundaries, gr, psib_with_phase, phase_correction = hm%phase%phase_corr(:, psib%ik), &
1706 buff_phase_corr = hm%phase%buff_phase_corr, offset=k_offset*n_boundary_points, async=.true.)
1707 else
1708 call psib%copy_data_to(gr%np, psib_with_phase)
1709 call boundaries_set(gr%der%boundaries, gr, psib_with_phase)
1710 end if
1711
1712 call psib_with_phase%do_pack(copy = .true.)
1713
1716
1717 ! ---------------------------------------------------------
1718 subroutine hamiltonian_elec_diagonal (hm, mesh, diag, ik)
1719 type(hamiltonian_elec_t), intent(in) :: hm
1720 class(mesh_t), intent(in) :: mesh
1721 real(real64), contiguous, intent(out) :: diag(:,:)
1722 integer, intent(in) :: ik
1723
1724 integer :: idim, ip, ispin
1725
1726 real(real64), allocatable :: ldiag(:)
1727
1729
1730 safe_allocate(ldiag(1:mesh%np))
1731
1732 diag = m_zero
1733
1734 call derivatives_lapl_diag(hm%der, ldiag)
1735
1736 do idim = 1, hm%d%dim
1737 diag(1:mesh%np, idim) = -m_half/hm%mass*ldiag(1:mesh%np)
1738 end do
1739
1740 select case (hm%d%ispin)
1741
1742 case (unpolarized, spin_polarized)
1743 ispin = hm%d%get_spin_index(ik)
1744 diag(1:mesh%np, 1) = diag(1:mesh%np, 1) + hm%ep%vpsl(1:mesh%np)
1745
1746 case (spinors)
1747 do ip = 1, mesh%np
1748 diag(ip, 1) = diag(ip, 1) + hm%ep%vpsl(ip)
1749 diag(ip, 2) = diag(ip, 2) + hm%ep%vpsl(ip)
1750 end do
1751
1752 end select
1753
1754 call hm%ks_pot%add_vhxc(diag)
1755
1757 end subroutine hamiltonian_elec_diagonal
1759
1760
1761#include "undef.F90"
1762#include "real.F90"
1763#include "hamiltonian_elec_inc.F90"
1764
1765#include "undef.F90"
1766#include "complex.F90"
1767#include "hamiltonian_elec_inc.F90"
1768
1769end module hamiltonian_elec_oct_m
1770
1771!! Local Variables:
1772!! mode: f90
1773!! coding: utf-8
1774!! 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:429
subroutine, public accel_finish()
Definition: accel.F90:985
pure logical function, public accel_is_enabled()
Definition: accel.F90:419
integer, parameter, public accel_mem_read_only
Definition: accel.F90:196
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:193
real(real64), parameter, public m_zero
Definition: global.F90:191
integer, parameter, public rdmft
Definition: global.F90:237
integer, parameter, public hartree_fock
Definition: global.F90:237
integer, parameter, public independent_particles
Theory level.
Definition: global.F90:237
integer, parameter, public generalized_kohn_sham_dft
Definition: global.F90:237
integer, parameter, public kohn_sham_dft
Definition: global.F90:237
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
Definition: global.F90:229
real(real64), parameter, public m_one
Definition: global.F90:192
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_init(this, namespace, space, level, gr, ions, st, mc, kpoints)
Definition: lda_u.F90:284
subroutine, public lda_u_update_basis(this, space, gr, ions, st, psolver, namespace, kpoints, has_phase)
Definition: lda_u.F90:697
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:824
subroutine, public lda_u_end(this)
Definition: lda_u.F90:655
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:1091
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1063
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:455
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:240
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:599
integer pure function, public symmetries_number(this)
Definition: symmetries.F90:553
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:641
logical pure function, public family_is_hybrid(xcs)
Returns true if the functional is an hybrid functional.
Definition: xc.F90:656
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)