Octopus
berry.F90
Go to the documentation of this file.
1!! Copyright (C) 2010 D. Strubbe
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module berry_oct_m
22 use debug_oct_m
25 use global_oct_m
26 use grid_oct_m
29 use ions_oct_m
30 use, intrinsic :: iso_fortran_env
33 use mesh_oct_m
38 use parser_oct_m
39 use smear_oct_m
40 use space_oct_m
43 use v_ks_oct_m
44
45 implicit none
46
47 private
48 public :: &
49 berry_t, &
50 berry_init, &
56
57 type berry_t
58 private
59
60 integer :: max_iter_berry
61 end type berry_t
62contains
63
64 ! ---------------------------------------------------------
65 subroutine berry_init(this, namespace)
66 type(berry_t), intent(inout) :: this
67 type(namespace_t), intent(in) :: namespace
68
69 push_sub(berry_init)
70
71 !%Variable MaximumIterBerry
72 !%Type integer
73 !%Default 10
74 !%Section SCF::Convergence
75 !%Description
76 !% Maximum number of iterations for the Berry potential, within each SCF iteration.
77 !% Only applies if a <tt>StaticElectricField</tt> is applied in a periodic direction.
78 !% The code will move on to the next SCF iteration even if convergence
79 !% has not been achieved. -1 means unlimited.
80 !%End
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)
83
84
85 pop_sub(berry_init)
86 end subroutine berry_init
87
88 ! ---------------------------------------------------------
89 subroutine berry_perform_internal_scf(this, namespace, space, eigensolver, gr, st, hm, &
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
102
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
107
109
110 assert(allocated(hm%vberry))
111
112 if (st%parallel_in_states) then
113 call messages_not_implemented("Berry phase parallel in states", namespace=namespace)
114 end if
115
116 call calc_dipole(dipole, space, gr, st, ions)
117
118 do iberry = 1, this%max_iter_berry
119 eigensolver%converged = 0
120 call eigensolver%run(namespace, gr, st, hm, space, ext_partners, iter)
121
122 !Calculation of the Berry potential
123 call berry_potential(st, namespace, space, gr, hm%ions%latt, hm%ep%e_field, hm%vberry)
124
125 !Calculation of the corresponding energy
126 hm%energy%berry = berry_energy_correction(st, space, gr, hm%ions%latt, &
127 hm%ep%e_field(1:space%periodic_dim), hm%vberry(1:gr%np, 1:hm%d%nspin))
128
129 !We recompute the KS potential
130 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=.false.)
131
132 dipole_prev = dipole
133 call calc_dipole(dipole, space, gr, st, ions)
134 write(message(1),'(a,9f12.6)') 'Dipole = ', dipole(1:space%dim)
135 call messages_info(1, namespace=namespace)
136
137 berry_conv = .true.
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))
142 else
143 berry_conv = berry_conv .and. (abs(dipole(idir) - dipole_prev(idir)) < tol)
144 end if
145 end do
146
147 if (berry_conv) exit
148 end do
149
151 end subroutine berry_perform_internal_scf
153 ! ---------------------------------------------------------
154 !TODO: This should be a method of the electronic system directly
155 ! that can be exposed
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
160 type(states_elec_t), intent(in) :: st
161 type(ions_t), intent(in) :: ions
162
163 integer :: ispin, idir
164 real(real64) :: e_dip(space%dim + 1, st%d%nspin), n_dip(space%dim), nquantumpol
165
166 push_sub(calc_dipole)
167
168 assert(.not. ions%latt%nonorthogonal)
169
170 dipole(1:space%dim) = m_zero
171
172 do ispin = 1, st%d%nspin
173 call dmf_multipoles(mesh, st%rho(:, ispin), 1, e_dip(:, ispin))
174 end do
175
176 n_dip = ions%dipole()
177
178 do idir = 1, space%dim
179 ! in periodic directions use single-point Berry`s phase calculation
180 if (idir <= space%periodic_dim) then
181 dipole(idir) = -n_dip(idir) - berry_dipole(st, mesh, ions%latt, space, idir)
182
183 ! use quantum of polarization to reduce to smallest possible magnitude
184 nquantumpol = nint(dipole(idir)/norm2(ions%latt%rlattice(:, idir)))
185 dipole(idir) = dipole(idir) - nquantumpol * norm2(ions%latt%rlattice(:, idir))
186 ! in aperiodic directions use normal dipole formula
187
188 else
189 e_dip(idir + 1, 1) = sum(e_dip(idir + 1, :))
190 dipole(idir) = -n_dip(idir) - e_dip(idir + 1, 1)
191 end if
192 end do
193
194 pop_sub(calc_dipole)
195 end subroutine calc_dipole
196
197
198 ! ---------------------------------------------------------
209 real(real64) function berry_dipole(st, mesh, latt, space, dir) result(dipole)
210 type(states_elec_t), intent(in) :: st
211 class(mesh_t), intent(in) :: mesh
212 type(lattice_vectors_t), intent(in) :: latt
213 class(space_t), intent(in) :: space
214 integer, intent(in) :: dir
215
216 integer :: ik
217 complex(real64) :: det
218
219 push_sub(berry_dipole)
220
221 !There is a missing reduction here.
222 !However, as we are limited to a single k-point at the moment, this is not a problem.
223 assert(st%nik == 1)
224
225 dipole = m_zero
226 do ik = st%d%kpt%start, st%d%kpt%end ! determinants for different spins and k-points multiply since matrix is block-diagonal
227 det = berry_phase_det(st, mesh, latt, space, dir, ik)
228 !Check in case det is zero. This should never occur in practical calculations
229 if (abs(det) > m_tiny) then
230 dipole = dipole + aimag(log(det))
231 end if
232 end do
233
234 dipole = -(norm2(latt%rlattice(:,dir)) / (m_two*m_pi)) * dipole
235
236 pop_sub(berry_dipole)
237 end function berry_dipole
238
239
240 ! ---------------------------------------------------------
241 !! \f[
242 !! det <\psi_i|exp(-i(2*\pi/L)x)|\psi_j>
243 !! \f]
244 !! E Yaschenko, L Fu, L Resca, R Resta, Phys. Rev. B 58, 1222-1229 (1998)
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
252
253 integer :: ist, noccst, gvector(3)
254 complex(real64), allocatable :: matrix(:, :), tmp(:), phase(:)
255 complex(real64), allocatable :: psi(:, :), psi2(:, :)
256
257 push_sub(berry_phase_det)
258
259 ! find how many states are occupied on this k-point. Formalism only works if semiconducting anyway.
260 noccst = 0
261 do ist = 1, st%nst
262 if (st%occ(ist, ik) > m_epsilon) noccst = ist
263 end do
264
265 safe_allocate(matrix(1:noccst, 1:noccst))
266
267 gvector(:) = 0
268 gvector(dir) = 1
269 call berry_phase_matrix(st, mesh, latt, space, noccst, ik, ik, gvector, matrix)
270
271 if (noccst > 0) then
272 det = lalg_determinant(noccst, matrix, preserve_mat=.false.) ** st%smear%el_per_state
273 else
274 det = m_one
275 end if
276
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)
282
283 pop_sub(berry_phase_det)
284 end function berry_phase_det
285
286
287 ! ---------------------------------------------------------
288 subroutine berry_phase_matrix(st, mesh, latt, space, nst, ik, ik2, gvector, matrix)
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(:,:)
298
299 integer :: ist, ist2, idim, ip, idir
300 complex(real64), allocatable :: tmp(:), phase(:)
301 complex(real64), allocatable :: psi(:, :), psi2(:, :)
302 real(real64) :: norm
303 ! FIXME: real/cplx versions
305 push_sub(berry_phase_matrix)
306
307 assert(.not. st%parallel_in_states)
308 assert(.not. latt%nonorthogonal)
309
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))
314
315 phase = m_zero
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)
320 do ip = 1, mesh%np
321 phase(ip) = phase(ip) + exp(-m_zi*gvector(idir)*(m_two*m_pi/norm)*mesh%x(ip, idir))
322 end do
323 end do
324
325 do ist = 1, nst
326 call states_elec_get_state(st, mesh, ist, ik, psi)
327 do ist2 = 1, nst
328 call states_elec_get_state(st, mesh, ist2, ik2, psi2)
329 matrix(ist, ist2) = m_z0
330 do idim = 1, st%d%dim ! spinor components
331
332 do ip = 1, mesh%np
333 tmp(ip) = conjg(psi(ip, idim))*phase(ip)*psi2(ip, idim)
334 end do
335
336 matrix(ist, ist2) = matrix(ist, ist2) + zmf_integrate(mesh, tmp)
337 end do
338 end do !ist2
339 end do !ist
341 safe_deallocate_a(tmp)
342 safe_deallocate_a(phase)
343 safe_deallocate_a(psi)
344 safe_deallocate_a(psi2)
345
346 pop_sub(berry_phase_matrix)
347 end subroutine berry_phase_matrix
348
349
350 ! ---------------------------------------------------------
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(:,:)
365
366 integer :: ispin, idir
367 complex(real64) :: factor, det
368
369 push_sub(berry_potential)
370
371 if (latt%nonorthogonal) then
372 call messages_not_implemented("Berry phase for non-orthogonal cells", namespace=namespace)
373 end if
374
375 if (st%nik > 1) then
376 call messages_not_implemented("Berry phase with k-points", namespace=namespace)
377 end if
378
379 pot(1:mesh%np, 1:st%d%nspin) = m_zero
380
381 do ispin = 1, st%d%nspin
382 do idir = 1, space%periodic_dim
383 if (abs(e_field(idir)) > m_epsilon) then
384 ! calculate the ip-independent part first
385 det = berry_phase_det(st, mesh, latt, space, idir, ispin)
386 if (abs(det) > m_epsilon) then
387 factor = e_field(idir) * (norm2(latt%rlattice(:,idir)) / (m_two*m_pi)) / det
388 else
389 ! If det = 0, mu = -infinity, so this condition should never be reached
390 ! if things are working properly.
391 write(message(1),*) "Divide by zero: dir = ", idir, " Berry-phase determinant = ", det
392 call messages_fatal(1, namespace=namespace)
393 end if
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))))
396 end if
397 end do
398 end do
399
400 pop_sub(berry_potential)
401 end subroutine berry_potential
402
403
404 ! ---------------------------------------------------------
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(:,:)
412
413 integer :: ispin, idir
414
416
417 ! first we calculate expectation value of Berry potential, to subtract off
418 delta = m_zero
419 do ispin = 1, st%d%nspin
420 delta = delta - dmf_dotp(mesh, st%rho(1:mesh%np, ispin), vberry(1:mesh%np, ispin))
421 end do
422
423 ! the real energy contribution is -mu.E
424 do idir = 1, space%periodic_dim
425 delta = delta - berry_dipole(st, mesh, latt, space, idir) * e_field(idir)
426 end do
427
429 end function berry_energy_correction
430
431end module berry_oct_m
432
433!! Local Variables:
434!! mode: f90
435!! coding: utf-8
436!! End:
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)
Definition: berry.F90:186
subroutine berry_phase_matrix(st, mesh, latt, space, nst, ik, ik2, gvector, matrix)
Definition: berry.F90:384
complex(real64) function berry_phase_det(st, mesh, latt, space, dir, ik)
Definition: berry.F90:341
real(real64) function, public berry_energy_correction(st, space, mesh, latt, e_field, vberry)
Definition: berry.F90:501
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
Definition: berry.F90:453
subroutine, public berry_init(this, namespace)
Definition: berry.F90:161
subroutine, public calc_dipole(dipole, space, mesh, st, ions)
Definition: berry.F90:252
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.
Definition: berry.F90:305
real(real64), parameter, public m_two
Definition: global.F90:193
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
real(real64), parameter, public m_tiny
Definition: global.F90:208
This module implements the underlying real-space grid.
Definition: grid.F90:119
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.
Definition: mesh.F90:120
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
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)
Definition: v_ks.F90:749
Describes mesh distribution to nodes.
Definition: mesh.F90:187
The states_elec_t class contains all electronic wave functions.
int true(void)