Commit f6ae46fd authored by Alexander Neukirchen's avatar Alexander Neukirchen

Pre merge commit. Commented out the non-funtional source-free code. Numerical...

Pre merge commit. Commented out the non-funtional source-free code. Numerical problems yet to be solved...
parent 192b38f8
......@@ -66,7 +66,6 @@ CONTAINS
USE m_ylm
USE m_metagga
USE m_plot
USE m_sfTests
#ifdef CPP_MPI
USE m_mpi_bc_potden
#endif
......@@ -513,9 +512,6 @@ CONTAINS
END IF
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)
......
......@@ -34,6 +34,7 @@ math/intgr.F90
math/ylm4.F90
math/lh_tofrom_lm.f90
math/gradYlm.f90
math/analysistests.f90
)
if (FLEUR_USE_FFTMKL)
set(fleur_F90 ${fleur_F90} math/mkl_dfti.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
......@@ -483,8 +483,16 @@ MODULE m_intgr
REAL, INTENT (OUT) :: z(jri)
INTEGER :: i
REAL :: alpha, h
z=0.0
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
......@@ -492,4 +500,89 @@ MODULE m_intgr
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
......@@ -164,13 +164,15 @@ CONTAINS
END SUBROUTINE divergence
SUBROUTINE divergence2(stars,atoms,sphhar,vacuum,sym,cell,noco,bxc,div)
USE m_lh_tofrom_lm
USE m_gradYlm
!--------------------------------------------------------------------------
! Use the two divergence subroutines above to now put the complete diver-
! gence of a field into a t_potden variable.
! Use the interstitial/vacuum divergence subroutine and an external MT-gra-
! dient routine from juPhon to assemble the divergence of a field into a
! t_potden variable. The MT-gradient is first calculated in sperical har-
! monics coefficients.
!--------------------------------------------------------------------------
USE m_lh_tofrom_lm
USE m_gradYlm
IMPLICIT NONE
......@@ -377,11 +379,15 @@ CONTAINS
USE m_gradYlm
USE m_constants
!--------------------------------------------------------------------------
! Use the interstitial/vacuum gradient subroutine and an external MT-gra-
! dient routine from juPhon to assemble the gradient of a potenital into a
! t_potden variable. The MT-gradient is first calculated in sperical har-
! monics coefficients.
!--------------------------------------------------------------------------
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
......
......@@ -209,6 +209,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
USE m_analysistests
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
......@@ -229,17 +230,20 @@ 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)
REAL :: radii(atoms%jmtd,atoms%ntype), funcsr(atoms%jmtd,atoms%ntype,3), trueder(atoms%jmtd,atoms%ntype,3), trueint(atoms%jmtd,atoms%ntype,3), &
testders(atoms%jmtd,atoms%ntype,3,5), testints(atoms%jmtd,atoms%ntype,3,5)
CALL buildfunctions(atoms,radii,funcsr,trueder,trueint)
CALL derivtest(atoms,radii,funcsr,testders)
CALL integtest(atoms,radii,funcsr,testints)
! 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,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
......
......@@ -18,6 +18,7 @@ CONTAINS
USE m_vmatgen
USE m_types
USE m_rotate_mt_den_tofrom_local
USE m_sfTests
IMPLICIT NONE
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
......@@ -80,6 +81,10 @@ CONTAINS
! ENDDO
ENDIF
ENDIF
! Source-free testwise
!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)
END SUBROUTINE vgen_finalize
END MODULE m_vgen_finalize
......@@ -33,7 +33,7 @@ contains
#include"cpp_double.h"
use m_constants
use m_types
use m_intgr, only : intgr2, intgrt
use m_intgr!, only : intgr2, intgrt, intgr5
use m_phasy1
use m_sphbes
use m_od_phasy
......@@ -170,10 +170,19 @@ 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
! Source-free testwise
!if (dosf) then
!if (l==5) THEN
! integrand_2(1:300)=integrand_2(1:300)*(atoms%rmsh(1:300,n)/atoms%rmsh(300,n))**l
!end if
!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))
!vtl(lh,n)=(0.0,0.0)
!call intgr4(integrand_1(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_1(1:imax))
!call intgr4(integrand_2(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_2(1:imax))
!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) )
!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) )
......
......@@ -94,6 +94,8 @@ CONTAINS
USE m_vgen_coulomb
USE m_gradYlm
USE m_grdchlh
USE m_sphpts
USE m_checkdop
! Takes a vectorial quantity, i.e. a t_potden variable of dimension 3, and
! makes it into a source free vector field as follows:
......@@ -125,11 +127,12 @@ CONTAINS
TYPE(t_potden) :: divloc
TYPE(t_atoms) :: atloc
INTEGER :: n, jr, lh, lhmax, jcut
INTEGER :: n, jr, lh, lhmax, jcut, nat
REAL :: xp(3,dimension%nspd)
DO i=1,3
aVec(i)%mt(:,atoms%lmaxd**2:,:,:)=0.0
END DO
!DO i=1,3
! aVec(i)%mt(:,atoms%lmaxd**2:,:,:)=0.0
!END DO
CALL div%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype, &
atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN, &
......@@ -143,6 +146,7 @@ CONTAINS
fcut=1.e-6
!div%mt(:300,25,:,1)=0.0
!DO n=1,atoms%ntype
!lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
!DO lh=0, lhmax
......@@ -167,9 +171,9 @@ CONTAINS
CALL divpotgrad2(stars,atoms,sphhar,vacuum,sym,cell,noco,phi,cvec)
DO i=1,3
cvec(i)%mt(:,atoms%lmaxd**2:,:,:)=0.0
END DO
!DO i=1,3
! cvec(i)%mt(:,atoms%lmaxd**2:,:,:)=0.0
!END DO
DO i=1,3
CALL corrB(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)
......@@ -184,15 +188,33 @@ CONTAINS
CALL divergence2(stars,atoms,sphhar,vacuum,sym,cell,noco,corrB,checkdiv)
!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
nat=1
DO n=1,atoms%ntype
lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
DO lh=0, lhmax
aVec(3)%mt(:,lh,n,1)=aVec(3)%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
div%mt(:,lh,n,1)=div%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
checkdiv%mt(:,lh,n,1)=checkdiv%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
cvec(1)%mt(:,lh,n,1)=cvec(1)%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
cvec(2)%mt(:,lh,n,1)=cvec(2)%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
cvec(3)%mt(:,lh,n,1)=cvec(3)%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
END DO
DO lh=0, lhmax
aVec(3)%mt(:,lh,n,1)=aVec(3)%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
div%mt(:,lh,n,1)=div%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
checkdiv%mt(:,lh,n,1)=checkdiv%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
cvec(1)%mt(:,lh,n,1)=cvec(1)%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
cvec(2)%mt(:,lh,n,1)=cvec(2)%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
cvec(3)%mt(:,lh,n,1)=cvec(3)%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
END DO
CALL sphpts(xp,dimension%nspd,atoms%rmt(n),atoms%pos(1,nat))
CALL checkdop(xp,dimension%nspd,n,nat,0,-1,1,dimension,atoms,sphhar,stars,sym,vacuum,cell,oneD,div)
CALL checkdop(xp,dimension%nspd,n,nat,0,-1,1,dimension,atoms,sphhar,stars,sym,vacuum,cell,oneD,phi)
nat = nat + atoms%neq(n)
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