30 use,
intrinsic :: iso_fortran_env
60 integer :: max_iter_berry
66 type(berry_t),
intent(inout) :: this
67 type(namespace_t),
intent(in) :: namespace
81 call parse_variable(namespace,
'MaximumIterBerry', 10, this%max_iter_berry)
82 if (this%max_iter_berry < 0) this%max_iter_berry = huge(this%max_iter_berry)
90 iter, ks, ions, ext_partners)
91 type(berry_t),
intent(in) :: this
92 type(namespace_t),
intent(in) :: namespace
93 type(electron_space_t),
intent(in) :: space
94 type(eigensolver_t),
intent(inout) :: eigensolver
95 type(grid_t),
intent(in) :: gr
96 type(states_elec_t),
intent(inout) :: st
97 type(hamiltonian_elec_t),
intent(inout) :: hm
98 integer,
intent(in) :: iter
99 type(v_ks_t),
intent(inout) :: ks
100 type(ions_t),
intent(in) :: ions
101 type(partner_list_t),
intent(in) :: ext_partners
103 integer :: iberry, idir
104 logical :: berry_conv
105 real(real64) :: dipole_prev(space%dim), dipole(space%dim)
106 real(real64),
parameter :: tol = 1e-5_real64
110 assert(
allocated(hm%vberry))
112 if (st%parallel_in_states)
then
118 do iberry = 1, this%max_iter_berry
119 eigensolver%converged = 0
120 call eigensolver%run(namespace, gr, st, hm, space, ext_partners, iter)
123 call berry_potential(st, namespace, space, gr, hm%ions%latt, hm%ep%e_field, hm%vberry)
127 hm%ep%e_field(1:space%periodic_dim), hm%vberry(1:gr%np, 1:hm%d%nspin))
130 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=.false.)
134 write(
message(1),
'(a,9f12.6)')
'Dipole = ', dipole(1:space%dim)
138 do idir = 1, space%periodic_dim
139 if (abs(dipole_prev(idir)) > 1e-10_real64)
then
140 berry_conv = berry_conv .and. (abs(dipole(idir) - dipole_prev(idir)) < tol &
141 .or.(abs((dipole(idir) - dipole_prev(idir)) / dipole_prev(idir)) < tol))
143 berry_conv = berry_conv .and. (abs(dipole(idir) - dipole_prev(idir)) < tol)
156 subroutine calc_dipole(dipole, space, mesh, st, ions)
157 real(real64),
intent(out) :: dipole(:)
158 class(
space_t),
intent(in) :: space
159 class(
mesh_t),
intent(in) :: mesh
161 type(
ions_t),
intent(in) :: ions
163 integer :: ispin, idir
164 real(real64) :: e_dip(space%dim + 1, st%d%nspin), n_dip(space%dim), nquantumpol
168 assert(.not. ions%latt%nonorthogonal)
170 dipole(1:space%dim) =
m_zero
172 do ispin = 1, st%d%nspin
176 n_dip = ions%dipole()
178 do idir = 1, space%dim
180 if (idir <= space%periodic_dim)
then
181 dipole(idir) = -n_dip(idir) -
berry_dipole(st, mesh, ions%latt, space, idir)
184 nquantumpol = nint(dipole(idir)/norm2(ions%latt%rlattice(:, idir)))
185 dipole(idir) = dipole(idir) - nquantumpol * norm2(ions%latt%rlattice(:, idir))
189 e_dip(idir + 1, 1) = sum(e_dip(idir + 1, :))
190 dipole(idir) = -n_dip(idir) - e_dip(idir + 1, 1)
209 real(real64) function
berry_dipole(st, mesh, latt, space, dir) result(dipole)
211 class(
mesh_t),
intent(in) :: mesh
213 class(
space_t),
intent(in) :: space
214 integer,
intent(in) :: dir
217 complex(real64) :: det
226 do ik = st%d%kpt%start, st%d%kpt%end
229 if (abs(det) >
m_tiny)
then
230 dipole = dipole + aimag(
log(det))
234 dipole = -(norm2(latt%rlattice(:,dir)) / (
m_two*
m_pi)) * dipole
245 complex(real64) function berry_phase_det(st, mesh, latt, space, dir, ik)
result(det)
246 type(states_elec_t),
intent(in) :: st
247 class(mesh_t),
intent(in) :: mesh
248 type(lattice_vectors_t),
intent(in) :: latt
249 class(space_t),
intent(in) :: space
250 integer,
intent(in) :: dir
251 integer,
intent(in) :: ik
253 integer :: ist, noccst, gvector(3)
254 complex(real64),
allocatable :: matrix(:, :), tmp(:), phase(:)
255 complex(real64),
allocatable :: psi(:, :), psi2(:, :)
262 if (st%occ(ist, ik) > m_epsilon) noccst = ist
265 safe_allocate(matrix(1:noccst, 1:noccst))
272 det = lalg_determinant(noccst, matrix, preserve_mat=.false.) ** st%smear%el_per_state
277 safe_deallocate_a(matrix)
278 safe_deallocate_a(tmp)
279 safe_deallocate_a(phase)
280 safe_deallocate_a(psi)
281 safe_deallocate_a(psi2)
289 type(states_elec_t),
intent(in) :: st
290 class(mesh_t),
intent(in) :: mesh
291 type(lattice_vectors_t),
intent(in) :: latt
292 class(space_t),
intent(in) :: space
293 integer,
intent(in) :: nst
294 integer,
intent(in) :: ik
295 integer,
intent(in) :: ik2
296 integer,
intent(in) :: gvector(:)
297 complex(real64),
intent(out) :: matrix(:,:)
299 integer :: ist, ist2, idim, ip, idir
300 complex(real64),
allocatable :: tmp(:), phase(:)
301 complex(real64),
allocatable :: psi(:, :), psi2(:, :)
307 assert(.not. st%parallel_in_states)
308 assert(.not. latt%nonorthogonal)
310 safe_allocate(tmp(1:mesh%np))
311 safe_allocate(phase(1:mesh%np))
312 safe_allocate(psi(1:mesh%np, 1:st%d%dim))
313 safe_allocate(psi2(1:mesh%np, 1:st%d%dim))
316 do idir = 1, space%dim
317 if (gvector(idir) == 0) cycle
318 norm = norm2(latt%rlattice(:,idir))
319 if (idir > space%periodic_dim) norm = norm * m_two * mesh%box%bounding_box_l(idir)
321 phase(ip) = phase(ip) +
exp(-m_zi*gvector(idir)*(m_two*m_pi/norm)*mesh%x(ip, idir))
326 call states_elec_get_state(st, mesh, ist, ik, psi)
328 call states_elec_get_state(st, mesh, ist2, ik2, psi2)
329 matrix(ist, ist2) = m_z0
330 do idim = 1, st%d%dim
333 tmp(ip) = conjg(psi(ip, idim))*phase(ip)*psi2(ip, idim)
336 matrix(ist, ist2) = matrix(ist, ist2) + zmf_integrate(mesh, tmp)
341 safe_deallocate_a(tmp)
342 safe_deallocate_a(phase)
343 safe_deallocate_a(psi)
344 safe_deallocate_a(psi2)
357 subroutine berry_potential(st, namespace, space, mesh, latt, e_field, pot)
358 type(states_elec_t),
intent(in) :: st
359 type(namespace_t),
intent(in) :: namespace
360 class(space_t),
intent(in) :: space
361 class(mesh_t),
intent(in) :: mesh
362 type(lattice_vectors_t),
intent(in) :: latt
363 real(real64),
intent(in) :: e_field(:)
364 real(real64),
intent(out) :: pot(:,:)
366 integer :: ispin, idir
367 complex(real64) :: factor, det
371 if (latt%nonorthogonal)
then
372 call messages_not_implemented(
"Berry phase for non-orthogonal cells", namespace=namespace)
376 call messages_not_implemented(
"Berry phase with k-points", namespace=namespace)
379 pot(1:mesh%np, 1:st%d%nspin) = m_zero
381 do ispin = 1, st%d%nspin
382 do idir = 1, space%periodic_dim
383 if (abs(e_field(idir)) > m_epsilon)
then
386 if (abs(det) > m_epsilon)
then
387 factor = e_field(idir) * (norm2(latt%rlattice(:,idir)) / (m_two*m_pi)) / det
391 write(message(1),*)
"Divide by zero: dir = ", idir,
" Berry-phase determinant = ", det
392 call messages_fatal(1, namespace=namespace)
394 pot(1:mesh%np, ispin) = pot(1:mesh%np, ispin) + &
395 aimag(factor *
exp(m_two * m_pi * m_zi * mesh%x(1:mesh%np, idir) / norm2(latt%rlattice(:,idir))))
405 real(real64) function berry_energy_correction(st, space, mesh, latt, e_field, vberry) result(delta)
406 type(states_elec_t),
intent(in) :: st
407 class(space_t),
intent(in) :: space
408 class(mesh_t),
intent(in) :: mesh
409 type(lattice_vectors_t),
intent(in) :: latt
410 real(real64),
intent(in) :: e_field(:)
411 real(real64),
intent(in) :: vberry(:,:)
413 integer :: ispin, idir
419 do ispin = 1, st%d%nspin
420 delta = delta - dmf_dotp(mesh, st%rho(1:mesh%np, ispin), vberry(1:mesh%np, ispin))
424 do idir = 1, space%periodic_dim
425 delta = delta -
berry_dipole(st, mesh, latt, space, idir) * e_field(idir)
double log(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
subroutine, public berry_perform_internal_scf(this, namespace, space, eigensolver, gr, st, hm, iter, ks, ions, ext_partners)
subroutine berry_phase_matrix(st, mesh, latt, space, nst, ik, ik2, gvector, matrix)
complex(real64) function berry_phase_det(st, mesh, latt, space, dir, ik)
real(real64) function, public berry_energy_correction(st, space, mesh, latt, e_field, vberry)
subroutine, public berry_potential(st, namespace, space, mesh, latt, e_field, pot)
local potential for electric enthalpy of uniform field in single-point Berry phase
subroutine, public berry_init(this, namespace)
subroutine, public calc_dipole(dipole, space, mesh, st, ions)
real(real64) function, public berry_dipole(st, mesh, latt, space, dir)
Uses the single-point Berry`s phase method to calculate dipole moment in a periodic system.
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
real(real64), parameter, public m_pi
some mathematical constants
real(real64), parameter, public m_tiny
This module implements the underlying real-space grid.
This module defines classes and functions for interaction partners.
This module defines various routines, operating on mesh functions.
subroutine, public dmf_multipoles(mesh, ff, lmax, multipole, mask)
This routine calculates the multipoles of a function ff.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_not_implemented(feature, namespace)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
This module handles spin dimensions of the states and the k-point distribution.
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Describes mesh distribution to nodes.
The states_elec_t class contains all electronic wave functions.