49 use,
intrinsic :: iso_fortran_env
91 use xc_functional_oct_m
94 use xc_functional_oct_m
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(:,:)
145 type(derivatives_t),
pointer,
private :: der
147 type(nonlocal_pseudopotential_t) :: vnl
149 type(ions_t),
pointer :: ions
150 real(real64) :: exx_coef
152 type(poisson_t) :: psolver
155 logical :: self_induced_magnetic
156 real(real64),
allocatable :: a_ind(:, :)
157 real(real64),
allocatable :: b_ind(:, :)
159 integer :: theory_level
160 type(xc_t),
pointer :: xc
161 type(xc_photons_t),
pointer :: xc_photons
167 logical,
private :: adjoint
170 real(real64),
private :: mass
173 logical,
private :: inh_term
174 type(states_elec_t) :: inh_st
178 type(oct_exchange_t) :: oct_exchange
180 type(scissor_t) :: scissor
182 real(real64) :: current_time
183 logical,
private :: is_applied_packed
186 type(lda_u_t) :: lda_u
187 integer :: lda_u_level
189 logical,
public :: time_zero
191 type(exchange_operator_t),
public :: exxop
193 type(kpoints_t),
pointer,
public :: kpoints => null()
195 type(partner_list_t) :: external_potentials
196 real(real64),
allocatable,
public :: v_ext_pot(:)
197 real(real64),
allocatable,
public :: v_static(:)
199 type(ion_electron_local_potential_t) :: v_ie_loc
202 type(magnetic_constrain_t) :: magnetic_constrain
208 type(mxll_coupling_t) :: mxll
209 type(zora_t),
pointer :: zora => null()
223 integer,
public,
parameter :: &
232 mc, kpoints, need_exchange, xc_photons)
236 type(
grid_t),
target,
intent(inout) :: gr
237 type(
ions_t),
target,
intent(inout) :: ions
240 integer,
intent(in) :: theory_level
241 type(xc_t),
target,
intent(in) :: xc
244 logical,
optional,
intent(in) :: need_exchange
248 logical :: need_exchange_
249 real(real64) :: rashba_coupling
256 hm%theory_level = theory_level
259 hm%kpoints => kpoints
285 if (space%dim /= 2)
then
286 write(
message(1),
'(a)')
'Rashba spin-orbit coupling can only be used for two-dimensional systems.'
292 call hm%hm_base%init(hm%d%nspin, hm%mass, rashba_coupling)
295 assert(
associated(gr%der%lapl))
296 hm%hm_base%kinetic => gr%der%lapl
298 safe_allocate(hm%energy)
305 if(
present(xc_photons))
then
306 hm%xc_photons => xc_photons
308 hm%xc_photons => null()
314 call hm%ks_pot%init(gr%np, gr%np_part, hm%d%nspin, hm%theory_level, family_is_mgga_with_exc(hm%xc))
317 call poisson_init(hm%psolver, namespace, space, gr%der, mc, gr%stencil, st%qtot)
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)
323 hm%zora =>
zora_t(namespace, hm%der, hm%d, hm%ep, hm%mass)
326 call hm%v_ie_loc%init(gr, hm%psolver, hm%ions, namespace)
328 call hm%nlcc%init(gr, hm%ions)
329 safe_allocate(st%rho_core(1:gr%np))
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))
373 hm%inh_term = .false.
378 call hm%phase%init(gr, hm%d%kpt, hm%kpoints, st%d, space)
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())
408 if(hm%lda_u_level /=
dft_u_none .and. hm%phase%is_allocated())
then
423 call parse_variable(namespace,
'HamiltonianApplyPacked', .
true., hm%is_applied_packed)
425 if (hm%theory_level ==
hartree_fock .and. st%parallel_in_states)
then
430 .and. st%parallel_in_states)
then
451 .or. hm%theory_level ==
rdmft .or. need_exchange_ .or. &
453 .or. hm%xc%functional(
func_x,1)%id == xc_oep_x_slater &
454 .or.
bitand(hm%xc%family, xc_family_oep) /= 0)
then
456 if (hm%xc%functional(
func_x,1)%id == xc_oep_x_slater)
then
459 else if (
bitand(hm%xc%family, xc_family_oep) /= 0 .or. hm%theory_level ==
rdmft)
then
461 hm%kpoints, hm%xc%cam_omega, hm%xc%cam_alpha, hm%xc%cam_beta)
470 if (gr%use_curvilinear)
then
472 hm%is_applied_packed = .false.
473 call messages_write(
'Cannot use CUDA or OpenCL as curvilinear coordinates are used.')
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.')
507 if (
associated(hm%xc_photons))
then
508 if (hm%xc_photons%wants_to_renormalize_mass())
then
510 call hm%set_mass(namespace, hm%xc_photons%get_renormalized_mass())
523 class(*),
pointer :: potential
528 safe_allocate(hm%v_ext_pot(1:gr%np))
529 hm%v_ext_pot(1:gr%np) =
m_zero
531 call iter%start(hm%external_potentials)
532 do while (iter%has_next())
533 potential => iter%get_next()
534 select type (potential)
537 call potential%allocate_memory(gr)
538 call potential%calculate(namespace, gr, hm%psolver)
541 select case (potential%type)
546 if (.not.
allocated(hm%ep%b_field))
then
547 safe_allocate(hm%ep%b_field(1:3))
548 hm%ep%b_field(1:3) =
m_zero
550 hm%ep%b_field(1:3) = hm%ep%b_field(1:3) + potential%b_field(1:3)
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
556 call lalg_axpy(gr%np, space%dim,
m_one, potential%a_static, hm%ep%a_static)
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
563 hm%ep%e_field(1:space%dim) = hm%ep%e_field(1:space%dim) + potential%e_field(1:space%dim)
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
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
579 if (hm%kpoints%use_symmetries)
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."
591 call potential%deallocate_memory()
604 class(*),
pointer :: potential
608 call iter%start(hm%external_potentials)
609 do while (iter%has_next())
610 potential => iter%get_next()
611 select type (potential)
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."
632 logical :: external_potentials_present
633 logical :: kick_present
637 if (
allocated(hm%ep%e_field) .and. space%is_periodic() .and. .not.
list_has_gauge_field(ext_partners))
then
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))
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
674 call hm%hm_base%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)
685 safe_deallocate_p(hm%zora)
704 safe_deallocate_a(hm%energy)
706 if (hm%pcm%run_pcm)
call pcm_end(hm%pcm)
708 call hm%v_ie_loc%end()
711 call iter%start(hm%external_potentials)
712 do while (iter%has_next())
713 potential => iter%get_next()
714 safe_deallocate_p(potential)
716 call hm%external_potentials%empty()
717 safe_deallocate_a(hm%v_static)
743 real(real64),
intent(in) :: delta(:)
744 real(real64),
intent(in) :: emin
755 hm%spectral_middle_point = (emax + emin) /
m_two
756 hm%spectral_half_span = (emax - emin) /
m_two
791 if (hm%inh_term)
then
793 hm%inh_term = .false.
805 if (.not. hm%adjoint)
then
808 hm%abs_boundaries%mf = -hm%abs_boundaries%mf
825 hm%abs_boundaries%mf = -hm%abs_boundaries%mf
837 class(
mesh_t),
intent(in) :: mesh
839 class(
space_t),
intent(in) :: space
841 real(real64),
optional,
intent(in) :: time
843 integer :: ispin, ip, idir, iatom, ilaser
844 real(real64) :: aa(space%dim), time_
845 real(real64),
allocatable :: vp(:,:)
848 real(real64) :: am(space%dim)
853 this%current_time =
m_zero
854 if (
present(time)) this%current_time = time
859 call this%hm_base%clear(mesh%np)
866 if (
present(time) .or. this%time_zero)
then
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
875 this%hm_base%potential(:, ispin), time_)
881 safe_allocate(vp(1:mesh%np, 1:space%dim))
882 vp(1:mesh%np, 1:space%dim) =
m_zero
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
891 call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
892 safe_deallocate_a(vp)
898 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/
p_c
903 assert(
allocated(this%hm_base%uniform_vector_potential))
905 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - am/
p_c
911 if (
associated(gfield))
then
914 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/
p_c
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)
924 if (
associated(this%xc_photons))
then
925 if(this%xc_photons%lpfmf .and.
allocated(this%xc_photons%mf_vector_potential))
then
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
936 if (
allocated(this%ep%a_static))
then
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)
951 call this%hm_base%accel_copy_pot(mesh)
954 if (
allocated(this%ep%b_field))
then
957 this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
962 call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
967 if (this%mxll%test_equad)
then
979 integer :: ik, imat, nmat, max_npoints, offset
981 integer :: iphase, nphase
985 if ((.not. this%kpoints%gamma_only()) .or.
allocated(this%hm_base%uniform_vector_potential))
then
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)
997 if (
allocated(this%hm_base%uniform_vector_potential))
then
999 call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
1004 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1008 max_npoints = this%vnl%max_npoints
1009 nmat = this%vnl%nprojector_matrices
1012 if (this%phase%is_allocated() .and.
allocated(this%vnl%projector_matrices))
then
1015 if (this%der%boundaries%spiralBC) nphase = 3
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))
1021 type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
1024 this%vnl%nphase = nphase
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
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)
1038 if (
accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0)
then
1040 this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
1041 offset = offset, async=.
true.)
1043 offset = offset + this%vnl%projector_matrices(imat)%npoints
1064 class(
mesh_t),
intent(in) :: mesh
1065 logical,
optional,
intent(in) :: accumulate
1067 integer :: ispin, ip
1074 do ispin = 1, this%d%nspin
1077 this%hm_base%potential(ip, ispin) =
m_zero
1084 do ispin = 1, this%d%nspin
1085 if (ispin <= 2)
then
1089 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%ep%vpsl(ip) + this%v_ext_pot(ip)
1093 if (this%pcm%run_pcm)
then
1094 if (this%pcm%solute)
then
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)
1102 if (this%pcm%localf)
then
1105 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
1106 this%pcm%v_ext_rs(ip)
1116 this%hm_base%Impotential(ip, ispin) = this%hm_base%Impotential(ip, ispin) + this%abs_boundaries%mf(ip)
1127 call this%zora%update(this%der, this%hm_base%potential)
1131 do ispin = 1, this%d%nspin
1133 if (
allocated(this%hm_base%zeeman_pot))
then
1136 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%hm_base%zeeman_pot(ip, ispin)
1142 if (this%mxll%test_equad)
then
1145 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%mxll%e_quadrupole_pot(ip)
1154 call this%ks_pot%add_vhxc(this%hm_base%potential)
1156 call this%hm_base%accel_copy_pot(mesh)
1166 type(
grid_t),
intent(in) :: gr
1167 type(
ions_t),
target,
intent(inout) :: ions
1170 real(real64),
optional,
intent(in) :: time
1175 call epot_generate(this%ep, namespace, gr, this%ions, this%d)
1180 if (
allocated(this%ep%e_field) .and. space%periodic_dim < space%dim)
then
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()
1191 call lalg_axpy(gr%np,
m_one, this%v_ie_loc%potential(:,1), this%ep%vpsl)
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)
1202 call this%vnl%build(space, gr, this%ep)
1207 if (this%ep%non_local .and. .not. this%vnl%apply_projector_matrices)
then
1209 call messages_write(
'Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.')
1212 call messages_write(
'Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.',&
1214 call messages_write(
'Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
1221 if (this%pcm%run_pcm)
then
1224 if (this%pcm%solute)
then
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))
1237 call lda_u_update_basis(this%lda_u, space, gr, ions, st, this%psolver, namespace, this%kpoints, &
1238 this%phase%is_allocated())
1248 time = this%current_time
1256 apply = this%is_applied_packed
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(:,:)
1275 real(real64),
allocatable :: vlocal(:)
1278 safe_allocate(vlocal(1:mesh%np_part))
1280 call epot_local_potential(hm%ep, namespace, space, latt, mesh, species, pos, ia, vlocal)
1282 do idim = 1, hm%d%dim
1283 vpsi(1:mesh%np, idim) = vlocal(1:mesh%np) * psi(1:mesh%np, idim)
1286 safe_deallocate_a(vlocal)
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)
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
1311 call profiling_in(
"HAMILTONIAN_ELEC_UPDATE_EXT_POT")
1313 this%current_time = m_zero
1314 this%current_time = time(1)
1317 call this%hm_base%clear(mesh%np)
1320 call this%hm_base%allocate_field(mesh, field_potential, &
1321 complex_potential = this%abs_boundaries%abtype == imaginary_absorbing)
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
1334 call laser_potential(lasers%lasers(ilaser), mesh, velectric, time_)
1337 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + mu(itime) * velectric(ip)
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.)
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
1350 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) &
1351 - mu(itime) * vp(ip, idir)/p_c
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.)
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
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
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)
1388 if (
allocated(this%ep%a_static))
then
1389 call this%hm_base%allocate_field(mesh, field_vector_potential, .false.)
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)
1400 call this%hm_base%accel_copy_pot(mesh)
1403 if (
allocated(this%ep%b_field))
then
1404 call this%hm_base%allocate_field(mesh, field_uniform_magnetic_field, .false.)
1406 this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
1410 call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
1416 call profiling_out(
"HAMILTONIAN_ELEC_UPDATE_EXT_POT")
1422 integer :: ik, imat, nmat, max_npoints, offset, iphase, nphase
1426 if ((.not. this%kpoints%gamma_only()) .or.
allocated(this%hm_base%uniform_vector_potential))
then
1428 call profiling_in(
'UPDATE_PHASES')
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)
1435 call profiling_out(
'UPDATE_PHASES')
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)
1442 max_npoints = this%vnl%max_npoints
1443 nmat = this%vnl%nprojector_matrices
1446 if (this%phase%is_allocated() .and.
allocated(this%vnl%projector_matrices))
then
1449 if (this%der%boundaries%spiralBC) nphase = 3
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)
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
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)
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), &
1474 offset = offset + this%vnl%projector_matrices(imat)%npoints
1488 logical,
intent(in) :: states_are_real
1492 if (hm%self_induced_magnetic)
then
1493 if (.not. states_are_real)
then
1496 message(1) =
'No current density for real states since it is identically zero.'
1497 call messages_warning(1)
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
1511 integer :: ik, ib, ist
1512 complex(real64),
allocatable :: psi(:, :)
1513 complex(real64),
allocatable :: psiall(:, :, :, :)
1517 do ik = st%d%kpt%start, st%d%kpt%end
1518 do ib = st%group%block_start, st%group%block_end
1523 if (oct_exchange_enabled(hm%oct_exchange))
then
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))
1527 call states_elec_get_state(st, mesh, psiall)
1529 call oct_exchange_prepare(hm%oct_exchange, mesh, psiall, hm%xc, hm%psolver, namespace)
1531 safe_deallocate_a(psiall)
1533 safe_allocate(psi(mesh%np_part, 1:hst%d%dim))
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)
1543 safe_deallocate_a(psi)
1553 subroutine magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase)
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
1563 complex(real64),
allocatable :: auxpsi(:, :), aux2psi(:, :)
1564 integer :: idim, ispin
1569 if (hm%d%dim /= 1)
then
1570 call messages_not_implemented(
"Magnus with spinors", namespace=namespace)
1573 safe_allocate( auxpsi(1:mesh%np_part, 1:hm%d%dim))
1574 safe_allocate(aux2psi(1:mesh%np, 1:hm%d%dim))
1576 ispin = hm%d%get_spin_index(ik)
1580 terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
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))
1588 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + vmagnus(1:mesh%np, ispin, 2)*psi(1:mesh%np, 1)
1591 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) - m_zi*vmagnus(1:mesh%np, ispin, 1)*auxpsi(1:mesh%np, 1)
1594 auxpsi(1:mesh%np, 1) = vmagnus(1:mesh%np, ispin, 1)*psi(1:mesh%np, 1)
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)
1599 safe_deallocate_a(auxpsi)
1600 safe_deallocate_a(aux2psi)
1605 subroutine vborders (mesh, hm, psi, hpsi)
1606 class(mesh_t),
intent(in) :: mesh
1608 complex(real64),
intent(in) :: psi(:)
1609 complex(real64),
intent(inout) :: hpsi(:)
1615 if (hm%abs_boundaries%abtype == imaginary_absorbing)
then
1617 hpsi(ip) = hpsi(ip) + m_zi*hm%abs_boundaries%mf(ip)*psi(ip)
1639 type(namespace_t),
intent(in) :: namespace
1640 real(real64),
intent(in) :: mass
1644 if (parse_is_defined(namespace,
'ParticleMass'))
then
1645 message(1) =
'Attempting to redefine a non-unit electron mass'
1646 call messages_fatal(1)
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
1668 integer :: k_offset, n_boundary_points
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.)
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.)
1682 call psib%copy_data_to(gr%np, psib_with_phase)
1683 call boundaries_set(gr%der%boundaries, gr, psib_with_phase)
1686 call psib_with_phase%do_pack(copy = .
true.)
1694 class(mesh_t),
intent(in) :: mesh
1695 real(real64),
intent(out) :: diag(:,:)
1696 integer,
intent(in) :: ik
1698 integer :: idim, ip, ispin
1700 real(real64),
allocatable :: ldiag(:)
1704 safe_allocate(ldiag(1:mesh%np))
1708 call derivatives_lapl_diag(hm%der, ldiag)
1710 do idim = 1, hm%d%dim
1711 diag(1:mesh%np, idim) = -m_half/hm%mass*ldiag(1:mesh%np)
1714 select case (hm%d%ispin)
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)
1722 diag(ip, 1) = diag(ip, 1) + hm%ep%vpsl(ip)
1723 diag(ip, 2) = diag(ip, 2) + hm%ep%vpsl(ip)
1728 call hm%ks_pot%add_vhxc(diag)
1737#include "hamiltonian_elec_inc.F90"
1740#include "complex.F90"
1741#include "hamiltonian_elec_inc.F90"
subroutine build_external_potentials()
subroutine external_potentials_checks()
subroutine build_interactions()
constant times a vector plus a vector
Copies a vector x, to a vector y.
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()
subroutine, public accel_finish()
pure logical function, public accel_is_enabled()
integer, parameter, public accel_mem_read_only
This module implements batches of mesh functions.
This module implements common operations on batches of mesh functions.
Module implementing boundary conditions in Octopus.
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)
integer, parameter, public scalar_relativistic_zora
subroutine, public epot_end(ep)
integer, parameter, public fully_relativistic_zora
subroutine, public epot_init(ep, namespace, gr, ions, psolver, ispin, xc_family, kpoints)
subroutine, public epot_generate(ep, namespace, mesh, ions, st_d)
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
real(real64), parameter, public m_zero
real(real64), parameter, public m_epsilon
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
real(real64), parameter, public m_one
This module implements the underlying real-space grid.
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.
integer, parameter, public kick_magnon_mode
subroutine, public kick_end(kick)
subroutine, public kick_init(kick, namespace, space, kpoints, nspin)
pure integer function, public kick_get_type(kick)
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)
subroutine, public lasers_nondipole_laser_field_step(this, field, time)
Retrieves the NDSFA vector_potential correction. The nondipole field is obtained for consecutive time...
logical function, public lasers_with_nondipole_field(lasers)
Check if a nondipole SFA correction should be computed for the given laser.
integer, parameter, public e_field_electric
integer, parameter, public e_field_vector_potential
subroutine, public laser_potential(laser, mesh, pot, time)
integer, parameter, public e_field_scalar_potential
integer pure elemental function, public laser_kind(laser)
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...
integer, parameter, public e_field_magnetic
integer, parameter, public dft_u_none
subroutine, public lda_u_update_basis(this, space, gr, ions, st, psolver, namespace, kpoints, has_phase)
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.
subroutine, public lda_u_end(this)
subroutine, public lda_u_init(this, namespace, space, level, gr, ions, st, mc, kpoints, has_phase)
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.
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_not_implemented(feature, namespace)
subroutine, public messages_warning(no_lines, all_nodes, namespace)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_experimental(name, namespace)
This module handles the communicators for the various parallelization strategies.
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
Some general things and nomenclature:
logical function, public parse_is_defined(namespace, name)
subroutine, public pcm_calc_pot_rs(pcm, mesh, psolver, ions, v_h, v_ext, kick, time_present, kick_time)
subroutine, public pcm_end(pcm)
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...
subroutine, public poisson_init(this, namespace, space, der, mc, stencil, qtot, label, solver, verbose, force_serial, force_cmplx)
subroutine, public poisson_end(this)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
subroutine, public projector_init_phases(this, dim, std, bnd, kpoints, vec_pot, vec_pot_var)
subroutine, public scissor_end(this)
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)
integer pure function, public symmetries_number(this)
type(type_t), public type_cmplx
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
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.
This module implements the ZORA terms for the Hamoiltonian.
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
The abstract Hamiltonian class defines a skeleton for specific implementations.
abstract class for general interaction partners
iterator for the list of partners
Describes mesh distribution to nodes.
Stores all communicators and groups.
The states_elec_t class contains all electronic wave functions.
This class described the 'photon-exchange' electron-photon xc functionals, based on QEDFT.
This class is responsible for calculating and applying the ZORA.