Skip to content
Snippets Groups Projects
Commit e1862c73 authored by Philipp Rüssmann's avatar Philipp Rüssmann
Browse files

Addded wronskian test writeout, works only wihth NEWSOSOL and without NOSOC... 

Addded wronskian test writeout, works only wihth NEWSOSOL and without NOSOC option (use SOCSCL=0 instead)
parent 2d79f4f7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
source/KKRhost/tmat_newsolver.F90
+ 18
− 12
View file @ e1862c73
......@@ -484,7 +484,7 @@ contains
! faster calculation of RLL.
! no irregular solutions are needed in self-consistent iterations
! because the t-matrix depends only on RLL
if (.not. set_cheby_nospeedup .and. .not. (calc_exchange_couplings .or. write_pkkr_operators)) then
if (.not. set_cheby_nospeedup .and. .not. (calc_exchange_couplings .or. write_pkkr_operators) .and. .not.calc_wronskian) then
call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), rll(:,:,:,ith), tmat0(:,:), ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, &
jlk_index, hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor, '1', use_sratrick, alpha0(:,:))
else
......@@ -585,10 +585,13 @@ contains
close (9999)
end if
! Calculate additional t-matrices for Jij-tensor calculation
if (t_dtmatjij_at%calculate .or. (t_wavefunctions%isave_wavefun(i1,ie)>0 .and.&
(t_wavefunctions%save_rllleft .or. t_wavefunctions%save_sllleft)) .or. &
((write_rhoq_input .and. ie==2) .and. (i1==mu0))) then ! rhoqtest
!----------------------------------------------------------------------------
! Calculate the left-hand side solution
!----------------------------------------------------------------------------
if ( t_dtmatjij_at%calculate .or. (t_wavefunctions%isave_wavefun(i1,ie)>0 .and. &
(t_wavefunctions%save_rllleft .or. t_wavefunctions%save_sllleft)) .or. &
((write_rhoq_input .and. ie==2) .and. (i1==mu0)) .or. & ! rhoqtest
calc_wronskian ) then
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Calculate the left-hand side solution this needs to be done for the
! calculation of t-matrices for Jij tensor or if wavefunctions should be saved
......@@ -642,7 +645,7 @@ contains
! faster calculation of RLL.
! no left solutions are needed in self-consistent iterations
! because the t-matrix depends only on RLL
if (.not. set_cheby_nospeedup .and. .not. ( calc_exchange_couplings .or. write_pkkr_operators)) then
if (.not. set_cheby_nospeedup .and. .not. ( calc_exchange_couplings .or. write_pkkr_operators) .and. .not.calc_wronskian) then
! do nothing
else
call rllsll(rpan_intervall, rnew, vnspll(:,:,:,ith), rllleft(:,:,:,ith), sllleft(:,:,:,ith), tmat0, ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, &
......@@ -688,11 +691,10 @@ contains
#endif
end if ! write_rhoq_input
!----------------------------------------------------------------------------
! Calculate the left-hand side solution
!----------------------------------------------------------------------------
end if ! t_dtmatJij_at%calculate .or. t_wavefunctions%Nwfsavemax>0
!----------------------------------------------------------------------------
! Calculate the left-hand side solution
!----------------------------------------------------------------------------
! save_wavefuncions
if (t_wavefunctions%nwfsavemax>0) then
......@@ -789,6 +791,9 @@ contains
! Calculation of the Wronskian. Just for nummerical checks
!#######################################################
if ( calc_wronskian ) then
open(9999, file='test_rmesh', form='formatted')
write(9999, '(50000E26.17)') rnew(1:irmdnew)
close(9999)
if (korbit==1) then
call calcwronskian(rll(:,:,:,ith), sll(:,:,:,ith), rllleft(:,:,:,ith), sllleft(:,:,:,ith), &
ncheb, npan_tot, ipan_intervall, rpan_intervall)
......@@ -796,6 +801,7 @@ contains
call calcwronskian(rll(:,:,:,ith), sll(:,:,:,ith), rll(:,:,:,ith), sll(:,:,:,ith), &
ncheb, npan_tot, ipan_intervall, rpan_intervall)
end if
stop
end if
end do ! IE loop
......@@ -829,7 +835,7 @@ contains
ipiv,tmat0,tmatll,alpha0,dtmatll,alphall,dalphall,jlk_index,nsra,lmmaxd,nth, &
lmax,vnspll,vnspll0,vnspll1,hlk,jlk,hlk2,jlk2,tmatsph,rll,sll,rllleft,sllleft)
use :: mod_datatypes, only: dp
use :: mod_runoptions, only: calc_exchange_couplings, write_rhoq_input
use :: mod_runoptions, only: calc_exchange_couplings, write_rhoq_input, calc_wronskian
use :: mod_constants, only: czero
use :: global_variables, only: korbit
use :: mod_profiling, only: memocc
......@@ -909,7 +915,7 @@ contains
sll = czero
! Left regular and irregular wavefunctions (used here only in case of XCPL or saving of left wavefunctions)
if (calc_exchange_couplings .or. (t_wavefunctions%save_rllleft .or. t_wavefunctions%save_sllleft .or. write_rhoq_input)) then
if (calc_exchange_couplings .or. (t_wavefunctions%save_rllleft .or. t_wavefunctions%save_sllleft .or. write_rhoq_input) .or. calc_wronskian) then
allocate (rllleft(nsra*lmmaxd,lmmaxd,irmdnew,0:nth-1), stat=i_stat)
call memocc(i_stat, product(shape(rllleft))*kind(rllleft), 'RLLLEFT', 'allocate_locals_tmat_newsolver')
rllleft = czero
......
source/common/wronskian.f90
+ 45
− 27
View file @ e1862c73
......@@ -26,28 +26,28 @@ contains
use mod_cheb, only: getCLambdaCinv
implicit none
! interface
complex (kind=dp), intent(in) :: rll(:,:,:) !! regular right wavefunction
complex (kind=dp), intent(in) :: sll(:,:,:) !! irregular right wavefunction
complex (kind=dp), intent(in) :: leftrll(:,:,:) !! regular left wavefunction
complex (kind=dp), intent(in) :: leftsll(:,:,:) !! irregular left wavefunction
complex (kind=dp), intent(in) :: rll(:,:,:) !! regular right wavefunction
complex (kind=dp), intent(in) :: sll(:,:,:) !! irregular right wavefunction
complex (kind=dp), intent(in) :: leftrll(:,:,:) !! regular left wavefunction
complex (kind=dp), intent(in) :: leftsll(:,:,:) !! irregular left wavefunction
integer, intent(in) :: ncheb !! number of Chebychev polynomials
integer, intent(in) :: npan_tot !! number of panels
integer, intent(in) :: ipan_intervall(:) !! index array for potential and shapefunction boundaries for each panel
real (kind=dp), intent(in) :: rpan_intervall(:) !! radial values of panel boundaries
integer, intent(in) :: ipan_intervall(0:npan_tot) !! index array for potential and shapefunction boundaries for each panel
real (kind=dp), intent(in) :: rpan_intervall(0:npan_tot) !! radial values of panel boundaries
! local arrays
complex (kind=dp), allocatable :: drlldr(:,:,:) !! derivative of rll
complex (kind=dp), allocatable :: dslldr(:,:,:) !! derivative of sll
complex (kind=dp), allocatable :: dleftrlldr(:,:,:) !! derivative of leftrll
complex (kind=dp), allocatable :: dleftslldr(:,:,:) !! derivative of leftsll
complex (kind=dp), allocatable :: fn(:),dfndr(:) !! working arrays for radial derivatives
complex (kind=dp), allocatable :: fn2(:),dfn2dr(:) !! working arrays for radial derivatives
complex (kind=dp), allocatable :: fn3(:),dfn3dr(:) !! working arrays for radial derivatives
complex (kind=dp), allocatable :: fn4(:),dfn4dr(:) !! working arrays for radial derivatives
complex (kind=dp), allocatable :: wronskian(:,:,:) !! wronskian of the first kind (eq. 4.44, PhD Bauer)
complex (kind=dp), allocatable :: wronskian2(:,:,:) !! wronskian of the second kind (eq. 4.37, PhD Bauer)
complex (kind=dp), allocatable :: drlldr(:,:,:) !! derivative of rll
complex (kind=dp), allocatable :: dslldr(:,:,:) !! derivative of sll
complex (kind=dp), allocatable :: dleftrlldr(:,:,:) !! derivative of leftrll
complex (kind=dp), allocatable :: dleftslldr(:,:,:) !! derivative of leftsll
complex (kind=dp), allocatable :: fn(:),dfndr(:) !! working arrays for radial derivatives
complex (kind=dp), allocatable :: fn2(:),dfn2dr(:) !! working arrays for radial derivatives
complex (kind=dp), allocatable :: fn3(:),dfn3dr(:) !! working arrays for radial derivatives
complex (kind=dp), allocatable :: fn4(:),dfn4dr(:) !! working arrays for radial derivatives
complex (kind=dp), allocatable :: wronskian(:,:,:) !! wronskian of the first kind (eq. 4.44, PhD Bauer)
complex (kind=dp), allocatable :: wronskian2(:,:,:) !! wronskian of the second kind (eq. 4.37, PhD Bauer)
complex (kind=dp) :: CLambdaCinv(0:ncheb,0:ncheb) !! complex version of integration matrix in Chebychev mesh (imaginary part is zero)
real (kind=dp) :: CLambdaCinv_temp(0:ncheb,0:ncheb) !! integration matrix in Chebychev mesh
real (kind=dp) :: widthfac !! width of panel = 1/2*(r(m)-r(m-1))
real (kind=dp) :: widthfac !! width of panel = 1/2*(r(m)-r(m-1))
! working indices and matrix sizes
integer :: ilm1,ilm2,ir,ipan,irstart,irstop
integer :: lmsize,lmsize2,nrmax
......@@ -56,21 +56,21 @@ contains
lmsize2 = ubound(rll,1)
nrmax = ubound(rll,3)
write(*,*) 'in calcwronskian:', lmsize,lmsize2,nrmax
allocate( drlldr(lmsize2,lmsize,nrmax),dslldr(lmsize2,lmsize,nrmax) )
allocate( dleftrlldr(lmsize2,lmsize,nrmax),dleftslldr(lmsize2,lmsize,nrmax) )
write(*,*) 'in calcwronskian:', lmsize, lmsize2, nrmax
allocate( drlldr(lmsize2,lmsize,nrmax), dslldr(lmsize2,lmsize,nrmax) )
allocate( dleftrlldr(lmsize2,lmsize,nrmax), dleftslldr(lmsize2,lmsize,nrmax) )
allocate( wronskian(lmsize,lmsize,nrmax) )
allocate( wronskian2(lmsize2,lmsize2,nrmax) )
allocate( fn(nrmax),fn2(nrmax),dfndr(nrmax),dfn2dr(nrmax) )
allocate( fn3(nrmax),fn4(nrmax),dfn3dr(nrmax),dfn4dr(nrmax) )
allocate( fn(nrmax), fn2(nrmax), dfndr(nrmax), dfn2dr(nrmax) )
allocate( fn3(nrmax), fn4(nrmax), dfn3dr(nrmax), dfn4dr(nrmax) )
! ############################################################
! Differentiate the wave functions
! ############################################################
! construct Chebychev intetgration matrix
! construct Chebychev integration matrix
call getCLambdaCinv(Ncheb, CLambdaCinv_temp(0:ncheb,0:ncheb))
! convert to complex numbers
CLambdaCinv(0:ncheb,0:ncheb) = cmplx(CLambdaCinv_temp(0:ncheb,0:ncheb), 0.0_dp)
! perform radial derivative of rll, rllleft, sll and sllleft
......@@ -84,7 +84,7 @@ contains
do ipan=1,npan_tot
irstart=ipan_intervall(ipan-1)+1
irstop = ipan_intervall(ipan)
widthfac = 2.0D0/(rpan_intervall(ipan)-rpan_intervall(ipan-1))
widthfac = 2.0_dp/(rpan_intervall(ipan)-rpan_intervall(ipan-1))
dfndr(irstart:irstop) = matvec_dzdz(CLambdaCinv,fn(irstart:irstop))
dfndr(irstart:irstop) = dfndr(irstart:irstop)*widthfac
......@@ -115,6 +115,15 @@ contains
! open output file for wronskian relations and write files out
! the files contail the left-hand site of eq. 4.44 and 4.37
open(unit=9246760,file='test_rpan_intervall')
write(9246760,'(50000F35.16)') rpan_intervall(:)
close(unit=9246760)
open(unit=9246761,file='test_rllleft')
open(unit=9246762,file='test_sllleft')
open(unit=9246763,file='test_rll')
open(unit=9246764,file='test_sll')
open(unit=3246762,file='test_wronskian')
open(unit=3246763,file='test_wronskian2')
write(3246762,'(a)') '# lm1, lm2, wronskian(lm1, lm2, 1:ir_max) (first kind)'
......@@ -122,18 +131,27 @@ contains
do ilm1=1, lmsize
do ilm2=1, lmsize
write(3246762,'(2i5,50000E25.14)') ilm2, ilm1, wronskian(ilm2,ilm1,:)
write(3246762,'(2i5,50000F35.16)') ilm2, ilm1, wronskian(ilm2,ilm1,:)
write(9246761,'(2i5,50000F35.16)') ilm2, ilm1, leftrll(ilm2,ilm1,:)
write(9246762,'(2i5,50000F35.16)') ilm2, ilm1, leftsll(ilm2,ilm1,:)
write(9246763,'(2i5,50000F35.16)') ilm2, ilm1, rll(ilm2,ilm1,:)
write(9246764,'(2i5,50000F35.16)') ilm2, ilm1, sll(ilm2,ilm1,:)
end do
end do
do ilm1=1, lmsize2
do ilm2=1, lmsize2
write(3246763,'(2i5,50000E25.14)') ilm2, ilm1, wronskian2(ilm2,ilm1,:)
write(3246763,'(2i5,50000F35.16)') ilm2, ilm1, wronskian2(ilm2,ilm1,:)
end do
end do
close(3246762)
close(3246763)
close(9246761)
close(9246762)
close(9246763)
close(9246764)
end subroutine calcwronskian
......
utils/check_wronskian.py
+ 31
− 1
View file @ e1862c73
#!/usr/bin/env python
from numpy import loadtxt, sqrt, sum
from matplotlib.pyplot import imshow, colorbar, show, figure, subplot, axis, subplots_adjust
from matplotlib.pyplot import imshow, colorbar, show, figure, subplot, axis, subplots_adjust, title
figure(figsize=(12,6))
subplots_adjust(left=0.04, bottom=0.25, right=0.99, top=0.75)
......@@ -24,4 +24,34 @@ imshow(d, interpolation='nearest')
axis('equal')
colorbar()
figure()
subplot(2,2,1)
title('rll')
d = sum(loadtxt('test_rll')[:,2:], axis=1)
d = d.reshape(int(sqrt(len(d))),-1)
imshow(d, interpolation='nearest')
axis('equal')
colorbar()
subplot(2,2,2)
title('rllleft')
d = sum(loadtxt('test_rllleft')[:,2:], axis=1)
d = d.reshape(int(sqrt(len(d))),-1)
imshow(d, interpolation='nearest')
axis('equal')
colorbar()
subplot(2,2,3)
title('sll')
d = sum(loadtxt('test_sll')[:,2:], axis=1)
d = d.reshape(int(sqrt(len(d))),-1)
imshow(d, interpolation='nearest')
axis('equal')
colorbar()
subplot(2,2,4)
title('sllleft')
d = sum(loadtxt('test_sllleft')[:,2:], axis=1)
d = d.reshape(int(sqrt(len(d))),-1)
imshow(d, interpolation='nearest')
axis('equal')
colorbar()
show()
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