Commit a4b6929e authored by Matthias Redies's avatar Matthias Redies

introduce fleur inputtype to coulomb matirx

parent f4e3d509
......@@ -104,7 +104,7 @@ CONTAINS
if(mpi%irank == 0) then
CALL open_hybinp_io2(mpdata, fi%hybinp, hybdat, fi%input, fi%atoms, fi%sym%invs)
CALL coulombmatrix(mpi, fi%atoms, fi%kpts, fi%cell, fi%sym, mpdata, fi%hybinp, hybdat, xcpot)
CALL coulombmatrix(mpi, fi, mpdata, hybdat, xcpot)
call close_hybinp_io2()
endif
......
......@@ -35,7 +35,7 @@ MODULE m_coulombmatrix
CONTAINS
SUBROUTINE coulombmatrix(mpi, atoms, kpts, cell, sym, mpdata, hybinp, hybdat, xcpot)
SUBROUTINE coulombmatrix(mpi, fi, mpdata, hybdat, xcpot)
USE m_types
USE m_types_hybdat
USE m_juDFT
......@@ -50,15 +50,11 @@ CONTAINS
use m_sphbes, only: sphbes
IMPLICIT NONE
TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_mpdata), intent(in) :: mpdata
TYPE(t_hybinp), INTENT(IN) :: hybinp
TYPE(t_hybdat), INTENT(IN) :: hybdat
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
TYPE(t_mpi), INTENT(IN) :: mpi
type(t_fleurinput), intent(in) :: fi
TYPE(t_mpdata), intent(in) :: mpdata
TYPE(t_hybdat), INTENT(IN) :: hybdat
! - local scalars -
......@@ -84,39 +80,39 @@ CONTAINS
! - local arrays -
INTEGER :: g(3)
INTEGER :: nbasm1(kpts%nkptf)
INTEGER :: nbasm1(fi%kpts%nkptf)
INTEGER, ALLOCATABLE :: pqnrm(:, :)
INTEGER :: rrot(3, 3, sym%nsym), invrrot(3, 3, sym%nsym)
INTEGER :: rrot(3, 3, fi%sym%nsym), invrrot(3, 3, fi%sym%nsym)
INTEGER, ALLOCATABLE :: iarr(:), POINTER(:, :, :, :)!,pointer(:,:,:)
INTEGER, ALLOCATABLE :: nsym_gpt(:, :), sym_gpt(:, :, :)
INTEGER :: nsym1(kpts%nkpt + 1), sym1(sym%nsym, kpts%nkpt + 1)
INTEGER :: nsym1(fi%kpts%nkpt + 1), sym1(fi%sym%nsym, fi%kpts%nkpt + 1)
INTEGER, ALLOCATABLE :: ngptm1(:)
INTEGER, ALLOCATABLE :: pgptm1(:,:)
LOGICAL :: calc_mt(kpts%nkpt)
LOGICAL :: calc_mt(fi%kpts%nkpt)
REAL :: q(3), q1(3), q2(3)
REAL :: integrand(atoms%jmtd), primf1(atoms%jmtd), primf2(atoms%jmtd)
REAL :: integrand(fi%atoms%jmtd), primf1(fi%atoms%jmtd), primf2(fi%atoms%jmtd)
REAL :: mat(maxval(mpdata%num_radbasfn)*(maxval(mpdata%num_radbasfn) + 1)/2)
REAL :: moment(maxval(mpdata%num_radbasfn), 0:maxval(hybinp%lcutm1), atoms%ntype), &
moment2(maxval(mpdata%num_radbasfn), atoms%ntype)
REAL :: sphbes_var(atoms%jmtd, 0:maxval(hybinp%lcutm1))
REAL :: sphbesmoment1(atoms%jmtd, 0:maxval(hybinp%lcutm1))
REAL :: rarr(0:hybinp%lexp + 1), rarr1(0:maxval(hybinp%lcutm1))
REAL :: moment(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype), &
moment2(maxval(mpdata%num_radbasfn), fi%atoms%ntype)
REAL :: sphbes_var(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
REAL :: sphbesmoment1(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
REAL :: rarr(0:fi%hybinp%lexp + 1), rarr1(0:maxval(fi%hybinp%lcutm1))
REAL, ALLOCATABLE :: gmat(:, :), qnrm(:)
REAL, ALLOCATABLE :: sphbesmoment(:, :, :)
REAL, ALLOCATABLE :: sphbes0(:, :, :)
REAL, ALLOCATABLE :: olap(:, :, :, :), integral(:, :, :, :)
REAL, ALLOCATABLE :: gridf(:, :)
REAL :: facA(0:MAX(2*atoms%lmaxd + maxval(hybinp%lcutm1) + 1, 4*MAX(maxval(hybinp%lcutm1), hybinp%lexp) + 1))
REAL :: facB(0:MAX(2*atoms%lmaxd + maxval(hybinp%lcutm1) + 1, 4*MAX(maxval(hybinp%lcutm1), hybinp%lexp) + 1))
REAL :: facC(-1:MAX(2*atoms%lmaxd + maxval(hybinp%lcutm1) + 1, 4*MAX(maxval(hybinp%lcutm1), hybinp%lexp) + 1))
COMPLEX :: cexp1(atoms%ntype), csumf(9)
COMPLEX :: structconst((2*hybinp%lexp + 1)**2, atoms%nat, atoms%nat, kpts%nkpt) ! nw = 1
COMPLEX :: y((hybinp%lexp + 1)**2), y1((hybinp%lexp + 1)**2), y2((hybinp%lexp + 1)**2)
COMPLEX :: dwgn(-maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), -maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), 0:maxval(hybinp%lcutm1), sym%nsym)
REAL :: facA(0:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
REAL :: facB(0:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
REAL :: facC(-1:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
COMPLEX :: cexp1(fi%atoms%ntype), csumf(9)
COMPLEX :: structconst((2*fi%hybinp%lexp + 1)**2, fi%atoms%nat, fi%atoms%nat, fi%kpts%nkpt) ! nw = 1
COMPLEX :: y((fi%hybinp%lexp + 1)**2), y1((fi%hybinp%lexp + 1)**2), y2((fi%hybinp%lexp + 1)**2)
COMPLEX :: dwgn(-maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), -maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), 0:maxval(fi%hybinp%lcutm1), fi%sym%nsym)
COMPLEX, ALLOCATABLE :: smat(:, :)
COMPLEX, ALLOCATABLE :: coulmat(:, :)
COMPLEX, ALLOCATABLE :: carr2(:, :), carr2a(:, :), carr2b(:, :)
......@@ -140,9 +136,9 @@ CONTAINS
call timestart("prep in coulomb")
if(mpi%is_root()) write (*,*) "start of coulomb calculation"
svol = SQRT(cell%vol)
fcoulfac = 4*pi_const/cell%vol
maxfac = MAX(2*atoms%lmaxd + maxval(hybinp%lcutm1) + 1, 4*MAX(maxval(hybinp%lcutm1), hybinp%lexp) + 1)
svol = SQRT(fi%cell%vol)
fcoulfac = 4*pi_const/fi%cell%vol
maxfac = MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1)
facA(0) = 1 !
facB(0) = 1 ! Define:
......@@ -153,12 +149,12 @@ CONTAINS
facC(i) = facC(i - 1)*(2*i + 1) !
END DO
CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
CALL intgrf_init(fi%atoms%ntype, fi%atoms%jmtd, fi%atoms%jri, fi%atoms%dx, fi%atoms%rmsh, gridf)
nbasm1 = hybdat%nbasp + mpdata%n_g(:)
! Calculate the structure constant
CALL structureconstant(structconst, cell, hybinp, atoms, kpts, mpi)
CALL structureconstant(structconst, fi%cell, fi%hybinp, fi%atoms, fi%kpts, mpi)
IF (mpi%irank == 0) WRITE (6, '(//A)') '### subroutine: coulombmatrix ###'
......@@ -169,7 +165,7 @@ CONTAINS
call timestart("coulomb allocation")
IF (ALLOCATED(coulomb)) deallocate(coulomb)
allocate(coulomb(hybdat%maxbasm1*(hybdat%maxbasm1 + 1)/2, kpts%nkpt), stat=ok, source=(0.0,0.0))
allocate(coulomb(hybdat%maxbasm1*(hybdat%maxbasm1 + 1)/2, fi%kpts%nkpt), stat=ok, source=(0.0,0.0))
IF (ok /= 0) call judft_error('coulombmatrix: failure allocation coulomb matrix')
call timestop("coulomb allocation")
......@@ -178,19 +174,19 @@ CONTAINS
! Generate Symmetry:
! Reduce list of g-Points so that only one of each symm-equivalent is calculated
! calculate rotations in reciprocal space
DO isym = 1, sym%nsym
IF (isym <= sym%nop) THEN
inviop = sym%invtab(isym)
rrot(:, :, isym) = TRANSPOSE(sym%mrot(:, :, inviop))
DO l = 0, maxval(hybinp%lcutm1)
dwgn(:, :, l, isym) = TRANSPOSE(hybinp%d_wgn2(-maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), &
-maxval(hybinp%lcutm1):maxval(hybinp%lcutm1), l, isym))
DO isym = 1, fi%sym%nsym
IF (isym <= fi%sym%nop) THEN
inviop = fi%sym%invtab(isym)
rrot(:, :, isym) = TRANSPOSE(fi%sym%mrot(:, :, inviop))
DO l = 0, maxval(fi%hybinp%lcutm1)
dwgn(:, :, l, isym) = TRANSPOSE(fi%hybinp%d_wgn2(-maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), &
-maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), l, isym))
END DO
ELSE
inviop = isym - sym%nop
inviop = isym - fi%sym%nop
rrot(:, :, isym) = -rrot(:, :, inviop)
dwgn(:, :, :, isym) = dwgn(:, :, :, inviop)
DO l = 0, maxval(hybinp%lcutm1)
DO l = 0, maxval(fi%hybinp%lcutm1)
DO m1 = -l, l
DO m2 = -l, -1
cdum = dwgn(m1, m2, l, isym)
......@@ -201,26 +197,26 @@ CONTAINS
END DO
END IF
END DO
invrrot(:, :, :sym%nop) = rrot(:, :, sym%invtab)
IF (sym%nsym > sym%nop) THEN
invrrot(:, :, sym%nop + 1:) = rrot(:, :, sym%invtab + sym%nop)
invrrot(:, :, :fi%sym%nop) = rrot(:, :, fi%sym%invtab)
IF (fi%sym%nsym > fi%sym%nop) THEN
invrrot(:, :, fi%sym%nop + 1:) = rrot(:, :, fi%sym%invtab + fi%sym%nop)
END IF
! Get symmetry operations that leave bk(:,ikpt) invariant -> sym1
nsym1 = 0
DO ikpt = 1, kpts%nkpt
DO ikpt = 1, fi%kpts%nkpt
isym1 = 0
DO isym = 1, sym%nsym
DO isym = 1, fi%sym%nsym
! temporary fix until bramat_trafo is correct
! for systems with symmetries including translations
IF (isym > sym%nop) THEN
isym2 = isym - sym%nop
IF (isym > fi%sym%nop) THEN
isym2 = isym - fi%sym%nop
ELSE
isym2 = isym
END IF
IF (ANY(abs(sym%tau(:, isym2)) > 1e-12)) CYCLE
IF (ANY(abs(fi%sym%tau(:, isym2)) > 1e-12)) CYCLE
IF (ALL(ABS(MATMUL(rrot(:, :, isym), kpts%bk(:, ikpt)) - kpts%bk(:, ikpt)) < 1e-12)) THEN
IF (ALL(ABS(MATMUL(rrot(:, :, isym), fi%kpts%bk(:, ikpt)) - fi%kpts%bk(:, ikpt)) < 1e-12)) THEN
isym1 = isym1 + 1
sym1(isym1, ikpt) = isym
END IF
......@@ -228,10 +224,10 @@ CONTAINS
nsym1(ikpt) = isym1
END DO
! Define reduced lists of G points -> pgptm1(:,ikpt), ikpt=1,..,nkpt
!if(allocated(pgptm1)) deallocate(hybinp%pgptm1)
allocate(pgptm1(maxval(mpdata%n_g),kpts%nkptf), source=0) !in mixedbasis
!if(allocated(pgptm1)) deallocate(fi%hybinp%pgptm1)
allocate(pgptm1(maxval(mpdata%n_g),fi%kpts%nkptf), source=0) !in mixedbasis
allocate(iarr(maxval(mpdata%n_g)), source=0)
allocate(POINTER(kpts%nkpt,&
allocate(POINTER(fi%kpts%nkpt,&
MINVAL(mpdata%g(1, :)) - 1:MAXVAL(mpdata%g(1, :)) + 1, &
MINVAL(mpdata%g(2, :)) - 1:MAXVAL(mpdata%g(2, :)) + 1, &
MINVAL(mpdata%g(3, :)) - 1:MAXVAL(mpdata%g(3, :)) + 1), &
......@@ -239,7 +235,7 @@ CONTAINS
allocate(ngptm1, mold=mpdata%n_g)
ngptm1 = 0
DO ikpt = 1, kpts%nkpt
DO ikpt = 1, fi%kpts%nkpt
DO igpt = 1, mpdata%n_g(ikpt)
g = mpdata%g(:, mpdata%gptm_ptr(igpt, ikpt))
POINTER(ikpt, g(1), g(2), g(3)) = igpt
......@@ -268,11 +264,11 @@ CONTAINS
call timestart("define gmat")
! Define gmat (symmetric)
i = (hybinp%lexp + 1)**2
i = (fi%hybinp%lexp + 1)**2
allocate(gmat(i, i))
gmat = 0
lm1 = 0
DO l1 = 0, hybinp%lexp
DO l1 = 0, fi%hybinp%lexp
DO m1 = -l1, l1
lm1 = lm1 + 1
lm2 = 0
......@@ -292,27 +288,27 @@ CONTAINS
! Calculate moments of MT functions
call timestart("calc moments of MT")
DO itype = 1, atoms%ntype
DO l = 0, hybinp%lcutm1(itype)
DO itype = 1, fi%atoms%ntype
DO l = 0, fi%hybinp%lcutm1(itype)
DO i = 1, mpdata%num_radbasfn(l, itype)
! note that mpdata%radbasfn_mt already contains the factor rgrid
moment(i, l, itype) = intgrf(atoms%rmsh(:, itype)**(l + 1)*mpdata%radbasfn_mt(:, i, l, itype), &
atoms, itype, gridf)
moment(i, l, itype) = intgrf(fi%atoms%rmsh(:, itype)**(l + 1)*mpdata%radbasfn_mt(:, i, l, itype), &
fi%atoms, itype, gridf)
END DO
END DO
DO i = 1, mpdata%num_radbasfn(0, itype)
moment2(i, itype) = intgrf(atoms%rmsh(:, itype)**3*mpdata%radbasfn_mt(:, i, 0, itype), &
atoms, itype, gridf)
moment2(i, itype) = intgrf(fi%atoms%rmsh(:, itype)**3*mpdata%radbasfn_mt(:, i, 0, itype), &
fi%atoms, itype, gridf)
END DO
END DO
call timestop("calc moments of MT")
call timestart("getnorm")
! Look for different qnorm = |k+G|, definition of qnrm and pqnrm.
CALL getnorm(kpts, mpdata%g, mpdata%n_g, mpdata%gptm_ptr, qnrm, nqnrm, pqnrm, cell)
allocate(sphbesmoment(0:hybinp%lexp, atoms%ntype, nqnrm), &
olap(maxval(mpdata%num_radbasfn), 0:maxval(hybinp%lcutm1), atoms%ntype, nqnrm), &
integral(maxval(mpdata%num_radbasfn), 0:maxval(hybinp%lcutm1), atoms%ntype, nqnrm))
CALL getnorm(fi%kpts, mpdata%g, mpdata%n_g, mpdata%gptm_ptr, qnrm, nqnrm, pqnrm, fi%cell)
allocate(sphbesmoment(0:fi%hybinp%lexp, fi%atoms%ntype, nqnrm), &
olap(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype, nqnrm), &
integral(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype, nqnrm))
sphbes_var = 0
sphbesmoment = 0
sphbesmoment1 = 0
......@@ -334,26 +330,26 @@ CONTAINS
!DO iqnrm = iqnrmstart, nqnrm, iqnrmstep
do iqnrm = 1,nqnrm
qnorm = qnrm(iqnrm)
DO itype = 1, atoms%ntype
ng = atoms%jri(itype)
rdum = atoms%rmt(itype)
DO itype = 1, fi%atoms%ntype
ng = fi%atoms%jri(itype)
rdum = fi%atoms%rmt(itype)
sphbes_var = 0
sphbesmoment1 = 0
IF (abs(qnorm) < 1e-12) THEN
sphbesmoment(0, itype, iqnrm) = rdum**3/3
DO i = 1, ng
sphbes_var(i, 0) = 1
sphbesmoment1(i, 0) = atoms%rmsh(i, itype)**2/3 + (rdum**2 - atoms%rmsh(i, itype)**2)/2
sphbesmoment1(i, 0) = fi%atoms%rmsh(i, itype)**2/3 + (rdum**2 - fi%atoms%rmsh(i, itype)**2)/2
END DO
ELSE
call sphbes(hybinp%lexp + 1, qnorm*rdum, rarr)
DO l = 0, hybinp%lexp
call sphbes(fi%hybinp%lexp + 1, qnorm*rdum, rarr)
DO l = 0, fi%hybinp%lexp
sphbesmoment(l, itype, iqnrm) = rdum**(l + 2)*rarr(l + 1)/qnorm
END DO
DO i = ng, 1, -1
rdum = atoms%rmsh(i, itype)
call sphbes(hybinp%lcutm1(itype) + 1, qnorm*rdum, rarr)
DO l = 0, hybinp%lcutm1(itype)
rdum = fi%atoms%rmsh(i, itype)
call sphbes(fi%hybinp%lcutm1(itype) + 1, qnorm*rdum, rarr)
DO l = 0, fi%hybinp%lcutm1(itype)
sphbes_var(i, l) = rarr(l)
IF (l /= 0) THEN; rdum1 = -rdum**(1 - l)*rarr(l - 1)
ELSE; rdum1 = -COS(qnorm*rdum)/qnorm
......@@ -364,16 +360,16 @@ CONTAINS
END DO
END DO
END IF
DO l = 0, hybinp%lcutm1(itype)
DO l = 0, fi%hybinp%lcutm1(itype)
DO n = 1, mpdata%num_radbasfn(l, itype)
! note that mpdata%radbasfn_mt already contains one factor rgrid
olap(n, l, itype, iqnrm) = &
intgrf(atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, n, l, itype)*sphbes_var(:, l), &
atoms, itype, gridf)
intgrf(fi%atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, n, l, itype)*sphbes_var(:, l), &
fi%atoms, itype, gridf)
integral(n, l, itype, iqnrm) = &
intgrf(atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, n, l, itype)*sphbesmoment1(:, l), &
atoms, itype, gridf)
intgrf(fi%atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, n, l, itype)*sphbesmoment1(:, l), &
fi%atoms, itype, gridf)
END DO
END DO
......@@ -392,28 +388,28 @@ CONTAINS
ix = 0
iy = 0
iy0 = 0
DO itype = 1, atoms%ntype
DO ineq = 1, atoms%neq(itype)
DO itype = 1, fi%atoms%ntype
DO ineq = 1, fi%atoms%neq(itype)
! Here the diagonal block matrices do not depend on ineq. In (1b) they do depend on ineq, though,
DO l = 0, hybinp%lcutm1(itype)
DO l = 0, fi%hybinp%lcutm1(itype)
DO n2 = 1, mpdata%num_radbasfn(l, itype)
! note that mpdata%radbasfn_mt already contains the factor rgrid
CALL primitivef(primf1, mpdata%radbasfn_mt(:, n2, l, itype) &
*atoms%rmsh(:, itype)**(l + 1), atoms%rmsh, atoms%dx, &
atoms%jri, atoms%jmtd, itype, atoms%ntype)
*fi%atoms%rmsh(:, itype)**(l + 1), fi%atoms%rmsh, fi%atoms%dx, &
fi%atoms%jri, fi%atoms%jmtd, itype, fi%atoms%ntype)
! -itype is to enforce inward integration
CALL primitivef(primf2, mpdata%radbasfn_mt(:atoms%jri(itype), n2, l, itype) &
/atoms%rmsh(:atoms%jri(itype), itype)**l, atoms%rmsh, atoms%dx, &
atoms%jri, atoms%jmtd, -itype, atoms%ntype)
CALL primitivef(primf2, mpdata%radbasfn_mt(:fi%atoms%jri(itype), n2, l, itype) &
/fi%atoms%rmsh(:fi%atoms%jri(itype), itype)**l, fi%atoms%rmsh, fi%atoms%dx, &
fi%atoms%jri, fi%atoms%jmtd, -itype, fi%atoms%ntype)
primf1(:atoms%jri(itype)) = primf1(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype), itype)**l
primf2 = primf2*atoms%rmsh(:, itype)**(l + 1)
primf1(:fi%atoms%jri(itype)) = primf1(:fi%atoms%jri(itype))/fi%atoms%rmsh(:fi%atoms%jri(itype), itype)**l
primf2 = primf2*fi%atoms%rmsh(:, itype)**(l + 1)
DO n1 = 1, n2
integrand = mpdata%radbasfn_mt(:, n1, l, itype)*(primf1 + primf2)
! call intgr0( (4*pimach())/(2*l+1)*integrand,rmsh(1,itype),dx(itype),jri(itype),mat(n2*(n2-1)/2+n1) )
mat(n2*(n2 - 1)/2 + n1) = (4*pi_const)/(2*l + 1) &
*intgrf(integrand, atoms, itype, gridf)
*intgrf(integrand, fi%atoms, itype, gridf)
END DO
END DO
......@@ -426,7 +422,7 @@ CONTAINS
iy = iy + 1
i = ix*(ix - 1)/2 + iy
j = n2*(n2 - 1)/2 + n1
coulomb(i, kpts%nkpt) = mat(j)
coulomb(i, fi%kpts%nkpt) = mat(j)
END DO
END DO
iy0 = ix
......@@ -445,7 +441,7 @@ CONTAINS
END IF
DO ikpt = 1, kpts%nkpt
DO ikpt = 1, fi%kpts%nkpt
! only the first rank handles the MT-MT part
call timestart("MT-MT part")
......@@ -453,11 +449,11 @@ CONTAINS
ix = 0
ic2 = 0
DO itype2 = 1, atoms%ntype
DO ineq2 = 1, atoms%neq(itype2)
DO itype2 = 1, fi%atoms%ntype
DO ineq2 = 1, fi%atoms%neq(itype2)
ic2 = ic2 + 1
lm2 = 0
DO l2 = 0, hybinp%lcutm1(itype2)
DO l2 = 0, fi%hybinp%lcutm1(itype2)
DO m2 = -l2, l2
lm2 = lm2 + 1
DO n2 = 1, mpdata%num_radbasfn(l2, itype2)
......@@ -466,10 +462,10 @@ CONTAINS
iy = 0
ic1 = 0
lp2: DO itype1 = 1, itype2
DO ineq1 = 1, atoms%neq(itype1)
DO ineq1 = 1, fi%atoms%neq(itype1)
ic1 = ic1 + 1
lm1 = 0
DO l1 = 0, hybinp%lcutm1(itype1)
DO l1 = 0, fi%hybinp%lcutm1(itype1)
DO m1 = -l1, l1
lm1 = lm1 + 1
DO n1 = 1, mpdata%num_radbasfn(l1, itype1)
......@@ -479,10 +475,10 @@ CONTAINS
l = l1 + l2
lm = l**2 + l + m1 - m2 + 1
idum = ix*(ix - 1)/2 + iy
coulmat(iy, ix) = coulomb(idum, kpts%nkpt) &
coulmat(iy, ix) = coulomb(idum, fi%kpts%nkpt) &
+ EXP(CMPLX(0.0, 1.0)*2*pi_const* &
dot_PRODUCT(kpts%bk(:, ikpt), &
atoms%taual(:, ic2) - atoms%taual(:, ic1))) &
dot_PRODUCT(fi%kpts%bk(:, ikpt), &
fi%atoms%taual(:, ic2) - fi%atoms%taual(:, ic1))) &
*rdum*structconst(lm, ic1, ic2, ikpt)
coulmat(ix, iy) = CONJG(coulmat(iy, ix))
END DO
......@@ -497,11 +493,11 @@ CONTAINS
END DO
END DO
IF (sym%invs) THEN
IF (fi%sym%invs) THEN
!symmetrize makes the Coulomb matrix real symmetric
CALL symmetrize(coulmat, hybdat%nbasp, hybdat%nbasp, 3, .FALSE., &
atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
mpdata%num_radbasfn, sym)
fi%atoms, fi%hybinp%lcutm1, maxval(fi%hybinp%lcutm1), &
mpdata%num_radbasfn, fi%sym)
ENDIF
coulomb(:hybdat%nbasp*(hybdat%nbasp + 1)/2, ikpt) = packmat(coulmat)
......@@ -527,7 +523,7 @@ CONTAINS
coulmat = 0
call timestart("loop over interst.")
DO ikpt = 1, kpts%nkpt !1,kpts%nkpt
DO ikpt = 1, fi%kpts%nkpt !1,fi%kpts%nkpt
coulmat = 0
! start to loop over interstitial plane waves
......@@ -535,25 +531,25 @@ CONTAINS
igpt = pgptm1(igpt0, ikpt)
igptp = mpdata%gptm_ptr(igpt, ikpt)
ix = hybdat%nbasp + igpt
q = MATMUL(kpts%bk(:, ikpt) + mpdata%g(:, igptp), cell%bmat)
q = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%cell%bmat)
qnorm = norm2(q)
iqnrm = pqnrm(igpt, ikpt)
IF (ABS(qnrm(iqnrm) - qnorm) > 1e-12) then
call judft_error('coulombmatrix: qnorm does not equal corresponding & element in qnrm (bug?)') ! We shouldn't stop here!
endif
call ylm4(2, MATMUL(kpts%bk(:, kpts%nkpt), cell%bmat), y1)
call ylm4(2, MATMUL(mpdata%g(:, igptp), cell%bmat), y2)
call ylm4(hybinp%lexp, q, y)
call ylm4(2, MATMUL(fi%kpts%bk(:, fi%kpts%nkpt), fi%cell%bmat), y1)
call ylm4(2, MATMUL(mpdata%g(:, igptp), fi%cell%bmat), y2)
call ylm4(fi%hybinp%lexp, q, y)
y1 = CONJG(y1); y2 = CONJG(y2); y = CONJG(y)
iy = 0
ic = 0
DO itype = 1, atoms%ntype
DO ineq = 1, atoms%neq(itype)
DO itype = 1, fi%atoms%ntype
DO ineq = 1, fi%atoms%neq(itype)
ic = ic + 1
lm = 0
DO l = 0, hybinp%lcutm1(itype)
DO l = 0, fi%hybinp%lcutm1(itype)
DO M = -l, l
lm = lm + 1
......@@ -561,15 +557,15 @@ CONTAINS
csum = 0
csumf = 0
ic1 = 0
DO itype1 = 1, atoms%ntype
DO ineq1 = 1, atoms%neq(itype1)
DO itype1 = 1, fi%atoms%ntype
DO ineq1 = 1, fi%atoms%neq(itype1)
ic1 = ic1 + 1
cexp = 4*pi_const*EXP(CMPLX(0.0, 1.0)*2*pi_const &
*(dot_PRODUCT(kpts%bk(:, ikpt) + mpdata%g(:, igptp), atoms%taual(:, ic1)) &
- dot_PRODUCT(kpts%bk(:, ikpt), atoms%taual(:, ic))))
*(dot_PRODUCT(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%atoms%taual(:, ic1)) &
- dot_PRODUCT(fi%kpts%bk(:, ikpt), fi%atoms%taual(:, ic))))
lm1 = 0
DO l1 = 0, hybinp%lexp
DO l1 = 0, fi%hybinp%lexp
l2 = l + l1 ! for structconst
idum = 1
cdum = sphbesmoment(l1, itype1, iqnrm)*CMPLX(0.0, 1.0)**(l1)*cexp
......@@ -584,16 +580,16 @@ CONTAINS
! add contribution of (2c) to csum and csumf coming from linear and quadratic orders of Y_lm*(G) / G * j_(l+1)(GS)
IF (ikpt == 1 .AND. l <= 2) THEN
cexp = EXP(CMPLX(0.0, 1.0)*2*pi_const*dot_PRODUCT(mpdata%g(:, igptp), atoms%taual(:, ic1))) &
*gmat(lm, 1)*4*pi_const/cell%vol
cexp = EXP(CMPLX(0.0, 1.0)*2*pi_const*dot_PRODUCT(mpdata%g(:, igptp), fi%atoms%taual(:, ic1))) &
*gmat(lm, 1)*4*pi_const/fi%cell%vol
csumf(lm) = csumf(lm) - cexp*SQRT(4*pi_const)* &
CMPLX(0.0, 1.0)**l*sphbesmoment(0, itype1, iqnrm)/facC(l - 1)
IF (l == 0) THEN
IF (igpt /= 1) THEN
csum = csum - cexp*(sphbesmoment(0, itype1, iqnrm)*atoms%rmt(itype1)**2 - &
csum = csum - cexp*(sphbesmoment(0, itype1, iqnrm)*fi%atoms%rmt(itype1)**2 - &
sphbesmoment(2, itype1, iqnrm)*2.0/3)/10
ELSE
csum = csum - cexp*atoms%rmt(itype1)**5/30
csum = csum - cexp*fi%atoms%rmt(itype1)**5/30
END IF
ELSE IF (l == 1) THEN
csum = csum + cexp*CMPLX(0.0, 1.0)*SQRT(4*pi_const) &
......@@ -613,7 +609,7 @@ CONTAINS
idum = ix*(ix - 1)/2
cdum = (4*pi_const)**2*CMPLX(0.0, 1.0)**(l)*y(lm) &
*EXP(CMPLX(0.0, 1.0)*2*pi_const &
*dot_PRODUCT(mpdata%g(:, igptp), atoms%taual(:, ic)))
*dot_PRODUCT(mpdata%g(:, igptp), fi%atoms%taual(:, ic)))
DO n = 1, mpdata%num_radbasfn(l, itype)
iy = iy + 1
......@@ -640,9 +636,9 @@ CONTAINS
END DO
END DO
IF (sym%invs) THEN
IF (