Commit 73087055 authored by Daniel Wortmann's avatar Daniel Wortmann

More bugfixes

Removed most of CPP_INVERSION
parent 2d811c6b
......@@ -76,7 +76,7 @@ CONTAINS
CALL setabc1locdn1(jspin, atoms,lapw, sym,usdus,kveclo,enough,nkvec,kvec,&
nbasf0,alo1,blo1,clo1)
nvmax=lapw%nv(jspin)
!---> loop over lapws
DO k = 1,nvmax
!calculate k+G
......
......@@ -53,7 +53,7 @@ CONTAINS
ALLOCATE ( lapw%k1(DIMENSION%nvd,DIMENSION%jspd),lapw%k2(DIMENSION%nvd,DIMENSION%jspd),&
lapw%k3(DIMENSION%nvd,DIMENSION%jspd),lapw%rk(DIMENSION%nvd,DIMENSION%jspd) )
ENDIF
lapw%rk = 0 ; lapw%k1 = 0 ; lapw%k2 = 0 ; lapw%k3 = 0
lapw%rk = 0 ; lapw%k1 = 0 ; lapw%k2 = 0 ; lapw%k3 = 0 ;lapw%nv=0
! ..
! ..
!---> in a spin-spiral calculation different basis sets are used for
......
......@@ -9,7 +9,7 @@ MODULE m_eigen
CONTAINS
SUBROUTINE eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
sym,kpts,DIMENSION, vacuum, input, cell, enpara_in,banddos, noco,jij, oneD,hybrid,&
it,eig_id,results,v)
it,eig_id,results,v,vx)
!*********************************************************************
! sets up and solves the eigenvalue problem for a basis of lapws.
!
......@@ -65,7 +65,7 @@ CONTAINS
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(INOUT) :: atoms!in u_setup n_u might be modified
TYPE(t_potden),INTENT(INOUT) :: v
TYPE(t_potden),INTENT(INOUT) :: v,vx
#ifdef CPP_MPI
INCLUDE 'mpif.h'
#endif
......@@ -291,6 +291,7 @@ CONTAINS
WRITE (*,*) 'the tlmplm%tuu, tlmplm%tdd etc.: ',err,' size: ',mlotot
CALL juDFT_error("eigen: Error during allocation of tlmplm, tdd etc.",calledby ="eigen")
ENDIF
lh0=1
CALL tlmplm(sphhar,atoms,DIMENSION,enpara, jsp,1,mpi, v%mt(1,0,1,jsp),lh0,input, td,ud)
IF (input%l_f) CALL write_tlmplm(td,vs_mmp,atoms%n_u>0,1,jsp,input%jspins)
CALL timestop("tlmplm")
......@@ -358,7 +359,7 @@ CONTAINS
IF( hybrid%l_hybrid ) THEN
!write overlap matrix b to direct access file olap
print *,"Wrong overlap matrix used, fix this later"
call olap%from_packed(dimension%nbasfcn,l_real,hamovlp%b_r,hamovlp%b_c)
call olap%from_packed(l_real,dimension%nbasfcn,hamovlp%b_r,hamovlp%b_c)
call write_olap(olap,nrec)
......@@ -371,7 +372,7 @@ CONTAINS
IF( hybrid%l_subvxc ) THEN
CALL subvxc(lapw,kpts%bk(:,nk),DIMENSION,input,jsp,vr0,atoms,ud,hybrid,enpara%el0,enpara%ello0,&
sym, atoms%nlotot,kveclo, cell,sphhar, stars, xcpot,mpi,&
oneD, hamovlp)
oneD, hamovlp,vx)
END IF
END IF ! hybrid%l_hybrid
......@@ -419,13 +420,13 @@ CONTAINS
zmat%nbands=ne_found
CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,lapw%nv(jsp),lapw%nmat,&
lapw%k1(:lapw%nv(jsp),jsp),lapw%k2 (:lapw%nv(jsp),jsp),lapw%k3(:lapw%nv(jsp),jsp),&
bkpt, kpts%wtkpt(nk),eig(:ne_found),enpara%el0(0:,:,jsp), enpara%ello0(:,:,jsp),enpara%evac0(:,jsp),&
atoms%nlotot,kveclo,mpi%n_size,mpi%n_rank,zMat)
bkpt, kpts%wtkpt(nk),eig(:ne_found),el=enpara%el0(0:,:,jsp),ello=enpara%ello0(:,:,jsp),evac=enpara%evac0(:,jsp),&
nlotot=atoms%nlotot,kveclo=kveclo,n_start=mpi%n_size,n_end=mpi%n_rank,zmat=zMat)
IF (noco%l_noco) THEN
CALL write_eig(eig_id, nk,2,ne_found,ne_all,lapw%nv(2),lapw%nmat,&
lapw%k1(:lapw%nv(2),2),lapw%k2 (:lapw%nv(2),2),lapw%k3(:lapw%nv(2),2),&
bkpt, kpts%wtkpt(nk),eig(:ne_found),enpara%el0(0:,:,2), enpara%ello0(:,:,2),enpara%evac0(:,2),&
atoms%nlotot,kveclo)
bkpt, kpts%wtkpt(nk),eig(:ne_found),el=enpara%el0(0:,:,2),ello= enpara%ello0(:,:,2),evac=enpara%evac0(:,2),&
nlotot=atoms%nlotot,kveclo=kveclo)
ENDIF
#if defined(CPP_MPI)
!RMA synchronization
......
......@@ -28,7 +28,7 @@ CONTAINS
!-----------------------------------------------------------------------
USE m_eig66_io, ONLY : read_eig
USE m_eig66_io, ONLY : read_eig,write_eig
#if defined(CPP_MPI)&&defined(CPP_NEVER)
USE m_mpi_col_eigJ
#endif
......@@ -235,6 +235,13 @@ CONTAINS
WRITE(attributes(2),'(a)') 'Htr'
IF (mpi%irank.EQ.0) CALL writeXMLElement('FermiEnergy',(/'value','units'/),attributes(1:2))
!Put w_iks into eig-file
DO jsp = 1,nspins
DO k = 1,kpts%nkpt
CALL write_eig(eig_id,k,jsp,w_iks=results%w_iks(:,k,jsp))
ENDDO
ENDDO
RETURN
8020 FORMAT (/,'FERMIE:',/,&
& 10x,'first approx. to ef (T=0) :',f10.6,' htr',&
......
......@@ -401,6 +401,7 @@ MODULE m_types
INTEGER :: lexp
INTEGER :: bands1 !Only read in
INTEGER :: bands2 !Only read in
INTEGER :: nbasp
INTEGER :: maxlcutm1
INTEGER :: maxindxm1
INTEGER :: maxbasm1
......@@ -452,7 +453,6 @@ MODULE m_types
INTEGER :: l1,l2,n1,n2
END TYPE prodtype
TYPE t_hybdat
INTEGER :: nbasp
INTEGER :: lmaxcd,maxindxc
REAL, ALLOCATABLE :: gridf(:,:) !alloc in util.F
INTEGER , ALLOCATABLE:: nindxc(:,:) !alloc in eigen_HF_init
......
......@@ -12,6 +12,8 @@ module m_types_rcmat
PROCEDURE :: multiply=>t_mat_multiply
PROCEDURE :: transpose=>t_mat_transpose
PROCEDURE :: from_packed=>t_mat_from_packed
PROCEDURE :: inverse =>t_mat_inverse
PROCEDURE :: to_packed=>t_mat_to_packed
END type t_mat
CONTAINS
......@@ -33,7 +35,9 @@ module m_types_rcmat
IF (mat%l_real) THEN
ALLOCATE(mat%data_r(mat%matsize1,mat%matsize2),STAT=err)
ALLOCATE(mat%data_c(0,0))
ELSE
ALLOCATE(mat%data_r(0,0))
ALLOCATE(mat%data_c(mat%matsize1,mat%matsize2),STAT=err)
ENDIF
......@@ -86,7 +90,7 @@ module m_types_rcmat
end IF
end SUBROUTINE t_mat_transpose
SUBROUTINE t_mat_from_packed(mat1,matsize,l_real,packed_r,packed_c)
SUBROUTINE t_mat_from_packed(mat1,l_real,matsize,packed_r,packed_c)
CLASS(t_mat),INTENT(INOUT) :: mat1
INTEGER,INTENT(IN) :: matsize
LOGICAL,INTENT(IN) :: l_real
......@@ -94,17 +98,68 @@ module m_types_rcmat
COMPLEX,INTENT(IN) :: packed_c(:)
INTEGER:: n,nn,i
call mat1%alloc(l_real,matsize)
call mat1%alloc(l_real,matsize,matsize)
i=1
DO n=1,matsize
DO nn=1,n
if (l_real) THEN
mat1%data_r(n,nn)=packed_r(i)
mat1%data_r(nn,n)=packed_r(i)
else
mat1%data_c(n,nn)=packed_c(i)
mat1%data_c(n,nn)=conjg(packed_c(i))
mat1%data_c(nn,n)=packed_c(i)
end if
i=i+1
end DO
end DO
end SUBROUTINE t_mat_from_packed
function t_mat_to_packed(mat)result(packed)
CLASS(t_mat),INTENT(IN) :: mat
COMPLEX :: packed(mat%matsize1*(mat%matsize1+1)/2)
integer :: n,nn,i
real,parameter :: tol=1e-5
if (mat%matsize1.ne.mat%matsize2) call judft_error("Could not pack no-square matrix",hint='This is a BUG, please report')
i=1
DO n=1,mat%matsize1
DO nn=1,n
if (mat%l_real) THEN
packed(i)=(mat%data_r(n,nn)+mat%data_r(nn,n))/2.
if (abs(mat%data_r(n,nn)-mat%data_r(nn,n))>tol) call judft_warn("Large unsymmetry in matrix packing")
else
packed(i)=(conjg(mat%data_c(n,nn))+mat%data_c(nn,n))/2.
if (abs(conjg(mat%data_c(n,nn))-mat%data_c(nn,n))>tol) call judft_warn("Large unsymmetry in matrix packing")
endif
i=i+1
end DO
end DO
end function t_mat_to_packed
subroutine t_mat_inverse(mat)
implicit none
CLASS(t_mat),INTENT(INOUT) :: mat
integer :: info
real, allocatable :: work_r(:)
integer, allocatable :: ipiv(:)
complex,allocatable :: work_c(:)
if (mat%matsize1.ne.mat%matsize2) call judft_error("Can only invert square matrices",hint="This is a BUG in FLEUR, please report")
ALLOCATE(ipiv(mat%matsize1))
if (mat%l_real) THEN
ALLOCATE(work_r(mat%matsize1))
call dgetrf(mat%matsize1,mat%matsize1,mat%data_r,size(mat%data_r,1),ipiv,info)
if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
call dgetri(mat%matsize1,mat%data_r,size(mat%data_r,1),ipiv,work_r,size(work_r),info)
if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
else
ALLOCATE(work_c(mat%matsize1))
call zgetrf(mat%matsize1,mat%matsize1,mat%data_c,size(mat%data_c,1),ipiv,info)
if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
call zgetri(mat%matsize1,mat%data_c,size(mat%data_c,1),ipiv,work_c,size(work_c),info)
if(info.ne.0) call judft_error("Failed to invert matrix: dpotrf failed.")
end if
end subroutine t_mat_inverse
end module m_types_rcmat
......@@ -41,7 +41,7 @@
! - local scalars -
INTEGER :: i,itype,iatom,ikpt,ineq,igpt,iband
INTEGER :: irecl_cmt,irecl_z
INTEGER :: irecl_cmt
INTEGER :: j,m
INTEGER :: l
INTEGER :: lm,lm1
......@@ -66,13 +66,7 @@
REAL :: rtaual(3)
REAL , ALLOCATABLE :: olapcb(:)
REAL , ALLOCATABLE :: olapcv_avg(:,:,:,:),olapcv_max(:,:,:,:)
#ifdef CPP_INVERSION
REAL , ALLOCATABLE :: olappp(:,:)
REAL :: z(dimension%nbasfcn,dimension%neigd,nkpti)
#else
COMPLEX , ALLOCATABLE :: olappp(:,:)
COMPLEX :: z(dimension%nbasfcn,dimension%neigd,nkpti)
#endif
TYPE(t_mat),ALLOCATABLE :: z(:)
COMPLEX :: cmt(dimension%neigd,hybrid%maxlmindx,atoms%nat,nkpti)
COMPLEX :: y((atoms%lmaxd+1)**2)
......@@ -85,7 +79,11 @@
& 'x','x','x','x','x','x','x','x','x','x','x','x','x' /)
LOGICAL :: l_mism = .true.
ALLOCATE(z(nkpti))
DO ikpt=1,nkpti
call z(ikpt)%alloc(sym%invs,dimension%nbasfcn,dimension%neigd)
ENDDO
IF ( mpi%irank == 0 ) WRITE(6,'(//A)') '### checkolap ###'
cmt = 0
......@@ -238,18 +236,9 @@
ALLOCATE (carr1(maxval(hybdat%nbands),(atoms%lmaxd+1)**2))
ALLOCATE (carr2(maxval(hybdat%nbands),(atoms%lmaxd+1)**2))
ALLOCATE (carr3(maxval(hybdat%nbands),(atoms%lmaxd+1)**2))
#ifdef CPP_INVERSION
irecl_z = dimension%nbasfcn*dimension%neigd*8
#else
irecl_z = dimension%nbasfcn*dimension%neigd*16
#endif
OPEN(unit=778,file='z',form='unformatted',access='direct',&
& recl=irecl_z)
DO ikpt = 1,nkpti
READ(778,rec=ikpt) z(:,:,ikpt)
call read_z(z(ikpt),ikpt)
END DO
CLOSE(778)
iatom = 0
DO itype = 1,atoms%ntype
......@@ -292,8 +281,12 @@
DO m = -l,l
lm = lm + 1
DO iband = 1,hybdat%nbands(ikpt)
carr2(iband,lm) = carr2(iband,lm) + cdum * z(igpt,iband,ikpt) * y(lm)
END DO
if (z(1)%l_real) THEN
carr2(iband,lm) = carr2(iband,lm) + cdum * z(ikpt)%data_r(igpt,iband) * y(lm)
Else
carr2(iband,lm) = carr2(iband,lm) + cdum * z(ikpt)%data_c(igpt,iband) * y(lm)
END if
end DO
END DO
END DO
END DO
......
......@@ -26,1662 +26,1245 @@
! In the code: coulomb(IJ,1) = v where only the upper triangle (I<=J) is stored,
! IJ
! (1)
! coulfac(IJ,lm) = v
! IJ,lm
! coulfac(IJ,lm) = v IJ,lm
!
! For the PW contribution we have to construct plane waves within the MT spheres with the help
! of spherical Bessel functions. The value lexp (LEXP in gwinp) is the corresponding cutoff.
!
MODULE m_coulomb
MODULE m_coulomb
CONTAINS
CONTAINS
SUBROUTINE coulombmatrix(mpi,atoms,kpts,cell, sym, hybrid, xcpot,l_restart)
USE m_constants , ONLY : pi_const
USE m_olap , ONLY : olap_pw,gptnorm
USE m_trafo , ONLY : symmetrize,bramat_trafo
USE m_util , ONLY : sphbessel,intgrf,intgrf_init, harmonicsr,primitivef
USE m_hsefunctional, ONLY : change_coulombmatrix
#ifdef CPP_MPI
USE m_mpi_work_dist, ONLY : work_dist_coulomb
#endif
USE m_wrapper
USE m_icorrkeys
SUBROUTINE coulombmatrix(mpi,atoms,kpts,cell, sym, hybrid, xcpot,l_restart)
USE m_constants , ONLY : pi_const
USE m_olap , ONLY : olap_pw,gptnorm
USE m_trafo , ONLY : symmetrize,bramat_trafo
USE m_util , ONLY : sphbessel,intgrf,intgrf_init, harmonicsr,primitivef
USE m_hsefunctional, ONLY : change_coulombmatrix
USE m_wrapper
USE m_icorrkeys
USE m_io_hybrid
USE m_types
IMPLICIT NONE
TYPE(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_hybrid),INTENT(INOUT) :: hybrid
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_atoms),INTENT(IN) :: atoms
IMPLICIT NONE
TYPE(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_hybrid),INTENT(INOUT) :: hybrid
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_atoms),INTENT(IN) :: atoms
! - scalars -
LOGICAL , INTENT(IN) :: l_restart
! - local scalars -
INTEGER :: maxfac
INTEGER :: inviop
INTEGER :: nqnrm,iqnrm,iqnrm1,iqnrm2, iqnrmstart,iqnrmstep
INTEGER :: itype,l ,ix,iy,iy0,i,j,lm,l1,l2,m1, m2,ineq,idum,ikpt,ikpt0,ikpt1
INTEGER :: lm1,lm2,itype1,itype2,ineq1,ineq2,n, n1,n2,ng
INTEGER :: ic,ic1,ic2,ic3,ic4,ic5,ic6,ic7,ic8
INTEGER :: igpt,igpt1,igpt2,igptp,igptp1,igptp2
INTEGER :: isym,isym1,isym2,igpt0
INTEGER :: maxlcut,ok
INTEGER :: m
INTEGER :: ikptmin,ikptmax,nkminmax
LOGICAL :: lsym
REAL :: rdum,rdum1,rdum2
REAL :: svol,qnorm,qnorm1,qnorm2,gnorm
REAL :: fcoulfac
REAL :: time1,time2
COMPLEX :: cdum,cdum1,cexp,csum
COMPLEX , PARAMETER :: img = (0d0,1d0)
! - local arrays -
INTEGER :: g(3)
INTEGER :: nbasm1(kpts%nkptf)
INTEGER , ALLOCATABLE :: pqnrm(:,:)
INTEGER :: rrot(3,3,sym%nsym),invrrot(3,3,sym%nsym)
INTEGER , ALLOCATABLE :: iarr(:),POINTER(:,:,:,:)!,pointer(:,:,:)
INTEGER :: igptmin(kpts%nkpt),igptmax(kpts%nkpt)
INTEGER , ALLOCATABLE :: nsym_gpt(:,:), sym_gpt(:,:,:)
INTEGER :: nsym1(kpts%nkpt+1), sym1(sym%nsym,kpts%nkpt+1)
LOGICAL :: calc_mt(kpts%nkpt)
REAL :: q(3),q1(3),q2(3)
REAL :: integrand(atoms%jmtd),primf1(atoms%jmtd), primf2(atoms%jmtd)
REAL :: mat(hybrid%maxindxm1*(hybrid%maxindxm1+1)/2)
REAL :: moment(hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype), moment2(hybrid%maxindxm1,atoms%ntype)
REAL :: sphbes(atoms%jmtd,0:hybrid%maxlcutm1)
REAL :: sphbesmoment1(atoms%jmtd,0:hybrid%maxlcutm1)
REAL :: rarr(0:hybrid%lexp+1),rarr1(0:hybrid%maxlcutm1)
REAL , ALLOCATABLE :: fac(:),sfac(:),facfac(:)
REAL , ALLOCATABLE :: gmat(:,:),qnrm(:)
REAL , ALLOCATABLE :: sphbesmoment(:,:,:)
REAL , ALLOCATABLE :: sphbes0(:,:,:)
REAL , ALLOCATABLE :: olap(:,:,:,:),integral(:,:,:,:)
REAL , ALLOCATABLE :: gridf(:,:)
COMPLEX :: cexp1(atoms%ntype),csumf(9)
COMPLEX :: structconst((2*hybrid%lexp+1)**2 ,atoms%nat,atoms%nat, kpts%nkpt) ! nw = 1
COMPLEX :: y((hybrid%lexp+1)**2),y1((hybrid%lexp+1)**2), y2((hybrid%lexp+1)**2)
COMPLEX :: dwgn(-hybrid%maxlcutm1:hybrid%maxlcutm1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1,sym%nsym)
COMPLEX , ALLOCATABLE :: smat(:,:)
COMPLEX , ALLOCATABLE :: coulmat(:,:)
COMPLEX , ALLOCATABLE :: carr2(:,:),carr2a(:,:), carr2b(:,:)
COMPLEX , ALLOCATABLE :: structconst1(:,:)
REAL , ALLOCATABLE :: coulomb_mt1(:,:,:,:,:)
!REAL , ALLOCATABLE :: coulomb(:,:) !At the moment we always calculate a complex coulomb matrix
REAL , ALLOCATABLE :: coulomb_mt2_r(:,:,:,:,:), coulomb_mt3_r(:,:,:,:)
REAL , ALLOCATABLE :: coulomb_mtir_r(:,:,:), coulombp_mtir_r(:,:)
COMPLEX , ALLOCATABLE :: coulomb(:,:)
COMPLEX , ALLOCATABLE :: coulomb_mt2_c(:,:,:,:,:), coulomb_mt3_c(:,:,:,:)
COMPLEX , ALLOCATABLE :: coulomb_mtir_c(:,:,:), coulombp_mtir_c(:,:)
INTEGER :: ishift,ishift1
INTEGER :: iatom,iatom1
INTEGER :: indx1,indx2,indx3,indx4
LOGICAL :: l_found,l_warn,l_warned, l_plot = .FALSE.!.true.!.false.
TYPE(t_mat) :: olapm,coulhlp
! - scalars -
LOGICAL , INTENT(IN) :: l_restart
! - local scalars -
INTEGER :: maxfac
! INTEGER :: irecl_coulomb
INTEGER :: inviop
INTEGER :: nqnrm,iqnrm,iqnrm1,iqnrm2, iqnrmstart,iqnrmstep
INTEGER :: itype,l ,ix,iy,iy0,i,j,lm,l1,l2,m1, m2,ineq,idum,ikpt,ikpt0,ikpt1
INTEGER :: lm1,lm2,itype1,itype2,ineq1,ineq2,n, n1,n2,ng
INTEGER :: ic,ic1,ic2,ic3,ic4,ic5,ic6,ic7,ic8
INTEGER :: igpt,igpt1,igpt2,igptp,igptp1,igptp2
INTEGER :: isym,isym1,isym2,igpt0
INTEGER :: maxlcut,ok
INTEGER :: nbasp,m
INTEGER :: ikptmin,ikptmax,nkminmax
#ifdef CPP_INVERSION
INTEGER , PARAMETER :: coul_size = 8
#else
INTEGER , PARAMETER :: coul_size = 16
#endif
LOGICAL :: lsym
REAL :: rdum,rdum1,rdum2
REAL :: svol,qnorm,qnorm1,qnorm2,gnorm
REAL :: fcoulfac
REAL :: time1,time2
COMPLEX :: cdum,cdum1,cexp,csum
COMPLEX , PARAMETER :: img = (0d0,1d0)
! - local arrays -
INTEGER :: g(3)
INTEGER :: nbasm1(kpts%nkptf)
INTEGER , ALLOCATABLE :: pqnrm(:,:)
INTEGER :: rrot(3,3,sym%nsym),invrrot(3,3,sym%nsym)
INTEGER , ALLOCATABLE :: iarr(:),pointer(:,:,:,:)!,pointer(:,:,:)
INTEGER :: igptmin(kpts%nkpt),igptmax(kpts%nkpt)
INTEGER , ALLOCATABLE :: nsym_gpt(:,:), sym_gpt(:,:,:)
INTEGER :: nsym1(kpts%nkpt+1), sym1(sym%nsym,kpts%nkpt+1)
LOGICAL :: calc_mt(kpts%nkpt)
REAL :: q(3),q1(3),q2(3)
REAL :: integrand(atoms%jmtd),primf1(atoms%jmtd),&
& primf2(atoms%jmtd)
REAL :: mat(hybrid%maxindxm1*(hybrid%maxindxm1+1)/2)
REAL :: moment(hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype),&
& moment2(hybrid%maxindxm1,atoms%ntype)
REAL :: sphbes(atoms%jmtd,0:hybrid%maxlcutm1)
REAL :: sphbesmoment1(atoms%jmtd,0:hybrid%maxlcutm1)
REAL :: rarr(0:hybrid%lexp+1),rarr1(0:hybrid%maxlcutm1)
REAL , ALLOCATABLE :: fac(:),sfac(:),facfac(:)
REAL , ALLOCATABLE :: gmat(:,:),qnrm(:)
REAL , ALLOCATABLE :: sphbesmoment(:,:,:)
REAL , ALLOCATABLE :: sphbes0(:,:,:)
REAL , ALLOCATABLE :: olap(:,:,:,:),integral(:,:,:,:)
REAL , ALLOCATABLE :: gridf(:,:)
COMPLEX :: cexp1(atoms%ntype),csumf(9)
COMPLEX :: structconst((2*hybrid%lexp+1)**2 ,atoms%nat,atoms%nat,&
& kpts%nkpt) ! nw = 1
COMPLEX :: y((hybrid%lexp+1)**2),y1((hybrid%lexp+1)**2),&
& y2((hybrid%lexp+1)**2)
COMPLEX :: dwgn(-hybrid%maxlcutm1:hybrid%maxlcutm1,&
& -hybrid%maxlcutm1:hybrid%maxlcutm1,&
& 0:hybrid%maxlcutm1,sym%nsym)
COMPLEX , ALLOCATABLE :: smat(:,:)
COMPLEX , ALLOCATABLE :: coulmat(:,:)
COMPLEX , ALLOCATABLE :: carr2(:,:),carr2a(:,:),&
& carr2b(:,:)
COMPLEX , ALLOCATABLE :: structconst1(:,:)
REAL , ALLOCATABLE :: coulomb_mt1(:,:,:,:,:)
#if ( defined(CPP_INVERSION) )
REAL , ALLOCATABLE :: coulomb(:,:),coulhlp(:,:),&
& coulhlp1(:,:)
REAL , ALLOCATABLE :: involapm(:,:),olapm(:,:),evec(:,:),&
& invevec(:,:)
REAL , ALLOCATABLE :: coulomb_mt2(:,:,:,:,:),&
& coulomb_mt3(:,:,:,:)
REAL , ALLOCATABLE :: coulomb_mtir(:,:,:),&
& coulombp_mtir(:,:)
#else
COMPLEX , ALLOCATABLE :: coulomb(:,:),coulhlp(:,:),&
& coulhlp1(:,:)
COMPLEX , ALLOCATABLE :: olapm(:,:)
COMPLEX , ALLOCATABLE :: coulomb_mt2(:,:,:,:,:),&
& coulomb_mt3(:,:,:,:)
COMPLEX , ALLOCATABLE :: coulomb_mtir(:,:,:),&
& coulombp_mtir(:,:)
#endif
INTEGER :: ishift,ishift1
INTEGER :: iatom,iatom1
INTEGER :: indx1,indx2,indx3,indx4
LOGICAL :: l_found,l_warn,l_warned,&
& l_plot = .false.!.true.!.false.
#ifdef CPP_NOSPMVEC
CHARACTER*8 , PARAMETER :: fname = 'coulomb'
#else
CHARACTER*8 , PARAMETER :: fname = 'coulomb1'
#endif
#ifndef CPP_MPI
INTEGER :: irecl_coulomb
#else
INCLUDE 'mpif.h'
INTEGER :: fh,ierr,reqd,length,irank2,isize2
INTEGER, ALLOCATABLE :: req(:)
INTEGER(KIND=MPI_OFFSET_KIND) :: filesize,offset0,offset,&
& irecl_coulomb
INTEGER(KIND=MPI_ADDRESS_KIND) :: iaddr
#ifdef CPP_INVERSION
INTEGER, PARAMETER :: datatype = MPI_REAL8
#else
INTEGER, PARAMETER :: datatype = MPI_COMPLEX16
#endif
CHARACTER(MPI_MAX_ERROR_STRING) :: errmsg
#endif
CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,gridf)
!hybrid%maxgptm = maxval(hybrid%ngptm(:))
nbasm1 = 0
nbasp = 0
DO itype=1,atoms%ntype
DO i=1,atoms%neq(itype)
DO l=0,hybrid%lcutm1(itype)
DO M=-l,l
DO j=1,hybrid%nindxm1(l,itype)
nbasp = nbasp + 1
END DO
END DO
END DO
END DO
END DO
hybrid%maxbasm1 = nbasp + hybrid%maxgptm
nbasm1 = nbasp + hybrid%ngptm(:)
svol = sqrt(cell%vol)
fcoulfac = 4*pi_const/cell%vol
maxlcut = maxval( atoms%lmax )
maxfac = max(2*maxlcut+hybrid%maxlcutm1+1,4*max(hybrid%maxlcutm1,hybrid%lexp)+1)
ALLOCATE ( fac( 0:maxfac),sfac( 0:maxfac),facfac(-1:maxfac) )
fac(0) = 1 !
sfac(0) = 1 ! Define:
facfac(-1:0) = 1 ! fac(i) = i!
DO i=1,maxfac ! sfac(i) = sqrt(i!)
fac(i) = fac(i-1)*i ! facfac(i) = (2i+1)!!
sfac(i) = sfac(i-1)*sqrt(i*1d0) !
facfac(i) = facfac(i-1)*(2*i+1) !
END DO
CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,gridf)