49 use,
intrinsic :: iso_fortran_env
137 type(space_t),
private :: space
138 type(states_elec_dim_t) :: d
139 type(hamiltonian_elec_base_t) :: hm_base
140 type(phase_t) :: phase
141 type(energy_t),
allocatable :: energy
142 type(absorbing_boundaries_t) :: abs_boundaries
143 type(ks_potential_t) :: ks_pot
144 real(real64),
allocatable :: vberry(:,:)
146 type(derivatives_t),
pointer,
private :: der
148 type(nonlocal_pseudopotential_t) :: vnl
150 type(ions_t),
pointer :: ions
151 real(real64) :: exx_coef
153 type(poisson_t) :: psolver
156 logical :: self_induced_magnetic
157 real(real64),
allocatable :: a_ind(:, :)
158 real(real64),
allocatable :: b_ind(:, :)
160 integer :: theory_level
161 type(xc_t),
pointer :: xc
162 type(xc_photons_t),
pointer :: xc_photons
168 logical,
private :: adjoint
171 real(real64),
private :: mass
174 logical,
private :: inh_term
175 type(states_elec_t) :: inh_st
179 type(oct_exchange_t) :: oct_exchange
181 type(scissor_t) :: scissor
183 real(real64) :: current_time
184 logical,
private :: is_applied_packed
187 type(lda_u_t) :: lda_u
188 integer :: lda_u_level
190 logical,
public :: time_zero
192 type(exchange_operator_t),
public :: exxop
194 type(kpoints_t),
pointer,
public :: kpoints => null()
196 type(partner_list_t) :: external_potentials
197 real(real64),
allocatable,
public :: v_ext_pot(:)
198 real(real64),
allocatable,
public :: v_static(:)
200 type(ion_electron_local_potential_t) :: v_ie_loc
203 type(magnetic_constrain_t) :: magnetic_constrain
209 type(mxll_coupling_t) :: mxll
210 type(zora_t),
pointer :: zora => null()
225 integer,
public,
parameter :: &
234 mc, kpoints, need_exchange, xc_photons)
238 type(
grid_t),
target,
intent(inout) :: gr
239 type(
ions_t),
target,
intent(inout) :: ions
242 integer,
intent(in) :: theory_level
243 type(
xc_t),
target,
intent(in) :: xc
246 logical,
optional,
intent(in) :: need_exchange
247 type(
xc_photons_t),
optional,
target,
intent(in) :: xc_photons
250 logical :: need_exchange_
251 real(real64) :: rashba_coupling
258 hm%theory_level = theory_level
261 hm%kpoints => kpoints
287 if (space%dim /= 2)
then
288 write(
message(1),
'(a)')
'Rashba spin-orbit coupling can only be used for two-dimensional systems.'
294 call hm%hm_base%init(hm%d%nspin, hm%mass, rashba_coupling)
297 assert(
associated(gr%der%lapl))
298 hm%hm_base%kinetic => gr%der%lapl
300 safe_allocate(hm%energy)
307 if(
present(xc_photons))
then
308 hm%xc_photons => xc_photons
310 hm%xc_photons => null()
319 call poisson_init(hm%psolver, namespace, space, gr%der, mc, gr%stencil, st%qtot)
322 call epot_init(hm%ep, namespace, gr, hm%ions, hm%psolver, hm%d%ispin, hm%xc%family, hm%kpoints)
323 call kick_init(hm%kick, namespace, space, hm%kpoints, hm%d%ispin)
325 hm%zora =>
zora_t(namespace, hm%der, hm%d, hm%ep, hm%mass)
328 call hm%v_ie_loc%init(gr, hm%psolver, hm%ions, namespace)
330 call hm%nlcc%init(gr, hm%ions)
331 safe_allocate(st%rho_core(1:gr%np))
364 call parse_variable(namespace,
'CalculateSelfInducedMagneticField', .false., hm%self_induced_magnetic)
365 if (hm%self_induced_magnetic)
then
366 safe_allocate(hm%a_ind(1:gr%np_part, 1:space%dim))
367 safe_allocate(hm%b_ind(1:gr%np_part, 1:space%dim))
375 hm%inh_term = .false.
380 call hm%phase%init(gr, hm%d%kpt, hm%kpoints, st%d, space)
400 call lda_u_init(hm%lda_u, namespace, space, hm%lda_u_level, gr, ions, st, mc, hm%kpoints)
409 if(hm%lda_u_level /=
dft_u_none .and. hm%phase%is_allocated())
then
424 call parse_variable(namespace,
'HamiltonianApplyPacked', .
true., hm%is_applied_packed)
426 if (hm%theory_level ==
hartree_fock .and. st%parallel_in_states)
then
431 .and. st%parallel_in_states)
then
450 if (hm%xc%compute_exchange(hm%theory_level) .or. need_exchange_)
then
455 else if (
bitand(hm%xc%family, xc_family_oep) /= 0 .or. hm%theory_level ==
rdmft)
then
457 hm%kpoints, hm%xc%cam)
458 if (hm%theory_level ==
rdmft) hm%exxop%useACE = .false.
467 if (gr%use_curvilinear)
then
469 hm%is_applied_packed = .false.
470 call messages_write(
'Cannot use CUDA or OpenCL as curvilinear coordinates are used.')
473 call messages_write(
'Cannot use CUDA or OpenCL as curvilinear coordinates are used.', new_line = .
true.)
474 call messages_write(
'Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
504 if (
associated(hm%xc_photons))
then
505 if (hm%xc_photons%wants_to_renormalize_mass())
then
507 call hm%set_mass(namespace, hm%xc_photons%get_renormalized_mass())
511 if (hm%xc%compute_exchange(hm%theory_level) .or. need_exchange_)
call hm%exxop%write_info(namespace)
521 class(*),
pointer :: potential
526 safe_allocate(hm%v_ext_pot(1:gr%np))
527 hm%v_ext_pot(1:gr%np) =
m_zero
529 call iter%start(hm%external_potentials)
530 do while (iter%has_next())
531 potential => iter%get_next()
532 select type (potential)
535 call potential%allocate_memory(gr)
536 call potential%calculate(namespace, gr, hm%psolver)
539 select case (potential%type)
545 message(1) =
"Cannot use static magnetic field with real wavefunctions"
549 if (.not.
allocated(hm%ep%b_field))
then
550 safe_allocate(hm%ep%b_field(1:3))
551 hm%ep%b_field(1:3) =
m_zero
553 hm%ep%b_field(1:3) = hm%ep%b_field(1:3) + potential%b_field(1:3)
555 if (.not.
allocated(hm%ep%a_static))
then
556 safe_allocate(hm%ep%a_static(1:gr%np, 1:space%dim))
557 hm%ep%a_static(1:gr%np, 1:space%dim) =
m_zero
559 call lalg_axpy(gr%np, space%dim,
m_one, potential%a_static, hm%ep%a_static)
562 if (.not.
allocated(hm%ep%e_field))
then
563 safe_allocate(hm%ep%e_field(1:space%dim))
564 hm%ep%e_field(1:space%dim) =
m_zero
566 hm%ep%e_field(1:space%dim) = hm%ep%e_field(1:space%dim) + potential%e_field(1:space%dim)
569 if (space%periodic_dim < space%dim)
then
570 if (.not.
allocated(hm%v_static))
then
571 safe_allocate(hm%v_static(1:gr%np))
572 hm%v_static(1:gr%np) =
m_zero
574 if (.not.
allocated(hm%ep%v_ext))
then
575 safe_allocate(hm%ep%v_ext(1:gr%np_part))
576 hm%ep%v_ext(1:gr%np_part) =
m_zero
582 if (hm%kpoints%use_symmetries)
then
586 message(1) =
"The StaticElectricField breaks (at least) one of the symmetries used to reduce the k-points."
587 message(2) =
"Set SymmetryBreakDir equal to StaticElectricField."
594 call potential%deallocate_memory()
607 class(*),
pointer :: potential
611 call iter%start(hm%external_potentials)
612 do while (iter%has_next())
613 potential => iter%get_next()
614 select type (potential)
618 message(1) =
"Applying StaticElectricField in a periodic direction is only accurate for large supercells."
619 message(2) =
"Single-point Berry phase is not appropriate when k-point sampling is needed."
635 logical :: external_potentials_present
636 logical :: kick_present
640 if (
allocated(hm%ep%e_field) .and. space%is_periodic() .and. .not.
list_has_gauge_field(ext_partners))
then
643 if (any(abs(hm%ep%e_field(1:space%periodic_dim)) >
m_epsilon))
then
644 safe_allocate(hm%vberry(1:gr%np, 1:hm%d%nspin))
655 call pcm_init(hm%pcm, namespace, space, ions, gr, st%qtot, st%val_charge, external_potentials_present, kick_present)
656 if (hm%pcm%run_pcm)
then
679 call hm%hm_base%end()
685 safe_deallocate_a(hm%vberry)
686 safe_deallocate_a(hm%a_ind)
687 safe_deallocate_a(hm%b_ind)
688 safe_deallocate_a(hm%v_ext_pot)
690 safe_deallocate_p(hm%zora)
709 safe_deallocate_a(hm%energy)
711 if (hm%pcm%run_pcm)
call pcm_end(hm%pcm)
713 call hm%v_ie_loc%end()
716 call iter%start(hm%external_potentials)
717 do while (iter%has_next())
718 potential => iter%get_next()
719 safe_deallocate_p(potential)
721 call hm%external_potentials%empty()
722 safe_deallocate_a(hm%v_static)
748 integer,
intent(in) :: terms
769 real(real64),
intent(in) ::
delta(:)
770 real(real64),
intent(in) :: emin
781 hm%spectral_middle_point = (emax + emin) /
m_two
782 hm%spectral_half_span = (emax - emin) /
m_two
817 if (hm%inh_term)
then
819 hm%inh_term = .false.
831 if (.not. hm%adjoint)
then
834 hm%abs_boundaries%mf = -hm%abs_boundaries%mf
851 hm%abs_boundaries%mf = -hm%abs_boundaries%mf
863 class(
mesh_t),
intent(in) :: mesh
865 class(
space_t),
intent(in) :: space
867 real(real64),
optional,
intent(in) :: time
869 integer :: ispin, ip, idir, iatom, ilaser
870 real(real64) :: aa(space%dim), time_
871 real(real64),
allocatable :: vp(:,:)
874 real(real64) :: am(space%dim)
879 this%current_time =
m_zero
880 if (
present(time)) this%current_time = time
885 call this%hm_base%clear(mesh%np)
892 if (
present(time) .or. this%time_zero)
then
895 if(
associated(lasers))
then
896 do ilaser = 1, lasers%no_lasers
897 select case (
laser_kind(lasers%lasers(ilaser)))
899 do ispin = 1, this%d%spin_channels
901 this%hm_base%potential(:, ispin), time_)
907 safe_allocate(vp(1:mesh%np, 1:space%dim))
908 vp(1:mesh%np, 1:space%dim) =
m_zero
912 do idir = 1, space%dim
913 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + vp(ip, idir)/
p_c
917 call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
918 safe_deallocate_a(vp)
924 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/
p_c
929 assert(
allocated(this%hm_base%uniform_vector_potential))
931 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - am/
p_c
937 if (
associated(gfield))
then
940 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/
p_c
944 if (
allocated(this%ep%e_field) .and.
associated(gfield))
then
945 this%hm_base%uniform_vector_potential(1:space%periodic_dim) = &
946 this%hm_base%uniform_vector_potential(1:space%periodic_dim) - time_*this%ep%e_field(1:space%periodic_dim)
950 if (
associated(this%xc_photons))
then
951 if(this%xc_photons%lpfmf .and.
allocated(this%xc_photons%mf_vector_potential))
then
954 this%hm_base%uniform_vector_potential(1:space%dim) = &
955 this%hm_base%uniform_vector_potential(1:space%dim) - this%xc_photons%mf_vector_potential(1:space%dim)/
p_c
962 if (
allocated(this%ep%a_static))
then
967 do idir = 1, space%dim
968 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + this%ep%a_static(ip, idir)
977 call this%hm_base%accel_copy_pot(mesh)
980 if (
allocated(this%ep%b_field))
then
983 this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
988 call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
993 if (this%mxll%test_equad)
then
1005 integer :: ik, imat, nmat, max_npoints, offset
1007 integer :: iphase, nphase
1011 if ((.not. this%kpoints%gamma_only()) .or.
allocated(this%hm_base%uniform_vector_potential))
then
1015 do iatom = 1, this%ep%natoms
1016 call projector_init_phases(this%ep%proj(iatom), space%dim, this%d, this%der%boundaries, this%kpoints, &
1017 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1023 if (
allocated(this%hm_base%uniform_vector_potential))
then
1025 call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
1030 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1034 max_npoints = this%vnl%max_npoints
1035 nmat = this%vnl%nprojector_matrices
1038 if (this%phase%is_allocated() .and.
allocated(this%vnl%projector_matrices))
then
1041 if (this%der%boundaries%spiralBC) nphase = 3
1043 if (.not.
allocated(this%vnl%projector_phases))
then
1044 safe_allocate(this%vnl%projector_phases(1:max_npoints, 1:nphase, nmat, this%d%kpt%start:this%d%kpt%end))
1047 type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
1050 this%vnl%nphase = nphase
1055 do ik = this%d%kpt%start, this%d%kpt%end
1056 do imat = 1, this%vnl%nprojector_matrices
1057 iatom = this%vnl%projector_to_atom(imat)
1058 do iphase = 1, nphase
1060 do ip = 1, this%vnl%projector_matrices(imat)%npoints
1061 this%vnl%projector_phases(ip, iphase, imat, ik) = this%ep%proj(iatom)%phase(ip, iphase, ik)
1064 if (
accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0)
then
1066 this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
1067 offset = offset, async=.
true.)
1069 offset = offset + this%vnl%projector_matrices(imat)%npoints
1090 class(
mesh_t),
intent(in) :: mesh
1091 logical,
optional,
intent(in) :: accumulate
1093 integer :: ispin, ip
1100 do ispin = 1, this%d%nspin
1103 this%hm_base%potential(ip, ispin) =
m_zero
1110 do ispin = 1, this%d%nspin
1111 if (ispin <= 2)
then
1115 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%ep%vpsl(ip) + this%v_ext_pot(ip)
1119 if (this%pcm%run_pcm)
then
1120 if (this%pcm%solute)
then
1123 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
1124 this%pcm%v_e_rs(ip) + this%pcm%v_n_rs(ip)
1128 if (this%pcm%localf)
then
1131 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + &
1132 this%pcm%v_ext_rs(ip)
1142 this%hm_base%Impotential(ip, ispin) = this%hm_base%Impotential(ip, ispin) + this%abs_boundaries%mf(ip)
1153 call this%zora%update(this%der, this%hm_base%potential)
1157 do ispin = 1, this%d%nspin
1159 if (
allocated(this%hm_base%zeeman_pot))
then
1162 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%hm_base%zeeman_pot(ip, ispin)
1168 if (this%mxll%test_equad)
then
1171 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + this%mxll%e_quadrupole_pot(ip)
1180 call this%ks_pot%add_vhxc(this%hm_base%potential)
1182 call this%hm_base%accel_copy_pot(mesh)
1192 type(
grid_t),
intent(in) :: gr
1193 type(
ions_t),
target,
intent(inout) :: ions
1196 real(real64),
optional,
intent(in) :: time
1201 call epot_generate(this%ep, namespace, gr, this%ions, this%d)
1206 if (
allocated(this%ep%e_field) .and. space%periodic_dim < space%dim)
then
1212 this%v_ie_loc%atoms_dist => ions%atoms_dist
1213 this%v_ie_loc%atom => ions%atom
1214 call this%v_ie_loc%calculate()
1217 call lalg_axpy(gr%np,
m_one, this%v_ie_loc%potential(:,1), this%ep%vpsl)
1221 if (this%ep%nlcc)
then
1222 call this%nlcc%end()
1223 call this%nlcc%init(gr, ions)
1224 call this%nlcc%calculate()
1225 call lalg_copy(gr%np, this%nlcc%density(:,1), st%rho_core)
1228 call this%vnl%build(space, gr, this%ep)
1233 if (this%ep%non_local .and. .not. this%vnl%apply_projector_matrices)
then
1235 call messages_write(
'Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.')
1238 call messages_write(
'Relativistic pseudopotentials have not been fully implemented for CUDA or OpenCL.',&
1240 call messages_write(
'Calculation will not be continued. To force execution, set AllowCPUonly = yes.')
1247 if (this%pcm%run_pcm)
then
1250 if (this%pcm%solute)
then
1257 if (this%pcm%localf .and.
allocated(this%v_static))
then
1258 call pcm_calc_pot_rs(this%pcm, gr, this%psolver, v_ext = this%ep%v_ext(1:gr%np_part))
1263 call lda_u_update_basis(this%lda_u, space, gr, ions, st, this%psolver, namespace, this%kpoints, &
1264 this%phase%is_allocated())
1274 time = this%current_time
1282 apply = this%is_applied_packed
1290 type(namespace_t),
intent(in) :: namespace
1291 class(space_t),
intent(in) :: space
1292 type(lattice_vectors_t),
intent(in) :: latt
1293 class(species_t),
intent(in) :: species
1294 real(real64),
intent(in) :: pos(1:space%dim)
1295 integer,
intent(in) :: ia
1296 class(mesh_t),
intent(in) :: mesh
1297 complex(real64),
intent(in) :: psi(:,:)
1298 complex(real64),
intent(out) :: vpsi(:,:)
1301 real(real64),
allocatable :: vlocal(:)
1304 safe_allocate(vlocal(1:mesh%np_part))
1306 call epot_local_potential(hm%ep, namespace, space, latt, mesh, species, pos, ia, vlocal)
1308 do idim = 1, hm%d%dim
1309 vpsi(1:mesh%np, idim) = vlocal(1:mesh%np) * psi(1:mesh%np, idim)
1312 safe_deallocate_a(vlocal)
1323 class(space_t),
intent(in) :: space
1324 class(mesh_t),
intent(in) :: mesh
1325 type(partner_list_t),
intent(in) :: ext_partners
1326 real(real64),
intent(in) :: time(1:2)
1327 real(real64),
intent(in) :: mu(1:2)
1329 integer :: ispin, ip, idir, iatom, ilaser, itime
1330 real(real64) :: aa(space%dim), time_
1331 real(real64),
allocatable :: vp(:,:)
1332 real(real64),
allocatable :: velectric(:)
1333 type(lasers_t),
pointer :: lasers
1334 type(gauge_field_t),
pointer :: gfield
1337 call profiling_in(
"HAMILTONIAN_ELEC_UPDATE_EXT_POT")
1339 this%current_time = m_zero
1340 this%current_time = time(1)
1343 call this%hm_base%clear(mesh%np)
1346 call this%hm_base%allocate_field(mesh, field_potential, &
1347 complex_potential = this%abs_boundaries%abtype == imaginary_absorbing)
1352 lasers => list_get_lasers(ext_partners)
1353 if(
associated(lasers))
then
1354 do ilaser = 1, lasers%no_lasers
1355 select case (laser_kind(lasers%lasers(ilaser)))
1356 case (e_field_scalar_potential, e_field_electric)
1357 safe_allocate(velectric(1:mesh%np))
1358 do ispin = 1, this%d%spin_channels
1360 call laser_potential(lasers%lasers(ilaser), mesh, velectric, time_)
1363 this%hm_base%potential(ip, ispin) = this%hm_base%potential(ip, ispin) + mu(itime) * velectric(ip)
1366 safe_deallocate_a(velectric)
1367 case (e_field_magnetic)
1368 call this%hm_base%allocate_field(mesh, field_vector_potential + field_uniform_magnetic_field, .false.)
1370 safe_allocate(vp(1:mesh%np, 1:space%dim))
1371 vp(1:mesh%np, 1:space%dim) = m_zero
1372 call laser_vector_potential(lasers%lasers(ilaser), mesh, vp, time_)
1373 do idir = 1, space%dim
1376 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) &
1377 - mu(itime) * vp(ip, idir)/p_c
1381 call laser_field(lasers%lasers(ilaser), this%hm_base%uniform_magnetic_field(1:space%dim), time_)
1382 safe_deallocate_a(vp)
1383 case (e_field_vector_potential)
1384 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
1387 call laser_field(lasers%lasers(ilaser), aa, time_)
1388 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) &
1389 - mu(itime) * aa/p_c
1395 gfield => list_get_gauge_field(ext_partners)
1396 if (
associated(gfield))
then
1397 call this%hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
1398 call gauge_field_get_vec_pot(gfield, aa)
1399 this%hm_base%uniform_vector_potential(1:space%dim) = this%hm_base%uniform_vector_potential(1:space%dim) - aa/p_c
1406 if (
allocated(this%ep%e_field) .and.
associated(gfield))
then
1407 this%hm_base%uniform_vector_potential(1:space%periodic_dim) = &
1408 this%hm_base%uniform_vector_potential(1:space%periodic_dim) - time_*this%ep%e_field(1:space%periodic_dim)
1414 if (
allocated(this%ep%a_static))
then
1415 call this%hm_base%allocate_field(mesh, field_vector_potential, .false.)
1419 do idir = 1, space%dim
1420 this%hm_base%vector_potential(idir, ip) = this%hm_base%vector_potential(idir, ip) + this%ep%a_static(ip, idir)
1426 call this%hm_base%accel_copy_pot(mesh)
1429 if (
allocated(this%ep%b_field))
then
1430 call this%hm_base%allocate_field(mesh, field_uniform_magnetic_field, .false.)
1432 this%hm_base%uniform_magnetic_field(idir) = this%hm_base%uniform_magnetic_field(idir) + this%ep%b_field(idir)
1436 call this%hm_base%update_magnetic_terms(mesh, this%ep%gyromagnetic_ratio, this%d%ispin)
1442 call profiling_out(
"HAMILTONIAN_ELEC_UPDATE_EXT_POT")
1448 integer :: ik, imat, nmat, max_npoints, offset, iphase, nphase
1452 if ((.not. this%kpoints%gamma_only()) .or.
allocated(this%hm_base%uniform_vector_potential))
then
1454 call profiling_in(
'UPDATE_PHASES')
1456 do iatom = 1, this%ep%natoms
1457 call projector_init_phases(this%ep%proj(iatom), space%dim, this%d, this%der%boundaries, this%kpoints, &
1458 vec_pot = this%hm_base%uniform_vector_potential, vec_pot_var = this%hm_base%vector_potential)
1461 call profiling_out(
'UPDATE_PHASES')
1464 if (
allocated(this%hm_base%uniform_vector_potential))
then
1465 call this%phase%update(mesh, this%d%kpt, this%kpoints, this%d, space, this%hm_base%uniform_vector_potential)
1468 max_npoints = this%vnl%max_npoints
1469 nmat = this%vnl%nprojector_matrices
1472 if (this%phase%is_allocated() .and.
allocated(this%vnl%projector_matrices))
then
1475 if (this%der%boundaries%spiralBC) nphase = 3
1477 if (.not.
allocated(this%vnl%projector_phases))
then
1478 safe_allocate(this%vnl%projector_phases(1:max_npoints, nphase, nmat, this%d%kpt%start:this%d%kpt%end))
1479 if (accel_is_enabled())
then
1480 call accel_create_buffer(this%vnl%buff_projector_phases, accel_mem_read_only, &
1481 type_cmplx, this%vnl%total_points*nphase*this%d%kpt%nlocal)
1486 do ik = this%d%kpt%start, this%d%kpt%end
1487 do imat = 1, this%vnl%nprojector_matrices
1488 iatom = this%vnl%projector_to_atom(imat)
1489 do iphase = 1, nphase
1491 do ip = 1, this%vnl%projector_matrices(imat)%npoints
1492 this%vnl%projector_phases(ip, imat, iphase, ik) = this%ep%proj(iatom)%phase(ip, iphase, ik)
1495 if (accel_is_enabled() .and. this%vnl%projector_matrices(imat)%npoints > 0)
then
1496 call accel_write_buffer(this%vnl%buff_projector_phases, &
1497 this%vnl%projector_matrices(imat)%npoints, this%vnl%projector_phases(1:, iphase, imat, ik), &
1500 offset = offset + this%vnl%projector_matrices(imat)%npoints
1514 logical,
intent(in) :: states_are_real
1518 if (hm%self_induced_magnetic)
then
1519 if (.not. states_are_real)
then
1522 message(1) =
'No current density for real states since it is identically zero.'
1523 call messages_warning(1)
1532 type(namespace_t),
intent(in) :: namespace
1533 type(grid_t),
intent(in) :: gr
1534 type(states_elec_t),
intent(inout) :: st
1535 type(states_elec_t),
intent(inout) :: hst
1537 integer :: ik, ib, ist
1538 complex(real64),
allocatable :: psi(:, :)
1539 complex(real64),
allocatable :: psiall(:, :, :, :)
1543 do ik = st%d%kpt%start, st%d%kpt%end
1544 do ib = st%group%block_start, st%group%block_end
1549 if (oct_exchange_enabled(hm%oct_exchange))
then
1551 safe_allocate(psiall(gr%np_part, 1:hst%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
1553 call states_elec_get_state(st, gr, psiall)
1555 call oct_exchange_prepare(hm%oct_exchange, gr, psiall, hm%xc, hm%psolver, namespace)
1557 safe_deallocate_a(psiall)
1559 safe_allocate(psi(gr%np_part, 1:hst%d%dim))
1563 call states_elec_get_state(hst, gr, ist, ik, psi)
1564 call oct_exchange_operator(hm%oct_exchange, namespace, gr, psi, ist, ik)
1565 call states_elec_set_state(hst, gr, ist, ik, psi)
1569 safe_deallocate_a(psi)
1579 subroutine magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase)
1581 type(namespace_t),
intent(in) :: namespace
1582 class(mesh_t),
intent(in) :: mesh
1583 complex(real64),
contiguous,
intent(inout) :: psi(:,:)
1584 complex(real64),
contiguous,
intent(out) :: hpsi(:,:)
1585 integer,
intent(in) :: ik
1586 real(real64),
intent(in) :: vmagnus(:, :, :)
1587 logical,
optional,
intent(in) :: set_phase
1589 complex(real64),
allocatable :: auxpsi(:, :), aux2psi(:, :)
1590 integer :: idim, ispin
1595 if (hm%d%dim /= 1)
then
1596 call messages_not_implemented(
"Magnus with spinors", namespace=namespace)
1599 safe_allocate( auxpsi(1:mesh%np_part, 1:hm%d%dim))
1600 safe_allocate(aux2psi(1:mesh%np, 1:hm%d%dim))
1602 ispin = hm%d%get_spin_index(ik)
1606 terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
1609 do idim = 1, hm%d%dim
1610 call lalg_copy(mesh%np, auxpsi(:, idim), hpsi(:, idim))
1611 hpsi(1:mesh%np, idim) = hpsi(1:mesh%np, idim) + hm%ep%Vpsl(1:mesh%np)*psi(1:mesh%np,idim)
1612 call vborders(mesh, hm, psi(:, idim), hpsi(:, idim))
1614 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + vmagnus(1:mesh%np, ispin, 2)*psi(1:mesh%np, 1)
1617 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) - m_zi*vmagnus(1:mesh%np, ispin, 1)*auxpsi(1:mesh%np, 1)
1620 auxpsi(1:mesh%np, 1) = vmagnus(1:mesh%np, ispin, 1)*psi(1:mesh%np, 1)
1622 terms = term_kinetic + term_non_local_potential, set_phase = set_phase)
1623 hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + m_zi*aux2psi(1:mesh%np, 1)
1625 safe_deallocate_a(auxpsi)
1626 safe_deallocate_a(aux2psi)
1631 subroutine vborders (mesh, hm, psi, hpsi)
1632 class(mesh_t),
intent(in) :: mesh
1634 complex(real64),
intent(in) :: psi(:)
1635 complex(real64),
intent(inout) :: hpsi(:)
1641 if (hm%abs_boundaries%abtype == imaginary_absorbing)
then
1643 hpsi(ip) = hpsi(ip) + m_zi*hm%abs_boundaries%mf(ip)*psi(ip)
1665 type(namespace_t),
intent(in) :: namespace
1666 real(real64),
intent(in) :: mass
1670 if (parse_is_defined(namespace,
'ParticleMass'))
then
1671 message(1) =
'Attempting to redefine a non-unit electron mass'
1672 call messages_fatal(1)
1689 type(grid_t),
intent(in) :: gr
1690 type(distributed_t),
intent(in) :: kpt
1691 type(wfs_elec_t),
intent(in) :: psib
1692 type(wfs_elec_t),
intent(out) :: psib_with_phase
1694 integer :: k_offset, n_boundary_points
1698 call psib%copy_to(psib_with_phase)
1699 if (hm%phase%is_allocated())
then
1700 call hm%phase%apply_to(gr, gr%np, conjugate = .false., psib = psib_with_phase, src = psib, async=.
true.)
1703 k_offset = psib%ik - kpt%start
1704 n_boundary_points = int(gr%np_part - gr%np)
1705 call boundaries_set(gr%der%boundaries, gr, psib_with_phase, phase_correction = hm%phase%phase_corr(:, psib%ik), &
1706 buff_phase_corr = hm%phase%buff_phase_corr, offset=k_offset*n_boundary_points, async=.
true.)
1708 call psib%copy_data_to(gr%np, psib_with_phase)
1709 call boundaries_set(gr%der%boundaries, gr, psib_with_phase)
1712 call psib_with_phase%do_pack(copy = .
true.)
1720 class(mesh_t),
intent(in) :: mesh
1721 real(real64),
contiguous,
intent(out) :: diag(:,:)
1722 integer,
intent(in) :: ik
1724 integer :: idim, ip, ispin
1726 real(real64),
allocatable :: ldiag(:)
1730 safe_allocate(ldiag(1:mesh%np))
1734 call derivatives_lapl_diag(hm%der, ldiag)
1736 do idim = 1, hm%d%dim
1737 diag(1:mesh%np, idim) = -m_half/hm%mass*ldiag(1:mesh%np)
1740 select case (hm%d%ispin)
1742 case (unpolarized, spin_polarized)
1743 ispin = hm%d%get_spin_index(ik)
1744 diag(1:mesh%np, 1) = diag(1:mesh%np, 1) + hm%ep%vpsl(1:mesh%np)
1748 diag(ip, 1) = diag(ip, 1) + hm%ep%vpsl(ip)
1749 diag(ip, 2) = diag(ip, 2) + hm%ep%vpsl(ip)
1754 call hm%ks_pot%add_vhxc(diag)
1763#include "hamiltonian_elec_inc.F90"
1766#include "complex.F90"
1767#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, cam)
subroutine, public exchange_operator_end(this)
logical function, public list_has_gauge_field(partners)
type(gauge_field_t) function, pointer, public list_get_gauge_field(partners)
logical function, public list_has_lasers(partners)
type(lasers_t) function, pointer, public list_get_lasers(partners)
integer, parameter, public external_pot_from_file
potential, defined in a file
subroutine, public load_external_potentials(external_potentials, namespace)
integer, parameter, public external_pot_charge_density
user-defined function for charge density
integer, parameter, public external_pot_usdef
user-defined function for local potential
integer, parameter, public external_pot_static_efield
Static electric field.
integer, parameter, public external_pot_static_bfield
Static magnetic field.
subroutine, public gauge_field_get_vec_pot(this, vec_pot)
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
integer, parameter, public rdmft
integer, parameter, public hartree_fock
integer, parameter, public independent_particles
Theory level.
integer, parameter, public generalized_kohn_sham_dft
integer, parameter, public kohn_sham_dft
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 term_mgga
integer, parameter, public field_potential
subroutine vborders(mesh, hm, psi, hpsi)
pure logical function, public hamiltonian_elec_apply_packed(this)
subroutine, public hamiltonian_elec_set_inh(hm, st)
subroutine, public zvmask(mesh, hm, st)
subroutine, public zhamiltonian_elec_apply_batch(hm, namespace, mesh, psib, hpsib, terms, set_bc)
subroutine, public hamiltonian_elec_adjoint(hm)
subroutine, public hamiltonian_elec_end(hm)
subroutine, public hamiltonian_elec_copy_and_set_phase(hm, gr, kpt, psib, psib_with_phase)
Copy a batch to another batch and apply the Bloch phase to it.
subroutine, public zhamiltonian_elec_apply_single(hm, namespace, mesh, psi, hpsi, ist, ik, terms, set_bc, set_phase)
pure logical function hamiltonian_elec_needs_mgga_term(hm, terms)
logical function, public hamiltonian_elec_has_kick(hm)
subroutine, public magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase)
logical function hamiltonian_elec_hermitian(hm)
subroutine, public dhamiltonian_elec_apply_single(hm, namespace, mesh, psi, hpsi, ist, ik, terms, set_bc, set_phase)
subroutine, public hamiltonian_elec_epot_generate(this, namespace, space, gr, ions, ext_partners, st, time)
real(real64) function, public hamiltonian_elec_get_time(this)
subroutine, public dvmask(mesh, hm, st)
logical function, public hamiltonian_elec_needs_current(hm, states_are_real)
subroutine, public hamiltonian_elec_remove_inh(hm)
subroutine zhamiltonian_elec_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
subroutine, public zhamiltonian_elec_apply_atom(hm, namespace, space, latt, species, pos, ia, mesh, psi, vpsi)
subroutine, public zhamiltonian_elec_apply_all(hm, namespace, gr, st, hst)
subroutine, public hamiltonian_elec_diagonal(hm, mesh, diag, ik)
subroutine, public hamiltonian_elec_update_pot(this, mesh, accumulate)
Update the KS potential of the electronic Hamiltonian.
integer, parameter, public velocity
subroutine dhamiltonian_elec_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus)
subroutine hamiltonian_elec_update(this, mesh, namespace, space, ext_partners, time)
(re-)build the Hamiltonian for the next application:
subroutine hamiltonian_elec_span(hm, delta, emin, namespace)
subroutine, public hamiltonian_elec_init(hm, namespace, space, gr, ions, ext_partners, st, theory_level, xc, mc, kpoints, need_exchange, xc_photons)
subroutine dhamiltonian_elec_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
pure logical function, public hamiltonian_elec_inh_term(hm)
subroutine hamiltonian_elec_set_mass(this, namespace, mass)
set the effective electron mass, checking whether it was previously redefined.
subroutine, public dhamiltonian_elec_apply_batch(hm, namespace, mesh, psib, hpsib, terms, set_bc)
subroutine, public hamiltonian_elec_update_with_ext_pot(this, mesh, space, ext_partners, time, mu)
This is an extension of "hamiltonian_elec_update_pot" to be used by the CFM4 propagator....
subroutine, public hamiltonian_elec_not_adjoint(hm)
subroutine zhamiltonian_elec_apply(hm, namespace, mesh, psib, hpsib, terms, set_bc)
This module defines classes and functions for interaction partners.
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.
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_init(this, namespace, space, level, gr, ions, st, mc, kpoints)
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)
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)
real(real64), dimension(:,:), allocatable delta
D_E matrix in JCP 139, 024105 (2013).
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)
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public states_elec_dim_copy(dout, din)
subroutine, public states_elec_dim_end(dim)
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
This module provides routines for communicating states when using states parallelization.
integer pure function, public symmetries_identity_index(this)
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
type(xc_cam_t), parameter, public cam_exact_exchange
Use only Hartree Fock exact exchange.
integer, parameter, public xc_oep_x_slater
Slater approximation to the exact exchange.
integer, parameter, public func_x
logical pure function, public family_is_mgga_with_exc(xcs)
Is the xc function part of the mGGA family with an energy functional.
logical pure function, public family_is_hybrid(xcs)
Returns true if the functional is an hybrid functional.
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.