Octopus
lda_u.F90
Go to the documentation of this file.
1!! Copyright (C) 2016-2020 N. Tancogne-Dejean
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 lda_u_oct_m
22 use accel_oct_m
25 use batch_oct_m
27 use comm_oct_m
28 use debug_oct_m
32 use energy_oct_m
33 use global_oct_m
34 use grid_oct_m
36 use ions_oct_m
40 use loct_oct_m
43 use math_oct_m
44 use mesh_oct_m
47 use mpi_oct_m
53 use parser_oct_m
55 use phase_oct_m
58 use space_oct_m
65 use types_oct_m
68
69 implicit none
70
71 private
72
73 public :: &
74 lda_u_t, &
75 lda_u_init, &
80 lda_u_end, &
100 dlda_u_rvu, &
101 zlda_u_rvu, &
106
107
108 integer, public, parameter :: &
109 DFT_U_NONE = 0, &
110 dft_u_empirical = 1, &
111 dft_u_acbn0 = 2
112
113 integer, public, parameter :: &
114 DFT_U_FLL = 0, &
115 dft_u_amf = 1, &
117
118
121 type lda_u_t
122 private
123 integer, public :: level = dft_u_none
124
125 ! DFT+U basic variables
126 real(real64), allocatable, public :: dn(:,:,:,:)
127 real(real64), allocatable :: dV(:,:,:,:)
128
129 ! ACBN0 variables
130 complex(real64), allocatable, public :: zn(:,:,:,:)
131 complex(real64), allocatable :: zV(:,:,:,:)
132 real(real64), allocatable, public :: dn_alt(:,:,:,:)
133 complex(real64), allocatable, public :: zn_alt(:,:,:,:)
134
135 real(real64), allocatable :: renorm_occ(:,:,:,:,:)
136
137 ! Coulomb integrales
138 real(real64), allocatable, public :: coulomb(:,:,:,:,:)
139 ! !<(for the ACBN0 functional)
140 complex(real64), allocatable, public :: zcoulomb(:,:,:,:,:,:,:)
141 ! !< (for the ACBN0 functional with spinors)
142
143 type(orbitalbasis_t), public :: basis
144 type(orbitalset_t), pointer, public :: orbsets(:) => null()
145 integer, public :: norbsets = 0
146
147 integer, public :: nspins = 0
148 integer, public :: spin_channels = 0
149 integer :: nspecies = 0
150 integer, public :: maxnorbs = 0
151 integer :: max_np = 0
152
153 logical :: useAllOrbitals = .false.
154 logical :: skipSOrbitals = .true.
155 logical :: freeze_occ = .false.
156 logical :: freeze_u = .false.
157 logical, public :: intersite = .false.
158 real(real64) :: intersite_radius = m_zero
159 logical, public :: basisfromstates = .false.
160 real(real64) :: acbn0_screening = m_one
161 integer, allocatable :: basisstates(:)
162 integer, allocatable :: basisstates_os(:)
163 ! !! the state specified in basisstate(:)
164 logical :: rot_inv = .false.
165 integer :: double_couting = dft_u_fll
166 integer :: sm_poisson = sm_poisson_direct
167 real(real64), allocatable :: dc_alpha(:)
168
169 type(lattice_vectors_t), pointer :: latt
170
171 type(distributed_t) :: orbs_dist
172
173 ! Intersite interaction variables
174 integer, public :: maxneighbors = 0
175 real(real64), allocatable, public :: dn_ij(:,:,:,:,:), dn_alt_ij(:,:,:,:,:), dn_alt_ii(:,:,:,:,:)
176 complex(real64), allocatable, public :: zn_ij(:,:,:,:,:), zn_alt_ij(:,:,:,:,:), zn_alt_ii(:,:,:,:,:)
177
178 ! Symmetrization-related variables
179 logical :: symmetrize_occ_matrices
180 integer, allocatable :: inv_map_symm(:,:)
181 integer :: nsym
182 real(real64), allocatable :: symm_weight(:,:,:,:)
183 end type lda_u_t
184
185contains
186
187 ! ---------------------------------------------------------
188 subroutine lda_u_init(this, namespace, space, level, gr, ions, st, mc, kpoints)
189 type(lda_u_t), target, intent(inout) :: this
190 type(namespace_t), intent(in) :: namespace
191 class(space_t), intent(in) :: space
192 integer, intent(in) :: level
193 type(grid_t), intent(in) :: gr
194 type(ions_t), target, intent(in) :: ions
195 type(states_elec_t), intent(in) :: st
196 type(multicomm_t), intent(in) :: mc
197 type(kpoints_t), intent(in) :: kpoints
198
199 integer :: is, ierr
200 type(block_t) :: blk
201
202 push_sub(lda_u_init)
204 assert(.not. (level == dft_u_none))
205
206 call messages_print_with_emphasis(msg="DFT+U", namespace=namespace)
207 if (gr%parallel_in_domains) call messages_experimental("dft+u parallel in domains", namespace=namespace)
208 this%level = level
209
210 this%latt => ions%latt
211
212 !%Variable DFTUBasisFromStates
213 !%Type logical
214 !%Default no
215 !%Section Hamiltonian::DFT+U
216 !%Description
217 !% If set to yes, Octopus will construct the localized basis from
218 !% user-defined states. The states are taken at the Gamma point (or the first k-point of the
219 !% states in the restart_proj folder.
220 !% The states are defined via the block DFTUBasisStates
221 !%End
222 call parse_variable(namespace, 'DFTUBasisFromStates', .false., this%basisfromstates)
223 if (this%basisfromstates) call messages_experimental("DFTUBasisFromStates", namespace=namespace)
224
225 !%Variable DFTUDoubleCounting
226 !%Type integer
227 !%Default dft_u_fll
228 !%Section Hamiltonian::DFT+U
229 !%Description
230 !% This variable selects which DFT+U
231 !% double counting term is used.
232 !%Option dft_u_fll 0
233 !% (Default) The Fully Localized Limit (FLL)
234 !%Option dft_u_amf 1
235 !% (Experimental) Around mean field double counting, as defined in PRB 44, 943 (1991) and PRB 49, 14211 (1994).
236 !%Option dft_u_mix 2
237 !% (Experimental) Mixed double countind term as introduced by Petukhov et al., PRB 67, 153106 (2003).
238 !% This recovers the FLL and AMF as limiting cases.
239 !%End
240 call parse_variable(namespace, 'DFTUDoubleCounting', dft_u_fll, this%double_couting)
241 call messages_print_var_option('DFTUDoubleCounting', this%double_couting, namespace=namespace)
242 if (this%double_couting /= dft_u_fll) call messages_experimental("DFTUDoubleCounting /= dft_u_ffl", namespace=namespace)
243 if (st%d%ispin == spinors .and. this%double_couting /= dft_u_fll) then
244 call messages_not_implemented("AMF and MIX double counting with spinors", namespace=namespace)
245 end if
247 !%Variable DFTUPoissonSolver
248 !%Type integer
249 !%Section Hamiltonian::DFT+U
250 !%Description
251 !% This variable selects which Poisson solver
252 !% is used to compute the Coulomb integrals over a submesh.
253 !% These are non-periodic Poisson solvers.
254 !% The FFT Poisson solver with spherical cutoff is used by default.
255 !%
256 !%Option dft_u_poisson_direct 0
257 !% Direct Poisson solver. Slow but working in all cases.
258 !%Option dft_u_poisson_isf 1
259 !% (Experimental) ISF Poisson solver on a submesh.
260 !% This does not work for non-orthogonal cells nor domain parallelization.
261 !%Option dft_u_poisson_psolver 2
262 !% (Experimental) PSolver Poisson solver on a submesh.
263 !% This does not work for non-orthogonal cells nor domain parallelization.
264 !% Requires the PSolver external library.
265 !%Option dft_u_poisson_fft 3
266 !% (Default) FFT Poisson solver on a submesh.
267 !% This uses the 0D periodic version of the FFT kernels.
268 !%End
269 call parse_variable(namespace, 'DFTUPoissonSolver', sm_poisson_fft, this%sm_poisson)
270 call messages_print_var_option('DFTUPoissonSolver', this%sm_poisson, namespace=namespace)
271 if (this%sm_poisson /= sm_poisson_direct .and. this%sm_poisson /= sm_poisson_fft) then
272 call messages_experimental("DFTUPoissonSolver different from dft_u_poisson_direct", namespace=namespace)
273 call messages_experimental("and dft_u_poisson_fft", namespace=namespace)
274 end if
275 if (this%sm_poisson == sm_poisson_isf) then
276 if (gr%parallel_in_domains) then
277 call messages_not_implemented("DFTUPoissonSolver=dft_u_poisson_isf with domain parallelization", namespace=namespace)
278 end if
279 if (ions%latt%nonorthogonal) then
280 call messages_not_implemented("DFTUPoissonSolver=dft_u_poisson_isf with non-orthogonal cells", namespace=namespace)
281 end if
282 end if
283 if (this%sm_poisson == sm_poisson_psolver) then
284#if !(defined HAVE_PSOLVER)
285 message(1) = "The PSolver Poisson solver cannot be used since the code was not compiled with the PSolver library."
286 call messages_fatal(1, namespace=namespace)
287#endif
288 if (gr%parallel_in_domains) then
289 call messages_not_implemented("DFTUPoissonSolver=dft_u_poisson_psolver with domain parallelization", namespace=namespace)
290 end if
291 if (ions%latt%nonorthogonal) then
292 call messages_not_implemented("DFTUPoissonSolver=dft_u_poisson_psolver with non-orthogonal cells", namespace=namespace)
293 end if
294 end if
295
296 if (this%level == dft_u_acbn0) then
297 !%Variable UseAllAtomicOrbitals
298 !%Type logical
299 !%Default no
300 !%Section Hamiltonian::DFT+U
301 !%Description
302 !% If set to yes, Octopus will determine the effective U for all atomic orbitals
303 !% from the peusopotential. Only available with ACBN0 functional.
304 !% It is strongly recommended to set AOLoewdin=yes when using the option.
305 !%End
306 call parse_variable(namespace, 'UseAllAtomicOrbitals', .false., this%useAllOrbitals)
307 if (this%useAllOrbitals) call messages_experimental("UseAllAtomicOrbitals", namespace=namespace)
308
309 !%Variable SkipSOrbitals
310 !%Type logical
311 !%Default no
312 !%Section Hamiltonian::DFT+U
313 !%Description
314 !% If set to yes, Octopus will determine the effective U for all atomic orbitals
315 !% from the peusopotential but s orbitals. Only available with ACBN0 functional.
316 !%End
317 call parse_variable(namespace, 'SkipSOrbitals', .true., this%skipSOrbitals)
318 if (.not. this%SkipSOrbitals) call messages_experimental("SkipSOrbitals", namespace=namespace)
319
320 !%Variable ACBN0Screening
321 !%Type float
322 !%Default 1.0
323 !%Section Hamiltonian::DFT+U
324 !%Description
325 !% If set to 0, no screening will be included in the ACBN0 functional, and the U
326 !% will be estimated from bare Hartree-Fock. If set to 1 (default), the full screening
327 !% of the U, as defined in the ACBN0 functional, is used.
328 !%End
329 call parse_variable(namespace, 'ACBN0Screening', m_one, this%acbn0_screening)
330 call messages_print_var_value('ACBN0Screening', this%acbn0_screening, namespace=namespace)
331
332 !%Variable ACBN0RotationallyInvariant
333 !%Type logical
334 !%Section Hamiltonian::DFT+U
335 !%Description
336 !% If set to yes, Octopus will use for U and J a formula which is rotationally invariant.
337 !% This is different from the original formula for U and J.
338 !% This is activated by default, except in the case of spinors, as this is not yet implemented in this case.
339 !%End
340 call parse_variable(namespace, 'ACBN0RotationallyInvariant', st%d%ispin /= spinors, this%rot_inv)
341 call messages_print_var_value('ACBN0RotationallyInvariant', this%rot_inv, namespace=namespace)
342 if (this%rot_inv .and. st%d%ispin == spinors) then
343 call messages_not_implemented("Rotationally invariant ACBN0 with spinors", namespace=namespace)
344 end if
345
346 !%Variable ACBN0IntersiteInteraction
347 !%Type logical
348 !%Default no
349 !%Section Hamiltonian::DFT+U
350 !%Description
351 !% If set to yes, Octopus will determine the effective intersite interaction V
352 !% Only available with ACBN0 functional.
353 !% It is strongly recommended to set AOLoewdin=yes when using the option.
354 !%End
355 call parse_variable(namespace, 'ACBN0IntersiteInteraction', .false., this%intersite)
356 call messages_print_var_value('ACBN0IntersiteInteraction', this%intersite, namespace=namespace)
357 if (this%intersite) call messages_experimental("ACBN0IntersiteInteraction", namespace=namespace)
358
359 if (this%intersite) then
360
361 !This is a non local operator. To make this working, one probably needs to apply the
362 ! symmetries to the generalized occupation matrices
363 if (kpoints%use_symmetries) then
364 call messages_not_implemented("Intersite interaction with kpoint symmetries", namespace=namespace)
365 end if
366
367 !%Variable ACBN0IntersiteCutoff
368 !%Type float
369 !%Section Hamiltonian::DFT+U
370 !%Description
371 !% The cutoff radius defining the maximal intersite distance considered.
372 !% Only available with ACBN0 functional with intersite interaction.
373 !%End
374 call parse_variable(namespace, 'ACBN0IntersiteCutoff', m_zero, this%intersite_radius, unit = units_inp%length)
375 if (abs(this%intersite_radius) < m_epsilon) then
376 call messages_write("ACBN0IntersiteCutoff must be greater than 0")
377 call messages_fatal(1, namespace=namespace)
378 end if
379
380 end if
381
382 end if
383
384 call lda_u_write_info(this, namespace=namespace)
385
386 if (.not. this%basisfromstates) then
387
388 call orbitalbasis_init(this%basis, namespace, space%periodic_dim)
389
390 if (states_are_real(st)) then
391 call dorbitalbasis_build(this%basis, namespace, ions, gr, st%d%kpt, st%d%dim, &
392 this%skipSOrbitals, this%useAllOrbitals)
393 else
394 call zorbitalbasis_build(this%basis, namespace, ions, gr, st%d%kpt, st%d%dim, &
395 this%skipSOrbitals, this%useAllOrbitals)
396 end if
397 this%orbsets => this%basis%orbsets
398 this%norbsets = this%basis%norbsets
399 this%maxnorbs = this%basis%maxnorbs
400 this%max_np = this%basis%max_np
401 this%nspins = st%d%nspin
402 this%spin_channels = st%d%spin_channels
403 this%nspecies = ions%nspecies
404
405 !We allocate the necessary ressources
406 if (states_are_real(st)) then
407 call dlda_u_allocate(this, st)
408 else
409 call zlda_u_allocate(this, st)
410 end if
411
412 call distributed_nullify(this%orbs_dist, this%norbsets)
413#ifdef HAVE_MPI
414 if (.not. gr%parallel_in_domains) then
415 call distributed_init(this%orbs_dist, this%norbsets, mpi_comm_world, "orbsets")
416 end if
417#endif
418
419 else
420
421 !%Variable DFTUBasisStates
422 !%Type block
423 !%Default none
424 !%Section Hamiltonian::DFT+U
425 !%Description
426 !% This block starts by a line containing a single integer describing the number of
427 !% orbital sets. One orbital set is a group of orbitals on which one adds a Hubbard U.
428 !% Each following line of this block contains the index of a state to be used to construct the
429 !% localized basis, followed by the index of the corresponding orbital set.
430 !% See DFTUBasisFromStates for details.
431 !%End
432 if (parse_block(namespace, 'DFTUBasisStates', blk) == 0) then
433 call parse_block_integer(blk, 0, 0, this%norbsets)
434 this%maxnorbs = parse_block_n(blk)-1
435 if (this%maxnorbs <1) then
436 write(message(1),'(a,i3,a,i3)') 'DFTUBasisStates must contains at least one state.'
437 call messages_fatal(1, namespace=namespace)
438 end if
439 safe_allocate(this%basisstates(1:this%maxnorbs))
440 safe_allocate(this%basisstates_os(1:this%maxnorbs))
441 do is = 1, this%maxnorbs
442 call parse_block_integer(blk, is, 0, this%basisstates(is))
443 call parse_block_integer(blk, is, 1, this%basisstates_os(is))
444 end do
445 call parse_block_end(blk)
446 else
447 write(message(1),'(a,i3,a,i3)') 'DFTUBasisStates must be specified if DFTUBasisFromStates=yes'
448 call messages_fatal(1, namespace=namespace)
449 end if
450
451 if (states_are_real(st)) then
452 call dorbitalbasis_build_empty(this%basis, gr, st%d%dim, this%norbsets, this%basisstates_os)
453 else
454 call zorbitalbasis_build_empty(this%basis, gr, st%d%dim, this%norbsets, this%basisstates_os)
455 end if
456
457 this%max_np = gr%np
458 this%nspins = st%d%nspin
459 this%spin_channels = st%d%spin_channels
460 this%nspecies = 1
461
462 this%orbsets => this%basis%orbsets
463
464 call distributed_nullify(this%orbs_dist, this%norbsets)
465
466 !We allocate the necessary ressources
467 if (states_are_real(st)) then
468 call dlda_u_allocate(this, st)
469 else
470 call zlda_u_allocate(this, st)
471 end if
472
473 call lda_u_loadbasis(this, namespace, space, st, gr, mc, ierr)
474 if (ierr /= 0) then
475 message(1) = "Unable to load DFT+U basis from selected states."
476 call messages_fatal(1)
477 end if
478
479 end if
480
481 ! Symmetrization of the occupation matrices
482 this%symmetrize_occ_matrices = st%symmetrize_density .or. kpoints%use_symmetries
483 if (this%basisfromstates) this%symmetrize_occ_matrices = .false.
484 if (this%symmetrize_occ_matrices) then
485 call build_symmetrization_map(this, ions, gr)
486 end if
487
488 safe_allocate(this%dc_alpha(this%norbsets))
489 this%dc_alpha = m_one
490
491 call messages_print_with_emphasis(namespace=namespace)
492
493 pop_sub(lda_u_init)
494 end subroutine lda_u_init
495
496
497 ! ---------------------------------------------------------
498 subroutine lda_u_init_coulomb_integrals(this, namespace, space, gr, st, psolver)
499 type(lda_u_t), target, intent(inout) :: this
500 type(namespace_t), intent(in) :: namespace
501 class(space_t), intent(in) :: space
502 type(grid_t), intent(in) :: gr
503 type(states_elec_t), intent(in) :: st
504 type(poisson_t), intent(in) :: psolver
505
506 logical :: complex_coulomb_integrals
507 integer :: ios
508 integer :: norbs
509
511
512 norbs = this%maxnorbs
513
514 if (.not. this%basisfromstates) then
515
516 if (this%level == dft_u_acbn0) then
517
518 complex_coulomb_integrals = .false.
519 do ios = 1, this%norbsets
520 if (this%orbsets(ios)%ndim > 1) complex_coulomb_integrals = .true.
521 end do
522
523 if (.not. complex_coulomb_integrals) then
524 write(message(1),'(a)') 'Computing the Coulomb integrals of the localized basis.'
525 call messages_info(1, namespace=namespace)
526 safe_allocate(this%coulomb(1:norbs,1:norbs,1:norbs,1:norbs, 1:this%norbsets))
527 if (states_are_real(st)) then
528 call dcompute_coulomb_integrals(this, namespace, space, gr, psolver)
529 else
530 call zcompute_coulomb_integrals(this, namespace, space, gr, psolver)
531 end if
532 else
533 assert(.not. states_are_real(st))
534 write(message(1),'(a)') 'Computing complex Coulomb integrals of the localized basis.'
535 call messages_info(1, namespace=namespace)
536 safe_allocate(this%zcoulomb(1:norbs, 1:norbs, 1:norbs, 1:norbs, 1:st%d%dim, 1:st%d%dim, 1:this%norbsets))
537 call compute_complex_coulomb_integrals(this, gr, st, psolver, namespace, space)
538 end if
539 end if
540
541 else
542 write(message(1),'(a)') 'Computing the Coulomb integrals of the localized basis.'
543 call messages_info(1, namespace=namespace)
544 safe_allocate(this%coulomb(1:norbs, 1:norbs, 1:norbs, 1:norbs, 1:this%norbsets))
545 if (states_are_real(st)) then
546 call dcompute_coulomb_integrals(this, namespace, space, gr, psolver)
547 else
548 call zcompute_coulomb_integrals(this, namespace, space, gr, psolver)
549 end if
550
551 end if
552
554
555 end subroutine lda_u_init_coulomb_integrals
556
557
558 ! ---------------------------------------------------------
559 subroutine lda_u_end(this)
560 implicit none
561 type(lda_u_t), intent(inout) :: this
562
563 push_sub(lda_u_end)
564
565 this%level = dft_u_none
566
567 safe_deallocate_a(this%dn)
568 safe_deallocate_a(this%zn)
569 safe_deallocate_a(this%dn_alt)
570 safe_deallocate_a(this%zn_alt)
571 safe_deallocate_a(this%dV)
572 safe_deallocate_a(this%zV)
573 safe_deallocate_a(this%coulomb)
574 safe_deallocate_a(this%zcoulomb)
575 safe_deallocate_a(this%renorm_occ)
576 safe_deallocate_a(this%dn_ij)
577 safe_deallocate_a(this%zn_ij)
578 safe_deallocate_a(this%dn_alt_ij)
579 safe_deallocate_a(this%zn_alt_ij)
580 safe_deallocate_a(this%dn_alt_ii)
581 safe_deallocate_a(this%zn_alt_ii)
582 safe_deallocate_a(this%basisstates)
583 safe_deallocate_a(this%basisstates_os)
584 safe_deallocate_a(this%dc_alpha)
585 safe_deallocate_a(this%inv_map_symm)
586 safe_deallocate_a(this%symm_weight)
587
588 nullify(this%orbsets)
589 call orbitalbasis_end(this%basis)
590
591 this%max_np = 0
592
593 if (.not. this%basisfromstates) then
594 call distributed_end(this%orbs_dist)
595 end if
596
597 pop_sub(lda_u_end)
598 end subroutine lda_u_end
599
600 ! When moving the ions, the basis must be reconstructed
601 subroutine lda_u_update_basis(this, space, gr, ions, st, psolver, namespace, kpoints, has_phase)
602 type(lda_u_t), target, intent(inout) :: this
603 class(space_t), intent(in) :: space
604 type(grid_t), intent(in) :: gr
605 type(ions_t), target, intent(in) :: ions
606 type(states_elec_t), intent(in) :: st
607 type(poisson_t), intent(in) :: psolver
608 type(namespace_t), intent(in) :: namespace
609 type(kpoints_t), intent(in) :: kpoints
610 logical, intent(in) :: has_phase
611
612 integer :: ios, maxorbs, nspin
613
614 if(this%level == dft_u_none) return
615
616 push_sub(lda_u_update_basis)
617
618 if(.not. this%basisfromstates) then
619 !We clean the orbital basis, to be able to reconstruct it
620 call orbitalbasis_end(this%basis)
621 nullify(this%orbsets)
622
623 !We now reconstruct the basis
624 if (states_are_real(st)) then
625 call dorbitalbasis_build(this%basis, namespace, ions, gr, st%d%kpt, st%d%dim, &
626 this%skipSOrbitals, this%useAllOrbitals, verbose = .false.)
627 else
628 call zorbitalbasis_build(this%basis, namespace, ions, gr, st%d%kpt, st%d%dim, &
629 this%skipSOrbitals, this%useAllOrbitals, verbose = .false.)
630 end if
631 this%orbsets => this%basis%orbsets
632 this%max_np = this%basis%max_np
633 end if
634
635 !In case of intersite interaction we need to reconstruct the basis
636 if (this%intersite) then
637 this%maxneighbors = 0
638 do ios = 1, this%norbsets
639 call orbitalset_init_intersite(this%orbsets(ios), namespace, space, ios, ions, gr%der, psolver, &
640 this%orbsets, this%norbsets, this%maxnorbs, this%intersite_radius, st%d%kpt, has_phase, &
641 this%sm_poisson, this%basisfromstates, this%basis%combine_j_orbitals)
642 this%maxneighbors = max(this%maxneighbors, this%orbsets(ios)%nneighbors)
643 end do
644
645 maxorbs = this%maxnorbs
646 nspin = this%nspins
647
648 if (states_are_real(st)) then
649 safe_deallocate_a(this%dn_ij)
650 safe_allocate(this%dn_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
651 this%dn_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_zero
652 safe_deallocate_a(this%dn_alt_ij)
653 safe_allocate(this%dn_alt_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
654 this%dn_alt_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_zero
655 safe_deallocate_a(this%dn_alt_ii)
656 safe_allocate(this%dn_alt_ii(1:2,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
657 this%dn_alt_ii(1:2,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_zero
658 else
659 safe_deallocate_a(this%zn_ij)
660 safe_allocate(this%zn_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
661 this%zn_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_z0
662 safe_deallocate_a(this%zn_alt_ij)
663 safe_allocate(this%zn_alt_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
664 this%zn_alt_ij(1:maxorbs,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_z0
665 safe_deallocate_a(this%zn_alt_ii)
666 safe_allocate(this%zn_alt_ii(1:2,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors))
667 this%zn_alt_ii(1:2,1:maxorbs,1:nspin,1:this%norbsets,1:this%maxneighbors) = m_z0
668 end if
669 end if
670
671 ! We rebuild the phase for the orbital projection, similarly to the one of the pseudopotentials
672 ! In case of a laser field, the phase is recomputed in hamiltonian_elec_update
673 if (has_phase) then
674 call lda_u_build_phase_correction(this, space, st%d, gr%der%boundaries, namespace, kpoints)
675 else
676 if(.not. this%basisfromstates) then
677 !In case there is no phase, we perform the orthogonalization here
678 if(this%basis%orthogonalization) then
679 call dloewdin_orthogonalize(this%basis, st%d%kpt, namespace)
680 else
681 if(debug%info .and. space%is_periodic()) then
682 call dloewdin_info(this%basis, st%d%kpt, namespace)
683 end if
684 end if
685 end if
686 end if
687
688 ! Rebuild the Coulomb integrals
689 if (allocated(this%coulomb)) then
690 safe_deallocate_a(this%coulomb)
691 end if
692 if (allocated(this%zcoulomb)) then
693 safe_deallocate_a(this%zcoulomb)
694 end if
695 call lda_u_init_coulomb_integrals(this, namespace, space, gr, st, psolver)
697 pop_sub(lda_u_update_basis)
698
699 end subroutine lda_u_update_basis
700
701 ! Interface for the X(update_occ_matrices) routines
702 subroutine lda_u_update_occ_matrices(this, namespace, mesh, st, phase, energy)
703 type(lda_u_t), intent(inout) :: this
704 type(namespace_t), intent(in) :: namespace
705 class(mesh_t), intent(in) :: mesh
706 type(states_elec_t), intent(inout) :: st
707 type(phase_t), intent(in) :: phase
708 type(energy_t), intent(inout) :: energy
709
710 if (this%level == dft_u_none .or. this%freeze_occ) return
712
713 if (states_are_real(st)) then
714 call dupdate_occ_matrices(this, namespace, mesh, st, energy%dft_u)
715 else
716 if (phase%is_allocated()) then
717 call zupdate_occ_matrices(this, namespace, mesh, st, energy%dft_u, phase)
718 else
719 call zupdate_occ_matrices(this, namespace, mesh, st, energy%dft_u)
720 end if
721 end if
722
724 end subroutine lda_u_update_occ_matrices
725
726
728 subroutine lda_u_build_phase_correction(this, space, std, boundaries, namespace, kpoints, vec_pot, vec_pot_var)
729 type(lda_u_t), intent(inout) :: this
730 class(space_t), intent(in) :: space
731 type(states_elec_dim_t), intent(in) :: std
732 type(boundaries_t), intent(in) :: boundaries
733 type(namespace_t), intent(in) :: namespace
734 type(kpoints_t), intent(in) :: kpoints
735 real(real64), optional, allocatable, intent(in) :: vec_pot(:)
736 real(real64), optional, allocatable, intent(in) :: vec_pot_var(:, :)
737
738 integer :: ios
739
740 if (boundaries%spiralBC) call messages_not_implemented("DFT+U with spiral boundary conditions", &
741 namespace=namespace)
742
744
745 write(message(1), '(a)') 'Debug: Building the phase correction for DFT+U orbitals.'
746 call messages_info(1, namespace=namespace, debug_only=.true.)
747
748 do ios = 1, this%norbsets
749 call orbitalset_update_phase(this%orbsets(ios), space%dim, std%kpt, kpoints, &
750 (std%ispin==spin_polarized), vec_pot, vec_pot_var)
751 call orbitalset_update_phase_shift(this%orbsets(ios), space%dim, std%kpt, kpoints, &
752 (std%ispin==spin_polarized), vec_pot, vec_pot_var)
753 end do
754
755 if (.not. this%basisfromstates) then
756 if (this%basis%orthogonalization) then
757 call zloewdin_orthogonalize(this%basis, std%kpt, namespace)
758 else
759 if (debug%info .and. space%is_periodic()) call zloewdin_info(this%basis, std%kpt, namespace)
760 end if
761 end if
762
764
765 end subroutine lda_u_build_phase_correction
766
767 ! ---------------------------------------------------------
768 subroutine compute_acbno_u_kanamori(this, st, kanamori)
769 type(lda_u_t), intent(in) :: this
770 type(states_elec_t), intent(in) :: st
771 real(real64), intent(out) :: kanamori(:,:)
772
773 if (this%nspins == 1) then
774 if (states_are_real(st)) then
775 call dcompute_acbno_u_kanamori_restricted(this, kanamori)
776 else
777 call zcompute_acbno_u_kanamori_restricted(this, kanamori)
778 end if
779 else
780 if (states_are_real(st)) then
781 call dcompute_acbno_u_kanamori(this, kanamori)
782 else
783 call zcompute_acbno_u_kanamori(this, kanamori)
784 end if
785 end if
786
787
788 end subroutine compute_acbno_u_kanamori
789
790 ! ---------------------------------------------------------
791 subroutine lda_u_freeze_occ(this)
792 type(lda_u_t), intent(inout) :: this
793
794 this%freeze_occ = .true.
795 end subroutine lda_u_freeze_occ
796
797 ! ---------------------------------------------------------
798 subroutine lda_u_freeze_u(this)
799 type(lda_u_t), intent(inout) :: this
800
801 this%freeze_u = .true.
802 end subroutine lda_u_freeze_u
803
804 ! ---------------------------------------------------------
805 subroutine lda_u_set_effectiveu(this, Ueff)
806 type(lda_u_t), intent(inout) :: this
807 real(real64), intent(in) :: ueff(:)
808
809 integer :: ios
810
811 push_sub(lda_u_set_effectiveu)
812
813 do ios = 1,this%norbsets
814 this%orbsets(ios)%Ueff = ueff(ios)
815 end do
816
817 pop_sub(lda_u_set_effectiveu)
818 end subroutine lda_u_set_effectiveu
819
820 ! ---------------------------------------------------------
821 subroutine lda_u_get_effectiveu(this, Ueff)
822 type(lda_u_t), intent(in) :: this
823 real(real64), intent(inout) :: ueff(:)
824
825 integer :: ios
826
827 push_sub(lda_u_get_effectiveu)
828
829 do ios = 1,this%norbsets
830 ueff(ios) = this%orbsets(ios)%Ueff
831 end do
832
833 pop_sub(lda_u_get_effectiveu)
834 end subroutine lda_u_get_effectiveu
835
836 ! ---------------------------------------------------------
837 subroutine lda_u_set_effectivev(this, Veff)
838 type(lda_u_t), intent(inout) :: this
839 real(real64), intent(in) :: veff(:)
840
841 integer :: ios, ncount
842
843 push_sub(lda_u_set_effectivev)
844
845 ncount = 0
846 do ios = 1, this%norbsets
847 this%orbsets(ios)%V_ij(1:this%orbsets(ios)%nneighbors,0) = veff(ncount+1:ncount+this%orbsets(ios)%nneighbors)
848 ncount = ncount + this%orbsets(ios)%nneighbors
849 end do
850
851 pop_sub(lda_u_set_effectivev)
852 end subroutine lda_u_set_effectivev
853
854 ! ---------------------------------------------------------
855 subroutine lda_u_get_effectivev(this, Veff)
856 type(lda_u_t), intent(in) :: this
857 real(real64), intent(inout) :: veff(:)
858
859 integer :: ios, ncount
860
861 push_sub(lda_u_get_effectivev)
862
863 ncount = 0
864 do ios = 1, this%norbsets
865 veff(ncount+1:ncount+this%orbsets(ios)%nneighbors) = this%orbsets(ios)%V_ij(1:this%orbsets(ios)%nneighbors,0)
866 ncount = ncount + this%orbsets(ios)%nneighbors
867 end do
868
869 pop_sub(lda_u_get_effectivev)
870 end subroutine lda_u_get_effectivev
871
872 ! ---------------------------------------------------------
873 subroutine lda_u_write_info(this, iunit, namespace)
874 type(lda_u_t), intent(in) :: this
875 integer, optional, intent(in) :: iunit
876 type(namespace_t), optional, intent(in) :: namespace
877
878 push_sub(lda_u_write_info)
879
880 write(message(1), '(1x)')
881 call messages_info(1, iunit=iunit, namespace=namespace)
882 if (this%level == dft_u_empirical) then
883 write(message(1), '(a)') "Method:"
884 write(message(2), '(a)') " [1] Dudarev et al., Phys. Rev. B 57, 1505 (1998)"
885 call messages_info(2, iunit=iunit, namespace=namespace)
886 else
887 if (.not. this%intersite) then
888 write(message(1), '(a)') "Method:"
889 write(message(2), '(a)') " [1] Agapito et al., Phys. Rev. X 5, 011006 (2015)"
890 else
891 write(message(1), '(a)') "Method:"
892 write(message(2), '(a)') " [1] Tancogne-Dejean, and Rubio, Phys. Rev. B 102, 155117 (2020)"
893 end if
894 call messages_info(2, iunit=iunit, namespace=namespace)
895 end if
896 write(message(1), '(a)') "Implementation:"
897 write(message(2), '(a)') " [1] Tancogne-Dejean, Oliveira, and Rubio, Phys. Rev. B 69, 245133 (2017)"
898 write(message(3), '(1x)')
899 call messages_info(3, iunit=iunit, namespace=namespace)
901 pop_sub(lda_u_write_info)
902
903 end subroutine lda_u_write_info
904
905 ! ---------------------------------------------------------
906 subroutine lda_u_loadbasis(this, namespace, space, st, mesh, mc, ierr)
907 type(lda_u_t), intent(inout) :: this
908 type(namespace_t), intent(in) :: namespace
909 class(space_t), intent(in) :: space
910 type(states_elec_t), intent(in) :: st
911 class(mesh_t), intent(in) :: mesh
912 type(multicomm_t), intent(in) :: mc
913 integer, intent(out) :: ierr
914
915 integer :: err, wfns_file, is, ist, idim, ik, ios, iorb
917 character(len=256) :: lines(3)
918 character(len=256), allocatable :: restart_file(:, :)
919 logical, allocatable :: restart_file_present(:, :)
920 character(len=12) :: filename
921 character(len=1) :: char
922 character(len=50) :: str
923 type(orbitalset_t), pointer :: os
924 integer, allocatable :: count(:)
925 real(real64) :: norm, center(space%dim)
926 real(real64), allocatable :: dpsi(:,:,:)
927 complex(real64), allocatable :: zpsi(:,:,:)
928
929
930 push_sub(lda_u_loadbasis)
931
932 ierr = 0
933
934 message(1) = "Debug: Loading DFT+U basis from states."
935 call messages_info(1, debug_only=.true.)
936
937 call restart_gs%init(namespace, restart_proj, restart_type_load, mc, err, mesh=mesh)
938
939 ! If any error occured up to this point then it is not worth continuing,
940 ! as there something fundamentally wrong with the restart files
941 if (err /= 0) then
942 call restart_gs%end()
943 message(1) = "Error loading DFT+U basis from states, cannot proceed with the calculation"
944 call messages_fatal(1)
945 pop_sub(lda_u_loadbasis)
946 return
947 end if
948
949 ! open files to read
950 wfns_file = restart_gs%open('wfns')
951 call restart_gs%read(wfns_file, lines, 2, err)
952 if (err /= 0) then
953 ierr = ierr - 2**5
954 else if (states_are_real(st)) then
955 read(lines(2), '(a)') str
956 if (str(2:8) == 'Complex') then
957 message(1) = "Cannot read real states from complex wavefunctions."
958 call messages_fatal(1, namespace=namespace)
959 else if (str(2:5) /= 'Real') then
960 message(1) = "Restart file 'wfns' does not specify real/complex; cannot check compatibility."
961 call messages_warning(1, namespace=namespace)
962 end if
963 end if
964 ! complex can be restarted from real, so there is no problem.
965
966 ! If any error occured up to this point then it is not worth continuing,
967 ! as there something fundamentally wrong with the restart files
968 if (err /= 0) then
969 call restart_gs%close(wfns_file)
970 call restart_gs%end()
971 message(1) = "Error loading DFT+U basis from states, cannot proceed with the calculation"
972 call messages_fatal(1)
973 pop_sub(lda_u_loadbasis)
974 return
975 end if
976
977 safe_allocate(restart_file(1:st%d%dim, 1:st%nst))
978 safe_allocate(restart_file_present(1:st%d%dim, 1:st%nst))
979 restart_file_present = .false.
980
981 ! Next we read the list of states from the files.
982 ! Errors in reading the information of a specific state from the files are ignored
983 ! at this point, because later we will skip reading the wavefunction of that state.
984 do
985 call restart_gs%read(wfns_file, lines, 1, err)
986 if (err == 0) then
987 read(lines(1), '(a)') char
988 if (char == '%') then
989 !We reached the end of the file
990 exit
991 else
992 read(lines(1), *) ik, char, ist, char, idim, char, filename
993 end if
994 end if
995
996 if (any(this%basisstates==ist) .and. ik == 1) then
997 restart_file(idim, ist) = trim(filename)
998 restart_file_present(idim, ist) = .true.
999 end if
1000 end do
1001 call restart_gs%close(wfns_file)
1002
1003 !We loop over the states we need
1004 safe_allocate(count(1:this%norbsets))
1005 count = 0
1006 do is = 1, this%maxnorbs
1007 ist = this%basisstates(is)
1008 ios = this%basisstates_os(is)
1009 count(ios) = count(ios)+1
1010 do idim = 1, st%d%dim
1011
1012 if (.not. restart_file_present(idim, ist)) then
1013 write(message(1), '(a,i3,a)') "Cannot read states ", ist, "from the projection folder"
1014 call messages_fatal(1, namespace=namespace)
1015 end if
1016
1017 if (states_are_real(st)) then
1018 call restart_gs%read_mesh_function(restart_file(idim, ist), mesh, &
1019 this%orbsets(ios)%dorb(:,idim,count(ios)), err)
1020 else
1021 call restart_gs%read_mesh_function(restart_file(idim, ist), mesh, &
1022 this%orbsets(ios)%zorb(:,idim,count(ios)), err)
1023 end if
1024
1025 end do
1026 end do
1027 safe_deallocate_a(count)
1028 safe_deallocate_a(restart_file)
1029 safe_deallocate_a(restart_file_present)
1030 call restart_gs%end()
1031
1032 ! Normalize the orbitals. This is important if we use Wannier orbitals instead of KS states
1033 if(this%basis%normalize) then
1034 do ios = 1, this%norbsets
1035 do iorb = 1, this%orbsets(ios)%norbs
1036 if (states_are_real(st)) then
1037 norm = dmf_nrm2(mesh, st%d%dim, this%orbsets(ios)%dorb(:,:,iorb))
1038 call lalg_scal(mesh%np, st%d%dim, m_one/norm, this%orbsets(ios)%dorb(:,:,iorb))
1039 else
1040 norm = zmf_nrm2(mesh, st%d%dim, this%orbsets(ios)%zorb(:,:,iorb))
1041 call lalg_scal(mesh%np, st%d%dim, m_one/norm, this%orbsets(ios)%zorb(:,:,iorb))
1042 end if
1043 end do
1044 end do
1045 end if
1046
1047 ! We rotate the orbitals in the complex plane to have them as close as possible to real functions
1048 if(states_are_complex(st) .and. st%d%dim == 1) then
1049 do ios = 1, this%norbsets
1050 do iorb = 1, this%orbsets(ios)%norbs
1051 call zmf_fix_phase(mesh, this%orbsets(ios)%zorb(:,1,iorb))
1052 end do
1053 end do
1054 end if
1055
1056 ! We determine the center of charge by computing <w|r|w>
1057 ! We could also determine the spread by \Omega = <w|r^2|w> - <w|r|w>^2
1058 do ios = 1, this%norbsets
1059 if (states_are_real(st)) then
1060 call dorbitalset_get_center_of_mass(this%orbsets(ios), space, mesh, this%latt)
1061 else
1062 call zorbitalset_get_center_of_mass(this%orbsets(ios), space, mesh, this%latt)
1063 end if
1064 end do
1065
1066 message(1) = "Debug: Converting the Wannier states to submeshes."
1067 call messages_info(1, debug_only=.true.)
1068
1069 ! We now transfer the states to a submesh centered on the center of mass of the Wannier orbitals
1070 this%max_np = 0
1071 do ios = 1, this%norbsets
1072 os => this%orbsets(ios)
1073 center = os%sphere%center
1074 safe_deallocate_a(os%sphere%center)
1075 if (states_are_real(st)) then
1076 safe_allocate(dpsi(1:mesh%np, 1:os%ndim, 1:os%norbs))
1077 dpsi(1:mesh%np, 1:os%ndim, 1:os%norbs) = os%dorb(1:mesh%np, 1:os%ndim, 1:os%norbs)
1078
1079 safe_deallocate_a(os%dorb)
1080 !We initialise the submesh corresponding to the orbital
1081 call submesh_init(os%sphere, space, mesh, this%latt, center, os%radius)
1082 safe_allocate(os%dorb(1:os%sphere%np, 1:os%ndim, 1:os%norbs))
1083 do iorb = 1, os%norbs
1084 do idim = 1, os%ndim
1085 call dsubmesh_copy_from_mesh(os%sphere, dpsi(:,idim,iorb), os%dorb(:,idim, iorb))
1086 end do
1087 end do
1088 safe_deallocate_a(dpsi)
1089 else
1090 safe_allocate(zpsi(1:mesh%np, 1:os%ndim, 1:os%norbs))
1091 zpsi(1:mesh%np, 1:os%ndim, 1:os%norbs) = os%zorb(1:mesh%np, 1:os%ndim, 1:os%norbs)
1092 safe_deallocate_a(os%zorb)
1093 !We initialise the submesh corresponding to the orbital
1094 call submesh_init(os%sphere, space, mesh, this%latt, center, os%radius)
1095 safe_allocate(os%zorb(1:os%sphere%np, 1:os%ndim, 1:os%norbs))
1096 do iorb = 1, os%norbs
1097 do idim = 1, os%ndim
1098 call zsubmesh_copy_from_mesh(os%sphere, zpsi(:,idim,iorb), os%zorb(:,idim, iorb))
1099 end do
1100 end do
1101 safe_deallocate_a(zpsi)
1102
1103 safe_allocate(os%phase(1:os%sphere%np, st%d%kpt%start:st%d%kpt%end))
1104 safe_allocate(os%eorb_submesh(1:os%sphere%np, 1:os%ndim, 1:os%norbs, st%d%kpt%start:st%d%kpt%end))
1105 end if
1106 os%use_submesh = .true. ! We are now on a submesh
1107 this%max_np = max(this%max_np, os%sphere%np)
1108 end do
1109
1110 this%basis%use_submesh = .true.
1111
1112 ! If we use GPUs, we need to transfert the orbitals on the device
1113 if (accel_is_enabled() .and. st%d%dim == 1) then
1114 do ios = 1, this%norbsets
1115 os => this%orbsets(ios)
1116
1117 os%ldorbs = max(accel_padded_size(os%sphere%np), 1)
1118 if (states_are_real(st)) then
1119 call accel_create_buffer(os%dbuff_orb, accel_mem_read_only, type_float, os%ldorbs*os%norbs)
1120 else
1121 call accel_create_buffer(os%zbuff_orb, accel_mem_read_only, type_cmplx, os%ldorbs*os%norbs)
1122 safe_allocate(os%buff_eorb(st%d%kpt%start:st%d%kpt%end))
1123
1124 do ik= st%d%kpt%start, st%d%kpt%end
1125 call accel_create_buffer(os%buff_eorb(ik), accel_mem_read_only, type_cmplx, os%ldorbs*os%norbs)
1126 end do
1127 end if
1128
1129 call accel_create_buffer(os%sphere%buff_map, accel_mem_read_only, type_integer, max(os%sphere%np, 1))
1130 call accel_write_buffer(os%sphere%buff_map, os%sphere%np, os%sphere%map)
1131
1132 do iorb = 1, os%norbs
1133 if(states_are_complex(st)) then
1134 call accel_write_buffer(os%zbuff_orb, os%sphere%np, os%zorb(:, 1, iorb), &
1135 offset = (iorb - 1)*os%ldorbs)
1136 else
1137 call accel_write_buffer(os%dbuff_orb, os%sphere%np, os%dorb(:, 1, iorb), &
1138 offset = (iorb - 1)*os%ldorbs)
1139 end if
1140 end do
1141 end do
1142 end if
1143
1144
1145 message(1) = "Debug: Loading DFT+U basis from states done."
1146 call messages_info(1, debug_only=.true.)
1147
1148 pop_sub(lda_u_loadbasis)
1149 end subroutine lda_u_loadbasis
1150
1152 subroutine build_symmetrization_map(this, ions, gr)
1153 type(lda_u_t), intent(inout) :: this
1154 type(ions_t), intent(in) :: ions
1155 type(grid_t), intent(in) :: gr
1156
1157 integer :: nsym, iop, ios, iatom, iatom_sym, ios_sym
1158
1159 push_sub(build_symmetrization_map)
1160
1161 nsym = ions%symm%nops
1162 safe_allocate(this%inv_map_symm(1:this%norbsets, 1:nsym))
1163 this%inv_map_symm = -1
1164
1165 this%nsym = nsym
1166
1167 do ios = 1, this%norbsets
1168 iatom = this%orbsets(ios)%iatom
1169 do iop = 1, nsym
1170 iatom_sym = ions%inv_map_symm_atoms(iatom, iop)
1171
1172 do ios_sym = 1, this%norbsets
1173 if (this%orbsets(ios_sym)%iatom == iatom_sym .and. this%orbsets(ios_sym)%norbs == this%orbsets(ios)%norbs &
1174 .and. this%orbsets(ios_sym)%nn == this%orbsets(ios)%nn .and. this%orbsets(ios_sym)%ll == this%orbsets(ios)%ll &
1175 .and. is_close(this%orbsets(ios_sym)%jj, this%orbsets(ios)%jj)) then
1176 this%inv_map_symm(ios, iop) = ios_sym
1177 exit
1178 end if
1179 end do
1180 assert(this%inv_map_symm(ios, iop) > 0)
1181 end do
1182 end do
1183
1184 safe_allocate(this%symm_weight(1:this%maxnorbs, 1:this%maxnorbs, 1:this%nsym, 1:this%norbsets))
1185 this%symm_weight = m_zero
1186
1187 do ios = 1, this%norbsets
1188 ! s-orbitals
1189 if (this%orbsets(ios)%norbs == 1) then
1190 this%symm_weight(1,1, 1:this%nsym, ios) = m_one
1191 cycle
1192 end if
1193
1194 ! Not implemented yet
1195 if (this%orbsets(ios)%ndim > 1) cycle
1196
1197 call orbitals_get_symm_weight(this%orbsets(ios), ions%space, ions%latt, gr, ions%symm, this%symm_weight(:,:,:,ios))
1198 end do
1199
1201 end subroutine build_symmetrization_map
1202
1204 subroutine orbitals_get_symm_weight(os, space, latt, gr, symm, weight)
1205 type(orbitalset_t), intent(in) :: os
1206 type(space_t), intent(in) :: space
1207 type(lattice_vectors_t), intent(in) :: latt
1208 type(grid_t), intent(in) :: gr
1209 type(symmetries_t), intent(in) :: symm
1210 real(real64), intent(inout) :: weight(:,:,:)
1211
1212 integer :: im, imp, iop, mm
1213 real(real64), allocatable :: orb(:,:), orb_sym(:), ylm(:)
1214 type(submesh_t) :: sphere
1215 real(real64) :: rc, norm, origin(space%dim)
1216
1217 push_sub(orbitals_get_symm_weight)
1218
1219 assert(os%ndim == 1)
1220
1221 safe_allocate(orb_sym(1:gr%np))
1222 safe_allocate(orb(1:gr%np, 1:os%norbs))
1223
1224 assert(2*os%ll+1 == os%norbs)
1225
1226 ! We generate an artificial submesh to compute the symmetries on it
1227 ! The radius is such that we fit in 20 points
1228 rc = (50.0_real64 * m_three/m_four/m_pi*product(gr%spacing))**m_third
1229 origin = m_zero
1230 call submesh_init(sphere, space, gr, latt, origin, rc)
1231
1232 safe_allocate(ylm(1:sphere%np))
1233
1234 ! We then compute the spherical harmonics in this submesh
1235 do im = 1, os%norbs
1236 mm = im-1-os%ll
1237 call loct_ylm(sphere%np, sphere%rel_x(1,1), sphere%r(1), os%ll, mm, ylm(1))
1238 orb(:,im) = m_zero
1239 call submesh_add_to_mesh(sphere, ylm, orb(:,im))
1240 norm = dmf_nrm2(gr, orb(:,im))
1241 call lalg_scal(gr%np, m_one / norm, orb(:,im))
1242 end do
1243 safe_deallocate_a(ylm)
1244
1245 ! Then we put them on the grid, rotate them
1246 do im = 1, os%norbs
1247 do iop = 1, symm%nops
1248 call dgrid_symmetrize_single(gr, iop, orb(:,im), orb_sym)
1249
1250 do imp = 1, os%norbs
1251 weight(im, imp, iop) = dmf_dotp(gr, orb(:,imp), orb_sym, reduce=.false.)
1252 end do
1253 end do
1254 end do
1255
1256 call gr%allreduce(weight)
1257
1258 safe_deallocate_a(orb)
1259 safe_deallocate_a(orb_sym)
1260
1262 end subroutine orbitals_get_symm_weight
1263
1264#include "dft_u_noncollinear_inc.F90"
1265
1266#include "undef.F90"
1267#include "real.F90"
1268#include "lda_u_inc.F90"
1269
1270#include "undef.F90"
1271#include "complex.F90"
1272#include "lda_u_inc.F90"
1273end module lda_u_oct_m
scales a vector by a constant
Definition: lalg_basic.F90:159
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:182
pure logical function, public accel_is_enabled()
Definition: accel.F90:419
integer, parameter, public accel_mem_read_only
Definition: accel.F90:196
This module implements batches of mesh functions.
Definition: batch.F90:135
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:118
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:124
type(debug_t), save, public debug
Definition: debug.F90:158
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public distributed_end(this)
subroutine, public distributed_nullify(this, total)
subroutine, public distributed_init(this, total, comm, tag, scalapack_compat)
Distribute N instances across M processes of communicator comm
integer, parameter, public spinors
integer, parameter, public spin_polarized
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_four
Definition: global.F90:195
real(real64), parameter, public m_third
Definition: global.F90:198
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
complex(real64), parameter, public m_z0
Definition: global.F90:201
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_one
Definition: global.F90:192
real(real64), parameter, public m_three
Definition: global.F90:194
This module implements the underlying real-space grid.
Definition: grid.F90:119
subroutine, public dgrid_symmetrize_single(gr, iop, field, symm_field, suppress_warning)
Definition: grid.F90:724
integer, parameter, public dft_u_amf
Definition: lda_u.F90:208
subroutine, public zlda_u_commute_r(this, mesh, space, d, namespace, psib, gpsib)
This routine computes [r,V_lda+u] .
Definition: lda_u.F90:4984
subroutine zlda_u_allocate(this, st)
Definition: lda_u.F90:5554
subroutine, public lda_u_get_effectiveu(this, Ueff)
Definition: lda_u.F90:917
subroutine compute_complex_coulomb_integrals(this, gr, st, psolver, namespace, space)
Definition: lda_u.F90:1381
subroutine, public dlda_u_force(this, namespace, space, mesh, st, iq, psib, grad_psib, force)
Definition: lda_u.F90:3169
subroutine, public dlda_u_apply(this, d, mesh, psib, hpsib)
Apply the +U nonlocal potential to psib and adds the result to hpsib.
Definition: lda_u.F90:1749
subroutine, public lda_u_set_effectiveu(this, Ueff)
Definition: lda_u.F90:901
subroutine build_symmetrization_map(this, ions, gr)
Builds a mapping between the orbital sets based on symmetries.
Definition: lda_u.F90:1248
subroutine dcompute_acbno_u_kanamori_restricted(this, kanamori)
This routine computes the Kanamori U, Up, and J.
Definition: lda_u.F90:2770
subroutine, public lda_u_init(this, namespace, space, level, gr, ions, st, mc, kpoints)
Definition: lda_u.F90:284
subroutine, public dlda_u_get_occupations(this, occ)
Definition: lda_u.F90:3474
subroutine, public dlda_u_rvu(this, mesh, space, d, namespace, psib, gpsib)
This routine computes .
Definition: lda_u.F90:3307
subroutine, public lda_u_update_basis(this, space, gr, ions, st, psolver, namespace, kpoints, has_phase)
Definition: lda_u.F90:697
subroutine orbitals_get_symm_weight(os, space, latt, gr, symm, weight)
Computes the weight of each rotated orbitals in the basis of the same localized subspace.
Definition: lda_u.F90:1300
subroutine dupdate_occ_matrices(this, namespace, mesh, st, lda_u_energy, phase)
This routine computes the values of the occupation matrices.
Definition: lda_u.F90:1809
subroutine, public zlda_u_apply(this, d, mesh, psib, hpsib)
Apply the +U nonlocal potential to psib and adds the result to hpsib.
Definition: lda_u.F90:3721
subroutine, public zcompute_dftu_energy(this, energy, st)
This routine computes the value of the double counting term in the DFT+U energy.
Definition: lda_u.F90:4182
integer, parameter, public dft_u_empirical
Definition: lda_u.F90:203
subroutine zcompute_acbno_u_kanamori_restricted(this, kanamori)
This routine computes the Kanamori U, Up, and J.
Definition: lda_u.F90:4742
subroutine dcompute_acbno_u_kanamori(this, kanamori)
This routine computes the Kanamori U, Up, and J.
Definition: lda_u.F90:2671
subroutine, public zlda_u_commute_r_single(this, mesh, space, d, namespace, ist, ik, psi, gpsi, has_phase)
Definition: lda_u.F90:4933
subroutine, public lda_u_freeze_occ(this)
Definition: lda_u.F90:887
integer, parameter, public dft_u_mix
Definition: lda_u.F90:208
subroutine, public lda_u_freeze_u(this)
Definition: lda_u.F90:894
subroutine, public zlda_u_rvu(this, mesh, space, d, namespace, psib, gpsib)
This routine computes .
Definition: lda_u.F90:5309
subroutine dlda_u_allocate(this, st)
Definition: lda_u.F90:3531
subroutine, public compute_acbno_u_kanamori(this, st, kanamori)
Definition: lda_u.F90:864
subroutine, public zlda_u_update_potential(this, st)
This routine computes the potential that, once multiplied by the projector Pmm' and summed over m and...
Definition: lda_u.F90:4256
subroutine lda_u_init_coulomb_integrals(this, namespace, space, gr, st, psolver)
Definition: lda_u.F90:594
subroutine, public lda_u_write_info(this, iunit, namespace)
Definition: lda_u.F90:969
subroutine, public dlda_u_commute_r(this, mesh, space, d, namespace, psib, gpsib)
This routine computes [r,V_lda+u] .
Definition: lda_u.F90:3012
subroutine, public zlda_u_force(this, namespace, space, mesh, st, iq, psib, grad_psib, force)
Definition: lda_u.F90:5171
subroutine dcompute_coulomb_integrals(this, namespace, space, gr, psolver)
Definition: lda_u.F90:2845
subroutine, public dlda_u_set_occupations(this, occ)
Definition: lda_u.F90:3417
subroutine, public lda_u_get_effectivev(this, Veff)
Definition: lda_u.F90:951
subroutine lda_u_loadbasis(this, namespace, space, st, mesh, mc, ierr)
Definition: lda_u.F90:1002
subroutine, public dlda_u_commute_r_single(this, mesh, space, d, namespace, ist, ik, psi, gpsi, has_phase)
Definition: lda_u.F90:2961
subroutine, public lda_u_build_phase_correction(this, space, std, boundaries, namespace, kpoints, vec_pot, vec_pot_var)
Build the phase correction to the global phase for all orbitals.
Definition: lda_u.F90:824
subroutine, public lda_u_end(this)
Definition: lda_u.F90:655
subroutine, public lda_u_set_effectivev(this, Veff)
Definition: lda_u.F90:933
subroutine, public lda_u_update_occ_matrices(this, namespace, mesh, st, phase, energy)
Definition: lda_u.F90:798
integer, parameter, public dft_u_acbn0
Definition: lda_u.F90:203
subroutine zcompute_acbno_u_kanamori(this, kanamori)
This routine computes the Kanamori U, Up, and J.
Definition: lda_u.F90:4643
subroutine zcompute_coulomb_integrals(this, namespace, space, gr, psolver)
Definition: lda_u.F90:4817
subroutine, public zlda_u_get_occupations(this, occ)
Definition: lda_u.F90:5497
subroutine zupdate_occ_matrices(this, namespace, mesh, st, lda_u_energy, phase)
This routine computes the values of the occupation matrices.
Definition: lda_u.F90:3781
subroutine, public dcompute_dftu_energy(this, energy, st)
This routine computes the value of the double counting term in the DFT+U energy.
Definition: lda_u.F90:2210
subroutine, public dlda_u_update_potential(this, st)
This routine computes the potential that, once multiplied by the projector Pmm' and summed over m and...
Definition: lda_u.F90:2284
subroutine, public zlda_u_set_occupations(this, occ)
Definition: lda_u.F90:5440
subroutine, public dloewdin_orthogonalize(basis, kpt, namespace)
Definition: loewdin.F90:214
subroutine, public zloewdin_info(basis, kpt, namespace)
Definition: loewdin.F90:747
subroutine, public zloewdin_orthogonalize(basis, kpt, namespace)
Definition: loewdin.F90:539
subroutine, public dloewdin_info(basis, kpt, namespace)
Definition: loewdin.F90:415
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
This module defines various routines, operating on mesh functions.
subroutine, public zmf_fix_phase(mesh, ff)
Fix the phase of complex function.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:898
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1063
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:594
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
subroutine, public orbitalbasis_end(this)
subroutine, public zorbitalbasis_build(this, namespace, ions, mesh, kpt, ndim, skip_s_orb, use_all_orb, verbose)
This routine is an interface for constructing the orbital basis.
subroutine, public dorbitalbasis_build(this, namespace, ions, mesh, kpt, ndim, skip_s_orb, use_all_orb, verbose)
This routine is an interface for constructing the orbital basis.
subroutine, public dorbitalbasis_build_empty(this, mesh, ndim, norbsets, map_os, verbose)
This routine constructd an empty orbital basis.
subroutine, public zorbitalbasis_build_empty(this, mesh, ndim, norbsets, map_os, verbose)
This routine constructd an empty orbital basis.
subroutine, public orbitalbasis_init(this, namespace, periodic_dim)
subroutine, public orbitalset_update_phase_shift(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
Build the phase shift for the intersite interaction.
Definition: orbitalset.F90:380
subroutine, public orbitalset_update_phase(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
Build the phase correction to the global phase in case the orbital crosses the border of the simulato...
Definition: orbitalset.F90:284
integer, parameter, public sm_poisson_psolver
integer, parameter, public sm_poisson_isf
subroutine, public zorbitalset_get_center_of_mass(os, space, mesh, latt)
subroutine, public orbitalset_init_intersite(this, namespace, space, ind, ions, der, psolver, os, nos, maxnorbs, rcut, kpt, has_phase, sm_poisson, basis_from_states, combine_j_orbitals)
integer, parameter, public sm_poisson_direct
subroutine, public dorbitalset_get_center_of_mass(os, space, mesh, latt)
integer, parameter, public sm_poisson_fft
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:615
integer, parameter, public restart_gs
Definition: restart.F90:156
integer, parameter, public restart_proj
Definition: restart.F90:156
integer, parameter, public restart_type_load
Definition: restart.F90:183
pure logical function, public states_are_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 zsubmesh_copy_from_mesh(this, phi, sphi, conjugate)
Definition: submesh.F90:1826
subroutine, public dsubmesh_copy_from_mesh(this, phi, sphi, conjugate)
Definition: submesh.F90:1289
subroutine, public submesh_init(this, space, mesh, latt, center, rc)
Definition: submesh.F90:282
type(type_t), public type_float
Definition: types.F90:135
type(type_t), public type_cmplx
Definition: types.F90:136
type(type_t), public type_integer
Definition: types.F90:137
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
This class contains information about the boundary conditions.
Definition: boundaries.F90:159
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
Class to describe DFT+U parameters.
Definition: lda_u.F90:216
Describes mesh distribution to nodes.
Definition: mesh.F90:187
Stores all communicators and groups.
Definition: multicomm.F90:208
A container for the phase.
Definition: phase.F90:179
class for organizing spins and k-points
The states_elec_t class contains all electronic wave functions.
A submesh is a type of mesh, used for the projectors in the pseudopotentials It contains points on a ...
Definition: submesh.F90:177
int true(void)