Skip to content

GitLab

fleur
fleur
Commit 89eed34b authored by Alexander Neukirchen's avatar Alexander Neukirchen
Browse files
Options

Various changes. Sourcefree sort of works now.

parent 40c0881e
Hide whitespace changes
Inline Side-by-side
Showing with 447 additions and 272 deletions
cdn/m_perp.f90
View file @ 89eed34b
......@@ -18,7 +18,7 @@ CONTAINS
USE m_constants
USE m_intgr, ONLY : intgr3
USE m_sphcoord
USE m_polangle
USE m_rotdenmat
USE m_types
IMPLICIT NONE
......@@ -59,7 +59,7 @@ CONTAINS
absmag=SQRT(mx*mx+my*my+mz*mz)
WRITE (6,8025) mx,my,mz,absmag
!---> determine the polar angles of the moment vector in the local frame
CALL sphcoord(mx,my,mz,betah,alphh)
CALL pol_angle(mx,my,mz,betah,alphh)
WRITE (6,8026) betah,alphh
8025 FORMAT(2x,'--> local frame: ','mx=',f9.5,' my=',f9.5,' mz=',f9.5,' |m|=',f9.5)
8026 FORMAT(2x,'-->',10x,' delta beta=',f9.5,&
......@@ -76,7 +76,7 @@ CONTAINS
mx = 2*REAL(rho21)
my = 2*AIMAG(rho21)
mz = rho11 - rho22
CALL sphcoord(mx,my,mz,betah,alphh)
CALL pol_angle(mx,my,mz,betah,alphh)
WRITE (6,8027) noco%beta(itype),noco%alph(itype)-alphdiff
WRITE (6,8028) betah,alphh-alphdiff
8027 FORMAT(2x,'-->',10x,' input noco%beta=',f9.5, ' input noco%alpha=',f9.5)
......@@ -103,7 +103,7 @@ CONTAINS
mz_mix = (-1.0) * mz_mix
ENDIF
! calculate angles alpha and beta in global frame
CALL sphcoord(mx_mix,my_mix,mz_mix,betah,alphh)
CALL pol_angle(mx_mix,my_mix,mz_mix,betah,alphh)
WRITE (6,8029) betah,alphh-alphdiff
8029 FORMAT(2x,'-->',10x,' new noco%beta =',f9.5, ' new noco%alpha =',f9.5)
noco%alph(itype) = alphh
......
main/fleur.F90
View file @ 89eed34b
......@@ -522,12 +522,21 @@ CONTAINS
!CALL builddivtest(stars,atoms,sphhar,vacuum,sym,cell,1,testDen)
!CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,inDen,1.0,testDen)
!CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,vtot,2.0,testDen)
CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,vtot,2.0,testDen)
!testDen(3)%mt(:,1,:,1)=testDen(3)%mt(:,0,:,1)*atoms%rmsh
testDen(3)%mt(:,1:,:,:)=0.0
!testDen(3)%mt(:,2:,:,:)=0.0
!testDen(3)%mt(:,0,:,:)=0.0
testDen(2)%mt(:,:,:,:)=0.0
testDen(1)%mt(:,:,:,:)=0.0
!testDen(3)%mt(:,0,:,1)*atoms%rmsh
!testDen(3)%mt(:,0,:,1)=0.0
!CALL checkplotinp()
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testDen ', testDen(1), 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)
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)
......@@ -535,18 +544,18 @@ CONTAINS
!CALL divpotgrad(stars,atoms,sphhar,vacuum,sym,cell,noco,testDen(3),testGrad)
!CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, .FALSE., .FALSE., 'testGrad ', testGrad(1), testGrad(1), testGrad(2), testGrad(3))
ALLOCATE (flh(atoms%jri(1),0:sphhar%nlh(atoms%ntypsy(1))),flm(atoms%jri(1),sphhar%nlh(atoms%ntypsy(1))+1),flh2(atoms%jri(1),0:sphhar%nlh(atoms%ntypsy(1))))
flh=inDen%mt(:,:,1,1)
flh(:,1)=-flh(:,0)
flh(:,2)=0*flh(:,0)
flh(:,3)=flh(:,0)
flh(:,4)=flh(:,0)
flh(:,5)=2*flh(:,0)
flh(:,6)=3*flh(:,0)
flh(:,7)=4*flh(:,0)
flh(:,8)=5*flh(:,0)
CALL lh_to_lm(atoms, sphhar, 1, flh, flm)
CALL lh_from_lm(atoms, sphhar, 1, flm, flh2)
!ALLOCATE (flh(atoms%jri(1),0:sphhar%nlh(atoms%ntypsy(1))),flm(atoms%jri(1),sphhar%nlh(atoms%ntypsy(1))+1),flh2(atoms%jri(1),0:sphhar%nlh(atoms%ntypsy(1))))
!flh=inDen%mt(:,:,1,1)
!flh(:,1)=-flh(:,0)
!flh(:,2)=0*flh(:,0)
!flh(:,3)=flh(:,0)
!flh(:,4)=flh(:,0)
!flh(:,5)=2*flh(:,0)
!flh(:,6)=3*flh(:,0)
!flh(:,7)=4*flh(:,0)
!flh(:,8)=5*flh(:,0)
!CALL lh_to_lm(atoms, sphhar, 1, flh, flm)
!CALL lh_from_lm(atoms, sphhar, 1, flm, flh2)
CALL add_usage_data("Iterations",iter)
......
math/CMakeLists.txt
View file @ 89eed34b
......@@ -33,6 +33,7 @@ math/ExpSave.f90
math/intgr.F90
math/ylm4.F90
math/lh_tofrom_lm.f90
math/gradYlm.f90
)
if (FLEUR_USE_FFTMKL)
set(fleur_F90 ${fleur_F90} math/mkl_dfti.f90)
......
math/gradYlm.f90 0 → 100644
View file @ 89eed34b
MODULE m_gradYlm
CONTAINS
! Not from me, taken from Spex, ask Christoph! Calculates derivative of a function with scalar argument lying on a muffin-tin mesh
subroutine Derivative(f, itype, atoms, df)
use m_types
implicit none
type(t_atoms), intent(in) :: atoms
integer, intent(in) :: itype
real, intent(in) :: f(atoms%jri(itype))
real, intent(out) :: df(atoms%jri(itype))
real :: h, r, d21, d32, d43, d31, d42, d41, df1, df2, s
real :: y0, y1, y2
integer :: i, n
n = atoms%jri(itype)
h = atoms%dx(itype)
r = atoms%rmsh(1, itype)
! use Lagrange interpolation of 3rd order (and averaging) for points 3 to n
d21 = r * (exp(h)-1) ; d32 = d21 * exp(h) ; d43 = d32 * exp(h)
d31 = d21 + d32 ; d42 = d32 + d43
d41 = d31 + d43
df(1) = d31*d41 / (d21*d32*d42) * f(2) + ( -1d0/d21 - 1d0/d31 - 1d0/d41) * f(1)&
& - d21*d41 / (d31*d32*d43) * f(3) + d21*d31 / (d41*d42*d43) * f(4)
df(2) = - d32*d42 / (d21*d31*d41) * f(1) + ( 1d0/d21 - 1d0/d32 - 1d0/d42) * f(2)&
& + d21*d42 / (d31*d32*d43) * f(3) - d21*d32 / (d41*d42*d43) * f(4)
df1 = d32*d43 / (d21*d31*d41) * f(1) - d31*d43 / (d21*d32*d42) * f(2) +&
& ( 1d0/d31 + 1d0/d32 - 1d0/d43 ) * f(3) + d31*d32 / (d41*d42*d43) * f(4)
do i = 3, n - 2
d21 = d32 ; d32 = d43 ; d43 = d43 * exp(h)
d31 = d42 ; d42 = d42 * exp(h)
d41 = d41 * exp(h)
df2 = - d32*d42 / (d21*d31*d41) * f(i-1) + ( 1d0/d21 - 1d0/d32 - 1d0/d42) * f(i) + &
& d21*d42 / (d31*d32*d43) * f(i+1) - d21*d32 / (d41*d42*d43) * f(i+2)
df(i) = ( df1 + df2 ) / 2
df1 = d32*d43 / (d21*d31*d41) * f(i-1) - d31*d43 / (d21*d32*d42) * f(i) +&
& ( 1d0/d31 + 1d0/d32 - 1d0/d43 ) * f(i+1) + d31*d32 / (d41*d42*d43) * f(i+2)
enddo
df(n-1) = df1
df(n) = - d42*d43 / (d21*d31*d41) * f(n-3) + d41*d43 / (d21*d32*d42) * f(n-2) -&
& d41*d42 / (d31*d32*d43) * f(n-1) + ( 1d0/d41 + 1d0/d42 + 1d0/d43 ) * f(n)
! for first two points use Lagrange interpolation of second order for log(f(i))
! or, as a fall-back, Lagrange interpolation with the conditions f(1), f(2), f(3), f'(3).
s = sign(1d0,f(1))
if(sign(1d0,f(2)) /= s .or. sign(1d0,f(3)) /= s .or. any(f(:3) == 0)) then
d21 = r * (exp(h)-1)
d32 = d21 * exp(h)
d31 = d21 + d32
s = df(3) / (d31*d32) - f(1) / (d21*d31**2) + f(2) / (d21*d32**2) - f(3) / (d31**2*d32) - f(3) / (d31*d32**2)
df(1) = - (d21+d31) / (d21*d31) * f(1) + d31 / (d21*d32) * f(2) - d21 / (d31*d32) * f(3) + d21*d31 * s
df(2) = - (d21-d32) / (d21*d32) * f(2) - d32 / (d21*d31) * f(1) + d21 / (d31*d32) * f(3) - d21*d32 * s
else
y0 = log(abs(f(1)))
y1 = log(abs(f(2)))
y2 = log(abs(f(3)))
df(1) = ( - 3*y0/2 + 2*y1 - y2/2 ) * f(1) / (h*r)
df(2) = (y2-y0)/2 * f(2) / (h*r*exp(h))
endif
end subroutine Derivative
SUBROUTINE gradYlm(atoms, r2FshMt, r2GrFshMt)
! Based on work for juphon by C. Gerhorst.
!---------------------------------------------------------------------------------------------------------------------------------
!> @author
!> Christian-Roman Gerhorst, Forschungszentrum Jülich: IAS1 / PGI1
!>
!> @brief
!> Calculates the spherical harmonic expansion coefficients of the muffin-tin gradient applied to an arbitrary function multiplied
!> by $r^2$. The resulting gradient expansion coefficients are multiplied by a factor $r^2$.
!>
!> @note
!> The ingoing function is assumed to be multiplied with $r^2$. The outgoing resulting function is also multiplied by $r^2$.
!>
!> @param[in] atoms : Contains atoms-related quantities; definition of its members in types.F90 file.
!> @param[in] r2FshMt : Lattice harmonic coefficients of muffin-tin quantity multiplied by a factor of r**2.
!> @param[out] r2GrFshMt : Spherical harmonic coefficients of muffin-tin quantity's gradient multiplied by a factor of r**2
!---------------------------------------------------------------------------------------------------------------------------------
use m_constants, only : fpi_const, ImagUnit
use m_gaunt, only : gaunt1
use m_types
USE m_grdchlh
implicit none
! Type parameter
! ***************
type(t_atoms), intent(in) :: atoms
! Array parameters
! ****************
complex, allocatable, intent(in) :: r2FshMt(:, :, :)
complex, allocatable, intent(out) :: r2GrFshMt(:, :, :, :)
! Local Scalar Variables
! **********************
! pfac : Prefactor
! tGaunt : Gaunt coefficient
! itype : Loop index for atom types
! ieqat : Loop index for equivalent atoms
! iatom : Loop index for all atoms
! imesh : Loop index for radial mesh point
! mqn_m : Magnetic quantum number m
! oqn_l : Orbital quantum number l
! mqn_mpp : Magnetic quantum number double primed to index the natural coordinates
! lmInput : Collective index for orbital and magnetic quantum number used for input function
! lmOutput : Collective index for orbital and magnetic quantum number used for output function
real :: pfac
real :: tGaunt
integer :: itype
integer :: ieqat
integer :: iatom
integer :: imesh
integer :: mqn_m
integer :: oqn_l
integer :: mqn_mpp
integer :: lmOutput
integer :: lmInput
! Local Array Variables
! *********************
! rDerFlhMtre : Real part of the radial derrivative applied to the incoming fuction coefficients
! rDerFlhMtim : Imaginary part of the radial derrivative applied to the incoming fuction coefficients
! rDerFlhMt : Radial derrivative of the incoming fuction
! r2GrFshMtNat : Expansion coefficients of the muffin-tin gradient applied to the incoming function. The coefficients are given
! in natural coordinates and multiplied by $r^2$
real, allocatable :: rDerFshMtre(:)
real, allocatable :: rDerFshMtim(:)
real, allocatable :: rDerJunk(:)
complex, allocatable :: rDerFshMt(:)
complex, allocatable :: r2GrFshMtNat(:, :, :, :)
!Matrix syntax idea from http://stackoverflow.com/questions/3708307/how-to-initialize-two-dimensional-arrays-in-fortran
complex, parameter, dimension(3, 3) :: Tmatrix = transpose(reshape([ &
& cmplx(1 / sqrt(2.), 0), cmplx(0, 0), cmplx(-1 / sqrt(2.), 0), &
& -ImagUnit / sqrt(2.), cmplx(0, 0), -ImagUnit / sqrt(2.), &
& cmplx(0, 0), cmplx(1, 0), cmplx(0, 0) &
& ], [3, 3] )) !< see my notes
! Initialization of additionaly required arrays.
allocate( r2GrFshMt(atoms%jmtd, ( atoms%lmaxd + 2 )**2, atoms%nat, 3) )
allocate( r2GrFshMtNat(atoms%jmtd, ( atoms%lmaxd + 2 )**2, atoms%nat, 3) )
allocate( rDerFshMtre(atoms%jmtd), rDerFshMtim(atoms%jmtd), rDerJunk(atoms%jmtd), rDerFshMt(atoms%jmtd) )
rDerFshMtre(:) = 0.
rDerFshMtim(:) = 0.
rDerFshMt(:) = 0.
r2GrFshMt = cmplx(0., 0.)
r2GrFshMtNat = cmplx(0., 0.)
pfac = sqrt( fpi_const / 3 )
do mqn_mpp = -1, 1
iatom = 0
do itype = 1, atoms%ntype
do ieqat = 1, atoms%neq(itype)
iatom = iatom + 1
do oqn_l = 0, atoms%lmax(itype)
do mqn_m = -oqn_l, oqn_l
! l + 1 block
if ( abs(mqn_m - mqn_mpp) <= oqn_l + 1 ) then
lmOutput = ( oqn_l + 1 ) * ( oqn_l + 2 ) + 1 + mqn_m - mqn_mpp
lmInput = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
rDerFshMtre(:) = 0.
rDerFshMtim(:) = 0.
! This is also a trade of between storage and performance, because derivative is called redundantly, maybe store it?
!CALL grdchlh(1,1,atoms%jri(itype),atoms%dx(itype),atoms%rmsh(1, itype),REAL(r2FshMt(:, lmInput, itype)),6,rDerFshMtre,rDerJunk)
!CALL grdchlh(1,1,atoms%jri(itype),atoms%dx(itype),atoms%rmsh(1, itype),AIMAG(r2FshMt(:, lmInput, itype)),6,rDerFshMtim,rDerJunk)
call Derivative( real(r2FshMt(:, lmInput, itype)), itype, atoms, rDerFshMtre )
call Derivative( aimag(r2FshMt(:, lmInput, itype)), itype, atoms, rDerFshMtim )
do imesh = 1, atoms%jri(itype)
rDerFshMt(imesh) = cmplx(rDerFshMtre(imesh), rDerFshMtim(imesh))
end do ! imesh
tGaunt = Gaunt1( oqn_l + 1, oqn_l, 1, mqn_m - mqn_mpp, mqn_m, -mqn_mpp, atoms%lmaxd + 1 )
do imesh = 1, atoms%jri(itype)
r2GrFshMtNat(imesh, lmOutput, iatom, mqn_mpp + 2) = r2GrFshMtNat(imesh, lmOutput, iatom, mqn_mpp + 2) + pfac * &
& (-1)**mqn_mpp * tGaunt &
& * (rDerFshMt(imesh) - ((oqn_l + 2)* r2FshMt(imesh, lmInput, iatom) / atoms%rmsh(imesh, itype)))
end do ! imesh
end if ! ( abs(mqn_m - mqn_mpp) <= oqn_l + 1 )
! l - 1 block
if ( abs(mqn_m - mqn_mpp) <= oqn_l - 1 ) then
lmInput = oqn_l * ( oqn_l + 1 ) + 1 + mqn_m
lmOutput = (oqn_l - 1) * oqn_l + 1 + mqn_m - mqn_mpp
rDerFshMtre(:) = 0.
rDerFshMtim(:) = 0.
! This is also a trade of between storage and performance, because derivative is called redundantly, maybe store it?
!CALL grdchlh(1,1,atoms%jri(itype),atoms%dx(itype),atoms%rmsh(1, itype),REAL(r2FshMt(:, lmInput, itype))/(atoms%rmsh(:,itype)**2),6,rDerFshMtre,rDerJunk)
!CALL grdchlh(1,1,atoms%jri(itype),atoms%dx(itype),atoms%rmsh(1, itype),AIMAG(r2FshMt(:, lmInput, itype))/(atoms%rmsh(:,itype)**2),6,rDerFshMtim,rDerJunk)
call Derivative( real(r2FshMt(:, lmInput, itype)), itype, atoms, rDerFshMtre )
call Derivative( aimag(r2FshMt(:, lmInput, itype)), itype, atoms, rDerFshMtim )
do imesh = 1, atoms%jri(itype)
rDerFshMt(imesh) = cmplx(rDerFshMtre(imesh), rDerFshMtim(imesh))
end do ! imesh
tGaunt = Gaunt1( oqn_l - 1, oqn_l, 1, mqn_m - mqn_mpp, mqn_m, -mqn_mpp, atoms%lmaxd + 1 )
do imesh = 1, atoms%jri(itype)
r2GrFshMtNat(imesh, lmOutput, iatom, mqn_mpp + 2) = r2GrFshMtNat(imesh, lmOutput, iatom, mqn_mpp + 2) + pfac * &
& (-1)**mqn_mpp * tGaunt * ( rDerFshMt(imesh) + ( (oqn_l - 1) * r2FshMt(imesh, lmInput, iatom)&
& / atoms%rmsh(imesh, itype) ) )
enddo ! imesh
end if ! ( abs(mqn_m - mqn_mpp) <= oqn_l - 1 )
end do ! mqn_m
end do ! oqn_l
end do ! ieqat
end do ! itype
end do ! mqn_mpp
! Conversion from natural to cartesian coordinates
iatom = 0
do itype = 1, atoms%ntype
do ieqat = 1, atoms%neq(itype)
iatom = iatom + 1
do oqn_l = 0, atoms%lmax(itype) + 1
do mqn_m = -oqn_l, oqn_l
lmOutput = oqn_l * (oqn_l + 1) + 1 + mqn_m
do imesh = 1, atoms%jri(itype)
r2GrFshMt(imesh, lmOutput, iatom, 1:3) = matmul( Tmatrix(1:3, 1:3), r2GrFshMtNat(imesh, lmOutput, iatom, 1:3) )
end do ! imesh
end do ! mqn_m
end do ! oqn_l
end do ! ieqat
end do ! itype
END SUBROUTINE gradYlm
SUBROUTINE divYlm(gradMtx, gradMty, gradMtz, divMt)
COMPLEX, INTENT(IN) :: gradMtx(:,:,:,:), gradMty(:,:,:,:), gradMtz(:,:,:,:) !r,lm,n,x
COMPLEX, INTENT(OUT) :: divMt(:,:,:)
!ALLOCATE (divMT(SIZE(gradMtx,1),SIZE(gradMtx,2),SIZE(gradMtx,3)))
divMT(:,:,:)=gradMtx(:,:,:,1)+gradMty(:,:,:,2)+gradMtz(:,:,:,3)
END SUBROUTINE divYlm
END MODULE m_gradYlm
optional/plot.f90
View file @ 89eed34b
......@@ -438,15 +438,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 +727,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
......
vgen/CMakeLists.txt
View file @ 89eed34b
......@@ -6,7 +6,7 @@ vgen/modcyli.f
vgen/modcylk.f
vgen/od_mkgxyz3.f
vgen/od_mkgz.f
vgen/sphcoord.f
vgen/pol_angle.f
vgen/vacp5_0.f
vgen/vacp5_z.f
vgen/visp5_0.f
......
vgen/divergence.f90
View file @ 89eed34b
......@@ -6,7 +6,7 @@
MODULE m_divergence
USE m_types
PRIVATE
PUBLIC :: mt_div, pw_div, divergence, mt_grad, pw_grad, divpotgrad
PUBLIC :: mt_div, pw_div, divergence, divergence2, mt_grad, pw_grad, divpotgrad, divpotgrad2
CONTAINS
SUBROUTINE mt_div(atoms,sphhar,sym,bxc,div)
......@@ -44,10 +44,6 @@ CONTAINS
CALL init_mt_grid(1, atoms, sphhar, .TRUE., sym, thet, phi)
DO i=1,3
bxc(i)%mt(:,1:,:,:)=0.0
END DO
DO n = 1, atoms%ntype
ALLOCATE (gradx%gr(3,atoms%jri(n)*nsp,1),grady%gr(3,atoms%jri(n)*nsp,1), &
gradz%gr(3,atoms%jri(n)*nsp,1))
......@@ -167,6 +163,62 @@ CONTAINS
END SUBROUTINE divergence
SUBROUTINE divergence2(stars,atoms,sphhar,vacuum,sym,cell,noco,bxc,div)
!--------------------------------------------------------------------------
! Use the two divergence subroutines above to now put the complete diver-
! gence of a field into a t_potden variable.
!--------------------------------------------------------------------------
USE m_lh_tofrom_lm
USE m_gradYlm
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_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: bxc
TYPE(t_potden), INTENT(INOUT) :: div
INTEGER :: i,iType,indmax
COMPLEX, ALLOCATABLE :: flm(:,:,:),grsflm1(:,:,:,:),grsflm2(:,:,:,:),grsflm3(:,:,:,:),divflm(:,:,:) ! (iR,lm,n[,x,i])
indmax=(atoms%lmaxd+1)**2
ALLOCATE(flm(atoms%jmtd,indmax,atoms%ntype))
ALLOCATE(divflm(atoms%jmtd,indmax,atoms%ntype))
DO i=1,3
DO iType=1, atoms%ntype
CALL lh_to_lm(atoms, sphhar, iType, bxc(i)%mt(:,:,iType,1), flm(:,:,iType))
END DO
IF (i==1) THEN
CALL gradYlm(atoms,flm,grsflm1)
ELSE IF (i==2) THEN
CALL gradYlm(atoms,flm,grsflm2)
ELSE
CALL gradYlm(atoms,flm,grsflm3)
END IF
END DO
DEALLOCATE(flm)
CALL divYlm(grsflm1(:,:,:,:),grsflm2(:,:,:,:),grsflm3(:,:,:,:), divflm)
DO iType=1, atoms%ntype
CALL lh_from_lm(atoms, sphhar, iType, divflm(:,1:indmax,iType), div%mt(:,0:,iType,1))
END DO
DEALLOCATE(divflm,grsflm1,grsflm2,grsflm3)
CALL pw_div(stars,sym,cell,noco,bxc,div)
END SUBROUTINE divergence2
SUBROUTINE mt_grad(n,atoms,sphhar,sym,den,gradphi)
!-----------------------------------------------------------------------------
!By use of the cartesian components of a field, its radial/angular derivati-
......@@ -232,9 +284,9 @@ CONTAINS
CALL mt_from_grid(atoms, sphhar, n, 1, grad_temp(:,:,i), gradphi(i)%mt(:,0:,n,:))
DO lh=0, lhmax
gradphi(i)%mt(:,lh,n,1) = gradphi(i)%mt(:,lh,n,1)*atoms%rmsh(:, n)**2
IF ((sphhar%llh(lh,1)/=0).AND.(sphhar%llh(lh,1)/=2)) THEN
gradphi(i)%mt(:,lh,n,1) = 0.0
END IF
!IF ((sphhar%llh(lh,1)/=0).AND.(sphhar%llh(lh,1)/=2)) THEN
! gradphi(i)%mt(:,lh,n,1) = 0.0
!END IF
END DO ! lh
END DO ! i
......@@ -318,4 +370,59 @@ CONTAINS
END SUBROUTINE divpotgrad
SUBROUTINE divpotgrad2(stars,atoms,sphhar,vacuum,sym,cell,noco,pot,grad)
USE m_types
USE m_lh_tofrom_lm
USE m_gradYlm
USE m_constants
IMPLICIT NONE
!-----------------------------------------------------------------------------
!Use the two gradient subroutines above to now put the complete gradient
!of a potential into a t_potden variable.
!-----------------------------------------------------------------------------
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_vacuum),INTENT(IN)