Skip to content

GitLab

fleur
fleur
Commit d2438580 authored by Robin Hilgers's avatar Robin Hilgers
Browse files
Options

Merge branch 'revert-9ade4527' into 'develop' 

iplot takes integers as input again

See merge request fleur/fleur!30
parents 8ca35011 3ffc2456
Hide whitespace changes
Inline Side-by-side
Showing with 567 additions and 66 deletions
init/old_inp/fleur_init_old.F90
View file @ d2438580
......@@ -184,7 +184,7 @@ CONTAINS
INQUIRE(file="cdn1",exist=l_opti)
IF (noco%l_noco) INQUIRE(file="rhomat_inp",exist=l_opti)
l_opti=.NOT.l_opti
IF ((sliceplot%iplot).OR.(input%strho).OR.(input%swsp).OR.&
IF ((sliceplot%iplot.NE.0).OR.(input%strho).OR.(input%swsp).OR.&
& (input%lflip).OR.(input%l_bmt)) l_opti = .TRUE.
!
......
init/old_inp/setup.f90
View file @ d2438580
......@@ -173,7 +173,7 @@
!
ENDIF ! (mpi%irank == 0)
CALL stepf(sym,stars,atoms,oneD, input,cell, vacuum,mpi)
IF (.NOT.sliceplot%iplot) THEN
IF (sliceplot%iplot.EQ.0) THEN
IF ( mpi%irank == 0 ) THEN
CALL convn(DIMENSION,atoms,stars)
......
init/postprocessInput.F90
View file @ d2438580
......@@ -523,7 +523,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
INQUIRE(file="cdn1",exist=l_opti)
if (noco%l_noco) INQUIRE(file="rhomat_inp",exist=l_opti)
l_opti=.not.l_opti
IF ((sliceplot%iplot).OR.(input%strho).OR.(input%swsp).OR.&
IF ((sliceplot%iplot.NE.0).OR.(input%strho).OR.(input%swsp).OR.&
(input%lflip).OR.(input%l_bmt)) l_opti = .TRUE.
IF (.NOT.l_opti) THEN
......@@ -546,7 +546,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts
CALL timestart("stepf")
CALL stepf(sym,stars,atoms,oneD,input,cell,vacuum,mpi)
CALL timestop("stepf")
IF (.NOT.sliceplot%iplot) THEN
IF (sliceplot%iplot.EQ.0) THEN
IF (mpi%irank.EQ.0) THEN
CALL convn(DIMENSION,atoms,stars)
CALL e_field(atoms,DIMENSION,stars,sym,vacuum,cell,input,field%efield)
......
inpgen/set_inp.f90
View file @ d2438580
......@@ -140,7 +140,7 @@
input%gauss= .false. ; input%tria = .false.
sliceplot%slice= .false. ; input%swsp = .false.
input%lflip= .false. ; banddos%vacdos= .false. ; input%integ = .false.
sliceplot%iplot= .false. ; input%score = .false. ; sliceplot%plpot = .false.
sliceplot%iplot= 0 ; input%score = .false. ; sliceplot%plpot = .false.
input%pallst = .false. ; obsolete%lwb = .false. ; vacuum%starcoeff = .false.
input%strho = .false. ; input%l_f = .false. ; atoms%l_geo(:) = .true.
noco%l_noco = noco%l_ss ; input%jspins = 1
......
io/r_inpXML.F90
View file @ d2438580
......@@ -1924,7 +1924,7 @@ CONTAINS
input%vchk = .FALSE.
input%cdinf = .FALSE.
sliceplot%iplot = .FALSE.
sliceplot%iplot = 0
input%score = .FALSE.
sliceplot%plpot = .FALSE.
......@@ -1962,7 +1962,7 @@ CONTAINS
numberNodes = xmlGetNumberOfNodes(xPathA)
IF (numberNodes.EQ.1) THEN
sliceplot%iplot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@iplot'))
sliceplot%iplot = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@iplot'))
input%score = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@score'))
sliceplot%plpot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plplot'))
END IF
......
io/rw_inp.f90
View file @ d2438580
......@@ -619,15 +619,17 @@
END IF
!
band = .false.
READ (UNIT=5,FMT=8050,END=992,ERR=992) sliceplot%iplot,input%score,sliceplot%plpot,band
WRITE (6,9240) sliceplot%iplot,input%score,sliceplot%plpot,band
READ (UNIT=5,FMT=8050,END=992,ERR=992) ldum,input%score,sliceplot%plpot,band
WRITE (6,9240) ldum,input%score,sliceplot%plpot,band
sliceplot%iplot=MERGE(1,0,ldum)
IF (band) THEN
banddos%dos=.true. ; banddos%ndir = -4
ENDIF
GOTO 993
992 BACKSPACE(5)
READ (UNIT=5,FMT=8050,END=99,ERR=99) sliceplot%iplot,input%score,sliceplot%plpot
WRITE (6,9240) sliceplot%iplot,input%score,sliceplot%plpot,band
READ (UNIT=5,FMT=8050,END=99,ERR=99) ldum,input%score,sliceplot%plpot
WRITE (6,9240) ldum,input%score,sliceplot%plpot,band
sliceplot%iplot=MERGE(1,0,ldum)
!
993 READ (UNIT=5,FMT='(i3,2f10.6,6x,i3,8x,l1)',END=99,ERR=99)&
& sliceplot%kk,sliceplot%e1s,sliceplot%e2s,sliceplot%nnne,input%pallst
......@@ -930,7 +932,7 @@
WRITE (5,*)
END IF
band = .false.
WRITE (5,9240) sliceplot%iplot,input%score,sliceplot%plpot,band
WRITE (5,9240) ldum,input%score,sliceplot%plpot,band
9240 FORMAT ('iplot=',l1,',score=',l1,',plpot=',l1,',band=',l1)
WRITE (5,9250) sliceplot%kk,sliceplot%e1s,sliceplot%e2s,sliceplot%nnne,input%pallst
9250 FORMAT (i3,2f10.6,',nnne=',i3,',pallst=',l1)
......
io/w_inpXML.f90
View file @ d2438580
......@@ -679,8 +679,8 @@ SUBROUTINE w_inpXML(&
395 FORMAT(' <unfoldingBand unfoldBand="',l1,'" supercellX="',i0,'" supercellY="',i0,'" supercellZ="',i0,'"/>')
WRITE (fileNum,395) banddos%unfoldband, banddos%s_cell_x, banddos%s_cell_y, banddos%s_cell_z
! <plotting iplot="F" score="F" plplot="F"/>
400 FORMAT(' <plotting iplot="',l1,'" score="',l1,'" plplot="',l1,'"/>')
! <plotting iplot="0" score="F" plplot="F"/>
400 FORMAT(' <plotting iplot="',i0,'" score="',l1,'" plplot="',l1,'"/>')
WRITE (fileNum,400) sliceplot%iplot,input%score,sliceplot%plpot
! <chargeDensitySlicing numkpt="0" minEigenval="0.000000" maxEigenval="0.000000" nnne="0" pallst="F"/>
......
io/xml/FleurInputSchema.xsd
View file @ d2438580
......@@ -580,7 +580,7 @@
</xsd:complexType>
<xsd:complexType name="PlottingType">
<xsd:attribute default="F" name="iplot" type="FleurBool" use="optional"/>
<xsd:attribute default="0" name="iplot" type="xsd:nonNegativeInteger" use="optional"/>
<xsd:attribute default="F" name="score" type="FleurBool" use="optional"/>
<xsd:attribute default="F" name="plplot" type="FleurBool" use="optional"/>
</xsd:complexType>
......
io/xml/FleurOutputSchema.xsd
View file @ d2438580
......@@ -820,7 +820,7 @@
</xsd:complexType>
<xsd:complexType name="PlottingType">
<xsd:attribute default="F" name="iplot" type="FleurBool" use="optional"/>
<xsd:attribute default="F" name="iplot" type="xsd:nonNegativeInteger" use="optional"/>
<xsd:attribute default="F" name="score" type="FleurBool" use="optional"/>
<xsd:attribute default="F" name="plplot" type="FleurBool" use="optional"/>
</xsd:complexType>
......
main/fleur_init.F90
View file @ d2438580
......@@ -156,7 +156,7 @@
kpts%ntet = 1
kpts%numSpecialPoints = 1
sliceplot%iplot=.FALSE.
sliceplot%iplot=0
sliceplot%kk = 0
sliceplot%e1s = 0.0
sliceplot%e2s = 0.0
......
main/optional.F90
View file @ d2438580
......@@ -94,7 +94,7 @@ CONTAINS
! ..
it = 1
IF (sliceplot%iplot .AND. (mpi%irank==0) ) THEN
IF ((sliceplot%iplot.NE.0 ).AND. (mpi%irank==0) ) THEN
IF (noco%l_noco) THEN
CALL pldngen(mpi,sym,stars,atoms,sphhar,vacuum,&
cell,input,noco,oneD,sliceplot)
......@@ -104,9 +104,9 @@ CONTAINS
IF (mpi%irank == 0) THEN
IF (sliceplot%plpot) input%score = .FALSE.
IF (sliceplot%iplot) THEN
IF (sliceplot%iplot.NE.0) THEN
CALL timestart("Plotting")
IF (input%strho) CALL juDFT_error("strho = T and iplot=T",calledby = "optional")
IF (input%strho) CALL juDFT_error("strho = T and iplot=/=0",calledby = "optional")
CALL plotdop(oneD,dimension,stars,vacuum,sphhar,atoms,&
input,sym,cell,sliceplot,noco)
CALL timestop("Plotting")
......@@ -116,7 +116,7 @@ CONTAINS
! --->generate starting charge density
!
strho=input%strho
IF (.NOT.(strho.OR.sliceplot%iplot)) THEN
IF (.NOT.(strho.OR.(sliceplot%iplot.NE.0))) THEN
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
IF (noco%l_noco) THEN
archiveType = CDN_ARCHIVE_TYPE_NOCO_const
......@@ -168,7 +168,7 @@ CONTAINS
ENDIF ! mpi%irank == 0
IF (sliceplot%iplot) CALL juDFT_end("density plot o.k.",mpi%irank)
IF (sliceplot%iplot.NE.0) CALL juDFT_end("density plot o.k.",mpi%irank)
IF (input%strho) CALL juDFT_end("starting density generated",mpi%irank)
IF (input%swsp) CALL juDFT_end("spin polarised density generated",mpi%irank)
IF (input%lflip) CALL juDFT_end("magnetic moments flipped",mpi%irank)
......
mpi/mpi_bc_all.F90
View file @ d2438580
......@@ -38,9 +38,9 @@ CONTAINS
INTEGER n
REAL rdum
! .. Local Arrays ..
INTEGER i(42),ierr(3)
INTEGER i(43),ierr(3)
REAL r(34)
LOGICAL l(46)
LOGICAL l(45)
! ..
! .. External Subroutines..
#ifdef CPP_MPI
......@@ -57,7 +57,7 @@ CONTAINS
i(27)=vacuum%nstars ; i(28)=vacuum%nstm ; i(29)=oneD%odd%nq2 ; i(30)=oneD%odd%nop
i(31)=input%gw ; i(32)=input%gw_neigd ; i(33)=hybrid%ewaldlambda ; i(34)=hybrid%lexp
i(35)=hybrid%bands1 ; i(36)=input%maxiter ; i(37)=input%imix ; i(38)=banddos%orbCompAtom
i(39)=input%kcrel;i(40)=banddos%s_cell_x;i(41)=banddos%s_cell_y;i(42)=banddos%s_cell_z
i(39)=input%kcrel;i(40)=banddos%s_cell_x;i(41)=banddos%s_cell_y;i(42)=banddos%s_cell_z; i(43)=sliceplot%iplot
r(1)=cell%omtil ; r(2)=cell%area ; r(3)=vacuum%delz ; r(4)=cell%z1 ; r(5)=input%alpha
r(6)=sliceplot%e1s ; r(7)=sliceplot%e2s ; r(8)=noco%theta; r(9)=noco%phi; r(10)=vacuum%tworkf
......@@ -71,7 +71,7 @@ CONTAINS
l(1)=input%eonly ; l(2)=input%l_useapw ; l(3)=input%secvar ; l(4)=sym%zrfs ; l(5)=input%film
l(6)=sym%invs ; l(7)=sym%invs2 ; l(8)=input%l_bmt ; l(9)=input%l_f ; l(10)=input%cdinf
l(11)=banddos%dos ; l(12) = hybrid%l_hybrid ; l(13)=banddos%vacdos ; l(14)=input%integ ; l(15)=sliceplot%iplot
l(11)=banddos%dos ; l(12) = hybrid%l_hybrid ; l(13)=banddos%vacdos ; l(14)=input%integ; l(15)=noco%l_spav
l(16)=input%strho ; l(17)=input%swsp ; l(18)=input%lflip
l(21)=input%pallst ; l(22)=sliceplot%slice ; l(23)=noco%l_soc ; l(24)=vacuum%starcoeff
l(25)=noco%l_noco ; l(26)=noco%l_ss; l(27)=noco%l_mperp; l(28)=noco%l_constr
......@@ -80,7 +80,7 @@ CONTAINS
l(38)=field%efield%l_segmented
l(39)=sym%symor ; l(40)=input%frcor ; l(41)=input%tria ; l(42)=field%efield%dirichlet
l(43)=field%efield%l_dirichlet_coeff ; l(44)=input%l_coreSpec ; l(45)=input%ldauLinMix
l(46)=noco%l_spav
ENDIF
!
CALL MPI_BCAST(i,SIZE(i),MPI_INTEGER,0,mpi%mpi_comm,ierr)
......@@ -91,7 +91,7 @@ CONTAINS
sliceplot%nnne=i(17) ; banddos%ndir=i(18) ; stars%mx1=i(19) ; stars%mx2=i(20) ; stars%mx3=i(21)
input%jspins=i(12) ; vacuum%nvac=i(13) ; input%itmax=i(14) ; sliceplot%kk=i(15) ; vacuum%layers=i(16)
stars%ng2=i(7) ; stars%ng3=i(8) ; vacuum%nmz=i(9) ; vacuum%nmzxy=i(10) ; obsolete%lepr=i(11)
atoms%ntype=i(3) ; banddos%orbCompAtom=i(38);banddos%s_cell_x=i(40);banddos%s_cell_y=i(41);banddos%s_cell_z=i(42)
atoms%ntype=i(3) ; banddos%orbCompAtom=i(38);banddos%s_cell_x=i(40);banddos%s_cell_y=i(41);banddos%s_cell_z=i(42) ;sliceplot%iplot=i(43)
input%coretail_lmax=i(2) ; input%kcrel=i(39)
stars%kimax=i(25);stars%kimax2=i(26)
!
......@@ -114,14 +114,14 @@ CONTAINS
input%pallst=l(21) ; sliceplot%slice=l(22) ; noco%l_soc=l(23) ; vacuum%starcoeff=l(24)
input%strho=l(16) ; input%swsp=l(17) ; input%lflip=l(18)
banddos%dos=l(11) ; hybrid%l_hybrid=l(12) ; banddos%vacdos=l(13) ; banddos%l_orb=l(33) ; banddos%l_mcd=l(34)
input%integ=l(14) ; sliceplot%iplot=l(15)
input%integ=l(14)
sym%invs=l(6) ; sym%invs2=l(7) ; input%l_bmt=l(8) ; input%l_f=l(9) ; input%cdinf=l(10)
input%eonly=l(1) ; input%secvar=l(3) ; sym%zrfs=l(4) ; input%film=l(5)
field%efield%l_segmented = l(38) ; sym%symor=l(39); field%efield%dirichlet = l(40)
field%efield%l_dirichlet_coeff = l(41) ; input%l_coreSpec=l(44) ; input%ldauLinMix=l(45)
banddos%unfoldband=l(35)
noco%l_mtNocoPot=l(36)
noco%l_spav=l(46)
noco%l_spav=l(15)
!
! -> Broadcast the arrays:
IF (field%efield%l_segmented) THEN
......
optional/plot.f90 0 → 100644
View file @ d2438580
!--------------------------------------------------------------------------------
! Copyright (c) 2019 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_plot
USE m_types
USE m_juDFT
USE m_constants
USE m_cdn_io
USE m_loddop
USE m_wrtdop
USE m_qfix
USE m_fft2d
USE m_fft3d
USE m_types
USE m_rotdenmat
PRIVATE
!------------------------------------------------
! A general purpose plotting routine for FLEUR.
!
! Based on older plotting routines in pldngen.f90
! and plotdop.f90 called by optional.F90 and now
! called within a scf-loop instead of as a post
! process functionality.
!
! A. Neukirchen & R. Hilgers, September 2019
!------------------------------------------------
INTEGER, PARAMETER :: PLOT_INPDEN_const=2 !ind_plot= 1
INTEGER, PARAMETER :: PLOT_OUTDEN_Y_CORE_const=4 !ind_plot= 2
INTEGER, PARAMETER :: PLOT_INPDEN_N_CORE_const=8 !ind_plot= 4
INTEGER, PARAMETER :: PLOT_POT_TOT_const=128 !ind_plot= 7
INTEGER, PARAMETER :: PLOT_POT_EXT_const=256 !ind_plot= 8
INTEGER, PARAMETER :: PLOT_POT_COU_const=512 !ind_plot= 9
INTEGER, PARAMETER :: PLOT_POT_VXC_const=1024 !ind_plot=10
PUBLIC :: checkplotinp, makeplots, procplot, vectorsplit, matrixsplit, scalarplot, vectorplot, matrixplot
CONTAINS
SUBROUTINE checkplotinp()
! Checks for existing plot input. If an ancient plotin file is used, an
! error is called. If no usable plot_inp exists, a new one is generated.
oldform = .false.
INQUIRE(file = "plotin", exist = oldform)
IF (oldform) THEN
CALL juDFT_error("Use of plotin file no longer supported",calledby = "plotdop")
END IF
INQUIRE(file = "plot_inp", exist = newform)
IF (.NOT.newform) THEN
OPEN(20,file ="plot_inp")
WRITE(20,'(i2,a5,l1)') 2,",xsf=",.true.
WRITE(20,*) "&PLOT twodim=t,cartesian=t"
WRITE(20,*) " vec1(1)=10.0 vec2(2)=10.0"
WRITE(20,*) " filename='plot1' /"
WRITE(20,*) "&PLOT twodim=f,cartesian=f"
WRITE(20,*) " vec1(1)=1.0 vec1(2)=0.0 vec1(3)=0.0 "
WRITE(20,*) " vec2(1)=0.0 vec2(2)=1.0 vec2(3)=0.0 "
WRITE(20,*) " vec3(1)=0.0 vec3(2)=0.0 vec3(3)=1.0 "
WRITE(20,*) " grid(1)=30 grid(2)=30 grid(3)=30 "
WRITE(20,*) " zero(1)=0.0 zero(2)=0.0 zero(3)=0.5 "
WRITE(20,*) " filename ='plot2' /"
CLOSE(20)
END IF
END SUBROUTINE checkplotinp
!--------------------------------------------------------------------------------------------
SUBROUTINE vectorsplit(stars,vacuum,atoms,sphhar,input,noco,denmat,cden,mden)
! Takes a 2D potential/density vector and rearanges it into two plottable
! seperate ones (e.g. [rho_up, rho_down] ---> n, m).
IMPLICIT NONE
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_input), INTENT(IN) :: input
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), INTENT(INOUT) :: denmat
TYPE(t_potden), INTENT(OUT) :: cden, mden
CALL cden%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
CALL mden%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
CALL SpinsToChargeAndMagnetisation(denmat)
cden%mt(:,0:,1:,1) = denmat%mt(:,0:,1:,1)
cden%pw(1:,1) = denmat%pw(1:,1)
cden%vacz(1:,1:,1) = denmat%vacz(1:,1:,1)
cden%vacxy(1:,1:,1:,1) = denmat%vacxy(1:,1:,1:,1)
mden%mt(:,0:,1:,1) = denmat%mt(:,0:,1:,2)
mden%pw(1:,1) = denmat%pw(1:,2)
mden%vacz(1:,1:,1) = denmat%vacz(1:,1:,2)
mden%vacxy(1:,1:,1:,1) = denmat%vacxy(1:,1:,1:,2)
CALL ChargeAndMagnetisationToSpins(denmat)
END SUBROUTINE vectorsplit
!--------------------------------------------------------------------------------------------
SUBROUTINE matrixsplit(mpi,sym,stars,atoms,sphhar,vacuum,cell,input,noco,oneD,sliceplot,factor,denmat,cden,mxden,myden,mzden)
! Takes a 2x2 potential/density matrix and rearanges it into four plottable
! seperate ones (e.g. rho_mat ---> n, mx, my, mz).
!
! This is basically 1:1 the old pldngen.f90 routine, courtesy of Philipp Kurz.
IMPLICIT NONE
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_input), INTENT(IN) :: input
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_sliceplot), INTENT(IN) :: sliceplot
REAL, INTENT(IN) :: factor
TYPE(t_potden), INTENT(INOUT) :: denmat
TYPE(t_potden), INTENT(OUT) :: cden, mxden, myden, mzden
! Local type instances
TYPE(t_input) :: inp
TYPE(t_potden) :: den
! Local scalars
INTEGER iden,ivac,ifft2,ifft3
INTEGER imz,ityp,iri,ilh,imesh,iter
REAL cdnup,cdndown,chden,mgden,theta,phi,zero,rho_11,rziw,fermiEnergyTemp
REAL rho_22,rho_21r,rho_21i,rhotot,mx,my,mz,fix,vz_r,vz_i
COMPLEX czero
! Local arrays
!---> off-diagonal part of the density matrix
COMPLEX, ALLOCATABLE :: cdom(:),cdomvz(:,:),cdomvxy(:,:,:)
COMPLEX, ALLOCATABLE :: qpw(:,:),rhtxy(:,:,:,:)
REAL, ALLOCATABLE :: rht(:,:,:),rho(:,:,:,:)
REAL, ALLOCATABLE :: rvacxy(:,:,:,:),ris(:,:),fftwork(:)
!---> for testing: output of offdiag. output density matrix. to plot the
!---> offdiag. part of the output density matrix, that part has to be
!---> written the file rhomt21 in cdnmt.
LOGICAL :: l_qfix
REAL :: cdn11, cdn22
COMPLEX :: cdn21
!---> end of test part
iter = 0 ! This is not clean!
zero = 0.0 ; czero = CMPLX(0.0,0.0)
ifft3 = 27*stars%mx1*stars%mx2*stars%mx3
ifft2 = 9*stars%mx1*stars%mx2
ALLOCATE (qpw(stars%ng3,4),rhtxy(vacuum%nmzxyd,stars%ng2-1,2,4),&
cdom(stars%ng3),cdomvz(vacuum%nmzd,2),cdomvxy(vacuum%nmzxyd,stars%ng2-1,2),&
ris(0:27*stars%mx1*stars%mx2*stars%mx3-1,4),fftwork(0:27*stars%mx1*stars%mx2*stars%mx3-1),&
rvacxy(0:9*stars%mx1*stars%mx2-1,vacuum%nmzxyd,2,4),&
rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,4),rht(vacuum%nmzd,2,4) )
!---> initialize arrays for the density matrix
rho(:,:,:,:) = zero ; qpw(:,:) = czero ; cdom(:) = czero
IF (input%film) THEN
cdomvz(:,:) = czero ; rhtxy(:,:,:,:) = czero
cdomvxy(:,:,:) = czero ; rht(:,:,:) = zero
END IF
! Save the density matrix to a work density
CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
den=denmat
rho(:,0:,1:,:input%jspins) = factor*den%mt(:,0:,1:,:input%jspins)
qpw(1:,:input%jspins) = factor*den%pw(1:,:input%jspins)
rht(1:,1:,:input%jspins) = factor*den%vacz(1:,1:,:input%jspins)
rhtxy(1:,1:,1:,:input%jspins) = factor*den%vacxy(1:,1:,1:,:input%jspins)
IF(noco%l_noco) THEN
cdom = factor*den%pw(:,3)
cdomvz(:,:) = CMPLX(factor*den%vacz(:,:,3),factor*den%vacz(:,:,4))
cdomvxy = factor*den%vacxy(:,:,:,3)
END IF
IF (.NOT. sliceplot%slice) THEN
CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
den%iter = iter
den%mt(:,0:,1:,:input%jspins) = rho(:,0:,1:,:input%jspins)
den%pw(1:,:input%jspins) = qpw(1:,:input%jspins)
den%vacz(1:,1:,:input%jspins) = rht(1:,1:,:input%jspins)
den%vacxy(1:,1:,1:,:input%jspins) = rhtxy(1:,1:,1:,:input%jspins)
IF(noco%l_noco) THEN
den%pw(:,3) = cdom
den%vacz(:,:,3) = REAL(cdomvz(:,:))
den%vacz(:,:,4) = AIMAG(cdomvz(:,:))
den%vacxy(:,:,:,3) = cdomvxy
END IF
CALL qfix(mpi,stars,atoms,sym,vacuum,sphhar,input,cell,oneD,den,noco%l_noco,.FALSE.,.true.,fix)
rho(:,0:,1:,:input%jspins) = den%mt(:,0:,1:,:input%jspins)
qpw(1:,:input%jspins) = den%pw(1:,:input%jspins)
rht(1:,1:,:input%jspins) = den%vacz(1:,1:,:input%jspins)
rhtxy(1:,1:,1:,:input%jspins) = den%vacxy(1:,1:,1:,:input%jspins)
IF(noco%l_noco) THEN
cdom = den%pw(:,3)
cdomvz(:,:) = CMPLX(den%vacz(:,:,3),den%vacz(:,:,4))
cdomvxy = den%vacxy(:,:,:,3)
END IF
END IF
!---> calculate the charge and magnetization density in the muffin tins
DO ityp = 1,atoms%ntype
DO ilh = 0,sphhar%nlh(atoms%ntypsy(ityp))
DO iri = 1,atoms%jri(ityp)
IF (SIZE(den%mt,4).LE.2) THEN
cdnup = rho(iri,ilh,ityp,1)
cdndown = rho(iri,ilh,ityp,2)
theta = noco%beta(ityp)
phi = noco%alph(ityp)
chden = cdnup + cdndown
mgden = cdnup - cdndown
rho(iri,ilh,ityp,1) = chden
rho(iri,ilh,ityp,2) = mgden*COS(phi)*SIN(theta)
rho(iri,ilh,ityp,3) = mgden*SIN(phi)*SIN(theta)
rho(iri,ilh,ityp,4) = mgden*COS(theta)
ELSE
!---> for testing: output of offdiag. output density matrix
cdn11 = rho(iri,ilh,ityp,1)
cdn22 = rho(iri,ilh,ityp,2)
cdn21 = CMPLX(den%mt(iri,ilh,ityp,3),den%mt(iri,ilh,ityp,4))
CALL rot_den_mat(noco%alph(ityp),noco%beta(ityp),cdn11,cdn22,cdn21)
rho(iri,ilh,ityp,1) = cdn11 + cdn22
rho(iri,ilh,ityp,2) = 2.0*REAL(cdn21)
! Note: The minus sign in the following line is temporary to adjust for differences in the offdiagonal
! part of the density between this fleur version and ancient (v0.26) fleur.
rho(iri,ilh,ityp,3) = -2.0*AIMAG(cdn21)
rho(iri,ilh,ityp,4) = cdn11 - cdn22
!---> end of test part
END IF
END DO
END DO
END DO
!---> fouriertransform the diagonal part of the density matrix
!---> in the interstitial, qpw, to real space (ris)
DO iden = 1,2
CALL fft3d(ris(0,iden),fftwork,qpw(1,iden),stars,1)
END DO
!---> fouriertransform the off-diagonal part of the density matrix
CALL fft3d(ris(0,3),ris(0,4),cdom(1),stars,+1)
!---> calculate the charge and magnetization density on the
!---> real space mesh
DO imesh = 0,ifft3-1
rho_11 = ris(imesh,1)
rho_22 = ris(imesh,2)
rho_21r = ris(imesh,3)
rho_21i = ris(imesh,4)
rhotot = rho_11 + rho_22
mx = 2*rho_21r
my = -2*rho_21i
mz = (rho_11-rho_22)
ris(imesh,1) = rhotot
ris(imesh,2) = mx
ris(imesh,3) = my
ris(imesh,4) = mz
END DO
!---> Fouriertransform the density matrix back to reciprocal space
DO iden = 1,4
fftwork=zero
CALL fft3d(ris(0,iden),fftwork,qpw(1,iden),stars,-1)
END DO
!---> fouriertransform the diagonal part of the density matrix
!---> in the vacuum, rz & rxy, to real space (rvacxy)
IF (input%film) THEN
DO iden = 1,2
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
rziw = 0.0
CALL fft2d(stars,rvacxy(0,imz,ivac,iden),fftwork,rht(imz,ivac,iden),&
rziw,rhtxy(imz,1,ivac,iden),vacuum%nmzxyd,1)
END DO
END DO
END DO
!---> fouriertransform the off-diagonal part of the density matrix
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
rziw = 0.0
vz_r = REAL(cdomvz(imz,ivac))
vz_i = AIMAG(cdomvz(imz,ivac))
CALL fft2d(stars,rvacxy(0,imz,ivac,3),rvacxy(0,imz,ivac,4),&
vz_r,vz_i,cdomvxy(imz,1,ivac),vacuum%nmzxyd,1)
END DO
END DO
!---> calculate the four components of the matrix potential on
!---> real space mesh
DO ivac = 1,vacuum%nvac
DO imz = 1,vacuum%nmzxyd
DO imesh = 0,ifft2-1
rho_11 = rvacxy(imesh,imz,ivac,1)
rho_22 = rvacxy(imesh,imz,ivac,2)
rho_21r = rvacxy(imesh,imz,ivac,3)
rho_21i = rvacxy(imesh,imz,ivac,4)
rhotot = rho_11 + rho_22
mx = 2*rho_21r
my = -2*rho_21i
mz = (rho_11-rho_22)
rvacxy(imesh,imz,ivac,1) = rhotot
rvacxy(imesh,imz,ivac,2) = mx
rvacxy(imesh,imz,ivac,3) = my