Commit e4921cb5 authored by Daniel Wortmann's avatar Daniel Wortmann

Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

parents d3c897e3 01f64608
......@@ -188,7 +188,7 @@ CONTAINS
wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin) - usdus%us(l,n,jspin)*usdus%duds(l,n,jspin)
IF (apw(l,n)) THEN
fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin)
dfj(l) = 0.0d0
dfj(l) = 0.0
ELSE
dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk
fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk
......
......@@ -107,7 +107,7 @@ CONTAINS
wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin)-usdus%us(l,n,jspin)*usdus%duds(l,n,jspin) !Wronski determinante
IF (apw(l,n)) THEN
fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin)
dfj(l) = 0.0d0
dfj(l) = 0.0
ELSE
dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk
fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk
......
......@@ -16,8 +16,8 @@ if (${CMAKE_Fortran_COMPILER_ID} MATCHES "Intel")
else()
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -mkl -qopenmp -assume byterecl")
endif()
#set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -xHost -O2 -g")
set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -xMIC-AVX512 -O2")
set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -xHost -O2 -g")
#set(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE} -xMIC-AVX512 -O2")
if (${CMAKE_Fortran_COMPILER_VERSION} VERSION_LESS "19.0.0.0")
set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -C -traceback -O0 -g -ftrapuv -check uninit -check pointers -DCPP_DEBUG -warn all")
else()
......
......@@ -192,19 +192,19 @@ CONTAINS
n_col = indxl2g(i, nb, hmat%blacsdata%mycol, 0, hmat%blacsdata%npcol)
n_row = numroc (n_col, nb, hmat%blacsdata%myrow, 0, hmat%blacsdata%nprow)
IF (hmat%l_real) THEN
hmat%data_r(n_row+1:hmat%matsize1,i) = 0.d0
hmat%data_r(n_row+1:hmat%matsize1,i) = 0.0
ELSE
hmat%data_c(n_row+1:hmat%matsize1,i) = 0.d0
hmat%data_c(n_row+1:hmat%matsize1,i) = 0.0
ENDIF
ENDDO
! Use the ev_dist array to store the calculated values for the lower part.
IF (hmat%l_real) THEN
CALL pdtran(hmat%global_size1,hmat%global_size1,1.d0,hmat%data_r,1,1,&
hmat%blacsdata%blacs_desc,0.d0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc)
CALL pdtran(hmat%global_size1,hmat%global_size1,1.0,hmat%data_r,1,1,&
hmat%blacsdata%blacs_desc,0.0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc)
ELSE
CALL pztranc(hmat%global_size1,hmat%global_size2,cmplx(1.d0,0.d0),hmat%data_c,1,1,&
hmat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc)
CALL pztranc(hmat%global_size1,hmat%global_size2,cmplx(1.0,0.0),hmat%data_c,1,1,&
hmat%blacsdata%blacs_desc,cmplx(0.0,0.0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc)
ENDIF
! Copy the calculated values to the lower part of the H matrix
......@@ -259,11 +259,11 @@ CONTAINS
! 2b. tmp2 = ev_dist**T
IF (hmat%l_real) THEN
CALL pdtran(ev_dist%global_size1,ev_dist%global_size1,1.d0,ev_dist%data_r,1,1,&
ev_dist%blacsdata%blacs_desc,0.d0,tmp2_r,1,1,ev_dist%blacsdata%blacs_desc)
CALL pdtran(ev_dist%global_size1,ev_dist%global_size1,1.0,ev_dist%data_r,1,1,&
ev_dist%blacsdata%blacs_desc,0.0,tmp2_r,1,1,ev_dist%blacsdata%blacs_desc)
ELSE
CALL pztranc(ev_dist%global_size1,ev_dist%global_size1,cmplx(1.0,0.0),ev_dist%data_c,1,1,&
ev_dist%blacsdata%blacs_desc,cmplx(0.d0,0.d0),tmp2_c,1,1,ev_dist%blacsdata%blacs_desc)
ev_dist%blacsdata%blacs_desc,cmplx(0.0,0.0),tmp2_c,1,1,ev_dist%blacsdata%blacs_desc)
ENDIF
! 2c. A = U**-T * tmp2 ( = U**-T * Aorig * U**-1 )
......@@ -307,11 +307,11 @@ CONTAINS
! Set lower half from upper half
IF (hmat%l_real) THEN
CALL pdtran(hmat%global_size1,hmat%global_size1,1.d0,hmat%data_r,1,1,&
hmat%blacsdata%blacs_desc,0.d0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc)
CALL pdtran(hmat%global_size1,hmat%global_size1,1.0,hmat%data_r,1,1,&
hmat%blacsdata%blacs_desc,0.0,ev_dist%data_r,1,1,ev_dist%blacsdata%blacs_desc)
ELSE
CALL pztranc(hmat%global_size1,hmat%global_size1,cmplx(1.0,0.0),hmat%data_c,1,1,&
hmat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc)
hmat%blacsdata%blacs_desc,cmplx(0.0,0.0),ev_dist%data_c,1,1,ev_dist%blacsdata%blacs_desc)
ENDIF
......@@ -396,11 +396,11 @@ CONTAINS
! mult_ah_b_complex needs the transpose of U**-1, thus tmp2 = (U**-1)**T
IF (hmat%l_real) THEN
CALL pdtran(smat%global_size1,smat%global_size1,1.d0,smat%data_r,1,1,&
smat%blacsdata%blacs_desc,0.d0,tmp2_r,1,1,smat%blacsdata%blacs_desc)
CALL pdtran(smat%global_size1,smat%global_size1,1.0,smat%data_r,1,1,&
smat%blacsdata%blacs_desc,0.0,tmp2_r,1,1,smat%blacsdata%blacs_desc)
ELSE
CALL pztranc(smat%global_size1,smat%global_size1,cmplx(1.d0,0.d0),smat%data_c,1,1,&
smat%blacsdata%blacs_desc,cmplx(0.d0,0.d0),tmp2_c,1,1,smat%blacsdata%blacs_desc)
CALL pztranc(smat%global_size1,smat%global_size1,cmplx(1.0,0.0),smat%data_c,1,1,&
smat%blacsdata%blacs_desc,cmplx(0.0,0.0),tmp2_c,1,1,smat%blacsdata%blacs_desc)
ENDIF
#if defined (CPP_ELPA_201705003)
......
......@@ -119,8 +119,8 @@
ENDDO
IF ( l_mcd ) THEN ! create an energy grid for mcd-spectra
e_lo = 9.9d+9
e_up = -9.9d+9
e_lo = 9.9*10.0**9
e_up = -9.9*10.0**9
DO jspin = 1,input%jspins
DO n = 1,atoms%ntype
DO icore = 1 , mcd%ncore(n)
......
......@@ -12,11 +12,11 @@ module m_corespec
! PARAMETERS
complex, parameter :: cone = cmplx(1.d0,0.d0)
complex, parameter :: cimu = cmplx(0.d0,1.d0)
real, parameter :: alpha = 7.29735257d-3
real, parameter :: mec2 = 0.51099891d6
real, parameter :: ecoredeep = 0.5d0
complex, parameter :: cone = cmplx(1.0,0.0)
complex, parameter :: cimu = cmplx(0.0,1.0)
real, parameter :: alpha = 7.29735257e-3
real, parameter :: mec2 = 0.51099891e6
real, parameter :: ecoredeep = 0.5
integer, parameter :: edgel(11) = (/0,1,1,2,2,3,3,4,4,5,5/)
integer, parameter :: edgej(11) = (/1,1,3,3,5,5,7,7,9,9,11/)
......
This diff is collapsed.
......@@ -42,10 +42,10 @@ MODULE m_corespec_io
csi%edge = ""
csi%edgeidx = 0
csi%lx = -1
csi%ek0 = 0.d0
csi%emn = -2.d0
csi%emx = 20.d0
csi%ein = 0.1d0
csi%ek0 = 0.0
csi%emn = -2.0
csi%emx = 20.0
csi%ein = 0.1
csi%nqphi = 12
csi%nqr = 20
......@@ -182,12 +182,12 @@ MODULE m_corespec_io
&"maximum l: ","csi%lx = ",csi%lx,"will be used"
! csi%ek0
if(csi%ek0.le.0.d0) then
if(csi%ek0.le.0.0) then
write(*,csmsgs) trim(smeno),"found csi%ek0 <= 0.0 !"//csmsgerr ; stop
endif
csi%ek0 = csi%ek0*1000.d0 ! conversion from keV to eV
csv%gamma = 1.d0+csi%ek0/mec2
csv%beta = sqrt(1.d0-1.d0/(csv%gamma**2))
csi%ek0 = csi%ek0*1000.0 ! conversion from keV to eV
csv%gamma = 1.0+csi%ek0/mec2
csv%beta = sqrt(1.0-1.0/(csv%gamma**2))
if(csi%verb.eq.1) then
write(*,csmsgses) trim(smeno),&
&"kinetic energy of incoming electrons: ","csi%ek0 = ",csi%ek0,&
......@@ -204,8 +204,8 @@ MODULE m_corespec_io
if(csi%emn.gt.csi%emx) then
write(*,csmsgs) trim(smeno),"found csi%emn > csi%emx !"//csmsgerr ; stop
endif
if(csi%ein.le.0.d0) then
write(*,csmsgs) trim(smeno),"found csi%ein <= 0.d0 !"//csmsgerr ; stop
if(csi%ein.le.0.0) then
write(*,csmsgs) trim(smeno),"found csi%ein <= 0.0 !"//csmsgerr ; stop
endif
if(((csi%emx-csi%emn)/csi%ein)-int((csi%emx-csi%emn)/csi%ein).ne.0) then
write(*,csmsgs) trim(smeno),&
......@@ -216,7 +216,7 @@ MODULE m_corespec_io
csv%egrid = (/(csi%emn+csi%ein*dble(i), i = 0,csv%nex)/)
csv%nen = 0
!!$ do i = 0,csv%nex
!!$ if(csv%egrid(i).ge.0.d0) then
!!$ if(csv%egrid(i).ge.0.0) then
!!$ csv%nen = i
!!$ exit
!!$ endif
......@@ -240,9 +240,9 @@ MODULE m_corespec_io
&csv%nen,"will be used"
if(.not.allocated(csv%eedge)) allocate(csv%eedge(csv%nljc))
csv%eedge = 0.d0
csv%eedge = 0.0
if(.not.allocated(csv%occ)) allocate(csv%occ(csv%nljc))
csv%occ = 0.d0
csv%occ = 0.0
l_cs = .true.
......
......@@ -74,8 +74,8 @@ CONTAINS
INTEGER,INTENT(IN) :: eig_id
! Local Scalars
INTEGER jsp,nk,nred,ne_all,ne_found,nk_i
INTEGER ne,lh0
INTEGER jsp,nk,nred,ne_all,ne_found
INTEGER ne, nk_i
INTEGER isp,i,j,err
LOGICAL l_wu,l_file,l_real,l_zref
INTEGER :: solver=0
......@@ -85,7 +85,6 @@ CONTAINS
COMPLEX :: unfoldingBuffer(SIZE(results%unfolding_weights,1),kpts%nkpt,input%jspins) ! needed for unfolding bandstructure mpi case
INTEGER, PARAMETER :: lmaxb = 3
REAL, ALLOCATABLE :: bkpt(:)
REAL, ALLOCATABLE :: eig(:), eigBuffer(:,:,:)
......
......@@ -44,7 +44,7 @@ CONTAINS
gg = rk(k)*gb(l)
IF ( apw(l) ) THEN
fj(k,l,ispin) = 1.0*con1 * ff / us(l,ispin)
gj(k,l,ispin) = 0.0d0
gj(k,l,ispin) = 0.0
ELSE
IF (l_flag) THEN
DO jspin = 1, jspins
......@@ -167,7 +167,7 @@ CONTAINS
gg = lapw%rk(k,intspin)*gb(l)
IF ( apw(l) ) THEN
fj(k,l,ispin,intspin) = 1.0*con1 * ff / usdus%us(l,n,ispin)
gj(k,l,ispin,intspin) = 0.0d0
gj(k,l,ispin,intspin) = 0.0
ELSE
IF (noco%l_constr.or.l_socfirst) THEN
DO jspin = 1, input%jspins
......
......@@ -398,8 +398,8 @@ SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_s
!---> update overlap and l-diagonal hamiltonian matrix
kj_end = MIN(ki,lapw%nv(iintsp))
VecHelpS = 0.d0
VecHelpH = 0.d0
VecHelpS = 0.0
VecHelpH = 0.0
DO l = 0,atoms%lmax(n)
fjkiln = fj(ki,l,jintsp)
gjkiln = gj(ki,l,jintsp)
......
......@@ -88,7 +88,7 @@ CONTAINS
gg = lapw%rk(k,iintsp)*gb(l)
! IF ( apw(l) ) THEN
! fj(k,l,n,iintsp) = 1.0*con1 * ff / usdus%us(l,n,isp)
! gj(k,l,n,iintsp) = 0.0d0
! gj(k,l,n,iintsp) = 0.0
! ELSE
!---> in a spin-spiral calculation fj and gj are needed
!---> both interstitial spin directions at the same time
......
......@@ -139,7 +139,7 @@ CONTAINS
wronk = usdus%uds(l,n,jspin)*usdus%dus(l,n,jspin) - usdus%us(l,n,jspin)*usdus%duds(l,n,jspin)
IF (apw(l,n)) THEN
fj(l) = 1.0*const * fj(l)/usdus%us(l,n,jspin)
dfj(l) = 0.0d0
dfj(l) = 0.0
ELSE
dfj(l) = const* (usdus%dus(l,n,jspin)*fj(l)-df*usdus%us(l,n,jspin))/wronk
fj(l) = const* (df*usdus%uds(l,n,jspin)-fj(l)*usdus%duds(l,n,jspin))/wronk
......
......@@ -297,7 +297,7 @@ else
IF (i1.EQ.1) nn = 0
IF (i1.EQ.2) nn = nsz(1)
zhelp2(:,:) = 0.d0
zhelp2(:,:) = 0.0
DO j = 1,nsize
DO i = 1,nsz(jsp)
zhelp2(i,j) = CONJG(hso(i+nn,j))
......@@ -305,11 +305,11 @@ else
ENDDO ! j
if (l_real) THEN
CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.d0,0.d0),CMPLX(zmat(jsp)%data_r(:,:)),&
zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(1,1,jsp2),zmat(1)%matsize1)
CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd,CMPLX(1.0,0.0),CMPLX(zmat(jsp)%data_r(:,:)),&
zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(1,1,jsp2),zmat(1)%matsize1)
else
CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.d0,0.d0),zmat(jsp)%data_c(:,:),&
zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.d0,0.d0), zso(:,:,jsp2),zmat(1)%matsize1)
CALL CPP_BLAS_cgemm("N","N",zmat(1)%matsize1,2*dimension%neigd,dimension%neigd, CMPLX(1.0,0.0),zmat(jsp)%data_c(:,:),&
zmat(1)%matsize1, zhelp2,DIMENSION%neigd,CMPLX(0.0,0.0), zso(:,:,jsp2),zmat(1)%matsize1)
endif
ENDDO !isp
......
......@@ -287,10 +287,10 @@ c-------------------------------------------------------------------
w=0
nqvect=0
nshort=0
ReJq=0.d0
ImJq=0.d0
ReJq=0.0
ImJq=0.0
sqsin=(sin(thetaJ))**2
tpi = 2.d0 * pimach()
tpi = 2.0 * pimach()
limit=nmagn-1
IF (nmagn.gt.mtypes) limit=mtypes
......@@ -375,7 +375,7 @@ c...
DO ii=1,2
Dabsq(:)=ABS(q(:,nn,qcount)+((-1)**ii)*q(:,nnn,qcount))
IDabsq(:)=NINT(Dabsq(:))
divi(:)=ABS(Dabsq(:)/FLOAT(IDabsq(:))-1.d0)
divi(:)=ABS(Dabsq(:)/FLOAT(IDabsq(:))-1.0)
IF(((Dabsq(1).LT.tol).OR.(divi(1).LT.tol)).AND.
& ((Dabsq(2).LT.tol).OR.(divi(2).LT.tol)).AND.
& ((Dabsq(3).LT.tol).OR.(divi(3).LT.tol)))THEN
......@@ -397,12 +397,12 @@ c...
c... Now calculate Jq=Re(Jq)+i*Im(Jq)
mu=1
DO imt=1,mtypes
ReJq(mu,mu,qcount)=-2.d0*(seigv(mu,mu,qcount,1)
ReJq(mu,mu,qcount)=-2.0*(seigv(mu,mu,qcount,1)
& -seigv0(mu,mu,1))/(M(mu)*M(mu)*sqsin)
ImJq(mu,mu,qcount)=0.d0
ImJq(mu,mu,qcount)=0.0
DO remt=mu+1,mu+nmagtype(imt)-1
ReJq(remt,remt,qcount)=ReJq(mu,mu,qcount)
ImJq(remt,remt,qcount)=0.d0
ImJq(remt,remt,qcount)=0.0
ENDDO!remt
mu=mu+nmagtype(imt)
ENDDO !imt
......@@ -411,10 +411,10 @@ c... Now calculate Jq=Re(Jq)+i*Im(Jq)
DO nu=mu+1,nmagn
ReJq(mu,nu,qcount)=((seigv0(mu,nu,2)-
& seigv(mu,nu,qcount,1))/(M(mu)*M(nu)*sqsin))
& -(0.5d0*M(mu)*ReJq(mu,mu,qcount)/M(nu))-
& (0.5d0*M(nu)*ReJq(nu,nu,qcount)/M(mu))
& -(0.5*M(mu)*ReJq(mu,mu,qcount)/M(nu))-
& (0.5*M(nu)*ReJq(nu,nu,qcount)/M(mu))
IF(invs)THEN
ImJq(mu,nu,qcount)=0.d0
ImJq(mu,nu,qcount)=0.0
ELSE
ImJq(mu,nu,qcount)=((seigv(mu,nu,qcount,2)
& -seigv(mu,nu,qcount,1))/
......@@ -479,7 +479,7 @@ c ... for one magnetic atom per unit cell
qcount=nqpt-1
lwork=2*nshort
ALLOCATE (Cmat(qcount,nshort),DelE(qcount),work(lwork))
Cmat=0.d0
Cmat=0.0
IF (nshort.GE.nqpt)THEN
WRITE(*,*) ' Please supply the data for', nshort,
& 'q-points different from zero'
......@@ -492,7 +492,7 @@ c ... for one magnetic atom per unit cell
scp=(q(1,1,n)*R(1,atsh,nn)
& +q(2,1,n)*R(2,atsh,nn)
& +q(3,1,n)*R(3,atsh,nn))*tpi
Cmat(n,nn)=Cmat(n,nn)-1.d0+cos(scp)
Cmat(n,nn)=Cmat(n,nn)-1.0+cos(scp)
ENDDO
ENDDO
DelE(n)=ReJq(1,1,n)*2000 ! multiply by 2000 to get [mRy/muB**2]
......@@ -504,7 +504,7 @@ c ... for one magnetic atom per unit cell
& work,lwork,info)
c The routine dgels returns the solution, J(n), in the array DelE(n)
Tc=0.d0
Tc=0.0
DO n=1,nshort
Tc=Tc+nat(n)*DelE(n) !Mean-field Tc=1/3*(Sum_i(J_0,i))
WRITE(115,5005) n,lenR(n),DelE(n) ! J in units [mRy/muB**2]
......@@ -527,7 +527,7 @@ c... Perform the back-Fourier transform
wrJ=0
DO atsh=1,nat(nnn)
IF(atsh.gt.shmax) STOP 'jcoff2:increase shmax!'
J=0.d0
J=0.0
DO n=1,nqpt-1
DO nn=1,nop
IF(w(nn,n).EQ.1)THEN
......@@ -539,7 +539,7 @@ c... Perform the back-Fourier transform
ENDIF
ENDDO !nn
ENDDO !n (qpts)
J=(J/float(nqvect))*2000.d0 ! J in units [mRy/muB**2]
J=(J/float(nqvect))*2000.0 ! J in units [mRy/muB**2]
DO i=1,wrJ !A check for non-equivalent sub-shells
IF(ABS(J-Jw(i)).LE.(tol))GOTO 55
ENDDO
......@@ -565,7 +565,7 @@ c... In case of only one magnetic atom per unit cell, calculate the mean-field
mu=mu+nmagtype(imt)
ENDDO !imt
IF(nmagn.EQ.1) THEN
Tc=157.889*M(1)*M(1)*Tc/3.d0
Tc=157.889*M(1)*M(1)*Tc/3.0
WRITE(115,*) '# Tc(mean field)= ',Tc
ENDIF
5008 FORMAT(i4,i4,7(1x,f14.10))
......
MODULE m_hf_init
!
! preparations for HF and hybrid functional calculation
!
!
! preparations for HF and hybrid functional calculation
!
CONTAINS
SUBROUTINE hf_init(hybrid,kpts,atoms,input,DIMENSION,hybdat,l_real)
USE m_types
USE m_read_core
USE m_util
USE m_io_hybrid
IMPLICIT NONE
TYPE(t_hybrid),INTENT(INOUT) :: hybrid
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_input),INTENT(IN) :: input
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_hybdat),INTENT(OUT) :: hybdat
LOGICAL,INTENT(IN) :: l_real
INTEGER:: itype,ieq,l,m,i,nk,l1,l2,m1,m2,ok
SUBROUTINE hf_init(hybrid, kpts, atoms, input, DIMENSION, hybdat, l_real)
USE m_types
USE m_read_core
USE m_util
USE m_io_hybrid
IMPLICIT NONE
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_input), INTENT(IN) :: input
TYPE(t_dimension), INTENT(IN) :: DIMENSION
TYPE(t_hybdat), INTENT(OUT) :: hybdat
LOGICAL, INTENT(IN) :: l_real
!initialize hybdat%gridf for radial integration
CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
!Alloc variables
ALLOCATE(hybdat%lmaxc(atoms%ntype))
ALLOCATE(hybdat%bas1(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype))
ALLOCATE(hybdat%bas2(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype))
ALLOCATE(hybdat%bas1_MT(hybrid%maxindx,0:atoms%lmaxd,atoms%ntype))
ALLOCATE(hybdat%drbas1_MT(hybrid%maxindx,0:atoms%lmaxd,atoms%ntype))
!sym%tau = oneD%ods%tau
INTEGER:: itype, ieq, l, m, i, nk, l1, l2, m1, m2, ok
! preparations for core states
CALL core_init(dimension,input,atoms, hybdat%lmaxcd,hybdat%maxindxc)
ALLOCATE( hybdat%nindxc(0:hybdat%lmaxcd,atoms%ntype), stat = ok )
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation hybdat%nindxc'
ALLOCATE( hybdat%core1(atoms%jmtd,hybdat%maxindxc,0:hybdat%lmaxcd,atoms%ntype), stat=ok )
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation core1'
ALLOCATE( hybdat%core2(atoms%jmtd,hybdat%maxindxc,0:hybdat%lmaxcd,atoms%ntype), stat=ok )
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation core2'
ALLOCATE( hybdat%eig_c(hybdat%maxindxc,0:hybdat%lmaxcd,atoms%ntype), stat=ok )
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation hybdat%eig_c'
hybdat%nindxc = 0 ; hybdat%core1 = 0 ; hybdat%core2 = 0 ; hybdat%eig_c = 0
!initialize hybdat%gridf for radial integration
CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, hybdat%gridf)
! pre-calculate gaunt coefficients
!Alloc variables
ALLOCATE (hybdat%lmaxc(atoms%ntype))
ALLOCATE (hybdat%bas1(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype))
ALLOCATE (hybdat%bas2(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype))
ALLOCATE (hybdat%bas1_MT(hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype))
ALLOCATE (hybdat%drbas1_MT(hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype))
hybdat%maxfac = max(2*atoms%lmaxd+hybrid%maxlcutm1+1,2*hybdat%lmaxcd+2*atoms%lmaxd+1)
ALLOCATE ( hybdat%fac( 0:hybdat%maxfac),hybdat%sfac( 0:hybdat%maxfac),stat=ok)
IF( ok .ne. 0 ) STOP 'eigen_hf: failure allocation fac,hybdat%sfac'
hybdat%fac(0) = 1
hybdat%sfac(0) = 1
DO i=1,hybdat%maxfac
hybdat%fac(i) = hybdat%fac(i-1)*i ! hybdat%fac(i) = i!
hybdat%sfac(i) = hybdat%sfac(i-1)*sqrt(i*1d0) ! hybdat%sfac(i) = sqrt(i!)
END DO
!sym%tau = oneD%ods%tau
! preparations for core states
CALL core_init(dimension, input, atoms, hybdat%lmaxcd, hybdat%maxindxc)
ALLOCATE (hybdat%nindxc(0:hybdat%lmaxcd, atoms%ntype), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation hybdat%nindxc'
ALLOCATE (hybdat%core1(atoms%jmtd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation core1'
ALLOCATE (hybdat%core2(atoms%jmtd, hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation core2'
ALLOCATE (hybdat%eig_c(hybdat%maxindxc, 0:hybdat%lmaxcd, atoms%ntype), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation hybdat%eig_c'
hybdat%nindxc = 0; hybdat%core1 = 0; hybdat%core2 = 0; hybdat%eig_c = 0
ALLOCATE ( hybdat%gauntarr( 2, 0:atoms%lmaxd, 0:atoms%lmaxd, 0:hybrid%maxlcutm1, -atoms%lmaxd:atoms%lmaxd ,-hybrid%maxlcutm1:hybrid%maxlcutm1 ),stat=ok)
IF( ok .ne. 0 ) STOP 'eigen: failure allocation hybdat%gauntarr'
hybdat%gauntarr = 0
DO l2 = 0,atoms%lmaxd
DO l1 = 0,atoms%lmaxd
DO l = abs(l1-l2),min(l1+l2,hybrid%maxlcutm1)
DO m = -l,l
DO m1 = -l1,l1
m2 = m1 + m ! Gaunt condition -m1+m2-m = 0
IF(abs(m2).le.l2) hybdat%gauntarr(1,l1,l2,l,m1,m) = gaunt(l1,l2,l,m1,m2,m,hybdat%maxfac,hybdat%fac,hybdat%sfac)
m2 = m1 - m ! switch role of l2-index
IF(abs(m2).le.l2) hybdat%gauntarr(2,l1,l2,l,m1,m) = gaunt(l2,l1,l,m2,m1,m,hybdat%maxfac,hybdat%fac,hybdat%sfac)
END DO
END DO
END DO
END DO
END DO
! pre-calculate gaunt coefficients
!skip_kpt = .false.
hybdat%maxfac = max(2*atoms%lmaxd + hybrid%maxlcutm1 + 1, 2*hybdat%lmaxcd + 2*atoms%lmaxd + 1)
ALLOCATE (hybdat%fac(0:hybdat%maxfac), hybdat%sfac(0:hybdat%maxfac), stat=ok)
IF (ok /= 0) STOP 'eigen_hf: failure allocation fac,hybdat%sfac'
hybdat%fac(0) = 1
hybdat%sfac(0) = 1
DO i = 1, hybdat%maxfac
hybdat%fac(i) = hybdat%fac(i - 1)*i ! hybdat%fac(i) = i!
hybdat%sfac(i) = hybdat%sfac(i - 1)*sqrt(i*1.0) ! hybdat%sfac(i) = sqrt(i!)
END DO
END SUBROUTINE hf_init
ALLOCATE (hybdat%gauntarr(2, 0:atoms%lmaxd, 0:atoms%lmaxd, 0:hybrid%maxlcutm1, -atoms%lmaxd:atoms%lmaxd, -hybrid%maxlcutm1:hybrid%maxlcutm1), stat=ok)
IF (ok /= 0) STOP 'eigen: failure allocation hybdat%gauntarr'
hybdat%gauntarr = 0
DO l2 = 0, atoms%lmaxd
DO l1 = 0, atoms%lmaxd
DO l = abs(l1 - l2), min(l1 + l2, hybrid%maxlcutm1)
DO m = -l, l
DO m1 = -l1, l1
m2 = m1 + m ! Gaunt condition -m1+m2-m = 0
IF (abs(m2) <= l2) hybdat%gauntarr(1, l1, l2, l, m1, m) = gaunt(l1, l2, l, m1, m2, m, hybdat%maxfac, hybdat%fac, hybdat%sfac)
m2 = m1 - m ! switch role of l2-index
IF (abs(m2) <= l2) hybdat%gauntarr(2, l1, l2, l, m1, m) = gaunt(l2, l1, l, m2, m1, m, hybdat%maxfac, hybdat%fac, hybdat%sfac)
END DO
END DO
END DO
END DO
END DO
!skip_kpt = .false.
END SUBROUTINE hf_init
END MODULE m_hf_init
......@@ -40,95 +40,95 @@ MODULE m_add_vnonlocal
! c
! M.Betzinger (09/07) c
! c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c
CONTAINS
CONTAINS
SUBROUTINE add_vnonlocal(nk,lapw,atoms,hybrid,dimension,kpts,jsp,results,xcpot,noco,hmat)
SUBROUTINE add_vnonlocal(nk, lapw, atoms, hybrid, dimension, kpts, jsp, results, xcpot, noco, hmat)
USE m_symm_hf, ONLY: symm_hf
USE m_util, ONLY: intgrf,intgrf_init
USE m_symm_hf, ONLY: symm_hf
USE m_util, ONLY: intgrf, intgrf_init
USE m_exchange_valence_hf
USE m_exchange_core
USE m_symmetrizeh
USE m_wrapper
USE m_hsefunctional, ONLY: exchange_vccvHSE,exchange_ccccHSE
USE m_hsefunctional, ONLY: exchange_vccvHSE, exchange_ccccHSE
USE m_types
USE m_io_hybrid
IMPLICIT NONE
TYPE(t_results), INTENT(INOUT) :: results
CLASS(t_xcpot), INTENT(IN) :: xcpot
TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_lapw), INTENT(IN) :: lapw
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_mat), INTENT(INOUT) :: hmat
INTEGER, INTENT(IN) :: jsp
INTEGER, INTENT(IN) :: nk
TYPE(t_results), INTENT(INOUT) :: results
CLASS(t_xcpot), INTENT(IN) :: xcpot
TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_hybrid), INTENT(INOUT) :: hybrid
TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_lapw), INTENT(IN) :: lapw
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_mat), INTENT(INOUT) :: hmat
INTEGER, INTENT(IN) :: jsp
INTEGER, INTENT(IN) :: nk
! local scalars
INTEGER :: n,nn,iband,nbasfcn
INTEGER :: n, nn, iband, nbasfcn
REAL :: a_ex
TYPE(t_mat) :: olap,tmp,v_x,z
COMPLEX :: exch(dimension%neigd,dimension%neigd)
TYPE(t_mat) :: olap, tmp, v_x, z
COMPLEX :: exch(dimension%neigd, dimension%neigd)
! initialize weighting factor for HF exchange part
a_ex=xcpot%get_exchange_weight()
a_ex = xcpot%get_exchange_weight()
nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
CALL v_x%init(hmat%l_real,nbasfcn,nbasfcn)
nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco)
CALL v_x%init(hmat%l_real, nbasfcn, nbasfcn)
CALL read_v_x(v_x,kpts%nkpt*(jsp-1)+nk)
CALL read_v_x(v_x, kpts%nkpt*(jsp - 1) + nk)
! add non-local x-potential to the hamiltonian hmat
DO n = 1, v_x%matsize1
DO nn = 1, n
DO nn = 1, n
IF (hmat%l_real) THEN
hmat%data_r(nn,n) = hmat%data_r(nn,n) - a_ex*v_x%data_r(nn,n)
hmat%data_r(nn, n) = hmat%data_r(nn, n) - a_ex*v_x%data_r(nn, n)
ELSE
hmat%data_c(nn,n) = hmat%data_c(nn,n) - a_ex*v_x%