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