63 type(projector_matrix_t),
allocatable,
public :: projector_matrices(:)
64 integer,
public :: nprojector_matrices
65 logical,
public :: apply_projector_matrices
66 logical,
public :: has_non_local_potential
67 integer :: full_projection_size
68 integer,
public :: max_npoints
69 integer,
public :: total_points
71 logical :: projector_mix
72 complex(real64),
allocatable,
public :: projector_phases(:, :, :, :)
73 integer,
allocatable,
public :: projector_to_atom(:)
75 integer,
allocatable :: regions(:)
76 integer,
public :: nphase
81 type(accel_mem_t) :: buff_offsets
82 type(accel_mem_t) :: buff_matrices
83 type(accel_mem_t) :: buff_maps
84 type(accel_mem_t) :: buff_scals
85 type(accel_mem_t) :: buff_position
86 type(accel_mem_t) :: buff_pos
87 type(accel_mem_t) :: buff_invmap
88 type(accel_mem_t),
public :: buff_projector_phases
89 type(accel_mem_t) :: buff_mix
90 logical :: projector_self_overlap
91 real(real64),
pointer,
public :: spin(:,:,:) => null()
128 real(real64),
allocatable :: dprojection(:, :)
129 complex(real64),
allocatable :: zprojection(:, :)
130 type(accel_mem_t) :: buff_projection
131 type(accel_mem_t) :: buff_spin_to_phase
140 class(nonlocal_pseudopotential_t),
intent(inout) :: this
144 this%apply_projector_matrices = .false.
145 this%has_non_local_potential = .false.
146 this%nprojector_matrices = 0
148 this%projector_self_overlap = .false.
163 if (
allocated(this%projector_matrices))
then
177 do iproj = 1, this%nprojector_matrices
180 safe_deallocate_a(this%regions)
181 safe_deallocate_a(this%projector_matrices)
182 safe_deallocate_a(this%projector_phases)
183 safe_deallocate_a(this%projector_to_atom)
198 class(
space_t),
intent(in) :: space
200 type(
epot_t),
target,
intent(in) :: epot
202 integer :: iatom, iproj, ll, lmax, lloc, mm, ic, jc
203 integer :: nmat, imat, ip, iorder
204 integer :: nregion, jatom, katom, iregion
205 integer,
allocatable :: order(:),
head(:), region_count(:)
206 logical,
allocatable :: atom_counted(:)
220 safe_allocate(order(1:epot%natoms))
221 safe_allocate(
head(1:epot%natoms + 1))
222 safe_allocate(region_count(1:epot%natoms))
223 safe_allocate(atom_counted(1:epot%natoms))
225 this%projector_self_overlap = .false.
226 atom_counted = .false.
232 nregion = nregion + 1
233 assert(nregion <= epot%natoms)
235 region_count(nregion) = 0
237 do iatom = 1, epot%natoms
238 if (atom_counted(iatom)) cycle
243 assert(
associated(epot%proj(iatom)%sphere%mesh))
244 do jatom = 1, region_count(nregion)
245 katom = order(
head(nregion) + jatom - 1)
247 overlap =
submesh_overlap(epot%proj(iatom)%sphere, epot%proj(katom)%sphere, space)
252 if (.not. overlap)
then
255 region_count(nregion) = region_count(nregion) + 1
256 order(
head(nregion) - 1 + region_count(nregion)) = iatom
257 atom_counted(iatom) = .
true.
262 head(nregion + 1) =
head(nregion) + region_count(nregion)
264 if (all(atom_counted))
exit
267 safe_deallocate_a(atom_counted)
268 safe_deallocate_a(region_count)
275 do iregion = 1, nregion
276 do iatom =
head(iregion),
head(iregion + 1) - 1
278 do jatom =
head(iregion), iatom - 1
280 assert(.not.
submesh_overlap(epot%proj(order(iatom))%sphere, epot%proj(order(jatom))%sphere, space))
291 this%nprojector_matrices = 0
292 this%apply_projector_matrices = .false.
293 this%has_non_local_potential = .false.
294 this%nregions = nregion
297 do iorder = 1, epot%natoms
298 iatom = order(iorder)
301 this%has_non_local_potential = .
true.
306 do iorder = 1, epot%natoms
307 iatom = order(iorder)
310 this%nprojector_matrices = this%nprojector_matrices + 1
311 this%apply_projector_matrices = .
true.
316 if (mesh%use_curvilinear) this%apply_projector_matrices = .false.
318 if (.not. this%apply_projector_matrices)
then
319 safe_deallocate_a(order)
320 safe_deallocate_a(
head)
327 safe_allocate(this%projector_matrices(1:this%nprojector_matrices))
328 safe_allocate(this%regions(1:this%nprojector_matrices + 1))
329 safe_allocate(this%projector_to_atom(1:epot%natoms))
331 this%full_projection_size = 0
332 this%regions(this%nregions + 1) = this%nprojector_matrices + 1
334 this%projector_mix = .false.
337 do iregion = 1, this%nregions
338 this%regions(iregion) = iproj + 1
339 do iorder =
head(iregion),
head(iregion + 1) - 1
341 iatom = order(iorder)
347 pmat => this%projector_matrices(iproj)
349 this%projector_to_atom(iproj) = iatom
351 lmax = epot%proj(iatom)%lmax
352 lloc = epot%proj(iatom)%lloc
359 if (ll == lloc) cycle
361 nmat = nmat + epot%proj(iatom)%kb_p(ll, mm)%n_c
371 if (ll == lloc) cycle
373 kb_p => epot%proj(iatom)%kb_p(ll, mm)
375 call lalg_copy(pmat%npoints, kb_p%p(:, ic), pmat%dprojectors(:, imat))
376 pmat%scal(imat) = kb_p%e(ic)*mesh%vol_pp(1)
382 this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
386 this%projector_mix = .
true.
391 if (ll == lloc) cycle
411 if (ll == lloc) cycle
413 hgh_p => epot%proj(iatom)%hgh_p(ll, mm)
420 pmat%zmix(imat - 1 + ic, imat - 1 + jc, 1) = hgh_p%h(ic, jc) +
m_half*mm*hgh_p%k(ic, jc)
421 pmat%zmix(imat - 1 + ic, imat - 1 + jc, 2) = hgh_p%h(ic, jc) -
m_half*mm*hgh_p%k(ic, jc)
424 pmat%zmix(imat - 1 + ic, imat + 3 - 1 + jc, 3) =
m_half*hgh_p%k(ic, jc) * &
425 sqrt(real(ll*(ll+1)-mm*(mm+1), real64))
429 pmat%zmix(imat - 1 + ic, imat - 3 - 1 + jc, 4) =
m_half*hgh_p%k(ic, jc) * &
430 sqrt(real(ll*(ll+1)-mm*(mm-1), real64))
437 pmat%dmix(imat - 1 + ic, imat - 1 + jc) = hgh_p%h(ic, jc)
444 call lalg_copy(pmat%npoints, hgh_p%zp(:, ic), pmat%zprojectors(:, imat))
446 call lalg_copy(pmat%npoints, hgh_p%dp(:, ic), pmat%dprojectors(:, imat))
448 pmat%scal(imat) = mesh%volume_element
455 this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
460 this%projector_mix = .
true.
464 if (lloc /= 0) nmat = nmat + epot%proj(iatom)%kb_p(1, 1)%n_c
467 if (ll == lloc) cycle
469 nmat = nmat + epot%proj(iatom)%rkb_p(ll, mm)%n_c
474 has_mix_matrix = .
true., is_cmplx = .
true.)
481 kb_p => epot%proj(iatom)%kb_p(1, 1)
484 pmat%zmix(ic, ic, 1:2) = kb_p%e(ic)
485 do ip = 1, pmat%npoints
486 pmat%zprojectors(ip, ic) = kb_p%p(ip, ic)
488 pmat%scal(ic) = mesh%volume_element
495 if (ll == lloc) cycle
497 rkb_p => epot%proj(iatom)%rkb_p(ll, mm)
500 do ic = 0, rkb_p%n_c/2-1
501 pmat%zmix(imat + ic*2, imat + ic*2, 1) = rkb_p%f(ic*2+1, 1, 1)
502 pmat%zmix(imat + ic*2, imat + ic*2, 2) = rkb_p%f(ic*2+1, 2, 2)
504 pmat%zmix(imat + ic*2+1, imat + ic*2+1, 1) = rkb_p%f(ic*2+2, 1, 1)
505 pmat%zmix(imat + ic*2+1, imat + ic*2+1, 2) = rkb_p%f(ic*2+2, 2, 2)
508 pmat%zmix(imat + ic*2+rkb_p%n_c, imat + ic*2, 4) = rkb_p%f(ic*2+1, 2, 1)
509 pmat%zmix(imat + ic*2+1+rkb_p%n_c, imat + ic*2+1, 4) = rkb_p%f(ic*2+2, 2, 1)
513 pmat%zmix(imat + ic*2-rkb_p%n_c, imat + ic*2, 3) = rkb_p%f(ic*2+1, 1, 2)
514 pmat%zmix(imat + ic*2+1-rkb_p%n_c, imat + ic*2+1, 3) = rkb_p%f(ic*2+2, 1, 2)
519 call lalg_copy(pmat%npoints, rkb_p%ket(:, ic, 1, 1), pmat%zprojectors(:, imat))
520 pmat%scal(imat) = mesh%volume_element
528 this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
534 pmat%map => epot%proj(iatom)%sphere%map
535 pmat%position => epot%proj(iatom)%sphere%rel_x
537 pmat%regions = epot%proj(iatom)%sphere%regions
539 this%full_projection_size = this%full_projection_size + pmat%nprojs
544 if (mesh%parallel_in_domains)
then
545 call mesh%mpi_grp%allreduce_inplace(this%projector_self_overlap, 1, mpi_logical, mpi_lor)
548 safe_deallocate_a(order)
549 safe_deallocate_a(
head)
551 this%total_points = 0
554 do imat = 1, this%nprojector_matrices
555 pmat => this%projector_matrices(imat)
557 this%max_npoints = max(this%max_npoints, pmat%npoints)
558 this%max_nprojs = max(this%max_nprojs, pmat%nprojs)
559 this%total_points = this%total_points + pmat%npoints
570 integer :: matrix_size, scal_size
571 integer,
allocatable :: cnt(:), invmap(:, :), invmap2(:), pos(:)
572 integer,
allocatable :: offsets(:, :)
573 integer,
parameter :: OFFSET_SIZE = 6
574 integer,
parameter :: POINTS = 1, projs = 2, matrix = 3, map = 4, scal = 5, mix = 6
575 integer :: ip, is, ii, ipos, mix_offset
579 safe_allocate(offsets(1:offset_size, 1:this%nprojector_matrices))
580 safe_allocate(cnt(1:mesh%np))
599 this%total_points = 0
604 do imat = 1, this%nprojector_matrices
605 pmat => this%projector_matrices(imat)
607 this%max_npoints = max(this%max_npoints, pmat%npoints)
608 this%max_nprojs = max(this%max_nprojs, pmat%nprojs)
610 offsets(points, imat) = pmat%npoints
611 offsets(projs, imat) = pmat%nprojs
613 offsets(matrix, imat) = matrix_size
614 matrix_size = matrix_size + pmat%npoints*pmat%nprojs
616 offsets(map, imat) = this%total_points
617 this%total_points = this%total_points + pmat%npoints
619 offsets(scal, imat) = scal_size
620 scal_size = scal_size + pmat%nprojs
622 offsets(mix, imat) = mix_offset
623 if (
allocated(pmat%dmix))
then
624 mix_offset = mix_offset + pmat%nprojs**2
625 else if (
allocated(pmat%zmix))
then
626 mix_offset = mix_offset + 4 * pmat%nprojs**2
628 offsets(mix, imat) = -1
631 do is = 1, pmat%npoints
633 cnt(ip) = cnt(ip) + 1
637 safe_allocate(invmap(1:max(maxval(cnt), 1), 1:mesh%np))
638 safe_allocate(invmap2(1:max(maxval(cnt)*mesh%np, 1)))
639 safe_allocate(pos(1:mesh%np + 1))
643 do imat = 1, this%nprojector_matrices
644 pmat => this%projector_matrices(imat)
645 do is = 1, pmat%npoints
647 cnt(ip) = cnt(ip) + 1
648 invmap(cnt(ip), ip) = ii
658 invmap2(ipos) = invmap(ii, ip)
664 if (this%projector_matrices(1)%is_cmplx)
then
673 if (mix_offset > 0)
then
674 if (
allocated(this%projector_matrices(1)%zmix))
then
682 do imat = 1, this%nprojector_matrices
683 pmat => this%projector_matrices(imat)
685 if (pmat%npoints > 0)
then
686 if (pmat%is_cmplx)
then
687 call accel_write_buffer(this%buff_matrices, pmat%nprojs*pmat%npoints, pmat%zprojectors, offset = offsets(matrix, imat))
689 call accel_write_buffer(this%buff_matrices, pmat%nprojs*pmat%npoints, pmat%dprojectors, offset = offsets(matrix, imat))
691 call accel_write_buffer(this%buff_maps, pmat%npoints, pmat%map, offset = offsets(map, imat))
692 call accel_write_buffer(this%buff_position, 3*pmat%npoints, pmat%position, offset = 3*offsets(map, imat))
694 call accel_write_buffer(this%buff_scals, pmat%nprojs, pmat%scal, offset = offsets(scal, imat))
695 if (offsets(mix, imat) /= -1)
then
696 if (
allocated(pmat%zmix))
then
697 call accel_write_buffer(this%buff_mix, 4*pmat%nprojs**2, pmat%zmix, offset = offsets(mix, imat))
699 call accel_write_buffer(this%buff_mix, pmat%nprojs**2, pmat%dmix, offset = offsets(mix, imat))
706 call accel_write_buffer(this%buff_offsets, offset_size*this%nprojector_matrices, offsets)
715 safe_deallocate_a(offsets)
716 safe_deallocate_a(cnt)
717 safe_deallocate_a(invmap)
718 safe_deallocate_a(invmap2)
719 safe_deallocate_a(pos)
729 logical pure function nonlocal_pseudopotential_self_overlap(this) result(projector_self_overlap)
732 projector_self_overlap = this%projector_self_overlap
737#include "nonlocal_pseudopotential_inc.F90"
740#include "complex.F90"
741#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)
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....