Commit 03299042 authored by Robin Hilgers's avatar Robin Hilgers

Merge branch 'source-free' into 'develop'

Source free

See merge request fleur/fleur!40
parents cd027d89 85fccf07
......@@ -67,6 +67,8 @@ CONTAINS
USE m_metagga
USE m_plot
#ifdef CPP_MPI
USE m_mpi_bc_potden
#endif
......@@ -102,8 +104,7 @@ CONTAINS
TYPE(t_coreSpecInput) :: coreSpecInput
TYPE(t_wann) :: wann
TYPE(t_potden) :: vTot, vx, vCoul, vTemp, vxcForPlotting
TYPE(t_potden) :: inDen, outDen, EnergyDen, dummyDen
TYPE(t_potden), DIMENSION(3) :: testDen, testGrad
TYPE(t_potden) :: inDen, outDen, EnergyDen
CLASS(t_xcpot), ALLOCATABLE :: xcpot
CLASS(t_forcetheo), ALLOCATABLE :: forcetheo
......@@ -119,6 +120,8 @@ CONTAINS
INCLUDE 'mpif.h'
INTEGER :: ierr(2),n
#endif
REAL, ALLOCATABLE :: flh(:,:),flh2(:,:)
COMPLEX, ALLOCATABLE :: flm(:,:)
mpi%mpi_comm = mpi_comm
......@@ -525,27 +528,8 @@ IF (.FALSE.) CALL rotateMagnetFromSpinAxis(noco,vacuum,sphhar,stars,sym,oneD,cel
CALL juDFT_end("Stopped self consistency loop after plots have been generated.")
END IF
END DO scfloop ! DO WHILE (l_cont)
! Test: Build a field, for which the theoretical divergence etc. are known and
! compare with the result of the routine.
!CALL builddivtest(stars,atoms,sphhar,vacuum,sym,cell,1,testDen)
!CALL makeBxc(stars,atoms,sphhar,vacuum,input,noco,vTot,testDen)
!CALL matrixsplit(stars, atoms, sphhar, vacuum, input, noco, 1.0, inDen, dummyDen, testDen(1), testDen(2), testDen(3))
!CALL checkplotinp()
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDen ', dummyDen, testDen(1), testDen(2), testDen(3))
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDeny ', testDen(2))
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDenz ', testDen(3))
!CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,testDen)
!DO i=1,3
! CALL testGrad(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
! ALLOCATE(testGrad(i)%pw_w,mold=testGrad(i)%pw)
!ENDDO
!CALL divpotgrad(stars,atoms,sphhar,vacuum,sym,cell,noco,testDen(2),testGrad)
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testGrad ', testGrad(1), testGrad(1), testGrad(2), testGrad(3))
CALL add_usage_data("Iterations",iter)
IF (mpi%irank.EQ.0) CALL closeXMLElement('scfLoop')
......
......@@ -80,7 +80,7 @@ CONTAINS
!sum up both spins in den into workden
CALL den%sum_both_spin(workden)
CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,workden,vCoul,results)
CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,.FALSE.,workden,vCoul,results)
CALL vCoul%copy_both_spin(vTot)
vCoul%mt(:,:,:,input%jspins)=vCoul%mt(:,:,:,1)
......
......@@ -27,4 +27,7 @@ math/ExpSave.f90
math/intgr.F90
math/ylm4.F90
math/grule.f90
math/lh_tofrom_lm.f90
math/gradYlm.f90
math/analysistests.f90
)
MODULE m_analysistests
USE m_types
INTEGER, PARAMETER :: numtests=3
INTEGER, PARAMETER :: numders=5
! 1 = grdchlh with ndvgrd=3
! 5 = SPEX Derivative
INTEGER, PARAMETER :: numints=5
CONTAINS
SUBROUTINE buildfunctions(atoms,radii,funcsr,trueder,trueint)
TYPE(t_atoms), INTENT(IN) :: atoms
REAL, INTENT(OUT) :: radii(atoms%jmtd,atoms%ntype), funcsr(atoms%jmtd,atoms%ntype,numtests), trueder(atoms%jmtd,atoms%ntype,numtests), trueint(atoms%jmtd,atoms%ntype,numtests)
REAL :: Rmt(atoms%jmtd,atoms%ntype)
INTEGER :: n
radii=atoms%rmsh
DO n=1,atoms%ntype
Rmt(:,n)=atoms%rmt(n)
END DO
funcsr(:,:,1)=radii(:,:)**2
trueder(:,:,1)=2*radii(:,:)
trueint(:,:,1)=(1./3.)*radii(:,:)**3
funcsr(:,:,2)=radii(:,:)**3
trueder(:,:,2)=3*radii(:,:)**2
trueint(:,:,2)=(1./4.)*radii(:,:)**4
funcsr(:,:,3)=radii(:,:)**2-(radii(:,:)**3)/Rmt(:,:)
trueder(:,:,3)=2*radii(:,:)-3*(radii(:,:)**2)/Rmt(:,:)
trueint(:,:,3)=(1./3.)*radii(:,:)**3-(1./4.)*(radii(:,:)**4)/Rmt(:,:)
END SUBROUTINE buildfunctions
SUBROUTINE derivtest(atoms,radii,funcsr,testders)
USE m_grdchlh
USE m_gradYlm
TYPE(t_atoms), INTENT(IN) :: atoms
REAL, INTENT(IN) :: radii(atoms%jmtd,atoms%ntype), funcsr(atoms%jmtd,atoms%ntype,numtests)
REAL, INTENT(OUT) :: testders(atoms%jmtd,atoms%ntype,numtests,numders)
INTEGER :: n
REAL :: junkder(atoms%jmtd)
DO n=1,atoms%ntype
CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),radii(:,n),funcsr(:,n,1),3, testders(:,n,1,1),junkder)
CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),radii(:,n),funcsr(:,n,1),4, testders(:,n,1,2),junkder)
CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),radii(:,n),funcsr(:,n,1),5, testders(:,n,1,3),junkder)
CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),radii(:,n),funcsr(:,n,1),6, testders(:,n,1,4),junkder)
CALL Derivative(funcsr(:,n,1), n, atoms, testders(:,n,1,5))
CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),radii(:,n),funcsr(:,n,2),3, testders(:,n,2,1),junkder)
CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),radii(:,n),funcsr(:,n,2),4, testders(:,n,2,2),junkder)
CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),radii(:,n),funcsr(:,n,2),5, testders(:,n,2,3),junkder)
CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),radii(:,n),funcsr(:,n,2),6, testders(:,n,2,4),junkder)
CALL Derivative(funcsr(:,n,2), n, atoms, testders(:,n,2,5))
CALL Derivative(funcsr(:,n,3), n, atoms, testders(:,n,3,5))
END DO
END SUBROUTINE derivtest
SUBROUTINE integtest(atoms,radii,funcsr,testints)
USE m_intgr
TYPE(t_atoms), INTENT(IN) :: atoms
REAL, INTENT(IN) :: radii(atoms%jmtd,atoms%ntype), funcsr(atoms%jmtd,atoms%ntype,numtests)
REAL, INTENT(OUT) :: testints(atoms%jmtd,atoms%ntype,numtests,numints)
INTEGER :: n
DO n=1,atoms%ntype
CALL intgr2(funcsr(:atoms%jri(n),n,1),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,1,1))
CALL intgr2(funcsr(:atoms%jri(n),n,2),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,2,1))
CALL intgr2(funcsr(:atoms%jri(n),n,3),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,3,1))
CALL intgrt(funcsr(:atoms%jri(n),n,1),radii(:,n),atoms%jri(n),testints(:atoms%jri(n),n,1,2))
CALL intgrt(funcsr(:atoms%jri(n),n,2),radii(:,n),atoms%jri(n),testints(:atoms%jri(n),n,2,2))
CALL intgrt(funcsr(:atoms%jri(n),n,3),radii(:,n),atoms%jri(n),testints(:atoms%jri(n),n,3,2))
CALL intgrtlog(funcsr(:atoms%jri(n),n,1),radii(:,n),atoms%jri(n),testints(:atoms%jri(n),n,1,3))
CALL intgrtlog(funcsr(:atoms%jri(n),n,2),radii(:,n),atoms%jri(n),testints(:atoms%jri(n),n,2,3))
CALL intgrtlog(funcsr(:atoms%jri(n),n,3),radii(:,n),atoms%jri(n),testints(:atoms%jri(n),n,3,3))
CALL intgr4(funcsr(:atoms%jri(n),n,1),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,1,4))
CALL intgr4(funcsr(:atoms%jri(n),n,2),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,2,4))
CALL intgr4(funcsr(:atoms%jri(n),n,3),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,3,4))
CALL intgr5(funcsr(:atoms%jri(n),n,1),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,1,5))
CALL intgr5(funcsr(:atoms%jri(n),n,2),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,2,5))
CALL intgr5(funcsr(:atoms%jri(n),n,3),radii(:,n),atoms%dx(n),atoms%jri(n),testints(:atoms%jri(n),n,3,5))
END DO
END SUBROUTINE integtest
END MODULE m_analysistests
This diff is collapsed.
......@@ -13,14 +13,19 @@ MODULE m_intgr
! decaying exponential between the first mesh point and infinity.
! y contains the nmz function values tabulated at a spacing of h.
!
! intgrt:
! Uses the basic trapezoidal rule of integration. Indefinite.
! - A. Neukirchen, '19
!
! integrator: ---- input ---- output
! intgr0: definite y r0 h jri | z
! intgr1: indefinite y r1 h jri | z(jri)
! intgr2: indefinite y rmsh(jri) h jri | z(jri)
! intgr3: definite y rmsh(jri) h jri | z
! intgz0: definite y tail n nmz | z
! intgz1: indefinite y tail n nmz | z(nmz)
!
! intgz0: definite y tail n nmz ! z
! intgz1: indefinite y tail n nmz ! z(nmz)
! intgrt: indefinite y x jri | z(jri)
!
! m. weinert
!**********************************************************************
......@@ -472,6 +477,112 @@ MODULE m_intgr
end subroutine intgz1ComplexReverse
SUBROUTINE intgrt(y,x,jri,z)
INTEGER, INTENT (IN) :: jri
REAL, INTENT (IN) :: x(jri), y(jri)
REAL, INTENT (OUT) :: z(jri)
INTEGER :: i
REAL :: alpha, h
z = zero
h=LOG(x(2))-LOG(x(1))
IF (y(1)*y(2).GT.zero) THEN
alpha = 1.0 + log(y(2)/y(1))/h
IF (alpha.GT.zero) z(1) = x(1)*y(1)/alpha
ENDIF
z=z(1)
DO i=2, jri
z(i:)=z(i:)+(y(i-1)+y(i))*(x(i)-x(i-1))/2.0
END DO
END SUBROUTINE intgrt
SUBROUTINE intgrtlog(y,x,jri,z)
INTEGER, INTENT (IN) :: jri
REAL, INTENT (IN) :: x(jri), y(jri)
REAL, INTENT (OUT) :: z(jri)
INTEGER :: i
REAL :: logr(jri),rf(jri)
REAL :: dr, alpha
logr=LOG(x)
dr=x(2)-x(1)
rf=x*y
DO i=2, jri
z(i:)=z(i:)+(rf(i-1)+rf(i))*dr/2.0
END DO
END SUBROUTINE intgrtlog
! Testwise: optional integrators for source-free purposes.
SUBROUTINE intgr4(y,rmsh,h,jri,z)
! Modified version of intgr2 with a different approach to the first few
! points. For point 1 through 6 we use the trapezoid integrator.
INTEGER, INTENT (IN) :: jri
REAL, INTENT (IN) :: h
REAL, INTENT (IN) :: rmsh(jri), y(jri)
REAL, INTENT (OUT):: z(jri)
REAL :: dr, r(7)
INTEGER :: i, j
REAL :: yr(7)
CALL intgrt(y(1:nr1),rmsh(1:nr1),nr1,z(1:nr1))
!---> simpson integration, j>nr-1
dr = exp(h)
DO i = 1,7
r(i) = rmsh(i)
yr(i) = rmsh(i)*y(i)
ENDDO
DO i = 1,nr
r(i) = h*ih(i)*r(i)/h0
ENDDO
DO j = nr,jri
z(j) = z(j-nr1) + CPP_BLAS_sdot(nr,r,1,y(j-nr1),1)
DO i = 1,7
r(i) = dr*r(i)
ENDDO
ENDDO
RETURN
END SUBROUTINE intgr4
SUBROUTINE intgr5(y,x,h,jri,z)
INTEGER, INTENT (IN) :: jri
REAL, INTENT (IN) :: h
REAL, INTENT (IN) :: x(jri), y(jri)
REAL, INTENT (OUT):: z(jri)
REAL :: dr, r(7), alpha
INTEGER :: i, j
REAL :: yr(7)
z(1) = zero
IF (y(1)*y(2).GT.zero) THEN
alpha = 1.0 + log(y(2)/y(1))/h
IF (alpha.GT.zero) z(1) = x(1)*y(1)/alpha
ENDIF
z(2)=z(1)+h*(x(1)*y(1)+x(2)*y(2))/2
z(3)=z(1)+h*(x(1)*y(1)+4*x(2)*y(2)+x(3)*y(3))/3.
z(4)=z(1)+3*h*(x(1)*y(1)+3*x(2)*y(2)+3*x(3)*y(3)+x(4)*y(4))/8.
DO j = 5,jri
z(j) = z(j-4) + 2*h*(7*x(j-4)*y(j-4)+32*x(j-3)*y(j-3)+12*x(j-2)*y(j-2)+32*x(j-1)*y(j-1)+7*x(j)*y(j))/45
ENDDO
RETURN
END SUBROUTINE intgr5
END MODULE m_intgr
MODULE m_lh_tofrom_lm
CONTAINS
SUBROUTINE lh_to_lm(atoms, lathar, iType, flh, flm)
USE m_types
IMPLICIT NONE
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: lathar
INTEGER, INTENT(IN) :: iType
REAL, INTENT(IN) :: flh(:,0:) ! (iR,iLH)
COMPLEX, INTENT(INOUT) :: flm(:,:) ! (iR,lm)
INTEGER :: iAtom, iLH, ns, l, iM, m, lm, iR, lh, imem
COMPLEX, ALLOCATABLE :: cMat(:,:), gVec(:,:)
REAL, ALLOCATABLE :: fVec(:,:)
iAtom = SUM(atoms%neq(:iType-1)) + 1
ns = atoms%ntypsy(iAtom)
flm = CMPLX(0.0,0.0)
DO l = 0, atoms%lmax(iType)
ALLOCATE (cMat(-l:l,-l:l))
ALLOCATE (fVec(-l:l,atoms%jri(iType)))
ALLOCATE (gVec(-l:l,atoms%jri(iType)))
cMat=CMPLX(0.0,0.0)
DO M = -l, l
lh = l*(l+1)+M
fVec(M,:)=flh(:,lh)
DO imem = 1, lathar%nmem(lh,ns)
cMat(M,lathar%mlh(imem,lh,ns))=lathar%clnu(imem,lh,ns)
END DO
END DO
DO iR = 1, atoms%jri(iType)
gVec(:,iR)=MATMUL(fVec(:,iR),cMat)
END DO
DO m = -l, l
flm(:,l*(l+1)+1+m)=gVec(m,:)
END DO
DEALLOCATE (cMat,fVec,gVec)
END DO
END SUBROUTINE lh_to_lm
SUBROUTINE lh_from_lm(atoms, lathar, iType, flm, flh)
USE m_types
IMPLICIT NONE
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: lathar
INTEGER, INTENT(IN) :: iType
COMPLEX, INTENT(IN) :: flm(:,:) ! (iR,lm)
REAL, INTENT(INOUT) :: flh(:,0:) ! (iR,iLH)
INTEGER :: iAtom, iLH, ns, l, iM, m, lm, iR, info, lh, imem, info2
INTEGER, ALLOCATABLE :: ipiv(:)
COMPLEX, ALLOCATABLE :: cMat(:,:),fVec(:,:)
EXTERNAL zgetrf, zgetrs
iAtom = SUM(atoms%neq(:iType-1)) + 1
ns = atoms%ntypsy(iAtom)
flh = 0.0
DO l = 0, atoms%lmax(iType)
ALLOCATE (cMat(-l:l,-l:l))
ALLOCATE (fVec(-l:l,atoms%jri(iType)))
ALLOCATE (ipiv(-l:l))
cMat=CMPLX(0.0,0.0)
DO M = -l, l
lh = l*(l+1)+M
fVec(M,:)=flm(:,lh+1)
DO imem = 1, lathar%nmem(lh,ns)
cMat(M,lathar%mlh(imem,lh,ns))=lathar%clnu(imem,lh,ns)
END DO
END DO
CALL zgetrf(2*l+1,2*l+1,cMat,2*l+1,ipiv,info)
DO iR = 1, atoms%jri(iType)
CALL zgetrs('T',2*l+1,1,cMat,2*l+1,ipiv,fVec(:,iR),2*l+1,info2)
END DO
DO M = -l, l
flh(:,l*(l+1)+M)=REAL(fVec(M,:))
END DO
DEALLOCATE (cMat,fVec,ipiv)
END DO
END SUBROUTINE lh_from_lm
END MODULE m_lh_tofrom_lm
......@@ -58,7 +58,7 @@ CONTAINS
#endif
IF ( .NOT. input%film ) THEN
CALL vgen_coulomb( 1, mpi, DIMENSION, oneD, input, field, vacuum, sym, stars, cell, &
sphhar, atoms, resDen, vYukawa )
sphhar, atoms, .FALSE., resDen, vYukawa )
ELSE
if( mpi%irank == 0 ) then
call resDenMod%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN )
......
......@@ -13,15 +13,11 @@ MODULE m_plot
!
! Based on the older plotting routines pldngen.f90 and plotdop.f90 that were
! originally called by optional.F90 and are now used in the scf-loop instead.
! At the cost of no reduced postprocess functionality, this allowed us to re-
! move I/O (using plot.hdf files) from the plotting routine completely.
! At the cost of reduced postprocess functionality, this allowed us to re-
! move I/O (using plot.hdf/text files) from the plotting routine completely.
!
! TODO:
! - plot_inp files are still in use and should be replaced by a plotting type.
! - Only the .xsf format for xcrysden is supported right now. Further exten-
! sion should include several plotting options, most importantly a way to
! neatly plot vectorial quantities like the magnetisation density as vectors
! on a grid.
!
! A. Neukirchen & R. Hilgers, September 2019
!------------------------------------------------
......@@ -438,15 +434,6 @@ CONTAINS
mzden%vacz(1:,1:,1) = rht(1:,1:,4)
mzden%vacxy(1:,1:,1:,1) = rhtxy(1:,1:,1:,4)
!DO ityp=1, atoms%ntype
! DO iri=1, atoms%jri(ityp)
! cden%mt(iri,0:,ityp,:) = cden%mt(iri,0:,ityp,:) * atoms%rmsh(iri,ityp)**2
! mxden%mt(iri,0:,ityp,:) = mxden%mt(iri,0:,ityp,:) * atoms%rmsh(iri,ityp)**2
! myden%mt(iri,0:,ityp,:) = myden%mt(iri,0:,ityp,:) * atoms%rmsh(iri,ityp)**2
! mzden%mt(iri,0:,ityp,:) = mzden%mt(iri,0:,ityp,:) * atoms%rmsh(iri,ityp)**2
! END DO
!END DO
END IF
DEALLOCATE (rho, qpw, rht, rhtxy, cdomvz, cdomvxy, &
......@@ -736,55 +723,44 @@ CONTAINS
xdnout(5) = SQRT(ABS(xdnout(2)**2+xdnout(3)**2+xdnout(4)**2))
IF (xdnout(5)<eps) THEN
xdnout(5)= 0.0
!xdnout(6)= -tpi_const
!xdnout(7)= -tpi_const
xdnout(6)= 0.0
xdnout(7)= 0.0
xdnout(6)= -tpi_const
xdnout(7)= -tpi_const
ELSE
DO j = 1, 3
help(j) = xdnout(1+j)/xdnout(5)
END DO
!IF (help(3)<0.5) THEN
IF (help(3)<0.5) THEN
xdnout(6)= ACOS(help(3))
!ELSE
! xdnout(6)= pi_const/2.0-ASIN(help(3))
!END IF
ELSE
xdnout(6)= pi_const/2.0-ASIN(help(3))
END IF
IF (SQRT(ABS(help(1)**2+help(2)**2)) < eps) THEN
!xdnout(7)= -tpi_const
xdnout(7)= 0.0
xdnout(7)= -tpi_const
ELSE
!IF ( ABS(help(1)) > ABS(help(2)) ) THEN
! xdnout(7)= ABS(ATAN(help(2)/help(1)))
!ELSE
! xdnout(7)= pi_const/2.0-ABS(ATAN(help(1)/help(2)))
!END IF
!IF (help(2)<0.0) THEN
! xdnout(7)= -xdnout(7)
!END IF
!IF (help(1)<0.0) THEN
! xdnout(7)= pi_const-xdnout(7)
!END IF
!phi0=0
!DO WHILE (xdnout(7)-pi_const*phi0 > +pi_const)
! xdnout(7)= xdnout(7)-tpi_const
!END DO
!DO WHILE (xdnout(7)-pi_const*phi0 < -pi_const)
! xdnout(7)= xdnout(7)+tpi_const
!END DO
!IF (ABS(xdnout(2)-xdnout(3))<eps) THEN
! IF (xdnout(2)>0) THEN
! xdnout(7)=pi_const/4.0
! ELSE
! xdnout(7)=-3*pi_const/4.0
! END IF
IF (xdnout(2)>eps) THEN
xdnout(7)=ATAN(xdnout(3)/xdnout(2))
ELSE IF (ABS(xdnout(2))<eps) THEN
xdnout(7)=SIGN(1.0, xdnout(3))*pi_const/2.0
ELSE IF ((xdnout(2)<-eps).AND.(xdnout(3)>=0.0)) THEN
xdnout(7)=ATAN(xdnout(3)/xdnout(2))+pi_const
IF ( ABS(help(1)) > ABS(help(2)) ) THEN
xdnout(7)= ABS(ATAN(help(2)/help(1)))
ELSE
xdnout(7)=ATAN(xdnout(3)/xdnout(2))-pi_const
xdnout(7)= pi_const/2.0-ABS(ATAN(help(1)/help(2)))
END IF
IF (help(2)<0.0) THEN
xdnout(7)= -xdnout(7)
END IF
IF (help(1)<0.0) THEN
xdnout(7)= pi_const-xdnout(7)
END IF
phi0=0
DO WHILE (xdnout(7)-pi_const*phi0 > +pi_const)
xdnout(7)= xdnout(7)-tpi_const
END DO
DO WHILE (xdnout(7)-pi_const*phi0 < -pi_const)
xdnout(7)= xdnout(7)+tpi_const
END DO
IF (ABS(xdnout(2)-xdnout(3))<eps) THEN
IF (xdnout(2)>0) THEN
xdnout(7)=pi_const/4.0
ELSE
xdnout(7)=-3*pi_const/4.0
END IF
END IF
END IF
END IF
......@@ -1071,52 +1047,6 @@ CONTAINS
END SUBROUTINE procplot
SUBROUTINE plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
noco, div, phi, divGrx, divGry, divGrz, &
xcBmodx, xcBmody, xcBmodz, div2)
IMPLICIT NONE
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_input), INTENT(IN) :: input
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), INTENT(IN) :: div, phi, divGrx, divGry, divGrz, &
xcBmodx, xcBmody, xcBmodz, div2
LOGICAL :: xsf
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .FALSE., 'div ', div)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'phiDiv ', phi)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'gradPhiDiv ', divGrx, divGrx, divGry, divGrz)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .TRUE., 'bCorrected ', xcBmodx, xcBmodx, xcBmody, xcBmodz)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .FALSE., 'divCorrected ', div2)
INQUIRE(file="div.xsf",exist=xsf)
IF (xsf) THEN
OPEN (120, FILE='gradPhiDiv_f.xsf', STATUS='OLD')
CLOSE (120, STATUS="DELETE")
OPEN (120, FILE='bCorrected_f.xsf', STATUS='OLD')
CLOSE (120, STATUS="DELETE")
END IF
END SUBROUTINE plotBtest
SUBROUTINE makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
noco, denmat, plot_const, sliceplot)
......
......@@ -451,7 +451,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars