64 type(projector_matrix_t),
allocatable,
public :: projector_matrices(:)
65 integer,
public :: nprojector_matrices
66 logical,
public :: apply_projector_matrices
67 logical,
public :: has_non_local_potential
68 integer :: full_projection_size
69 integer,
public :: max_npoints
70 integer,
public :: total_points
72 logical :: projector_mix
73 complex(real64),
allocatable,
public :: projector_phases(:, :, :, :)
74 integer,
allocatable,
public :: projector_to_atom(:)
76 integer,
allocatable :: regions(:)
77 integer,
public :: nphase
82 type(accel_mem_t) :: buff_offsets
83 type(accel_mem_t) :: buff_matrices
84 type(accel_mem_t) :: buff_maps
85 type(accel_mem_t) :: buff_scals
86 type(accel_mem_t) :: buff_position
87 type(accel_mem_t) :: buff_pos
88 type(accel_mem_t) :: buff_invmap
89 type(accel_mem_t),
public :: buff_projector_phases
90 type(accel_mem_t) :: buff_mix
91 logical :: projector_self_overlap
92 real(real64),
pointer,
public :: spin(:,:,:) => null()
129 real(real64),
allocatable :: dprojection(:, :)
130 complex(real64),
allocatable :: zprojection(:, :)
131 type(accel_mem_t) :: buff_projection
132 type(accel_mem_t) :: buff_spin_to_phase
133 type(accel_mem_t),
allocatable :: buff_phasepsi(:)
134 type(accel_mem_t),
allocatable :: buff_projection_temp(:)
143 class(nonlocal_pseudopotential_t),
intent(inout) :: this
147 this%apply_projector_matrices = .false.
148 this%has_non_local_potential = .false.
149 this%nprojector_matrices = 0
151 this%projector_self_overlap = .false.
166 if (
allocated(this%projector_matrices))
then
180 do iproj = 1, this%nprojector_matrices
183 safe_deallocate_a(this%regions)
184 safe_deallocate_a(this%projector_matrices)
185 safe_deallocate_a(this%projector_phases)
186 safe_deallocate_a(this%projector_to_atom)
201 class(
space_t),
intent(in) :: space
203 type(
epot_t),
target,
intent(in) :: epot
205 integer :: iatom, iproj, ll, lmax, lloc, mm, ic, jc
206 integer :: nmat, imat, ip, iorder
207 integer :: nregion, jatom, katom, iregion
208 integer,
allocatable :: order(:),
head(:), region_count(:)
209 logical,
allocatable :: atom_counted(:)
223 safe_allocate(order(1:epot%natoms))
224 safe_allocate(
head(1:epot%natoms + 1))
225 safe_allocate(region_count(1:epot%natoms))
226 safe_allocate(atom_counted(1:epot%natoms))
228 this%projector_self_overlap = .false.
229 atom_counted = .false.
235 nregion = nregion + 1
236 assert(nregion <= epot%natoms)
238 region_count(nregion) = 0
240 do iatom = 1, epot%natoms
241 if (atom_counted(iatom)) cycle
246 assert(
associated(epot%proj(iatom)%sphere%mesh))
247 do jatom = 1, region_count(nregion)
248 katom = order(
head(nregion) + jatom - 1)
250 overlap =
submesh_overlap(epot%proj(iatom)%sphere, epot%proj(katom)%sphere, space)
255 if (.not. overlap)
then
258 region_count(nregion) = region_count(nregion) + 1
259 order(
head(nregion) - 1 + region_count(nregion)) = iatom
260 atom_counted(iatom) = .
true.
265 head(nregion + 1) =
head(nregion) + region_count(nregion)
267 if (all(atom_counted))
exit
270 safe_deallocate_a(atom_counted)
271 safe_deallocate_a(region_count)
278 do iregion = 1, nregion
279 do iatom =
head(iregion),
head(iregion + 1) - 1
281 do jatom =
head(iregion), iatom - 1
283 assert(.not.
submesh_overlap(epot%proj(order(iatom))%sphere, epot%proj(order(jatom))%sphere, space))
294 this%nprojector_matrices = 0
295 this%apply_projector_matrices = .false.
296 this%has_non_local_potential = .false.
297 this%nregions = nregion
300 do iorder = 1, epot%natoms
301 iatom = order(iorder)
304 this%has_non_local_potential = .
true.
309 do iorder = 1, epot%natoms
310 iatom = order(iorder)
313 this%nprojector_matrices = this%nprojector_matrices + 1
314 this%apply_projector_matrices = .
true.
319 if (mesh%use_curvilinear) this%apply_projector_matrices = .false.
321 if (.not. this%apply_projector_matrices)
then
322 safe_deallocate_a(order)
323 safe_deallocate_a(
head)
330 safe_allocate(this%projector_matrices(1:this%nprojector_matrices))
331 safe_allocate(this%regions(1:this%nprojector_matrices + 1))
332 safe_allocate(this%projector_to_atom(1:epot%natoms))
334 this%full_projection_size = 0
335 this%regions(this%nregions + 1) = this%nprojector_matrices + 1
337 this%projector_mix = .false.
340 do iregion = 1, this%nregions
341 this%regions(iregion) = iproj + 1
342 do iorder =
head(iregion),
head(iregion + 1) - 1
344 iatom = order(iorder)
350 pmat => this%projector_matrices(iproj)
352 this%projector_to_atom(iproj) = iatom
354 lmax = epot%proj(iatom)%lmax
355 lloc = epot%proj(iatom)%lloc
362 if (ll == lloc) cycle
364 nmat = nmat + epot%proj(iatom)%kb_p(ll, mm)%n_c
374 if (ll == lloc) cycle
376 kb_p => epot%proj(iatom)%kb_p(ll, mm)
378 call lalg_copy(pmat%npoints, kb_p%p(:, ic), pmat%dprojectors(:, imat))
379 pmat%scal(imat) = kb_p%e(ic)*mesh%vol_pp(1)
385 this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
389 this%projector_mix = .
true.
394 if (ll == lloc) cycle
414 if (ll == lloc) cycle
416 hgh_p => epot%proj(iatom)%hgh_p(ll, mm)
423 pmat%zmix(imat - 1 + ic, imat - 1 + jc, 1) = hgh_p%h(ic, jc) +
m_half*mm*hgh_p%k(ic, jc)
424 pmat%zmix(imat - 1 + ic, imat - 1 + jc, 2) = hgh_p%h(ic, jc) -
m_half*mm*hgh_p%k(ic, jc)
427 pmat%zmix(imat - 1 + ic, imat + 3 - 1 + jc, 3) =
m_half*hgh_p%k(ic, jc) * &
428 sqrt(real(ll*(ll+1)-mm*(mm+1), real64))
432 pmat%zmix(imat - 1 + ic, imat - 3 - 1 + jc, 4) =
m_half*hgh_p%k(ic, jc) * &
433 sqrt(real(ll*(ll+1)-mm*(mm-1), real64))
440 pmat%dmix(imat - 1 + ic, imat - 1 + jc) = hgh_p%h(ic, jc)
447 call lalg_copy(pmat%npoints, hgh_p%zp(:, ic), pmat%zprojectors(:, imat))
449 call lalg_copy(pmat%npoints, hgh_p%dp(:, ic), pmat%dprojectors(:, imat))
451 pmat%scal(imat) = mesh%volume_element
458 this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
463 this%projector_mix = .
true.
467 if (lloc /= 0) nmat = nmat + epot%proj(iatom)%kb_p(1, 1)%n_c
470 if (ll == lloc) cycle
472 nmat = nmat + epot%proj(iatom)%rkb_p(ll, mm)%n_c
477 has_mix_matrix = .
true., is_cmplx = .
true.)
484 kb_p => epot%proj(iatom)%kb_p(1, 1)
487 pmat%zmix(ic, ic, 1:2) = kb_p%e(ic)
488 do ip = 1, pmat%npoints
489 pmat%zprojectors(ip, ic) = kb_p%p(ip, ic)
491 pmat%scal(ic) = mesh%volume_element
498 if (ll == lloc) cycle
500 rkb_p => epot%proj(iatom)%rkb_p(ll, mm)
503 do ic = 0, rkb_p%n_c/2-1
504 pmat%zmix(imat + ic*2, imat + ic*2, 1) = rkb_p%f(ic*2+1, 1, 1)
505 pmat%zmix(imat + ic*2, imat + ic*2, 2) = rkb_p%f(ic*2+1, 2, 2)
507 pmat%zmix(imat + ic*2+1, imat + ic*2+1, 1) = rkb_p%f(ic*2+2, 1, 1)
508 pmat%zmix(imat + ic*2+1, imat + ic*2+1, 2) = rkb_p%f(ic*2+2, 2, 2)
511 pmat%zmix(imat + ic*2+rkb_p%n_c, imat + ic*2, 4) = rkb_p%f(ic*2+1, 2, 1)
512 pmat%zmix(imat + ic*2+1+rkb_p%n_c, imat + ic*2+1, 4) = rkb_p%f(ic*2+2, 2, 1)
516 pmat%zmix(imat + ic*2-rkb_p%n_c, imat + ic*2, 3) = rkb_p%f(ic*2+1, 1, 2)
517 pmat%zmix(imat + ic*2+1-rkb_p%n_c, imat + ic*2+1, 3) = rkb_p%f(ic*2+2, 1, 2)
522 call lalg_copy(pmat%npoints, rkb_p%ket(:, ic, 1, 1), pmat%zprojectors(:, imat))
523 pmat%scal(imat) = mesh%volume_element
531 this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
537 pmat%map => epot%proj(iatom)%sphere%map
538 pmat%position => epot%proj(iatom)%sphere%rel_x
540 pmat%regions = epot%proj(iatom)%sphere%regions
542 this%full_projection_size = this%full_projection_size + pmat%nprojs
547 if (mesh%parallel_in_domains)
then
548 call mesh%mpi_grp%allreduce_inplace(this%projector_self_overlap, 1, mpi_logical, mpi_lor)
551 safe_deallocate_a(order)
552 safe_deallocate_a(
head)
554 this%total_points = 0
557 do imat = 1, this%nprojector_matrices
558 pmat => this%projector_matrices(imat)
560 this%max_npoints = max(this%max_npoints, pmat%npoints)
561 this%max_nprojs = max(this%max_nprojs, pmat%nprojs)
562 this%total_points = this%total_points + pmat%npoints
573 integer :: matrix_size, scal_size
574 integer,
allocatable :: cnt(:), invmap(:, :), invmap2(:), pos(:)
575 integer,
allocatable :: offsets(:, :)
576 integer,
parameter :: OFFSET_SIZE = 6
577 integer,
parameter :: POINTS = 1, projs = 2, matrix = 3, map = 4, scal = 5, mix = 6
578 integer :: ip, is, ii, ipos, mix_offset
582 safe_allocate(offsets(1:offset_size, 1:this%nprojector_matrices))
583 safe_allocate(cnt(1:mesh%np))
602 this%total_points = 0
607 do imat = 1, this%nprojector_matrices
608 pmat => this%projector_matrices(imat)
610 this%max_npoints = max(this%max_npoints, pmat%npoints)
611 this%max_nprojs = max(this%max_nprojs, pmat%nprojs)
613 offsets(points, imat) = pmat%npoints
614 offsets(projs, imat) = pmat%nprojs
616 offsets(matrix, imat) = matrix_size
617 matrix_size = matrix_size + pmat%npoints*pmat%nprojs
619 offsets(map, imat) = this%total_points
620 this%total_points = this%total_points + pmat%npoints
622 offsets(scal, imat) = scal_size
623 scal_size = scal_size + pmat%nprojs
625 offsets(mix, imat) = mix_offset
626 if (
allocated(pmat%dmix))
then
627 mix_offset = mix_offset + pmat%nprojs**2
628 else if (
allocated(pmat%zmix))
then
629 mix_offset = mix_offset + 4 * pmat%nprojs**2
631 offsets(mix, imat) = -1
634 do is = 1, pmat%npoints
636 cnt(ip) = cnt(ip) + 1
640 safe_allocate(invmap(1:max(maxval(cnt), 1), 1:mesh%np))
641 safe_allocate(invmap2(1:max(maxval(cnt)*mesh%np, 1)))
642 safe_allocate(pos(1:mesh%np + 1))
646 do imat = 1, this%nprojector_matrices
647 pmat => this%projector_matrices(imat)
648 do is = 1, pmat%npoints
650 cnt(ip) = cnt(ip) + 1
651 invmap(cnt(ip), ip) = ii
661 invmap2(ipos) = invmap(ii, ip)
667 if (this%projector_matrices(1)%is_cmplx)
then
676 if (mix_offset > 0)
then
677 if (
allocated(this%projector_matrices(1)%zmix))
then
685 do imat = 1, this%nprojector_matrices
686 pmat => this%projector_matrices(imat)
688 if (pmat%npoints > 0)
then
689 if (pmat%is_cmplx)
then
690 call accel_write_buffer(this%buff_matrices, pmat%nprojs*pmat%npoints, pmat%zprojectors, offset = offsets(matrix, imat))
692 call accel_write_buffer(this%buff_matrices, pmat%nprojs*pmat%npoints, pmat%dprojectors, offset = offsets(matrix, imat))
694 call accel_write_buffer(this%buff_maps, pmat%npoints, pmat%map, offset = offsets(map, imat))
695 call accel_write_buffer(this%buff_position, 3*pmat%npoints, pmat%position, offset = 3*offsets(map, imat))
697 call accel_write_buffer(this%buff_scals, pmat%nprojs, pmat%scal, offset = offsets(scal, imat))
698 if (offsets(mix, imat) /= -1)
then
699 if (
allocated(pmat%zmix))
then
700 call accel_write_buffer(this%buff_mix, 4*pmat%nprojs**2, pmat%zmix, offset = offsets(mix, imat))
702 call accel_write_buffer(this%buff_mix, pmat%nprojs**2, pmat%dmix, offset = offsets(mix, imat))
709 call accel_write_buffer(this%buff_offsets, offset_size*this%nprojector_matrices, offsets)
718 safe_deallocate_a(offsets)
719 safe_deallocate_a(cnt)
720 safe_deallocate_a(invmap)
721 safe_deallocate_a(invmap2)
722 safe_deallocate_a(pos)
732 logical pure function nonlocal_pseudopotential_self_overlap(this) result(projector_self_overlap)
735 projector_self_overlap = this%projector_self_overlap
740#include "nonlocal_pseudopotential_inc.F90"
743#include "complex.F90"
744#include "nonlocal_pseudopotential_inc.F90"
Copies a vector x, to a vector y.
double sqrt(double __x) __attribute__((__nothrow__
subroutine, public accel_release_buffer(this, async)
pure logical function, public accel_is_enabled()
integer, parameter, public accel_mem_read_only
type(accel_kernel_t), pointer head
This module implements batches of mesh functions.
This module implements common operations on batches of mesh functions.
This module contains interfaces for BLAS routines You should not use these routines directly....
integer, parameter, public spin_orbit
integer, parameter, public fully_relativistic_zora
real(real64), parameter, public m_zero
real(real64), parameter, public m_half
This module is intended to contain "only mathematical" functions and procedures.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
subroutine nonlocal_pseudopotential_destroy_proj(this)
Destroy the data of nonlocal_pseudopotential_t.
subroutine dnonlocal_pseudopotential_force(this, mesh, st, spiral_bnd, iqn, ndim, psi1b, psi2b, force)
calculate contribution to forces, from non-local potentials
subroutine znonlocal_pseudopotential_position_commutator(this, mesh, std, spiral_bnd, psib, commpsib, async)
apply the commutator between the non-local potential and the position to the wave functions.
subroutine znonlocal_pseudopotential_force(this, mesh, st, spiral_bnd, iqn, ndim, psi1b, psi2b, force)
calculate contribution to forces, from non-local potentials
subroutine dnonlocal_pseudopotential_start(this, mesh, std, spiral_bnd, psib, projection, async)
Start application of non-local potentials (stored in the Hamiltonian) to the wave functions.
subroutine dnonlocal_pseudopotential_finish(this, mesh, spiral_bnd, std, projection, vpsib)
finish the application of non-local potentials.
subroutine dnonlocal_pseudopotential_position_commutator(this, mesh, std, spiral_bnd, psib, commpsib, async)
apply the commutator between the non-local potential and the position to the wave functions.
subroutine znonlocal_pseudopotential_r_vnlocal(this, mesh, std, spiral_bnd, psib, commpsib)
Accumulates to commpsib the result of x V_{nl} | psib >
logical pure function nonlocal_pseudopotential_self_overlap(this)
Returns .true. if the Hamiltonian contains projectors, which overlap with themself.
subroutine nonlocal_pseudopotential_init(this)
initialize the nonlocal_pseudopotential_t object
subroutine dnonlocal_pseudopotential_r_vnlocal(this, mesh, std, spiral_bnd, psib, commpsib)
Accumulates to commpsib the result of x V_{nl} | psib >
subroutine znonlocal_pseudopotential_start(this, mesh, std, spiral_bnd, psib, projection, async)
Start application of non-local potentials (stored in the Hamiltonian) to the wave functions.
subroutine nonlocal_pseudopotential_build_proj(this, space, mesh, epot)
build the projectors for the application of pseudo-potentials
subroutine znonlocal_pseudopotential_finish(this, mesh, spiral_bnd, std, projection, vpsib)
finish the application of non-local potentials.
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_matrix_deallocate(this)
subroutine, public projector_matrix_allocate(this, nprojs, sphere, has_mix_matrix, is_cmplx)
logical elemental function, public projector_is(p, type)
logical elemental function, public projector_is_null(p)
integer, parameter, public proj_hgh
integer, parameter, public proj_rkb
integer, parameter, public proj_none
integer, parameter, public proj_kb
This module handles spin dimensions of the states and the k-point distribution.
logical function, public submesh_overlap(sm1, sm2, space)
type(type_t), public type_float
type(type_t), public type_cmplx
type(type_t), public type_integer
Describes mesh distribution to nodes.
nonlocal part of the pseudopotential
Class for projections of wave functions.
A set of projectors defined on a submesh.
The rkb_projector data type holds the KB projectors build with total angular momentum eigenfunctions....