Commit cbac3b75 authored by Matthias Redies's avatar Matthias Redies

put changes on branch

parent 55719b61
......@@ -9,7 +9,7 @@ MODULE m_cdncore
CONTAINS
SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
stars,cell,sphhar,atoms,vTot,outDen,moments,results)
stars,cell,sphhar,atoms,vTot,outDen,moments,results, EnergyDen)
USE m_constants
USE m_cdn_io
......@@ -27,21 +27,22 @@ SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
IMPLICIT NONE
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_input), INTENT(IN) :: input
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_potden), INTENT(IN) :: vTot
TYPE(t_potden), INTENT(INOUT) :: outDen
TYPE(t_moments), INTENT(INOUT) :: moments
TYPE(t_results), INTENT(INOUT) :: results
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_input), INTENT(IN) :: input
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_potden), INTENT(IN) :: vTot
TYPE(t_potden), INTENT(INOUT) :: outDen
TYPE(t_moments), INTENT(INOUT) :: moments
TYPE(t_results), INTENT(INOUT) :: results
TYPE(t_potden), INTENT(INOUT), OPTIONAL :: EnergyDen
INTEGER :: jspin, n, iType
REAL :: seig, rhoint, momint
......@@ -53,7 +54,7 @@ SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
REAL :: rhTemp(dimension%msh,atoms%ntype,dimension%jspd)
results%seigc = 0.0
IF (mpi%irank.EQ.0) THEN
IF (mpi%irank==0) THEN
DO jspin = 1,input%jspins
DO n = 1,atoms%ntype
moments%svdn(n,jspin) = outDen%mt(1,0,n,jspin) / (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
......@@ -61,10 +62,10 @@ SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
END DO
END IF
IF (input%kcrel.EQ.0) THEN
IF (input%kcrel==0) THEN
! Generate input file ecore for subsequent GW calculation
! 11.2.2004 Arno Schindlmayr
IF ((input%gw.eq.1 .or. input%gw.eq.3).AND.(mpi%irank.EQ.0)) THEN
IF ((input%gw==1 .or. input%gw==3).AND.(mpi%irank==0)) THEN
OPEN (15,file='ecore',status='unknown', action='write',form='unformatted')
END IF
......@@ -72,7 +73,7 @@ SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
tec = 0.0
qint = 0.0
IF (input%frcor) THEN
IF (mpi%irank.EQ.0) THEN
IF (mpi%irank==0) THEN
CALL readCoreDensity(input,atoms,dimension,rh,tec,qint)
END IF
#ifdef CPP_MPI
......@@ -82,51 +83,55 @@ SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
END IF
!add in core density
IF (mpi%irank.EQ.0) THEN
IF (input%kcrel.EQ.0) THEN
IF (mpi%irank==0) THEN
IF (input%kcrel==0) THEN
DO jspin = 1,input%jspins
CALL cored(input,jspin,atoms,outDen%mt,dimension,sphhar,vTot%mt(:,0,:,jspin), qint,rh,tec,seig)
CALL cored(input,jspin,atoms,outDen%mt,dimension,sphhar,vTot%mt(:,0,:,jspin), qint,rh ,tec,seig, EnergyDen%mt)
rhTemp(:,:,jspin) = rh(:,:,jspin)
results%seigc = results%seigc + seig
END DO
ELSE
IF(PRESENT(EnergyDen)) call juDFT_error("Energyden not implemented for relativistic")
CALL coredr(input,atoms,seig, outDen%mt,dimension,sphhar,vTot%mt(:,0,:,:),qint,rh)
results%seigc = results%seigc + seig
END IF
END IF
DO jspin = 1,input%jspins
IF (mpi%irank.EQ.0) THEN
IF (mpi%irank==0) THEN
DO n = 1,atoms%ntype
moments%stdn(n,jspin) = outDen%mt(1,0,n,jspin) / (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
END DO
END IF
IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
IF (jspin.EQ.2) THEN
IF ((noco%l_noco).AND.(mpi%irank==0)) THEN
IF (jspin==2) THEN
IF(PRESENT(EnergyDen)) call juDFT_error("Energyden not implemented for noco")
!pk non-collinear (start)
!add the coretail-charge to the constant interstitial
!charge (star 0), taking into account the direction of
!magnetisation of this atom
DO iType = 1,atoms%ntype
rhoint = (qint(iType,1) + qint(iType,2)) /cell%volint/input%jspins/2.0
momint = (qint(iType,1) - qint(iType,2)) /cell%volint/input%jspins/2.0
rhoint = (qint(iType,1) + qint(iType,2)) /(cell%volint * input%jspins * 2.0)
momint = (qint(iType,1) - qint(iType,2)) /(cell%volint * input%jspins * 2.0)
!rho_11
outDen%pw(1,1) = outDen%pw(1,1) + rhoint + momint*cos(noco%beta(iType))
!rho_22
outDen%pw(1,2) = outDen%pw(1,2) + rhoint - momint*cos(noco%beta(iType))
!real part rho_21
outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(iType))*sin(noco%beta(iType)),0.0)
outDen%pw(1,3) = outDen%pw(1,3) + cmplx( 0.5*momint *cos(noco%alph(iType))*sin(noco%beta(iType)),&
!imaginary part rho_21
outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(iType))*sin(noco%beta(iType)))
-0.5*momint *sin(noco%alph(iType))*sin(noco%beta(iType)))
END DO
!pk non-collinear (end)
END IF
ELSE
IF (input%ctail) THEN
IF(PRESENT(EnergyDen)) call juDFT_error("Energyden not implemented for ctail")
!+gu hope this works as well
CALL cdnovlp(mpi,sphhar,stars,atoms,sym,dimension,vacuum,&
cell,input,oneD,l_st,jspin,rh(:,:,jspin),&
outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz)
ELSE IF (mpi%irank.EQ.0) THEN
ELSE IF (mpi%irank==0) THEN
DO iType = 1,atoms%ntype
outDen%pw(1,jspin) = outDen%pw(1,jspin) + qint(iType,jspin) / (input%jspins * cell%volint)
END DO
......@@ -134,11 +139,11 @@ SUBROUTINE cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
END IF
END DO
IF (input%kcrel.EQ.0) THEN
IF (mpi%irank.EQ.0) THEN
IF (input%kcrel==0) THEN
IF (mpi%irank==0) THEN
CALL writeCoreDensity(input,atoms,dimension,rhTemp,tec,qint)
END IF
IF ((input%gw.eq.1 .or. input%gw.eq.3).AND.(mpi%irank.EQ.0)) CLOSE(15)
IF ((input%gw==1 .or. input%gw==3).AND.(mpi%irank==0)) CLOSE(15)
END IF
END SUBROUTINE cdncore
......
MODULE m_cored
CONTAINS
SUBROUTINE cored(&
& input,jspin,atoms,&
& rho,DIMENSION,&
& sphhar,&
& vr,&
& qint,rhc,tec,seig)
SUBROUTINE cored(input, jspin, atoms, rho, DIMENSION, sphhar, vr, qint, rhc, tec, seig, EnergyDen)
! *******************************************************
! ***** set up the core densities for compounds. *****
! ***** d.d.koelling *****
......@@ -29,20 +23,22 @@ CONTAINS
REAL, INTENT (OUT) :: seig
! ..
! .. Array Arguments ..
REAL, INTENT(IN) :: vr(atoms%jmtd,atoms%ntype)
REAL, INTENT(INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: rhc(DIMENSION%msh,atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: qint(atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: tec(atoms%ntype,DIMENSION%jspd)
REAL, INTENT(IN) :: vr(atoms%jmtd,atoms%ntype)
REAL, INTENT(INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: rhc(DIMENSION%msh,atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: qint(atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT) :: tec(atoms%ntype,DIMENSION%jspd)
REAL, INTENT(INOUT), OPTIONAL :: EnergyDen(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd)
! ..
! .. Local Scalars ..
REAL e,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2
REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight
REAL eig,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2
REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight, aux_weight
INTEGER i,j,jatom,korb,n,ncmsh,nm,nm1,nst ,l,ierr
! ..
! .. Local Arrays ..
REAL rhcs(DIMENSION%msh),rhoc(DIMENSION%msh),rhoss(DIMENSION%msh),vrd(DIMENSION%msh),f(0:3)
REAL rhcs_aux(DIMENSION%msh), rhoss_aux(DIMENSION%msh) !> quantities for energy density calculations
REAL occ(DIMENSION%nstd),a(DIMENSION%msh),b(DIMENSION%msh),ain(DIMENSION%msh),ahelp(DIMENSION%msh)
REAL occ_h(DIMENSION%nstd,2)
INTEGER kappa(DIMENSION%nstd),nprnc(DIMENSION%nstd)
......@@ -99,7 +95,7 @@ CONTAINS
WRITE (6,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
WRITE (16,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
DO j = 1,atoms%jri(jatom)
rhoss(j) = 0.
rhoss(j) = 0.0
vrd(j) = vr(j,jatom)
ENDDO
!
......@@ -137,23 +133,35 @@ CONTAINS
IF (occ(korb) /= 0.0) THEN
fn = nprnc(korb)
fj = iabs(kappa(korb)) - .5e0
weight = 2*fj + 1.e0
IF (bmu > 99.) weight = occ(korb)
fl = fj + (.5e0)*isign(1,kappa(korb))
e = -2* (z/ (fn+fl))**2
CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, e, a,b,ierr)
stateEnergies(korb) = e
WRITE (6,FMT=8010) fn,fl,fj,e,weight
WRITE (16,FMT=8010) fn,fl,fj,e,weight
eig = -2* (z/ (fn+fl))**2
CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, eig, a,b,ierr)
stateEnergies(korb) = eig
WRITE (6,FMT=8010) fn,fl,fj,eig,weight
WRITE (16,FMT=8010) fn,fl,fj,eig,weight
IF (ierr/=0) CALL juDFT_error("error in core-level routine" ,calledby ="cored")
IF (input%gw==1 .OR. input%gw==3) WRITE (15) NINT(fl),weight,e,&
IF (input%gw==1 .OR. input%gw==3) WRITE (15) NINT(fl),weight,eig,&
a(1:atoms%jri(jatom)),b(1:atoms%jri(jatom))
sume = sume + weight*e/input%jspins
sume = sume + weight*eig/input%jspins
DO j = 1,ncmsh
rhcs(j) = weight* (a(j)**2+b(j)**2)
rhcs(j) = weight* (a(j)**2+b(j)**2)
rhoss(j) = rhoss(j) + rhcs(j)
ENDDO
IF(present(EnergyDen)) THEN
rhoss_aux = rhoss
DO j = 1,ncmsh
! for energy density we want to multiply the weights
! with the eigenenergies
rhoss_aux(j) = rhoss_aux(j) + (rhcs(j) * eig)
ENDDO
ENDIF
ENDIF
ENDDO
......@@ -165,6 +173,14 @@ CONTAINS
rho(j,0,jatom,jspin) = rho(j,0,jatom,jspin) + rhoc(j)/sfp_const
ENDDO
IF(present(EnergyDen)) then
DO j = 1,nm
rhoc(j) = rhoss(j)/input%jspins
EnergyDen(j,0,jatom,jspin) = EnergyDen(j,0,jatom,jspin) &
+ rhoss_aux(j) /(input%jspins * sfp_const)
ENDDO
ENDIF
rhc(1:ncmsh,jatom,jspin) = rhoss(1:ncmsh) / input%jspins
rhc(ncmsh+1:DIMENSION%msh,jatom,jspin) = 0.0
......
This diff is collapsed.
......@@ -158,15 +158,15 @@ CONTAINS
veff = vTot
IF(xcpot%is_hybrid().AND.hybrid%l_subvxc) THEN
DO ispin = 1, input%jspins
DO ispin = 1, input%jspins
CALL convol(stars,vx%pw_w(:,ispin),vx%pw(:,ispin),stars%ufft)
END DO
veff%pw = vTot%pw - xcpot%get_exchange_weight() * vx%pw
veff%pw_w = vTot%pw_w - xcpot%get_exchange_weight() * vx%pw_w
veff%mt = vTot%mt - xcpot%get_exchange_weight() * vx%mt
exc%pw = exc%pw - xcpot%get_exchange_weight() * exc%pw
exc%pw_w = exc%pw_w - xcpot%get_exchange_weight() * exc%pw_w
exc%mt = exc%mt - xcpot%get_exchange_weight() * exc%mt
END DO
veff%pw = vTot%pw - xcpot%get_exchange_weight() * vx%pw
veff%pw_w = vTot%pw_w - xcpot%get_exchange_weight() * vx%pw_w
veff%mt = vTot%mt - xcpot%get_exchange_weight() * vx%mt
exc%pw = exc%pw - xcpot%get_exchange_weight() * exc%pw
exc%pw_w = exc%pw_w - xcpot%get_exchange_weight() * exc%pw_w
exc%mt = exc%mt - xcpot%get_exchange_weight() * exc%mt
END IF
results%te_veff = 0.0
......
set(fleur_F90 ${fleur_F90}
xc-pot/libxc_postprocess_gga.f90
xc-pot/metagga.F90
)
set(fleur_F77 ${fleur_F77}
......
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment