Commit 47573c88 authored by Gregor Michalicek's avatar Gregor Michalicek

Shorten kcrel branches in main/cdngen.F90

parent f4d1c6f1
......@@ -150,9 +150,17 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
CALL closeXMLElement('valenceDensity')
!---> changes
END IF ! mpi%irank = 0
IF (input%kcrel.EQ.0) THEN
results%seigc = 0.
results%seigc = 0.0
IF (mpi%irank.EQ.0) THEN
DO jspin = 1,input%jspins
DO n = 1,atoms%ntype
svdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
END DO
END DO
END IF
IF (input%kcrel.EQ.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
......@@ -170,129 +178,69 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
CALL mpi_bc_coreDen(mpi,atoms,input,dimension,rh,tec,qint)
#endif
END IF
END IF
IF (.NOT.sliceplot%slice) THEN
!add in core density
IF (mpi%irank.EQ.0) THEN
IF (input%kcrel.EQ.0) THEN
DO jspin = 1,input%jspins
CALL cored(input,jspin,atoms,outDen%mt,dimension,sphhar,vTot%mt(:,0,:,jspin), qint,rh,tec,seig)
rhTemp(:,:,jspin) = rh(:,:,jspin)
results%seigc = results%seigc + seig
END DO
ELSE
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 ((input%jspins.EQ.2).AND.(mpi%irank.EQ.0)) THEN
IF (mpi%irank.EQ.0) THEN
DO n = 1,atoms%ntype
svdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
stdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
END DO
END IF
!block 1 unnecessary for slicing: begin
IF (.NOT.sliceplot%slice) THEN
!add in core density
IF (mpi%irank.EQ.0) THEN
CALL cored(input,jspin,atoms,outDen%mt,dimension,sphhar,vTot%mt(:,0,:,jspin), qint,rh,tec,seig)
rhTemp(:,:,jspin) = rh(:,:,jspin)
results%seigc = results%seigc + seig
IF (input%jspins.EQ.2) THEN
DO n = 1,atoms%ntype
stdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
END DO
END IF
END IF ! mpi%irank = 0
!add core tail charge to outDen%pw
IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
IF (jspin.EQ.2) THEN
!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
IF (jspin .EQ. 2) THEN
DO ityp = 1,atoms%ntype
rhoint = (qint(ityp,1) + qint(ityp,2)) /cell%volint/input%jspins/2.0
momint = (qint(ityp,1) - qint(ityp,2)) /cell%volint/input%jspins/2.0
!rho_11
outDen%pw(1,1) = outDen%pw(1,1) + rhoint + momint*cos(noco%beta(ityp))
!rho_22
outDen%pw(1,2) = outDen%pw(1,2) + rhoint - momint*cos(noco%beta(ityp))
!real part rho_21
outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
!imaginary part rho_21
outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
END DO
END IF
DO ityp = 1,atoms%ntype
rhoint = (qint(ityp,1) + qint(ityp,2)) /cell%volint/input%jspins/2.0
momint = (qint(ityp,1) - qint(ityp,2)) /cell%volint/input%jspins/2.0
!rho_11
outDen%pw(1,1) = outDen%pw(1,1) + rhoint + momint*cos(noco%beta(ityp))
!rho_22
outDen%pw(1,2) = outDen%pw(1,2) + rhoint - momint*cos(noco%beta(ityp))
!real part rho_21
outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
!imaginary part rho_21
outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
END DO
!pk non-collinear (end)
ELSE IF (input%ctail) THEN
CALL cdnovlp(mpi,sphhar,stars,atoms,sym, dimension,vacuum,&
cell, input,oneD,l_st, jspin,rh(:,:,jspin),&
END IF
ELSE
IF (input%ctail) THEN
!+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
DO ityp = 1,atoms%ntype
outDen%pw(1,jspin) = outDen%pw(1,jspin) + qint(ityp,jspin)/input%jspins/cell%volint
END DO
END IF
!block 1 unnecessary for slicing: end
END IF ! .NOT.sliceplot%slice
END IF
END DO
END IF
END DO ! loop over spins
IF (input%kcrel.EQ.0) THEN
IF (mpi%irank.EQ.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)
ELSE ! input%kcrel.EQ.0
!relativistic core implementation : kcrel.eq.1
results%seigc = 0.
IF ((input%jspins.EQ.2).AND.(mpi%irank.EQ.0)) THEN
DO jspin = 1,input%jspins
DO n = 1,atoms%ntype
svdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
END DO
END DO
END IF
!block 1 unnecessary for slicing: begin
IF (.NOT.sliceplot%slice) THEN
!add in core density
IF (mpi%irank.EQ.0) THEN
CALL coredr(input,atoms,seig, outDen%mt,dimension,sphhar,vTot%mt(:,0,:,:),qint,rh)
results%seigc = results%seigc + seig
IF (input%jspins.EQ.2) THEN
DO jspin = 1,input%jspins
DO n = 1,atoms%ntype
stdn(n,jspin) = outDen%mt(1,0,n,jspin)/ (sfp_const*atoms%rmsh(1,n)*atoms%rmsh(1,n))
END DO
END DO
END IF
END IF
IF ((noco%l_noco).AND.(mpi%irank.EQ.0)) THEN
!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 ityp = 1,atoms%ntype
rhoint = (qint(ityp,1) + qint(ityp,2)) /cell%volint/input%jspins/2.0
momint = (qint(ityp,1) - qint(ityp,2)) /cell%volint/input%jspins/2.0
!rho_11
outDen%pw(1,1) = outDen%pw(1,1) + rhoint + momint*cos(noco%beta(ityp))
!rho_22
outDen%pw(1,2) = outDen%pw(1,2) + rhoint - momint*cos(noco%beta(ityp))
!real part rho_21
outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.5*momint *cos(noco%alph(ityp))*sin(noco%beta(ityp)),0.0)
!imaginary part rho_21
outDen%pw(1,3) = outDen%pw(1,3) + cmplx(0.0,-0.5*momint *sin(noco%alph(ityp))*sin(noco%beta(ityp)))
END DO
!pk non-collinear (end)
ELSE
DO jspin = 1,input%jspins
IF (input%ctail) THEN
!+gu hope this works as well
CALL cdnovlp(mpi,sphhar,stars,atoms,sym,dimension,vacuum,&
cell,input,oneD,l_st,jspin,rh(1,1,jspin),&
outDen%pw,outDen%vacxy,outDen%mt,outDen%vacz)
ELSE IF (mpi%irank.EQ.0) THEN
DO ityp = 1,atoms%ntype
outDen%pw(1,jspin) = outDen%pw(1,jspin) + qint(ityp,jspin)/input%jspins/cell%volint
END DO
END IF
END DO
END IF
!block 1 unnecessary for slicing: end
END IF ! .NOT.sliceplot%slice
! end relativistic core
END IF ! input%kcrel.EQ.0
END IF
IF (mpi%irank.EQ.0) THEN
!block 2 unnecessary for slicing: begin
......
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