Skip to content

GitLab

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

Moved all testing into one module. Tomorrow it is on.

parent c3bb3922
Hide whitespace changes
Inline Side-by-side
Showing with 283 additions and 255 deletions
main/fleur.F90
View file @ 6e5aa0d5
......@@ -66,8 +66,7 @@ CONTAINS
USE m_ylm
USE m_metagga
USE m_plot
USE m_xcBfield
USE m_lh_tofrom_lm
USE m_sfTests
#ifdef CPP_MPI
USE m_mpi_bc_potden
#endif
......@@ -101,8 +100,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
......@@ -514,48 +512,9 @@ CONTAINS
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 makeVectorField(stars,atoms,sphhar,vacuum,input,noco,inDen,1.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)
!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(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)
CALL sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,1,vTot,2.0)
CALL add_usage_data("Iterations",iter)
......
optional/plot.f90
View file @ 6e5aa0d5
......@@ -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
!------------------------------------------------
......@@ -1049,52 +1045,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., .FALSE., 'gradPhiDiv ', divGrx, divGrx, divGry, divGrz)
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .FALSE., '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)
......
vgen/CMakeLists.txt
View file @ 6e5aa0d5
......@@ -50,6 +50,7 @@ vgen/xy_av_den.f90
vgen/VYukawaFilm.f90
vgen/divergence.f90
vgen/xcBfield.f90
vgen/sfTests.f90
)
#vdW Stuff
set(fleur_F90 ${fleur_F90}
......
vgen/sfTests.f90 0 → 100644
View file @ 6e5aa0d5
MODULE m_sfTests
CONTAINS
SUBROUTINE plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
noco, xcB, div, phi, cvec, corrB, div2)
USE m_plot
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, div2
TYPE(t_potden), DIMENSION(3), INTENT(IN) :: xcB, cvec, corrB
LOGICAL :: xsf
CALL checkplotinp()
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .FALSE., 'bInitial ', xcB(1), xcB(1), xcB(2), xcB(3))
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., .FALSE., 'gradPhiDiv ', cvec(1), cvec(1), cvec(2), cvec(3))
CALL savxsf(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, &
.FALSE., .FALSE., 'bCorrected ', corrB(1), corrB(1), corrB(2), corrB(3))
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='bInitial_f.xsf', STATUS='OLD')
CLOSE (120, STATUS="DELETE")
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 buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,itest,Avec,denMat,factor)
USE m_mt_tofrom_grid
USE m_pw_tofrom_grid
USE m_xcBfield
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_noco), INTENT(IN) :: noco
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_cell), INTENT(IN) :: cell
INTEGER, INTENT(IN) :: itest
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: Avec
TYPE(t_potden), OPTIONAL, INTENT(IN) :: denMat
REAL, OPTIONAL, INTENT(IN) :: factor
INTEGER :: nsp, n, kt, kt2, ir, i, j, k, ifftxc3, ind
REAL :: r, th, ph, x, y, z, dx, dy, dz
REAL, ALLOCATABLE :: thet(:), phi(:), A_temp(:,:,:)!space grid, index
TYPE(t_gradients) :: grad
REAL :: vec1(3), vec2(3), vec3(3), zero(3), point(3)
INTEGER :: grid(3)
IF (itest.EQ.0) THEN
RETURN
END IF
IF (PRESENT(denMat)) THEN
CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denMat,factor,Avec)
RETURN
END IF
nsp = atoms%nsp()
ifftxc3=stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
DO i=1,3
CALL Avec(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(Avec(i)%pw_w,mold=Avec(i)%pw)
Avec(i)%mt = 0
! Temporary.
Avec(i)%pw = 0
Avec(i)%pw_w = CMPLX(0.0,0.0)
Avec(i)%vacxy = 0
Avec(i)%vacz = 0
END DO
ALLOCATE (thet(atoms%nsp()),phi(atoms%nsp()))
CALL init_mt_grid(1, atoms, sphhar, .TRUE., sym, thet, phi)
CALL init_pw_grid(.TRUE.,stars,sym,cell)
!--------------------------------------------------------------------------
! Test case 1.
! In MT: radial function f(r)=r*R_MT in direction e_r.
! In interstitial: f(r_vec)=r_vec.
! So by construction, the divergence for a set of atoms with R_MT=rmt(n)
! for every n should be an isotropic 3*R_MT.
! TODO: Add interstitial part and if-test (1---> this test).
DO n=1,atoms%ntype
ALLOCATE (A_temp(atoms%jri(n)*nsp,3,1))
kt = 0
DO ir = 1, atoms%jri(n)
r = atoms%rmsh(ir, n)
!print *, 'ir'
!print *, ir
DO k = 1, nsp
th = thet(k)
ph = phi(k)
!print *, k
!print *, th/pi_const
!print *, ph/pi_const
!(r/atoms%rmt(n))
A_temp(kt+k,1,1)=(r**2)*r**2!*SIN(th)*COS(ph)
A_temp(kt+k,2,1)=(r**2)*r**2!*SIN(th)*SIN(ph)
A_temp(kt+k,3,1)=(r**2)*r**2!*COS(th)
ENDDO ! k
kt = kt + nsp
END DO ! ir
CALL mt_from_grid(atoms, sphhar, n, 1, A_temp(:,1,:), Avec(1)%mt(:,0:,n,:))
CALL mt_from_grid(atoms, sphhar, n, 1, A_temp(:,2,:), Avec(2)%mt(:,0:,n,:))
CALL mt_from_grid(atoms, sphhar, n, 1, A_temp(:,3,:), Avec(3)%mt(:,0:,n,:))
!print *, 'A_z*r^2 3rd entry before fromto grid'
!print *, A_temp(3,3,1)
!CALL mt_to_grid(.FALSE., 1, atoms, sphhar, Avec(3)%mt(:,0:,n,:), n, grad, A_temp(:,3,:))
!kt=0
!DO ir = 1, atoms%jri(n)
! r = atoms%rmsh(ir, n)
! DO k = 1, nsp
! A_temp(kt+k,3,1)=A_temp(kt+k,3,1)*r**2
! ENDDO ! k
! kt = kt + nsp
!END DO ! ir
!print *, 'A_z*r^2 3rd entry after fromto grid'
!print *, A_temp(3,3,1)
DEALLOCATE (A_temp)
END DO ! n
ALLOCATE (A_temp(ifftxc3,3,1))
grid = (/stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft/)
vec1 = (/1.,0.,0./)
vec2 = (/0.,1.,0./)
vec3 = (/0.,0.,1./)
zero = (/0.,0.,0./)
vec1=matmul(cell%amat,vec1)
vec2=matmul(cell%amat,vec2)
vec3=matmul(cell%amat,vec3)
zero=matmul(cell%amat,zero)
DO k = 0, grid(3)-1
DO j = 0, grid(2)-1
DO i = 0, grid(1)-1
point = zero + vec1*REAL(i)/(grid(1)-1) +&
vec2*REAL(j)/(grid(2)-1) +&
vec3*REAL(k)/(grid(3)-1)
ind = k*grid(2)*grid(1) + j*grid(1) + i + 1
A_temp(ind,1,1)=SIN(i*2*pi_const/grid(1))
A_temp(ind,2,1)=SIN(j*2*pi_const/grid(2))
A_temp(ind,3,1)=SIN(k*2*pi_const/grid(3))
END DO
END DO
END DO !z-loop
DO i=1,3
CALL pw_from_grid(.TRUE.,stars,.TRUE.,A_temp(:,i,:),Avec(i)%pw,Avec(i)%pw_w)
END DO
! End test case 1.
!--------------------------------------------------------------------------
CALL finish_mt_grid
CALL finish_pw_grid()
END SUBROUTINE buildAtest
SUBROUTINE sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,itest,denMat,factor)
USE m_xcBfield
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_field), INTENT(INOUT) :: field
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
INTEGER, INTENT(IN) :: itest
TYPE(t_potden), OPTIONAL, INTENT(IN) :: denMat
REAL, OPTIONAL, INTENT(IN) :: factor
TYPE(t_potden), DIMENSION(3) :: aVec, cvec, corrB
TYPE(t_potden) :: div, phi, checkdiv
! 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 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)
RETURN
END IF
!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
END SUBROUTINE sftest
! 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)
! 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)
! END SUBROUTINE lhlmtest
END MODULE m_sfTests
vgen/vgen_coulomb.F90
View file @ 6e5aa0d5
......@@ -175,9 +175,9 @@ contains
call timestop("interstitial")
end if ! mpi%irank == 0
if (dosf) then
vCoul%pw(:,:)=0.0
end if
!if (dosf) then
! vCoul%pw(:,:)=0.0
!end if
! MUFFIN-TIN POTENTIAL
call timestart( "MT-spheres" )
......
vgen/xcBfield.f90
View file @ 6e5aa0d5
......@@ -17,8 +17,8 @@ MODULE m_xcBfield
! consistent source-free fields.
!-----------------------------------------------------------------------------
PUBLIC :: makeBxc, sourcefree, builddivtest ! TODO: Build a routine to pack A_vec
! back into the matrix correctly.
PUBLIC :: makeBxc, sourcefree ! TODO: Build a routine to pack A_vec
! back into the matrix correctly.
CONTAINS
SUBROUTINE makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denmat,factor,aVec)
......@@ -83,7 +83,7 @@ CONTAINS
END SUBROUTINE makeVectorField
SUBROUTINE sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,bxc)
SUBROUTINE sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,bxc,div,phi,cvec,corrB,checkdiv)
USE m_vgen_coulomb
! Takes a vectorial quantity, i.e. a t_potden variable of dimension 3, and
......@@ -110,24 +110,19 @@ CONTAINS
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: bxc
TYPE(t_potden), INTENT(OUT) :: div, phi, checkdiv
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: cvec, corrB
TYPE(t_atoms) :: atloc
TYPE(t_potden) :: div,div2,phi,checkdiv
TYPE(t_potden), DIMENSION(3) :: cvec, cvec2, corrB
INTEGER :: n, jr, lh, lhmax
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 div2%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(div2%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,div2)
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.
......@@ -137,7 +132,7 @@ CONTAINS
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.,div2,phi)
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))
......@@ -175,150 +170,6 @@ CONTAINS
!checkdiv%mt(:,2:,:,:)=0.0
!checkdiv%mt(:,0,:,:)=0.0
CALL plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
noco, div, phi, cvec(1), cvec(2), cvec(3), &
corrB(1), corrB(2), corrB(3), checkdiv)
END SUBROUTINE sourcefree
SUBROUTINE builddivtest(stars,atoms,sphhar,vacuum,sym,cell,itest,Avec)
USE m_mt_tofrom_grid
USE m_pw_tofrom_grid
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
INTEGER, INTENT(IN) :: itest
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: Avec
INTEGER :: nsp, n, kt, kt2, ir, i, j, k, ifftxc3, ind
REAL :: r, th, ph, x, y, z, dx, dy, dz
REAL, ALLOCATABLE :: thet(:), phi(:), A_temp(:,:,:)!space grid, index
TYPE(t_gradients) :: grad
REAL :: vec1(3), vec2(3), vec3(3), zero(3), point(3)
INTEGER :: grid(3)
IF (itest.EQ.0) THEN
RETURN
END IF
nsp = atoms%nsp()
ifftxc3=stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft
DO i=1,3
CALL Avec(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(Avec(i)%pw_w,mold=Avec(i)%pw)
Avec(i)%mt = 0
! Temporary.
Avec(i)%pw = 0
Avec(i)%pw_w = CMPLX(0.0,0.0)
Avec(i)%vacxy = 0
Avec(i)%vacz = 0
END DO
ALLOCATE (thet(atoms%nsp()),phi(atoms%nsp()))
CALL init_mt_grid(1, atoms, sphhar, .TRUE., sym, thet, phi)
CALL init_pw_grid(.TRUE.,stars,sym,cell)
!--------------------------------------------------------------------------
! Test case 1.
! In MT: radial function f(r)=r*R_MT in direction e_r.
! In interstitial: f(r_vec)=r_vec.
! So by construction, the divergence for a set of atoms with R_MT=rmt(n)
! for every n should be an isotropic 3*R_MT.
! TODO: Add interstitial part and if-test (1---> this test).
DO n=1,atoms%ntype
ALLOCATE (A_temp(atoms%jri(n)*nsp,3,1))
kt = 0
DO ir = 1, atoms%jri(n)
r = atoms%rmsh(ir, n)
!print *, 'ir'
!print *, ir
DO k = 1, nsp
th = thet(k)
ph = phi(k)
!print *, k
!print *, th/pi_const
!print *, ph/pi_const
!(r/atoms%rmt(n))
A_temp(kt+k,1,1)=(r**2)*r**2!*SIN(th)*COS(ph)
A_temp(kt+k,2,1)=(r**2)*r**2!*SIN(th)*SIN(ph)
A_temp(kt+k,3,1)=(r**2)*r**2!*COS(th)
ENDDO ! k
kt = kt + nsp
END DO ! ir
CALL mt_from_grid(atoms, sphhar, n, 1, A_temp(:,1,:), Avec(1)%mt(:,0:,n,:))
CALL mt_from_grid(atoms, sphhar, n, 1, A_temp(:,2,:), Avec(2)%mt(:,0:,n,:))
CALL mt_from_grid(atoms, sphhar, n, 1, A_temp(:,3,:), Avec(3)%mt(:,0:,n,:))
!print *, 'A_z*r^2 3rd entry before fromto grid'
!print *, A_temp(3,3,1)
!CALL mt_to_grid(.FALSE., 1, atoms, sphhar, Avec(3)%mt(:,0:,n,:), n, grad, A_temp(:,3,:))
!kt=0
!DO ir = 1, atoms%jri(n)
! r = atoms%rmsh(ir, n)
! DO k = 1, nsp
! A_temp(kt+k,3,1)=A_temp(kt+k,3,1)*r**2
! ENDDO ! k
! kt = kt + nsp
!END DO ! ir
!print *, 'A_z*r^2 3rd entry after fromto grid'
!print *, A_temp(3,3,1)
DEALLOCATE (A_temp)
END DO ! n
ALLOCATE (A_temp(ifftxc3,3,1))
grid = (/stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft/)
vec1 = (/1.,0.,0./)
vec2 = (/0.,1.,0./)
vec3 = (/0.,0.,1./)
zero = (/0.,0.,0./)
vec1=matmul(cell%amat,vec1)
vec2=matmul(cell%amat,vec2)
vec3=matmul(cell%amat,vec3)