Skip to content

GitLab

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

Sourcefree works fir arbitrary orders of l now, not only l=0 B.

parent 6e5aa0d5
Hide whitespace changes
Inline Side-by-side
Showing with 109 additions and 55 deletions
main/fleur.F90
View file @ a47f8850
......@@ -514,6 +514,7 @@ CONTAINS
END DO scfloop ! DO WHILE (l_cont)
!CALL sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,1,inDen,1.0)
CALL sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,1,vTot,2.0)
CALL add_usage_data("Iterations",iter)
......
math/gradYlm.f90
View file @ a47f8850
......@@ -194,8 +194,8 @@ CONTAINS
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 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)
......@@ -236,8 +236,6 @@ CONTAINS
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
......
math/intgr.F90
View file @ a47f8850
......@@ -13,6 +13,10 @@ 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)
......@@ -20,7 +24,8 @@ MODULE m_intgr
! intgr3: definite y rmsh(jri) h jri | z
! 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,19 @@ 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
z=0.0
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
END MODULE m_intgr
vgen/VYukawaFilm.f90
View file @ a47f8850
......@@ -110,7 +110,7 @@ module m_VYukawaFilm
! MUFFIN-TIN POTENTIAL
call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, &
call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, .FALSE., &
VYukawa%pw(:,1), den%mt(:,0:,:,1), VYukawa%potdenType, &
VYukawa%mt(:,0:,:,1) )
......@@ -951,7 +951,7 @@ module m_VYukawaFilm
! MUFFIN-TIN POTENTIAL
call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, &
call Vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, .FALSE., &
VYukawaModification%pw(:,1), den%mt(:,0:,:,1), VYukawaModification%potdenType, &
VYukawaModification%mt(:,0:,:,1) )
......
vgen/divergence.f90
View file @ a47f8850
......@@ -184,7 +184,7 @@ CONTAINS
TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: bxc
TYPE(t_potden), INTENT(INOUT) :: div
INTEGER :: i,iType,indmax
INTEGER :: i,iType,indmax, lh
COMPLEX, ALLOCATABLE :: flm(:,:,:),grsflm1(:,:,:,:),grsflm2(:,:,:,:),grsflm3(:,:,:,:),divflm(:,:,:) ! (iR,lm,n[,x,i])
indmax=(atoms%lmaxd+1)**2
......
vgen/sfTests.f90
View file @ a47f8850
......@@ -53,7 +53,7 @@ CONTAINS
END SUBROUTINE plotBtest
SUBROUTINE buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,itest,Avec,denMat,factor)
SUBROUTINE buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,itest,Avec,icut,denMat,factor)
USE m_mt_tofrom_grid
USE m_pw_tofrom_grid
USE m_xcBfield
......@@ -70,6 +70,7 @@ CONTAINS
TYPE(t_cell), INTENT(IN) :: cell
INTEGER, INTENT(IN) :: itest
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: Avec
INTEGER, INTENT(OUT) :: icut(atoms%ntype)
TYPE(t_potden), OPTIONAL, INTENT(IN) :: denMat
REAL, OPTIONAL, INTENT(IN) :: factor
......@@ -81,12 +82,14 @@ CONTAINS
REAL :: vec1(3), vec2(3), vec3(3), zero(3), point(3)
INTEGER :: grid(3)
icut=1
IF (itest.EQ.0) THEN
RETURN
END IF
IF (PRESENT(denMat)) THEN
CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denMat,factor,Avec)
CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denMat,factor,Avec,icut)
RETURN
END IF
......@@ -205,6 +208,7 @@ CONTAINS
SUBROUTINE sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,itest,denMat,factor)
USE m_xcBfield
USE m_divergence
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
......@@ -224,16 +228,22 @@ CONTAINS
TYPE(t_potden), DIMENSION(3) :: aVec, cvec, corrB
TYPE(t_potden) :: div, phi, checkdiv
INTEGER :: i, n, lh, l, icut(3,0:sphhar%nlhd,atoms%ntype)
REAL :: difftests(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
REAL :: g(atoms%jmtd)
! Test: Build a field and either compare with theoretical results or check,
! whether the sourcefree routine made it sourcefree.
IF (PRESENT(denMat)) THEN
CALL buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,1,aVec,denMat,factor)
CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,aVec,div,phi,cvec,corrB,checkdiv)
CALL buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,1,aVec,icut,denMat,factor)
DO i=1,3
aVec(i)%mt(:,atoms%lmaxd**2:,:,:)=0.0
END DO
CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,aVec,icut,div,phi,cvec,corrB,checkdiv)
CALL plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, aVec, div, phi, cvec, corrB, checkdiv)
ELSE
CALL buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,0,aVec)
CALL buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,0,aVec,icut)
RETURN
END IF
......@@ -248,6 +258,23 @@ CONTAINS
END SUBROUTINE sftest
SUBROUTINE difftester(atoms,n,l,f,g)
USE m_types
USE m_grdchlh
IMPLICIT NONE
TYPE(t_atoms), INTENT(IN) :: atoms
INTEGER, INTENT(IN) :: n,l
REAL, INTENT(IN) :: f(atoms%jri(n))
REAL, INTENT(OUT) :: g(atoms%jri(n))
REAL :: dfr(atoms%jri(n)), d2fr2(atoms%jri(n))
CALL grdchlh(1, 1, atoms%jri(n), atoms%dx(n), atoms%rmsh(1, n), f, 5, dfr, d2fr2)
g=(dfr-l*f/atoms%rmsh(:, n))!/(atoms%rmsh(1, n)**l) !must NOT be divergent towards the core
END SUBROUTINE difftester
! SUBROUTINE lhlmtest()
! 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)
......
vgen/vgen_coulomb.F90
View file @ a47f8850
......@@ -175,10 +175,6 @@ contains
call timestop("interstitial")
end if ! mpi%irank == 0
!if (dosf) then
! vCoul%pw(:,:)=0.0
!end if
! MUFFIN-TIN POTENTIAL
call timestart( "MT-spheres" )
#ifdef CPP_MPI
......@@ -186,16 +182,12 @@ contains
call MPI_BCAST( vcoul%pw, size(vcoul%pw), MPI_DOUBLE_COMPLEX, 0, mpi%mpi_comm, ierr )
CALL MPI_BARRIER(mpi%mpi_comm,ierr) !should be totally useless, but ...
#endif
call vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vCoul%pw(:,ispin), &
call vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, dosf, vCoul%pw(:,ispin), &
den%mt(:,0:,:,ispin), vCoul%potdenType, vCoul%mt(:,0:,:,ispin) )
call timestop( "MT-spheres" )
if( vCoul%potdenType == POTDEN_TYPE_POTYUK ) return
if ( mpi%irank == 0 ) then
CHECK_CONTINUITY: if ( input%vchk ) then
call timestart( "checking" )
......
vgen/vmts.F90
View file @ a47f8850
......@@ -2,7 +2,7 @@ module m_vmts
contains
subroutine vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, vpw, rho, potdenType, vr )
subroutine vmts( input, mpi, stars, sphhar, atoms, sym, cell, oneD, dosf, vpw, rho, potdenType, vr )
!-------------------------------------------------------------------------
! This subroutine calculates the lattice harmonics expansion coefficients
......@@ -33,7 +33,7 @@ contains
#include"cpp_double.h"
use m_constants
use m_types
use m_intgr, only : intgr2
use m_intgr, only : intgr2, intgrt
use m_phasy1
use m_sphbes
use m_od_phasy
......@@ -48,6 +48,7 @@ contains
type(t_sym), intent(in) :: sym
type(t_cell), intent(in) :: cell
type(t_oneD), intent(in) :: oneD
LOGICAL, INTENT(IN) :: dosf
complex, intent(in) :: vpw(:)!(stars%ng3,input%jspins)
real, intent(in) :: rho(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
integer, intent(in) :: potdenType
......@@ -59,8 +60,8 @@ contains
complex :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
real :: green_factor, termsR
real :: green_1 (1:atoms%jmtd), green_2 (1:atoms%jmtd)
real :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd)
real :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)
real :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd), integrand_3 (1:atoms%jmtd)
real :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd), integral_3 (1:atoms%jmtd)
real :: sbf(0:atoms%lmaxd)
real, allocatable, dimension(:,:) :: il, kl
......@@ -135,8 +136,6 @@ contains
deallocate( c_b )
#endif
! SPHERE INTERIOR CONTRIBUTION to the coefficients calculated from the
! values of the sphere Coulomb/Yukawa potential on the sphere boundary
......@@ -171,6 +170,10 @@ contains
integrand_2(1:imax) = green_2(1:imax) * rho(1:imax,lh,n)
call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
if (dosf) then
call intgrt(integrand_1(1:imax),atoms%rmsh(:,n),imax,integral_1(1:imax))
call intgrt(integrand_2(1:imax),atoms%rmsh(:,n),imax,integral_2(1:imax))
end if
termsR = integral_2(imax) + ( vtl(lh,n) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
vr(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
+ green_2(1:imax) * integral_1(1:imax) )
......
vgen/xcBfield.f90
View file @ a47f8850
......@@ -21,7 +21,7 @@ MODULE m_xcBfield
! back into the matrix correctly.
CONTAINS
SUBROUTINE makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denmat,factor,aVec)
SUBROUTINE makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denmat,factor,aVec,icut)
! Contructs the exchange-correlation magnetic field from the total poten-
! tial matrix or the magnetic density for the density matrix. Only used for
......@@ -43,10 +43,14 @@ CONTAINS
TYPE(t_potden), INTENT(IN) :: denmat
REAL, INTENT(IN) :: factor
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: aVec
INTEGER, INTENT(OUT) :: icut(3,0:sphhar%nlhd,atoms%ntype)
TYPE(t_potden), DIMENSION(4) :: dummyDen
INTEGER :: i, itype, ir
REAL :: r2
INTEGER :: i, itype, ir, lh
REAL :: r2(atoms%jmtd), fcut
fcut=1.e-12
icut=1
! Initialize and fill a dummy density array, that takes the initial result
! of matrixsplit.
......@@ -68,13 +72,15 @@ CONTAINS
atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE., &
POTDEN_TYPE_POTTOT,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(aVec(i)%pw_w,mold=aVec(i)%pw)
aVec(i)%mt(:,:,:,:) = dummyDen(i+1)%mt(:,:,:,:)
DO itype=1,atoms%ntype
DO ir=1,atoms%jri(itype)
DO lh=0, sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:itype - 1)) + 1))
IF (factor==2.0) THEN
r2=atoms%rmsh(ir,itype)**2
r2=atoms%rmsh(:,itype)**2
END IF
aVec(i)%mt(ir,0:,itype,:) = dummyDen(i+1)%mt(ir,0:,itype,:)*r2
END DO !ir
aVec(i)%mt(:,lh,itype,1) = aVec(i)%mt(:,lh,itype,1)*r2
!WHERE (ABS(aVec(i)%mt(ir,0:,itype,:)/r2) < fcut) aVec(i)%mt(ir,0:,itype,:) = 0.0
END DO !lh
END DO !itype
aVec(i)%pw(1:,:) = dummyDen(i+1)%pw(1:,:)
aVec(i)%vacz(1:,1:,:) = dummyDen(i+1)%vacz(1:,1:,:)
......@@ -83,8 +89,10 @@ CONTAINS
END SUBROUTINE makeVectorField
SUBROUTINE sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,bxc,div,phi,cvec,corrB,checkdiv)
SUBROUTINE sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,bxc,icut,div,phi,cvec,corrB,checkdiv)
USE m_vgen_coulomb
USE m_gradYlm
USE m_grdchlh
! Takes a vectorial quantity, i.e. a t_potden variable of dimension 3, and
! makes it into a source free vector field as follows:
......@@ -110,48 +118,48 @@ CONTAINS
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: bxc
INTEGER, INTENT(IN) :: icut(3,0:sphhar%nlhd,atoms%ntype)
TYPE(t_potden), INTENT(OUT) :: div, phi, checkdiv
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: cvec, corrB
TYPE(t_potden) :: divloc
TYPE(t_atoms) :: atloc
INTEGER :: n, jr, lh, lhmax
INTEGER :: n, jr, lh, lhmax, jcut
CALL div%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype, &
atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN, &
vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(div%pw_w,mold=div%pw)
!CALL divergence(stars,atoms,sphhar,vacuum,sym,cell,noco,bxc,div)
CALL divergence2(stars,atoms,sphhar,vacuum,sym,cell,noco,bxc,div)
atloc=atoms
atloc%zatom=0.0 !Local atoms variable with no charges; needed for the potential generation.
eps=1.e-10
fcut=1.e-6
!DO n=1,atoms%ntype
!lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
!DO lh=0, lhmax
!div%mt(:,lh,n,1)=div%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
!WHERE (ABS(div%mt(:,lh,n,1))<MAXVAL(ABS(div%mt(:,lh,n,1)))*fcut) div%mt(:,lh,n,1)=0.0
!END DO
!DO lh=0, lhmax
!div%mt(:,lh,n,1)=div%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
!END DO
!END DO
CALL phi%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTCOUL,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(phi%pw_w(SIZE(phi%pw,1),size(phi%pw,2)))
phi%pw_w = CMPLX(0.0,0.0)
CALL vgen_coulomb(1,mpi,dimension,oneD,input,field,vacuum,sym,stars,cell,sphhar,atloc,.TRUE.,div,phi)
!DO n=1,atoms%ntype
! lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
! DO jr = 1, atoms%jri(n)
! DO lh=0, lhmax
! !IF (ABS(phi%mt(jr,lh,n,1))<eps) THEN
! IF (lh/=1) THEN
! phi%mt(jr,lh,n,:)=0.0
! END IF
!END DO
! END DO
! END DO
DO i=1,3
CALL cvec(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(cvec(i)%pw_w,mold=cvec(i)%pw)
ENDDO
!CALL divpotgrad(stars,atoms,sphhar,vacuum,sym,cell,noco,phi,cvec)
CALL divpotgrad2(stars,atoms,sphhar,vacuum,sym,cell,noco,phi,cvec)
DO i=1,3
......@@ -167,8 +175,15 @@ CONTAINS
CALL divergence2(stars,atoms,sphhar,vacuum,sym,cell,noco,corrB,checkdiv)
!checkdiv%mt(:,2:,:,:)=0.0
!checkdiv%mt(:,0,:,:)=0.0
!DO n=1,atoms%ntype
!lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
!DO lh=0, lhmax
!checkdiv%mt(:,lh,n,1)=checkdiv%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
!END DO
!DO lh=0, lhmax
!checkdiv%mt(:,lh,n,1)=checkdiv%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
!END DO
!END DO
END SUBROUTINE sourcefree
......
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