Commit d6ea3890 authored by Matthias Redies's avatar Matthias Redies
Browse files

introduce mpbasis type

parent 2c951cf5
This diff is collapsed.
......@@ -56,7 +56,7 @@ MODULE m_exchange_valence_hf
CONTAINS
SUBROUTINE exchange_valence_hf(nk, kpts, nkpt_EIBZ, sym, atoms, hybrid, cell, dimension, input, jsp, hybdat, mnobd, lapw, &
SUBROUTINE exchange_valence_hf(nk, kpts, nkpt_EIBZ, sym, atoms, mpbasis, hybrid, cell, dimension, input, jsp, hybdat, mnobd, lapw, &
eig_irr, results, pointer_EIBZ, n_q, wl_iks, xcpot, noco, nsest, indx_sest, &
mpi, mat_ex)
......@@ -77,6 +77,7 @@ CONTAINS
TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_mpbasis), intent(in) :: mpbasis
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_input), INTENT(IN) :: input
TYPE(t_noco), INTENT(IN) :: noco
......@@ -139,10 +140,10 @@ CONTAINS
COMPLEX :: coulomb_mt2_c(maxval(hybrid%nindxm1) - 1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1 + 1, atoms%nat)
COMPLEX :: coulomb_mt3_c(maxval(hybrid%nindxm1) - 1, atoms%nat, atoms%nat)
REAL :: coulomb_mtir_r(((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm))* &
((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm) + 1)/2)
COMPLEX :: coulomb_mtir_c(((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm))* &
((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(hybrid%ngptm) + 1)/2)
REAL :: coulomb_mtir_r(((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(mpbasis%ngptm))* &
((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(mpbasis%ngptm) + 1)/2)
COMPLEX :: coulomb_mtir_c(((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(mpbasis%ngptm))* &
((hybrid%maxlcutm1 + 1)**2*atoms%nat + maxval(mpbasis%ngptm) + 1)/2)
LOGICAL :: occup(dimension%neigd)
CALL timestart("valence exchange calculation")
......@@ -203,7 +204,7 @@ CONTAINS
ikpt0 = pointer_EIBZ(ikpt)
n = hybrid%nbasp + hybrid%ngptm(ikpt0)
n = hybrid%nbasp + mpbasis%ngptm(ikpt0)
IF (hybrid%nbasm(ikpt0) /= n) call judft_error('error hybrid%nbasm')
nn = n*(n + 1)/2
......@@ -225,11 +226,11 @@ CONTAINS
IF (mat_ex%l_real) THEN
CALL wavefproducts_inv5(1, hybrid%nbands(nk), ibando, ibando + psize - 1, dimension, input, jsp, atoms, &
lapw, kpts, nk, ikpt0, hybdat, hybrid, cell, hybrid%nbasp, sym, &
lapw, kpts, nk, ikpt0, hybdat, mpbasis, hybrid, cell, hybrid%nbasp, sym, &
noco, nkqpt, cprod_vv_r)
ELSE
CALL wavefproducts_noinv5(1, hybrid%nbands(nk), ibando, ibando + psize - 1, nk, ikpt0, dimension, input, jsp, &!jsp,&
cell, atoms, hybrid, hybdat, kpts, lapw, sym, hybrid%nbasp, noco, nkqpt, cprod_vv_c)
cell, atoms, mpbasis, hybrid, hybdat, kpts, lapw, sym, hybrid%nbasp, noco, nkqpt, cprod_vv_c)
END IF
! The sparse matrix technique is not feasible for the HSE
......@@ -243,8 +244,8 @@ CONTAINS
exch_vv = exch_vv + &
dynamic_hse_adjustment(atoms%rmsh, atoms%rmt, atoms%dx, atoms%jri, atoms%jmtd, kpts%bkf(:, ikpt0), ikpt0, &
kpts%nkptf, cell%bmat, cell%omtil, atoms%ntype, atoms%neq, atoms%nat, atoms%taual, &
hybrid%lcutm1, hybrid%maxlcutm1, hybrid%nindxm1, maxval(hybrid%nindxm1), hybrid%gptm, &
hybrid%ngptm(ikpt0), hybrid%pgptm(:, ikpt0), hybrid%gptmd, hybrid%basm1, &
hybrid%lcutm1, hybrid%maxlcutm1, hybrid%nindxm1, maxval(hybrid%nindxm1), mpbasis%gptm, &
mpbasis%ngptm(ikpt0), hybrid%pgptm(:, ikpt0), mpbasis%gptmd, hybrid%basm1, &
hybrid%nbasm(ikpt0), iband1, hybrid%nbands(nk), nsest, ibando, psize, indx_sest, &
atoms%invsat, sym%invsatnr, mpi%irank, cprod_vv_r(:hybrid%nbasm(ikpt0), :, :), &
cprod_vv_c(:hybrid%nbasm(ikpt0), :, :), mat_ex%l_real, wl_iks(:iband1, nkqpt), n_q(ikpt))
......@@ -257,7 +258,7 @@ CONTAINS
CALL bra_trafo2(mat_ex%l_real, carr3_vv_r(:hybrid%nbasm(ikpt0), :, :), cprod_vv_r(:hybrid%nbasm(ikpt0), :, :), &
carr3_vv_c(:hybrid%nbasm(ikpt0), :, :), cprod_vv_c(:hybrid%nbasm(ikpt0), :, :), &
hybrid%nbasm(ikpt0), psize, hybrid%nbands(nk), kpts%bkp(ikpt0), ikpt0, kpts%bksym(ikpt0), sym, &
hybrid, kpts, cell, atoms, phase_vv)
mpbasis, hybrid, kpts, cell, atoms, phase_vv)
IF (mat_ex%l_real) THEN
cprod_vv_r(:hybrid%nbasm(ikpt0), :, :) = carr3_vv_r(:hybrid%nbasm(ikpt0), :, :)
ELSE
......@@ -278,11 +279,11 @@ CONTAINS
IF (mat_ex%l_real) THEN
carr1_v_r(:n) = 0
CALL spmvec_invs(atoms, hybrid, hybdat, ikpt0, kpts, cell, coulomb_mt1, coulomb_mt2_r, coulomb_mt3_r, &
CALL spmvec_invs(atoms, mpbasis, hybrid, hybdat, ikpt0, kpts, cell, coulomb_mt1, coulomb_mt2_r, coulomb_mt3_r, &
coulomb_mtir_r, cprod_vv_r(:n, iband, n1), carr1_v_r(:n))
ELSE
carr1_v_c(:n) = 0
CALL spmvec_noinvs(atoms, hybrid, hybdat, ikpt0, kpts, cell, coulomb_mt1, coulomb_mt2_c, coulomb_mt3_c, &
CALL spmvec_noinvs(atoms, mpbasis, hybrid, hybdat, ikpt0, kpts, cell, coulomb_mt1, coulomb_mt2_c, coulomb_mt3_c, &
coulomb_mtir_c, cprod_vv_c(:n, iband, n1), carr1_v_c(:n))
END IF
......
......@@ -44,7 +44,7 @@ MODULE m_hsfock
CONTAINS
SUBROUTINE hsfock(nk, atoms, hybrid, lapw, dimension, kpts, jsp, input, hybdat, eig_irr, sym, cell, noco, &
SUBROUTINE hsfock(nk, atoms, mpbasis, hybrid, lapw, dimension, kpts, jsp, input, hybdat, eig_irr, sym, cell, noco, &
results, it, mnobd, xcpot, mpi)
USE m_types
......@@ -69,6 +69,7 @@ CONTAINS
TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_lapw), INTENT(IN) :: lapw
TYPE(t_mpbasis), intent(in) :: mpbasis
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_hybdat), INTENT(INOUT) :: hybdat
TYPE(t_results), INTENT(INOUT) :: results
......@@ -165,7 +166,7 @@ CONTAINS
! calculate contribution from valence electrons to the
! HF exchange
ex%l_real = sym%invs
CALL exchange_valence_hf(nk, kpts, nkpt_EIBZ, sym, atoms, hybrid, cell, dimension, input, jsp, hybdat, mnobd, lapw, &
CALL exchange_valence_hf(nk, kpts, nkpt_EIBZ, sym, atoms, mpbasis, hybrid, cell, dimension, input, jsp, hybdat, mnobd, lapw, &
eig_irr, results, pointer_EIBZ, n_q, wl_iks, xcpot, noco, nsest, indx_sest, &
mpi, ex)
......
......@@ -9,7 +9,7 @@ MODULE m_calc_hybrid
CONTAINS
SUBROUTINE calc_hybrid(eig_id, hybrid, kpts, atoms, input, DIMENSION, mpi, noco, cell, oneD, &
SUBROUTINE calc_hybrid(eig_id, mpbasis, hybrid, kpts, atoms, input, DIMENSION, mpi, noco, cell, oneD, &
enpara, results, sym, xcpot, v, iterHF)
USE m_types
......@@ -27,6 +27,7 @@ CONTAINS
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: DIMENSION
TYPE(t_oneD), INTENT(IN) :: oneD
type(t_mpbasis), intent(inout) :: mpbasis
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_input), INTENT(IN) :: input
TYPE(t_noco), INTENT(IN) :: noco
......@@ -105,12 +106,12 @@ CONTAINS
!construct the mixed-basis
CALL timestart("generation of mixed basis")
CALL mixedbasis(atoms, kpts, input, cell, xcpot, hybrid, enpara, mpi, v)
CALL mixedbasis(atoms, kpts, input, cell, xcpot, mpbasis, hybrid, enpara, mpi, v)
CALL timestop("generation of mixed basis")
CALL open_hybrid_io2(hybrid, DIMENSION, atoms, sym%invs)
CALL open_hybrid_io2(mpbasis, hybrid, DIMENSION, atoms, sym%invs)
CALL coulombmatrix(mpi, atoms, kpts, cell, sym, hybrid, xcpot)
CALL coulombmatrix(mpi, atoms, kpts, cell, sym, mpbasis, hybrid, xcpot)
CALL hf_init(hybrid, atoms, input, DIMENSION, hybdat)
CALL timestop("Preparation for Hybrid functionals")
......@@ -124,7 +125,7 @@ CONTAINS
DO nk = 1, kpts%nkpt
!DO nk = mpi%n_start,kpts%nkpt,mpi%n_stride
CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, l_zref)
CALL hsfock(nk, atoms, hybrid, lapw, DIMENSION, kpts, jsp, input, hybdat, eig_irr, sym, cell, &
CALL hsfock(nk, atoms, mpbasis, hybrid, lapw, DIMENSION, kpts, jsp, input, hybdat, eig_irr, sym, cell, &
noco, results, iterHF, MAXVAL(hybrid%nobd), xcpot, mpi)
END DO
END DO
......
......@@ -38,7 +38,7 @@ MODULE m_mixedbasis
CONTAINS
SUBROUTINE mixedbasis(atoms, kpts, input, cell, xcpot, hybrid, enpara, mpi, v)
SUBROUTINE mixedbasis(atoms, kpts, input, cell, xcpot, mpbasis, hybrid, enpara, mpi, v)
USE m_judft
USE m_loddop, ONLY: loddop
......@@ -52,6 +52,7 @@ CONTAINS
TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_mpbasis), intent(inout) :: mpbasis
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_enpara), INTENT(IN) :: enpara
TYPE(t_input), INTENT(IN) :: input
......@@ -99,12 +100,12 @@ CONTAINS
IF (xcpot%is_name("exx")) CALL judft_error("EXX is not implemented in this version", calledby='mixedbasis')
! Deallocate arrays which might have been allocated in a previous run of this subroutine
IF (ALLOCATED(hybrid%ngptm)) deallocate(hybrid%ngptm)
IF (ALLOCATED(mpbasis%ngptm)) deallocate(mpbasis%ngptm)
IF (ALLOCATED(hybrid%ngptm1)) deallocate(hybrid%ngptm1)
IF (ALLOCATED(hybrid%nindxm1)) deallocate(hybrid%nindxm1)
IF (ALLOCATED(hybrid%pgptm)) deallocate(hybrid%pgptm)
IF (ALLOCATED(hybrid%pgptm1)) deallocate(hybrid%pgptm1)
IF (ALLOCATED(hybrid%gptm)) deallocate(hybrid%gptm)
IF (ALLOCATED(mpbasis%gptm)) deallocate(mpbasis%gptm)
IF (ALLOCATED(hybrid%basm1)) deallocate(hybrid%basm1)
CALL usdus%init(atoms, input%jspins)
......@@ -126,15 +127,15 @@ CONTAINS
! - - - - - - SETUP OF THE MIXED BASIS IN THE IR - - - - - - -
! construct G-vectors with cutoff smaller than gcutm
call gen_gvec(cell, kpts, hybrid)
call gen_gvec(cell, kpts, mpbasis, hybrid)
! construct IR mixed basis set for the representation of the non local exchange elements with cutoff gcutm
! first run to determine dimension of pgptm1
allocate(hybrid%ngptm1(kpts%nkptf))
hybrid%ngptm1 = 0
DO igpt = 1, hybrid%gptmd
g = hybrid%gptm(:,igpt)
DO igpt = 1, mpbasis%gptmd
g = mpbasis%gptm(:,igpt)
DO ikpt = 1, kpts%nkptf
kvec = kpts%bkf(:,ikpt)
rdum = norm2(MATMUL(kvec + g, cell%bmat))
......@@ -154,8 +155,8 @@ CONTAINS
allocate(length_kG(hybrid%maxgptm1, kpts%nkptf))
length_kG = 0
unsrt_pgptm = 0
DO igpt = 1, hybrid%gptmd
g = hybrid%gptm(:,igpt)
DO igpt = 1, mpbasis%gptmd
g = mpbasis%gptm(:,igpt)
DO ikpt = 1, kpts%nkptf
kvec = kpts%bkf(:,ikpt)
rdum = SUM(MATMUL(kvec + g, cell%bmat)**2)
......@@ -185,10 +186,10 @@ CONTAINS
IF (mpi%irank == 0) THEN
WRITE (6, '(/A)') 'Mixed basis'
WRITE (6, '(A,I5)') 'Number of unique G-vectors: ', hybrid%gptmd
WRITE (6, '(A,I5)') 'Number of unique G-vectors: ', mpbasis%gptmd
WRITE (6, *)
WRITE (6, '(3x,A)') 'IR Plane-wave basis with cutoff of gcutm (hybrid%gcutm/2*input%rkmax):'
WRITE (6, '(5x,A,I5)') 'Maximal number of G-vectors:', maxval(hybrid%ngptm)
WRITE (6, '(5x,A,I5)') 'Maximal number of G-vectors:', maxval(mpbasis%ngptm)
WRITE (6, *)
WRITE (6, *)
WRITE (6, '(3x,A)') 'IR Plane-wave basis for non-local exchange potential:'
......@@ -574,9 +575,9 @@ CONTAINS
END DO
END DO
END DO
hybrid%maxbasm1 = hybrid%nbasp + maxval(hybrid%ngptm)
hybrid%maxbasm1 = hybrid%nbasp + maxval(mpbasis%ngptm)
DO nk = 1, kpts%nkptf
hybrid%nbasm(nk) = hybrid%nbasp + hybrid%ngptm(nk)
hybrid%nbasm(nk) = hybrid%nbasp + mpbasis%ngptm(nk)
END DO
hybrid%maxlmindx = MAXVAL([(SUM([(hybrid%num_radfun_per_l(l, itype)*(2*l + 1), l=0, atoms%lmax(itype))]), itype=1, atoms%ntype)])
......@@ -665,12 +666,13 @@ CONTAINS
END DO
end subroutine gen_bas_fun
subroutine gen_gvec(cell, kpts, hybrid)
subroutine gen_gvec(cell, kpts, mpbasis, hybrid)
use m_types
USE m_util, ONLY: intgrf_init, intgrf, rorderpf
implicit NONE
type(t_cell), intent(in) :: cell
type(t_kpts), intent(in) :: kpts
TYPE(t_mpbasis), intent(inout) :: mpbasis
type(t_hybrid), intent(inout) :: hybrid
......@@ -684,9 +686,9 @@ CONTAINS
REAL, ALLOCATABLE :: length_kg(:,:) ! length of the vectors k + G
INTEGER, ALLOCATABLE :: ptr(:)
allocate(hybrid%ngptm(kpts%nkptf))
allocate(mpbasis%ngptm(kpts%nkptf))
hybrid%ngptm = 0
mpbasis%ngptm = 0
i = 0
n = -1
......@@ -712,7 +714,7 @@ CONTAINS
l_found_kg_in_sphere = .TRUE.
END IF
hybrid%ngptm(ikpt) = hybrid%ngptm(ikpt) + 1
mpbasis%ngptm(ikpt) = mpbasis%ngptm(ikpt) + 1
l_found_new_gpt = .TRUE.
END IF
END DO ! k-loop
......@@ -722,17 +724,17 @@ CONTAINS
IF (.NOT. l_found_new_gpt) EXIT
END DO
hybrid%gptmd = i
allocate(hybrid%gptm(3, hybrid%gptmd))
allocate(hybrid%pgptm(maxval(hybrid%ngptm), kpts%nkptf))
mpbasis%gptmd = i
allocate(mpbasis%gptm(3, mpbasis%gptmd))
allocate(hybrid%pgptm(maxval(mpbasis%ngptm), kpts%nkptf))
! Allocate and initialize arrays needed for G vector ordering
allocate(unsrt_pgptm(maxval(hybrid%ngptm), kpts%nkptf))
allocate(length_kG(maxval(hybrid%ngptm), kpts%nkptf))
allocate(unsrt_pgptm(maxval(mpbasis%ngptm), kpts%nkptf))
allocate(length_kG(maxval(mpbasis%ngptm), kpts%nkptf))
hybrid%gptm = 0
mpbasis%gptm = 0
hybrid%pgptm = 0
hybrid%ngptm = 0
mpbasis%ngptm = 0
i = 0
n = -1
......@@ -757,17 +759,17 @@ CONTAINS
IF (norm2(MATMUL(kvec + g, cell%bmat)) <= hybrid%gcutm) THEN
IF (.NOT. l_found_kg_in_sphere) THEN
i = i + 1
hybrid%gptm(:,i) = g
mpbasis%gptm(:,i) = g
l_found_kg_in_sphere = .TRUE.
END IF
hybrid%ngptm(ikpt) = hybrid%ngptm(ikpt) + 1
mpbasis%ngptm(ikpt) = mpbasis%ngptm(ikpt) + 1
l_found_new_gpt = .TRUE.
! Position of the vector is saved as pointer
unsrt_pgptm(hybrid%ngptm(ikpt), ikpt) = i
unsrt_pgptm(mpbasis%ngptm(ikpt), ikpt) = i
! Save length of vector k + G for array sorting
length_kG(hybrid%ngptm(ikpt), ikpt) = norm2(MATMUL(kvec + g, cell%bmat))
length_kG(mpbasis%ngptm(ikpt), ikpt) = norm2(MATMUL(kvec + g, cell%bmat))
END IF
END DO
END DO
......@@ -778,13 +780,13 @@ CONTAINS
! Sort pointers in array, so that shortest |k+G| comes first
DO ikpt = 1, kpts%nkptf
allocate(ptr(hybrid%ngptm(ikpt)))
allocate(ptr(mpbasis%ngptm(ikpt)))
! Divide and conquer algorithm for arrays > 1000 entries
divconq = MAX(0, INT(1.443*LOG(0.001*hybrid%ngptm(ikpt))))
divconq = MAX(0, INT(1.443*LOG(0.001*mpbasis%ngptm(ikpt))))
! create pointers which correspond to a sorted array
CALL rorderpf(ptr, length_kG(1:hybrid%ngptm(ikpt), ikpt), hybrid%ngptm(ikpt), divconq)
CALL rorderpf(ptr, length_kG(1:mpbasis%ngptm(ikpt), ikpt), mpbasis%ngptm(ikpt), divconq)
! rearrange old pointers
DO igpt = 1, hybrid%ngptm(ikpt)
DO igpt = 1, mpbasis%ngptm(ikpt)
hybrid%pgptm(igpt, ikpt) = unsrt_pgptm(ptr(igpt), ikpt)
END DO
deallocate(ptr)
......
......@@ -192,7 +192,7 @@
END DO
END DO
DO igpt = 1, hybrid%ngptm(ikpt)
DO igpt = 1, mpbasis%ngptm(ikpt)
coulhlp1(ic2, ic4) = coulhlp(ic1, nbasp + igpt)
IF (sym%invs) THEN
coulhlp1(ic4, ic2) = coulhlp1(ic2, ic4)
......
......@@ -4,7 +4,7 @@ CONTAINS
!Note this module contains a real/complex version of spmvec
SUBROUTINE spmvec_invs(&
& atoms, hybrid,&
& atoms, mpbasis, hybrid,&
& hybdat, ikpt, kpts, cell,&
& coulomb_mt1, coulomb_mt2, coulomb_mt3,&
& coulomb_mtir, vecin,&
......@@ -16,6 +16,7 @@ CONTAINS
USE m_juDFT
IMPLICIT NONE
TYPE(t_hybdat), INTENT(IN) :: hybdat
TYPE(t_mpbasis), intent(in) :: mpbasis
TYPE(t_hybrid), INTENT(IN) :: hybrid
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_kpts), INTENT(IN) :: kpts
......@@ -134,7 +135,7 @@ CONTAINS
! compute vecout for the index-range from ibasm+1:nbasm
indx1 = sum((/(((2*l + 1)*atoms%neq(itype), l=0, hybrid%lcutm1(itype)),&
& itype=1, atoms%ntype)/)) + hybrid%ngptm(ikpt)
& itype=1, atoms%ntype)/)) + mpbasis%ngptm(ikpt)
CALL dspmv('U', indx1, 1.0, coulomb_mtir, vecinhlp(ibasm + 1:), 1, 0.0, vecout(ibasm + 1:), 1)
iatom = 0
......@@ -209,7 +210,7 @@ CONTAINS
END SUBROUTINE spmvec_invs
SUBROUTINE spmvec_noinvs(&
& atoms, hybrid,&
& atoms, mpbasis, hybrid,&
& hybdat, ikpt, kpts, cell,&
& coulomb_mt1, coulomb_mt2, coulomb_mt3,&
& coulomb_mtir, vecin,&
......@@ -221,6 +222,7 @@ CONTAINS
USE m_juDFT
IMPLICIT NONE
TYPE(t_hybdat), INTENT(IN) :: hybdat
TYPE(t_mpbasis), INTENT(IN) :: mpbasis
TYPE(t_hybrid), INTENT(IN) :: hybrid
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_kpts), INTENT(IN) :: kpts
......@@ -341,7 +343,7 @@ CONTAINS
! compute vecout for the index-range from ibasm+1:nbasm
indx1 = sum((/(((2*l + 1)*atoms%neq(itype), l=0, hybrid%lcutm1(itype)),&
& itype=1, atoms%ntype)/)) + hybrid%ngptm(ikpt)
& itype=1, atoms%ntype)/)) + mpbasis%ngptm(ikpt)
call zhpmv('U', indx1, (1.0, 0.0), coulomb_mtir, vecinhlp(ibasm + 1), 1, (0.0, 0.0), vecout(ibasm + 1), 1)
iatom = 0
......
......@@ -536,7 +536,7 @@ CONTAINS
SUBROUTINE bra_trafo2( &
l_real, vecout_r, vecin_r, vecout_c, vecin_c, &
dim, nobd, nbands, ikpt0, ikpt1, iop, sym, &
hybrid, kpts, cell, atoms, &
mpbasis, hybrid, kpts, cell, atoms, &
phase)
! ikpt0 :: parent of ikpt1
......@@ -547,6 +547,7 @@ CONTAINS
USE m_util
USE m_types
IMPLICIT NONE
type(t_mpbasis), intent(in) :: mpbasis
TYPE(t_hybrid), INTENT(IN) :: hybrid
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
......@@ -704,26 +705,26 @@ CONTAINS
END DO
! PW
DO igptm = 1, hybrid%ngptm(ikpt0)
DO igptm = 1, mpbasis%ngptm(ikpt0)
igptp = hybrid%pgptm(igptm, ikpt0)
g1 = matmul(rrot, hybrid%gptm(:, igptp)) + g
g1 = matmul(rrot, mpbasis%gptm(:, igptp)) + g
igptm2 = 0
DO i = 1, hybrid%ngptm(ikpt1)
IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt1)))) <= 1E-06) THEN
DO i = 1, mpbasis%ngptm(ikpt1)
IF (maxval(abs(g1 - mpbasis%gptm(:, hybrid%pgptm(i, ikpt1)))) <= 1E-06) THEN
igptm2 = i
EXIT
END IF
END DO
IF (igptm2 == 0) THEN
WRITE (*, *) ikpt0, ikpt1, g1
WRITE (*, *) hybrid%ngptm(ikpt0), hybrid%ngptm(ikpt1)
WRITE (*, *) mpbasis%ngptm(ikpt0), mpbasis%ngptm(ikpt1)
WRITE (*, *)
WRITE (*, *) igptp, hybrid%gptm(:, igptp)
WRITE (*, *) igptp, mpbasis%gptm(:, igptp)
WRITE (*, *) g
WRITE (*, *) rrot
WRITE (*, *) "Failed tests:", g1
DO i = 1, hybrid%ngptm(ikpt1)
WRITE (*, *) hybrid%gptm(:, hybrid%pgptm(i, ikpt1))
DO i = 1, mpbasis%ngptm(ikpt1)
WRITE (*, *) mpbasis%gptm(:, hybrid%pgptm(i, ikpt1))
ENDDO
call judft_error('bra_trafo2: G-point not found in G-point set.')
END IF
......@@ -765,7 +766,7 @@ CONTAINS
! This routine is not very fast at the moment.
SUBROUTINE matrixtrafo(matout, matin, ikpt0, isym, lsymmetrize, atoms, &
kpts, sym, hybrid, cell, maxindxm, nindxm, nbasm, ngptmall, nbasp, &
kpts, sym, mpbasis, hybrid, cell, maxindxm, nindxm, nbasm, ngptmall, nbasp, &
lcutm, maxlcutm)
USE m_wrapper
......@@ -774,6 +775,7 @@ CONTAINS
USE m_types
USE m_juDFT
IMPLICIT NONE
type(t_mpbasis), intent(in) :: mpbasis
TYPE(t_hybrid), INTENT(IN) :: hybrid
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
......@@ -806,7 +808,7 @@ CONTAINS
! - local arrays -
INTEGER :: pnt(maxindxm, 0:maxlcutm, atoms%nat),&
& g(3), g1(3), iarr(hybrid%ngptm(ikpt0))
& g(3), g1(3), iarr(mpbasis%ngptm(ikpt0))
INTEGER :: rot(3, 3), invrot(3, 3), rrot(3, 3), invrrot(3, 3)
REAL :: rkpt(3), rkpthlp(3), rtaual(3)
......@@ -820,7 +822,7 @@ CONTAINS
COMPLEX :: dwgninv(-maxlcutm:maxlcutm,&
& -maxlcutm:maxlcutm,&
& 0:maxlcutm)
COMPLEX :: carr(hybrid%ngptm(ikpt0))
COMPLEX :: carr(mpbasis%ngptm(ikpt0))
! Transform back to unsymmetrized product basis in case of inversion symmetry.
matin1 = matin
......@@ -943,12 +945,12 @@ CONTAINS
END DO
! PW
DO igptm = 1, hybrid%ngptm(ikpt1)
DO igptm = 1, mpbasis%ngptm(ikpt1)
igptp = hybrid%pgptm(igptm, ikpt1)
g1 = matmul(invrrot, hybrid%gptm(:, igptp) - g)
g1 = matmul(invrrot, mpbasis%gptm(:, igptp) - g)
igptm2 = 0
DO i = 1, hybrid%ngptm(ikpt0)
IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt0)))) <= 1E-06) THEN
DO i = 1, mpbasis%ngptm(ikpt0)
IF (maxval(abs(g1 - mpbasis%gptm(:, hybrid%pgptm(i, ikpt0)))) <= 1E-06) THEN
igptm2 = i
EXIT
END IF
......@@ -956,7 +958,7 @@ CONTAINS
! igptm2 = pntgptm(g1(1),g1(2),g1(3),ikpt0)
IF (igptm2 == 0) call judft_error('matrixtrafo: G point not found in G-point set.')
cdum = exp(img*tpi_const*dot_product(kpts%bkf(:, ikpt1) + hybrid%gptm(:, igptp), sym%tau(:, iisym)))
cdum = exp(img*tpi_const*dot_product(kpts%bkf(:, ikpt1) + mpbasis%gptm(:, igptp), sym%tau(:, iisym)))
matout1(:, nbasp + igptm) = cdum*matin1(:, nbasp + igptm2)
......@@ -993,12 +995,12 @@ CONTAINS
END DO
! PW
DO igptm = 1, hybrid%ngptm(ikpt1)
DO igptm = 1, mpbasis%ngptm(ikpt1)
igptp = hybrid%pgptm(igptm, ikpt1)
g1 = matmul(invrrot, hybrid%gptm(:, igptp) - g)
g1 = matmul(invrrot, mpbasis%gptm(:, igptp) - g)
igptm2 = 0
DO i = 1, hybrid%ngptm(ikpt0)
IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt0)))) <= 1E-06) THEN
DO i = 1, mpbasis%ngptm(ikpt0)
IF (maxval(abs(g1 - mpbasis%gptm(:, hybrid%pgptm(i, ikpt0)))) <= 1E-06) THEN
igptm2 = i
EXIT
END IF
......@@ -1007,12 +1009,12 @@ CONTAINS
IF (igptm2 == 0) call judft_error('matrixtrafo: G point not found in G-point set.')
iarr(igptm) = igptm2
carr(igptm) = exp(-img*tpi_const* &
& dot_product(kpts%bkf(:, ikpt1) + hybrid%gptm(:, igptp), sym%tau(:, iisym)))
& dot_product(kpts%bkf(:, ikpt1) + mpbasis%gptm(:, igptp), sym%tau(:, iisym)))
! cdum = exp(-img * 2*pi * dot_product(bk(:,ikpt1)+gptm(:,igptp),tau(:,isym)))
! matout1(nbasp+igptm,:) = cdum * matin1(nbasp+igptm2,:)
END DO
DO i2 = 1, nbasm(ikpt1)
DO i1 = 1, hybrid%ngptm(ikpt1)
DO i1 = 1, mpbasis%ngptm(ikpt1)
matout1(nbasp + i1, i2) = carr(i1)*matin1(nbasp + iarr(i1), i2)
END DO
END DO
......@@ -1029,7 +1031,7 @@ CONTAINS
SUBROUTINE matrixtrafo1( &
matout, matin, ikpt0, isym, lsymmetrize, atoms, kpts, sym, &
hybrid, cell, maxindxm, nindxm, nbasm, ngptmall, nbasp, lcutm, maxlcutm)
mpbasis, hybrid, cell, maxindxm, nindxm, nbasm, ngptmall, nbasp, lcutm, maxlcutm)
USE m_wrapper
USE m_dwigner
......@@ -1039,6 +1041,7 @@ CONTAINS
USE m_juDFT
IMPLICIT NONE
type(t_mpbasis), intent(in) :: mpbasis
TYPE(t_hybrid), INTENT(IN) :: hybrid
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
......@@ -1071,7 +1074,7 @@ CONTAINS
! - local arrays -
INTEGER :: pnt(maxindxm, 0:maxlcutm, atoms%nat), g(3),&
& g1(3), iarr(hybrid%ngptm(ikpt0))