mirror of https://gitlab.com/QEF/q-e.git
Band interpolator added
This commit is contained in:
Ivan Carnimeo
parent
ee02eb80a5
commit
6fdb134bc5
|
|
@ -5,6 +5,7 @@ include ../../make.inc
|
||||||
# location of needed modules and included files (if any)
|
# location of needed modules and included files (if any)
|
||||||
MODFLAGS= $(BASEMOD_FLAGS) \
|
MODFLAGS= $(BASEMOD_FLAGS) \
|
||||||
$(MOD_FLAG)../../PW/src \
|
$(MOD_FLAG)../../PW/src \
|
||||||
|
$(MOD_FLAG)../../Modules/ \
|
||||||
$(MOD_FLAG)../../dft-d3/
|
$(MOD_FLAG)../../dft-d3/
|
||||||
|
|
||||||
PPOBJS = \
|
PPOBJS = \
|
||||||
|
|
@ -51,6 +52,10 @@ write_p_avg.o \
|
||||||
write_proj.o \
|
write_proj.o \
|
||||||
write_io_header.o \
|
write_io_header.o \
|
||||||
write_hamiltonians.o \
|
write_hamiltonians.o \
|
||||||
|
fouriermod.o \
|
||||||
|
globalmod.o \
|
||||||
|
shepardmod.o \
|
||||||
|
fourierdiffmod.o \
|
||||||
xc_vdW_scale_mod.o # added by Yang Jiao
|
xc_vdW_scale_mod.o # added by Yang Jiao
|
||||||
|
|
||||||
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
|
PWOBJS = ../../PW/src/libpw.a ../../KS_Solvers/libks_solvers.a ../../dft-d3/libdftd3qe.a
|
||||||
|
|
@ -60,7 +65,7 @@ MODULES = $(PWOBJS) $(QEMODS)
|
||||||
|
|
||||||
TLDEPS= pwlibs
|
TLDEPS= pwlibs
|
||||||
|
|
||||||
all : tldeps open_grid.x average.x bands.x dos.x epsilon.x initial_state.x fs.x \
|
all : tldeps open_grid.x average.x bands.x band_interpolation.x dos.x epsilon.x initial_state.x fs.x \
|
||||||
plan_avg.x plotband.x plotproj.x plotrho.x pmw.x pp.x projwfc.x \
|
plan_avg.x plotband.x plotproj.x plotrho.x pmw.x pp.x projwfc.x \
|
||||||
pawplot.x sumpdos.x pw2wannier90.x pw2critic.x pw2gw.x pw2gt.x \
|
pawplot.x sumpdos.x pw2wannier90.x pw2critic.x pw2gw.x pw2gt.x \
|
||||||
wannier_ham.x wannier_plot.x molecularpdos.x \
|
wannier_ham.x wannier_plot.x molecularpdos.x \
|
||||||
|
|
@ -89,6 +94,11 @@ bands.x : bands.o libpp.a $(MODULES)
|
||||||
bands.o libpp.a $(MODULES) $(QELIBS)
|
bands.o libpp.a $(MODULES) $(QELIBS)
|
||||||
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
|
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
|
||||||
|
|
||||||
|
band_interpolation.x : band_interpolation.o libpp.a $(MODULES)
|
||||||
|
$(LD) $(LDFLAGS) -o $@ \
|
||||||
|
band_interpolation.o libpp.a $(MODULES) $(QELIBS)
|
||||||
|
- ( cd ../../bin ; ln -fs ../PP/src/$@ . )
|
||||||
|
|
||||||
dos.x : dos.o libpp.a $(MODULES)
|
dos.x : dos.o libpp.a $(MODULES)
|
||||||
$(LD) $(LDFLAGS) -o $@ \
|
$(LD) $(LDFLAGS) -o $@ \
|
||||||
dos.o libpp.a $(MODULES) $(QELIBS)
|
dos.o libpp.a $(MODULES) $(QELIBS)
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,51 @@
|
||||||
|
!
|
||||||
|
! Copyright (C) 2001-2022 Quantum ESPRESSO group
|
||||||
|
! This file is distributed under the terms of the
|
||||||
|
! GNU General Public License. See the file `License'
|
||||||
|
! in the root directory of the present distribution,
|
||||||
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
! Author: Ivan Carnimeo (September 2021)
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
program band_interpolation
|
||||||
|
use globalmod
|
||||||
|
use shepardmod
|
||||||
|
use fouriermod
|
||||||
|
use fourierdiffmod
|
||||||
|
implicit none
|
||||||
|
!
|
||||||
|
Call read_input ()
|
||||||
|
!
|
||||||
|
ek = 0.0d0
|
||||||
|
!
|
||||||
|
if(TRIM(method).eq.'shepard') then
|
||||||
|
!
|
||||||
|
Call shepard (1, Nb, Nq, q, eq, Nk, k, ek)
|
||||||
|
!
|
||||||
|
elseif(TRIM(method).eq.'shepard-sphere') then
|
||||||
|
!
|
||||||
|
Call shepard (2, Nb, Nq, q, eq, Nk, k, ek)
|
||||||
|
!
|
||||||
|
elseif(TRIM(method).eq.'fourier') then
|
||||||
|
!
|
||||||
|
Call fourier (Nb, Nq, q, eq, Nk, k, ek, Nsym, at, bg, Op)
|
||||||
|
!
|
||||||
|
elseif(TRIM(method).eq.'fourier-diff') then
|
||||||
|
!
|
||||||
|
Call fourierdiff (Nb, Nq, q, eq, Nk, k, ek, Nsym, at, bg, Op)
|
||||||
|
!
|
||||||
|
else
|
||||||
|
!
|
||||||
|
write(*,*) 'ERROR: Wrong method ', TRIM(method)
|
||||||
|
stop
|
||||||
|
!
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
Call print_bands (TRIM(method))
|
||||||
|
!
|
||||||
|
Call deallocate_global ()
|
||||||
|
!
|
||||||
|
end program
|
||||||
|
|
@ -0,0 +1,141 @@
|
||||||
|
!
|
||||||
|
! Copyright (C) 2001-2022 Quantum ESPRESSO group
|
||||||
|
! This file is distributed under the terms of the
|
||||||
|
! GNU General Public License. See the file `License'
|
||||||
|
! in the root directory of the present distribution,
|
||||||
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
! Author: Ivan Carnimeo (September 2021)
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
MODULE fourierdiffmod
|
||||||
|
use fouriermod
|
||||||
|
implicit none
|
||||||
|
CONTAINS
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine fourierdiff(Nb, Nq, q, eq, Nk, k, ek, Nsym, at, bg, Op)
|
||||||
|
!
|
||||||
|
! compute the band structure with Fourier interpolation from
|
||||||
|
! Pickett W. E., Krakauer H., Allen P. B., Phys. Rev. B, vol. 38, issue 4, page 2721, 1988
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: Nq, Nk, Nb, Nsym
|
||||||
|
real(dp), intent(in) :: q(3,Nq), k(3,Nk), eq(Nq,Nb)
|
||||||
|
real(dp), intent(in) :: at(3,3)
|
||||||
|
real(dp), intent(in) :: bg(3,3)
|
||||||
|
real(dp), intent(in) :: Op(1:3,1:3,1:Nsym)
|
||||||
|
real(dp), intent(out) :: ek(Nk,Nb)
|
||||||
|
!
|
||||||
|
! local variables
|
||||||
|
!
|
||||||
|
integer :: Na, ib, ik
|
||||||
|
complex(dp), allocatable :: fStarsOnQ(:,:) ! Star functions at uniform q-points
|
||||||
|
complex(dp), allocatable :: fStarsOnK(:,:) ! Star functions at path of k-points
|
||||||
|
complex(dp), allocatable :: matQQ(:,:) ! this is exactly H in the reference article
|
||||||
|
complex(dp), allocatable :: matKQ(:,:) ! this is an intermediate quantity S_m(q)*S_m(k)/rho(R_m) to construct J
|
||||||
|
complex(dp), allocatable :: matJ(:,:) ! this is exactly J in the reference article
|
||||||
|
complex(dp), allocatable :: ek_c(:,:), eq_c(:,:) ! complex band energies (for ZGEMM)
|
||||||
|
real(dp) :: vec(3)
|
||||||
|
complex(dp) :: fStar
|
||||||
|
!
|
||||||
|
! for matrix inversion
|
||||||
|
complex(dp), allocatable :: matQQ_(:,:)
|
||||||
|
real(dp), allocatable :: rmatQQ(:,:)
|
||||||
|
real(dp), allocatable :: rmatQQ_(:,:)
|
||||||
|
real(dp), allocatable :: rmatJ(:,:)
|
||||||
|
!
|
||||||
|
! for linear system solution
|
||||||
|
integer :: INFO, istar, NCoeff
|
||||||
|
integer, allocatable :: IPIV(:)
|
||||||
|
real(dp) :: sqrtrho
|
||||||
|
complex(dp), allocatable :: matA(:,:), matX(:,:), matB(:,:)
|
||||||
|
complex(dp), allocatable :: matS(:)
|
||||||
|
complex(dp), allocatable :: matC(:,:) ! coefficients for m= 2 ,..., NStars
|
||||||
|
complex(dp), allocatable :: matC1(:) ! coefficient for m=1 are treated separately
|
||||||
|
!
|
||||||
|
write(*,'(A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,'(A)') 'Fourier difference interpolation method'
|
||||||
|
!
|
||||||
|
Na = Nq - 1 ! dimension of the linear system
|
||||||
|
!
|
||||||
|
write(*,'(A)') 'Creating b = e(q_i) - e(q_Nq) i = 1, ... , Nq-1 '
|
||||||
|
allocate( matB(Na,Nb), matX(Na,Nb) )
|
||||||
|
do ib = 1, Nb
|
||||||
|
matB(1:Na,ib) = (One, Zero) * ( eq(1:Na,ib) - eq(Nq,ib) )
|
||||||
|
end do
|
||||||
|
matX = matB ! matX will be overwritten with the solution of Ax=b
|
||||||
|
!
|
||||||
|
write(*,'(A)') 'Creating A'
|
||||||
|
write(*,'(A)') 'Creating Star functions...'
|
||||||
|
Call find_stars(NSym, Op, at, .true.)
|
||||||
|
!
|
||||||
|
write(*,*) 'Checking Star functions periodicity...'
|
||||||
|
Call check_stars(Nq, q, NSym, Op, bg)
|
||||||
|
!
|
||||||
|
! fStarsOnQ = [S_m(q_i)-S_m(q_Nq)] / sqrt(rho_m)
|
||||||
|
write(*,*) 'Computing fStarsOnQ...'
|
||||||
|
allocate( fStarsOnQ(Na,NStars), matS(NStars) )
|
||||||
|
fStarsOnQ = (Zero, Zero)
|
||||||
|
Call compute_stars(fStarsOnQ, Na, Nq, q, NSym, Op, 2, .true., matS)
|
||||||
|
!
|
||||||
|
! matQQ = fStarsOnQ * fStarsOnQ^T = sum_m [S_m(q_i)-S_m(q_Nq)]*[S_m(q_j)-S_m(q_Nq)] / rho_m
|
||||||
|
write(*,*) 'Computing fStarsOnQ * fStarsOnQ*...'
|
||||||
|
allocate(matQQ(Na,Na), matA(Na,Na))
|
||||||
|
matQQ = (Zero, Zero)
|
||||||
|
Call ZGEMM('N', 'C', Na, Na, NStars, (One,Zero), fStarsOnQ, Na, fStarsOnQ, Na, (Zero,Zero), matQQ, Na)
|
||||||
|
matA = matQQ
|
||||||
|
!
|
||||||
|
write(*,'(A)') 'Solving Ax = b'
|
||||||
|
allocate( IPIV(Na) )
|
||||||
|
Call ZGESV(Na, Nb, matQQ, Na, IPIV, matX, Na, INFO)
|
||||||
|
deallocate(IPIV)
|
||||||
|
write(*,'(A)') 'Checking A*x - b = 0...'
|
||||||
|
Call ZGEMM('N', 'N', Na, Nb, Na, (One,Zero), matA, Na, matX, Na, -(One,Zero), matB, Na)
|
||||||
|
Call MatCheck_k('A*x - b = 0', matB, Na, Nb)
|
||||||
|
!
|
||||||
|
! C_m,ib = rho^(-1)_m sum_m=2 lambda_iq,ib * [S_m(q_i)-S_m(q_Nq)] m = 2, ... NStars
|
||||||
|
write(*,*) 'Computing coefficients...'
|
||||||
|
allocate( matC(NStars,Nb), matC1(Nb) )
|
||||||
|
matC = (Zero, Zero)
|
||||||
|
Call ZGEMM('C', 'N', NStars, Nb, Na, (One, Zero), fStarsOnQ, Na, matX, Na, (Zero, Zero), matC, NStars)
|
||||||
|
do istar = 1, NStars
|
||||||
|
sqrtrho = sqrt_rho(VecStars(:,istar))
|
||||||
|
matC(istar,1:Nb) = matC(istar,1:Nb)/sqrtrho ! now matC has the right coefficients for m = 2, ...
|
||||||
|
end do
|
||||||
|
do ib = 1, Nb
|
||||||
|
matC1(ib) = eq(Nq,ib) - dot_product(matC(:,ib),matS(:))
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
write(*,*) 'Computing bands...'
|
||||||
|
! fStarsOnK = S_m(k)
|
||||||
|
write(*,*) 'Computing fStarsOnK...'
|
||||||
|
allocate( fStarsOnK(Nk,NStars), ek_c(Nk,Nb) )
|
||||||
|
fStarsOnK = (Zero, Zero)
|
||||||
|
Call compute_stars(fStarsOnK, Nk, Nk, k, NSym, Op, 0)
|
||||||
|
ek_c = (Zero, Zero)
|
||||||
|
Call ZGEMM('N', 'N', Nk, Nb, NStars, (One, Zero), fStarsOnK, Nk, matC, NStars, (Zero, Zero), ek_c, Nk)
|
||||||
|
!
|
||||||
|
! now add the C1 coefficient
|
||||||
|
vec(1:3) = Zero
|
||||||
|
do ik = 1, Nk
|
||||||
|
fStar = star_function(0, k(1:3,ik), vec, NSym, Op)
|
||||||
|
do ib = 1, Nb
|
||||||
|
ek_c(ik, ib) = ek_c(ik, ib) + matC1(ib) * fStar
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
ek(:,:) = dble(ek_c(:,:))
|
||||||
|
!
|
||||||
|
deallocate( matX, matB, matA, matC, matC1, ek_c )
|
||||||
|
deallocate( C )
|
||||||
|
deallocate( VecStars )
|
||||||
|
deallocate( fStarsOnQ )
|
||||||
|
deallocate( fStarsOnK )
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine fourierdiff
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
END MODULE
|
||||||
|
|
@ -0,0 +1,527 @@
|
||||||
|
!
|
||||||
|
! Copyright (C) 2001-2022 Quantum ESPRESSO group
|
||||||
|
! This file is distributed under the terms of the
|
||||||
|
! GNU General Public License. See the file `License'
|
||||||
|
! in the root directory of the present distribution,
|
||||||
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
! Author: Ivan Carnimeo (September 2021)
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
MODULE fouriermod
|
||||||
|
USE kinds, ONLY : dp
|
||||||
|
implicit none
|
||||||
|
save
|
||||||
|
real(dp), parameter :: eps = 0.000010d0, Zero = 0.0d0, One = 1.0d0, Two = 2.0d0, Four = 4.0d0
|
||||||
|
real(dp), parameter :: Pi = Four*atan(One)
|
||||||
|
!
|
||||||
|
! the largest Miller index used to generate all the lattice vectors inside an outer shell
|
||||||
|
integer :: NMax
|
||||||
|
integer :: NC
|
||||||
|
real(dp), allocatable :: C(:)
|
||||||
|
!
|
||||||
|
integer :: NStars ! total number of Star functions
|
||||||
|
real(dp), allocatable :: VecStars(:,:) ! symmetry inequivalent lattice vectors generating the Star functions (one per Star)
|
||||||
|
integer :: NUser ! (Optional) number of user-given star functions
|
||||||
|
real(dp), allocatable :: VecUser(:,:) ! (Optional) user-given star functions
|
||||||
|
!
|
||||||
|
CONTAINS
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine allocate_fourier( )
|
||||||
|
implicit none
|
||||||
|
allocate ( C(NC) )
|
||||||
|
return
|
||||||
|
end subroutine
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine fourier(Nb, Nq, q, eq, Nk, k, ek, Nsym, at, bg, Op)
|
||||||
|
!
|
||||||
|
! compute the band structure with Fourier interpolation from
|
||||||
|
! D. D. Koelling, J. H. Wood, J. Comput. Phys., 67, 253-262 (1986)
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: Nq, Nk, Nb, Nsym
|
||||||
|
real(dp), intent(in) :: q(3,Nq), k(3,Nk), eq(Nq,Nb)
|
||||||
|
real(dp), intent(in) :: at(3,3)
|
||||||
|
real(dp), intent(in) :: bg(3,3)
|
||||||
|
real(dp), intent(in) :: Op(1:3,1:3,1:Nsym)
|
||||||
|
real(dp), intent(out) :: ek(Nk,Nb)
|
||||||
|
!
|
||||||
|
! local variables
|
||||||
|
!
|
||||||
|
complex(dp), allocatable :: fStarsOnQ(:,:) ! Star functions at uniform q-points
|
||||||
|
complex(dp), allocatable :: fStarsOnK(:,:) ! Star functions at path of k-points
|
||||||
|
complex(dp), allocatable :: matQQ(:,:) ! this is exactly H in the reference article
|
||||||
|
complex(dp), allocatable :: matKQ(:,:) ! this is an intermediate quantity S_m(q)*S_m(k)/rho(R_m) to construct J
|
||||||
|
complex(dp), allocatable :: matJ(:,:) ! this is exactly J in the reference article
|
||||||
|
complex(dp), allocatable :: ek_c(:,:), eq_c(:,:) ! complex band energies (for ZGEMM)
|
||||||
|
!
|
||||||
|
! for matrix inversion
|
||||||
|
complex(dp), allocatable :: matQQ_(:,:)
|
||||||
|
real(dp), allocatable :: rmatQQ(:,:)
|
||||||
|
real(dp), allocatable :: rmatQQ_(:,:)
|
||||||
|
real(dp), allocatable :: rmatJ(:,:)
|
||||||
|
!
|
||||||
|
write(*,'(A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,'(A)') 'Fourier interpolation method'
|
||||||
|
write(*,*) 'Creating Star functions...'
|
||||||
|
Call find_stars(NSym, Op, at)
|
||||||
|
!
|
||||||
|
write(*,*) 'Checking Star functions periodicity...'
|
||||||
|
Call check_stars(Nq, q, NSym, Op, bg)
|
||||||
|
!
|
||||||
|
! fStarsOnQ = S_m(q) / sqrt(rho_m)
|
||||||
|
write(*,*) 'Computing fStarsOnQ...'
|
||||||
|
allocate( fStarsOnQ(Nq,NStars) )
|
||||||
|
fStarsOnQ = (Zero, Zero)
|
||||||
|
Call compute_stars(fStarsOnQ, Nq, Nq, q, NSym, Op, 2)
|
||||||
|
!
|
||||||
|
!!!!civn
|
||||||
|
!!! ! matQQ = fStarsOnQ * fStarsOnQ^T = sum_i S_m(q_i) S_n(q_i)
|
||||||
|
!!! write(*,*) 'Checking ortogonality fStarsOnQ^T * fStarsOnQ...'
|
||||||
|
!!! allocate(matQQ(NStars,NStars), rmatQQ(NStars,NStars))
|
||||||
|
!!! matQQ = (Zero, Zero)
|
||||||
|
!!! Call ZGEMM('C', 'N', NStars, NStars, Nq, (One,Zero), fStarsOnQ, Nq, fStarsOnQ, Nq, (Zero,Zero), matQQ, NStars)
|
||||||
|
!!! matQQ = matQQ / dble(Nq)
|
||||||
|
!!! Call matprt_k('matQQ', NStars, NStars, matQQ)
|
||||||
|
!!! rmatQQ = dble(matQQ)
|
||||||
|
!!! Call MatCheck('rmatQQ',rmatQQ,NStars,NStars)
|
||||||
|
!!! deallocate(matQQ, rmatQQ)
|
||||||
|
!!!stop
|
||||||
|
!
|
||||||
|
! matQQ = fStarsOnQ * fStarsOnQ^T = sum_m S_m(q_i) S_m(q_j) / rho_m
|
||||||
|
write(*,*) 'Computing fStarsOnQ * fStarsOnQ^T...'
|
||||||
|
allocate(matQQ(Nq,Nq))
|
||||||
|
matQQ = (Zero, Zero)
|
||||||
|
Call ZGEMM('N', 'C', Nq, Nq, NStars, (One,Zero), fStarsOnQ, Nq, fStarsOnQ, Nq, (Zero,Zero), matQQ, Nq)
|
||||||
|
!
|
||||||
|
! matQQ --> matQQ^(-1)
|
||||||
|
write(*,*) 'Inverting fStarsOnQ * fStarsOnQ^T...'
|
||||||
|
allocate( rmatQQ(Nq,Nq), rmatQQ_(Nq,Nq), rmatJ(Nq,Nq) )
|
||||||
|
rmatQQ(:,:) = dble(matQQ(:,:))
|
||||||
|
rmatQQ_(:,:) = rmatQQ(:,:)
|
||||||
|
rmatJ(:,:) = Zero
|
||||||
|
Call MatInv('G', Nq, rmatQQ)
|
||||||
|
write(*,*) 'Checking inverse...'
|
||||||
|
Call DGEMM("N", "N", Nq, Nq, Nq, One, rmatQQ, Nq, rmatQQ_, Nq, Zero, rmatJ, Nq)
|
||||||
|
Call MatCheck('rmatJ',rmatJ,Nq,Nq)
|
||||||
|
matQQ(:,:) = (One,Zero) * rmatQQ(:,:)
|
||||||
|
deallocate( rmatQQ, rmatQQ_, rmatJ )
|
||||||
|
!
|
||||||
|
! fStarsOnK = S_m(k) / sqrt(rho_m)
|
||||||
|
write(*,*) 'Computing fStarsOnK...'
|
||||||
|
allocate( fStarsOnK(Nk,NStars) )
|
||||||
|
fStarsOnK = (Zero, Zero)
|
||||||
|
Call compute_stars(fStarsOnK, Nk, Nk, k, NSym, Op, 2)
|
||||||
|
!
|
||||||
|
! matKQ = fStarsOnK * fStarsOnQ^T
|
||||||
|
allocate(matKQ(Nk,Nq))
|
||||||
|
matKQ = (Zero, Zero)
|
||||||
|
Call ZGEMM('N', 'C', Nk, Nq, NStars, (One,Zero), fStarsOnK, Nk, fStarsOnQ, Nq, (Zero,Zero), matKQ, Nk)
|
||||||
|
!
|
||||||
|
! matJ = matKQ * matQQ^(-1)
|
||||||
|
allocate(matJ(Nk,Nq))
|
||||||
|
matJ = (Zero, Zero)
|
||||||
|
Call ZGEMM('N', 'N', Nk, Nq, Nq, (One,Zero), matKQ, Nk, matQQ, Nq, (Zero,Zero), matJ, Nk)
|
||||||
|
!
|
||||||
|
! e(k) = matJ * e(q)
|
||||||
|
allocate(ek_c(Nk,Nb), eq_c(Nq,Nb))
|
||||||
|
eq_c(:,:) = (One,Zero) * eq(:,:)
|
||||||
|
ek_c(:,:) = (Zero,Zero)
|
||||||
|
Call ZGEMM('N', 'N', Nk, Nb, Nq, (One,Zero), matJ, Nk, eq_c, Nq, (Zero,Zero), ek_c, Nk)
|
||||||
|
!
|
||||||
|
ek(:,:) = dble(ek_c(:,:))
|
||||||
|
!
|
||||||
|
deallocate( C )
|
||||||
|
deallocate( ek_c, eq_c )
|
||||||
|
deallocate( VecStars )
|
||||||
|
deallocate( matQQ )
|
||||||
|
deallocate( matKQ )
|
||||||
|
deallocate( matJ )
|
||||||
|
deallocate( fStarsOnQ )
|
||||||
|
deallocate( fStarsOnK )
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine fourier
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine find_stars(NSym, Op, at, Skip000_)
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: NSym
|
||||||
|
real(dp), intent(in) :: at(3,3)
|
||||||
|
real(dp), intent(in) :: Op(1:3,1:3,1:Nsym)
|
||||||
|
logical, optional, intent(in) :: Skip000_ ! whether to skip Miller indeces (0 0 0)
|
||||||
|
!
|
||||||
|
! local variables
|
||||||
|
logical :: Skip000
|
||||||
|
!
|
||||||
|
! arrays containing the lattice vectors inside the outer shell defined by NMax
|
||||||
|
integer :: NAll, NPrint
|
||||||
|
real(dp), allocatable :: VecAll(:,:) ! all lattice vectors
|
||||||
|
real(dp), allocatable :: VecInq(:,:) ! symmetry inequivalent lattice vectors
|
||||||
|
real(dp), allocatable :: ModAll(:) ! modules of the lattice vectors
|
||||||
|
integer, allocatable :: MapAll(:) ! permutation map of the lattice vectors
|
||||||
|
!
|
||||||
|
integer :: ii, jj, kk, isym, ivec, jvec
|
||||||
|
real(dp) :: rdiff
|
||||||
|
real(dp) :: vec_i(3), mod_i, vecOp_i(3)
|
||||||
|
real(dp) :: vec_j(3), mod_j
|
||||||
|
!
|
||||||
|
if(present(Skip000_)) then
|
||||||
|
Skip000 = Skip000_
|
||||||
|
else
|
||||||
|
Skip000 = .false.
|
||||||
|
endif
|
||||||
|
!
|
||||||
|
NAll = (2 * NMax + 1 )**3 ! from -Nmax to Nmax is (2*Nmax + 1), for 3 space directions
|
||||||
|
if(Skip000) then
|
||||||
|
NAll = NAll - 1 ! remove the (0, 0, 0) lattice vector
|
||||||
|
write(*,*) 'Skipping the (0,0,0) lattice vector...'
|
||||||
|
else
|
||||||
|
write(*,*) 'Including the (0,0,0) lattice vector...'
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
if(NUser.gt.0) then
|
||||||
|
NAll = NAll + NUser
|
||||||
|
write(*,*) "Creating ", NAll, " vectors from ", NMax, " indexes and ", NUser, " user-given vectors"
|
||||||
|
else
|
||||||
|
write(*,*) "Creating ", NAll, " vectors from ", NMax, " indexes"
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
allocate ( VecAll(3,NAll), ModAll(NAll), MapAll(NAll) )
|
||||||
|
!
|
||||||
|
ivec = 0
|
||||||
|
do ii = -NMax, NMax
|
||||||
|
do jj= -NMax, NMax
|
||||||
|
do kk= -NMax, NMax
|
||||||
|
if(Skip000.and.(ii.eq.0.and.jj.eq.0.and.kk.eq.0)) cycle
|
||||||
|
ivec = ivec + 1
|
||||||
|
VecAll(:,ivec) = dble(ii)*at(:,1) + dble(jj)*at(:,2) + dble(kk)*at(:,3)
|
||||||
|
ModAll(ivec) = sqrt( VecAll(1,ivec)**2 + VecAll(2,ivec)**2 + VecAll(3,ivec)**2 )
|
||||||
|
MapAll(ivec) = ivec
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
if(NUser.gt.0) then
|
||||||
|
do ii = 1, NUser
|
||||||
|
ivec = ivec + 1
|
||||||
|
VecAll(:,ivec) = VecUser(:,ii)
|
||||||
|
ModAll(ivec) = sqrt( VecAll(1,ivec)**2 + VecAll(2,ivec)**2 + VecAll(3,ivec)**2 )
|
||||||
|
MapAll(ivec) = ivec
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
if(ivec.ne.NAll) then
|
||||||
|
write(*,*) "ERROR: wrong number of lattice vectors for a given NMax"
|
||||||
|
write(*,*) "NMax= ",NMax," ivec=",ivec," NAll=",NAll, " NUser=",NUser
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
!
|
||||||
|
write(*,*) "Sorting vectors in shells..."
|
||||||
|
Call hpsort_eps (NAll, ModAll, MapAll, eps)
|
||||||
|
!
|
||||||
|
write(*,*) "Removing symmetry equivalent lattice vectors..."
|
||||||
|
if(NUser.gt.0) write(*,*) "WARNING: user-given vectors will be removed they are symmetry equivalent"
|
||||||
|
!
|
||||||
|
allocate( VecInq(3,NAll) )
|
||||||
|
VecInq = Zero
|
||||||
|
!
|
||||||
|
NPrint = NAll / 10
|
||||||
|
NStars = 0
|
||||||
|
do jj = 1, NAll
|
||||||
|
!
|
||||||
|
if(jj.eq.1.or.mod(jj,NPrint).eq.0) write(*,'(5X,I10,A,f12.2,A,I10,A)') &
|
||||||
|
jj-1, ' (',dble(100*(jj-1))/dble(NAll), '%) vectors analysed ... ', NStars, " Stars found"
|
||||||
|
!
|
||||||
|
jvec = MapAll(jj)
|
||||||
|
vec_j(:) = VecAll(:,jvec) ! this is in the original ordering (jvec)
|
||||||
|
mod_j = ModAll(jj) ! sqrt(dot_product(vec_j,vec_j)) ! this has been sorted by hpsort_eps (jj)
|
||||||
|
!
|
||||||
|
! this loop is needed because in principle there could be two Stars
|
||||||
|
! generated by different vectors with the same module
|
||||||
|
do ivec = 1, NStars
|
||||||
|
vec_i(:) = VecInq(:,ivec)
|
||||||
|
mod_i = sqrt(dot_product(vec_i,vec_i))
|
||||||
|
!
|
||||||
|
if(abs(mod_i-mod_j).lt.eps) then ! if the modules are different they cannot be equivalent
|
||||||
|
do isym = 1, NSym
|
||||||
|
Call applyOp(isym, Op(1:3,1:3,isym), vec_i, vecOp_i)
|
||||||
|
rdiff=sqrt( (vecOp_i(1)-vec_j(1))**2 + (vecOp_i(2)-vec_j(2))**2 + (vecOp_i(3)-vec_j(3))**2 )
|
||||||
|
if(rdiff.lt.eps) go to 10 ! vec_j is equivalent to vec_j, skip it
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
! vec_j is not equivalent to any vec_i : it generates a new Star
|
||||||
|
NStars = NStars + 1
|
||||||
|
VecInq(:,NStars) = vec_j(:)
|
||||||
|
!
|
||||||
|
10 continue
|
||||||
|
!
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
write(*,*) NStars, " Stars of symmetry inequivalent lattice vectors found..."
|
||||||
|
!
|
||||||
|
! Knowing NStars we can now allocate VecStars with the right size
|
||||||
|
! and deallocate the over-sized VecInq
|
||||||
|
!
|
||||||
|
allocate( VecStars(3,NStars) )
|
||||||
|
!
|
||||||
|
VecStars = Zero
|
||||||
|
!
|
||||||
|
VecStars(:,:) = VecInq(:,1:NStars)
|
||||||
|
!
|
||||||
|
deallocate( VecAll, ModAll, MapAll, VecInq )
|
||||||
|
if(allocated(VecUser)) deallocate(VecUser)
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine find_stars
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine check_stars(Np, p, NSym, Op, bg)
|
||||||
|
! Check whether the Star functions have the periodicity
|
||||||
|
! of the reciprocal lattice
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: Np, NSym
|
||||||
|
real(dp), intent(in) :: p(1:3,1:Np), Op(1:3,1:3,1:NSym)
|
||||||
|
real(dp), intent(in) :: bg(3,3)
|
||||||
|
real(dp), parameter :: Thr=0.0010d0
|
||||||
|
real(dp) :: vec(3), gvec(3), pvec(3)
|
||||||
|
integer :: istar, ip, ig, jg, kg, isym
|
||||||
|
complex(dp) :: fp, fpg
|
||||||
|
!
|
||||||
|
do istar = 1, NStars
|
||||||
|
vec(:) = VecStars(:,istar)
|
||||||
|
do ip = 1, Np
|
||||||
|
fp = star_function(0, p(1:3,ip), vec, NSym, Op)
|
||||||
|
do ig = -1, 1
|
||||||
|
do jg = -1, 1
|
||||||
|
do kg = -1, 1
|
||||||
|
gvec(:) = dble(ig) * bg(:,1) + dble(jg) * bg(:,2) + dble(kg) * bg(:,3)
|
||||||
|
pvec(:) = p(:,ip) + gvec(:)
|
||||||
|
fpg = star_function(0, pvec(:), vec, NSym, Op)
|
||||||
|
if(abs(fp-fpg).gt.Thr) then
|
||||||
|
if(NUser.gt.0) then
|
||||||
|
write(*,'(A)') 'WARNING: A Star function does not have reciprocal lattice periodicity'
|
||||||
|
else
|
||||||
|
write(*,'(A)') 'ERROR: A Star function does not have reciprocal lattice periodicity'
|
||||||
|
end if
|
||||||
|
write(*,'(A,I5,A,3f12.4)') 'istar: ', istar, ' vec: ', vec(:)
|
||||||
|
write(*,'(A,I5,A,3f12.4,A,2f24.12)') 'ip: ', ip, ' P-vec: ', p(:,ip), ' fp: ', fp
|
||||||
|
write(*,'(3(A,I5),A,3f12.4)') 'ig: ', ig, ' jg: ', jg, ' kg: ', kg, ' G-vec: ', gvec(:)
|
||||||
|
write(*,'(A,3f12.4,A,2f24.12)') 'P+G-vec: ', pvec(:), ' fpg: ', fpg
|
||||||
|
if(NUser.gt.0) then
|
||||||
|
go to 30
|
||||||
|
else
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
isym = NSym/3 + 1
|
||||||
|
pvec = Zero
|
||||||
|
Call applyOp(isym, Op(1:3,1:3,isym), p(1:3,ip), pvec)
|
||||||
|
fpg = star_function(0, pvec(:), vec, NSym, Op)
|
||||||
|
if(abs(fp-fpg).gt.Thr) then
|
||||||
|
if(NUser.gt.0) then
|
||||||
|
write(*,'(A)') 'WARNING: A Star function does not have reciprocal lattice periodicity'
|
||||||
|
else
|
||||||
|
write(*,'(A)') 'ERROR: A Star function does not have reciprocal lattice periodicity'
|
||||||
|
end if
|
||||||
|
write(*,'(A,I5,A,3f12.4)') 'istar: ', istar, ' vec: ', vec(:)
|
||||||
|
write(*,'(A,I5,A,3f12.4,A,2f24.12)') 'ip: ', ip, ' P-vec: ', p(:,ip), ' fp: ', fp
|
||||||
|
write(*,'(A,I5,A,3f12.4,A,2f24.12)') 'isym: ', isym, ' OpP: ', pvec(:), ' fpg: ', fpg
|
||||||
|
if(NUser.le.0) stop
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
30 continue
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine check_stars
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine compute_stars(A, LDA, Np, p, NSym, Op, ialpha, DoDiff_, S)
|
||||||
|
!
|
||||||
|
! Computes the matrix A(LDA, NStars):
|
||||||
|
!
|
||||||
|
! A_im = alpha * S_m(p_i)
|
||||||
|
!
|
||||||
|
! where:
|
||||||
|
! m (=istar) is the index of a Star vector
|
||||||
|
! p is a set of Np k-points (either q or k)
|
||||||
|
! i = 1...LDA and LDA=Np
|
||||||
|
! alpha = One ... ialpha = 0
|
||||||
|
! 1/rho(R_m) ... ialpha = 1
|
||||||
|
! 1/sqrt(rho(R_m)) ... ialpha = 2
|
||||||
|
!
|
||||||
|
! if DoDiff_ is true, then the matrix is:
|
||||||
|
!
|
||||||
|
! A_im = alpha * [ S_m(p_i) - S_m(p_Np) ]
|
||||||
|
!
|
||||||
|
! where:
|
||||||
|
! i = 1...LDA and LDA=Np-1
|
||||||
|
! and the array S_m(p_Np) is also returned as output
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
complex(dp), intent(out) :: A(LDA,NStars)
|
||||||
|
complex(dp), optional, intent(out) :: S(NStars)
|
||||||
|
integer, intent(in) :: LDA, Np, NSym, ialpha
|
||||||
|
real(dp), intent(in) :: p(1:3,1:Np), Op(1:3,1:3,1:NSym)
|
||||||
|
logical, optional, intent(in) :: DoDiff_
|
||||||
|
!
|
||||||
|
logical :: DoDiff
|
||||||
|
!
|
||||||
|
real(dp) :: vec(3) ! R_m
|
||||||
|
real(dp) :: sqrtrho ! sqrt(rho(R_m))
|
||||||
|
integer :: istar, ip
|
||||||
|
complex(dp) :: alpha
|
||||||
|
complex(dp) :: fStar ! S_m(p_i)
|
||||||
|
complex(dp) :: fStarN ! S_m(p_Np)
|
||||||
|
!
|
||||||
|
if(ialpha.lt.0.or.ialpha.gt.2) then
|
||||||
|
write(*,*) 'ERROR: wrong ialpha in compute_stars'
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
if(present(DoDiff_)) then
|
||||||
|
DoDiff = DoDiff_
|
||||||
|
if(.not.present(S)) then
|
||||||
|
write(*,*) 'ERROR: please provide S with DoDiff=.true.'
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
else
|
||||||
|
DoDiff = .false.
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
write(*,'(A,L,3(A,I5))') 'compute_stars: DoDiff: ', DoDiff, ' LDA: ', LDA, ' Np:', Np, ' NStars:', NStars
|
||||||
|
if((DoDiff.and.(LDA.ne.(Np-1))).or.(.not.DoDiff.and.(LDA.ne.Np))) then
|
||||||
|
write(*,*) 'ERROR: Wrong dimensions in compute_stars'
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
do istar = 1, NStars
|
||||||
|
vec(:) = VecStars(:,istar) ! R_m
|
||||||
|
if(ialpha.eq.0) then
|
||||||
|
alpha = (One,Zero)
|
||||||
|
elseif(ialpha.eq.1) then
|
||||||
|
sqrtrho = sqrt_rho(vec)
|
||||||
|
alpha = (One,Zero) / sqrtrho / sqrtrho
|
||||||
|
elseif(ialpha.eq.2) then
|
||||||
|
alpha = (One,Zero) / sqrt_rho(vec)
|
||||||
|
end if
|
||||||
|
fStarN = (Zero, Zero)
|
||||||
|
if(DoDiff) then
|
||||||
|
fStarN = star_function(0, p(1:3,Np), vec, NSym, Op)
|
||||||
|
S(istar) = fStarN
|
||||||
|
end if
|
||||||
|
do ip = 1, LDA
|
||||||
|
fStar = star_function(0, p(1:3,ip), vec, NSym, Op)
|
||||||
|
A(ip,istar) = alpha * (fStar-fStarN)
|
||||||
|
end do
|
||||||
|
if(abs(aimag(A(ip,istar))).gt.eps) write(*,*) "Star function: ", ip, istar, A(ip,istar), " WARNING non zero imaginary part!!"
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine compute_stars
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
complex(dp) function star_function (iprint, p, vec, NSym, Op)
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: iprint
|
||||||
|
integer, intent(in) :: NSym
|
||||||
|
real(dp), intent(in) :: p(1:3), vec(3)
|
||||||
|
real(dp), intent(in) :: Op(1:3,1:3,1:NSym)
|
||||||
|
! local variables
|
||||||
|
real(dp) :: vecOp(3)
|
||||||
|
real(dp) :: diffMod
|
||||||
|
real(dp), allocatable :: vecInq(:,:) ! vecInq(NSym,3) a list of inequivalent vectors to remove duplicates
|
||||||
|
complex(dp) :: carg, cfunc
|
||||||
|
integer :: isym, jsym, NVec
|
||||||
|
!
|
||||||
|
cfunc = (Zero, Zero)
|
||||||
|
!
|
||||||
|
allocate( vecInq(3,NSym) )
|
||||||
|
!
|
||||||
|
vecInq(:,:) = Zero
|
||||||
|
NVec = 0
|
||||||
|
!
|
||||||
|
do isym = 1, NSym
|
||||||
|
Call applyOp(isym, Op(1:3,1:3,isym), vec, vecOp)
|
||||||
|
do jsym = 1, NVec
|
||||||
|
! for vecInq=Zero diffMod never .lt.eps
|
||||||
|
diffMod = sqrt((vecOp(1)-vecInq(1,jsym))**2 + (vecOp(2)-vecInq(2,jsym))**2 + (vecOp(3)-vecInq(3,jsym))**2 )
|
||||||
|
if(iprint.gt.0) write(*,'(2(I5,3f6.2,3x))') isym, vecOp(:), jsym, vecInq(:,jsym)
|
||||||
|
if(diffMod.lt.eps) go to 20 ! vecOp already included
|
||||||
|
end do
|
||||||
|
! vecOp is new
|
||||||
|
NVec = NVec + 1
|
||||||
|
vecInq(:,NVec) = vecOp(:)
|
||||||
|
carg = (Zero,One) * Two * Pi * dot_product(vecOp,p)
|
||||||
|
cfunc = cfunc + exp(carg)
|
||||||
|
20 continue
|
||||||
|
if(iprint.gt.0) write(*,'(I5,3x,2f12.4)') NVec, cfunc
|
||||||
|
end do
|
||||||
|
carg = (One,Zero) * sqrt(dble(NVec))
|
||||||
|
star_function = cfunc/carg
|
||||||
|
!
|
||||||
|
deallocate( vecInq )
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end function star_function
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
real(dp) function sqrt_rho(vec)
|
||||||
|
! computes:
|
||||||
|
! sqrt(C0 + C1 * |vec|^2 + C2 * |vec|^4)
|
||||||
|
implicit none
|
||||||
|
real(dp), intent(in) :: vec(1:3)
|
||||||
|
real(dp) :: rmod, rho
|
||||||
|
integer :: ic, iexp
|
||||||
|
!
|
||||||
|
rho = C(1)
|
||||||
|
if(NC.gt.1) then
|
||||||
|
rmod = sqrt(dot_product(vec,vec))
|
||||||
|
do ic = 2, NC
|
||||||
|
iexp = 2*(ic - 1 )
|
||||||
|
rho = rho + C(ic) * rmod**iexp
|
||||||
|
end do
|
||||||
|
end if
|
||||||
|
sqrt_rho = sqrt(rho)
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end function sqrt_rho
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine applyOp(isym, OpMat, vec, vecOp)
|
||||||
|
! apply symmetry operation in OpMat:
|
||||||
|
! vecOp = OpMat * vec
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: isym
|
||||||
|
real(dp), intent(in) :: OpMat(1:3,1:3), vec(1:3)
|
||||||
|
real(dp), intent(out) :: vecOp(1:3)
|
||||||
|
real(dp) :: vecErr
|
||||||
|
!
|
||||||
|
vecOp(:) = matmul(OpMat, vec)
|
||||||
|
!
|
||||||
|
vecErr = abs(sqrt(dot_product(vec,vec))-sqrt(dot_product(vecOp,vecOp)))
|
||||||
|
!
|
||||||
|
if (vecErr.gt.eps) then
|
||||||
|
write(*,*) "ERROR: non-unitary symmetry operation found"
|
||||||
|
write(*,*) "isym: ", isym
|
||||||
|
write(*,*) "vec: ", vec(:)
|
||||||
|
write(*,*) "vecOp: ", vecOp(:)
|
||||||
|
write(*,*) "vecErr: ", vecErr
|
||||||
|
Call MatPrt("OpMat", 3, 3, OpMat)
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine applyOp
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
END MODULE
|
||||||
|
|
@ -0,0 +1,298 @@
|
||||||
|
!
|
||||||
|
! Copyright (C) 2001-2022 Quantum ESPRESSO group
|
||||||
|
! This file is distributed under the terms of the
|
||||||
|
! GNU General Public License. See the file `License'
|
||||||
|
! in the root directory of the present distribution,
|
||||||
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
! Authors: Ivan Carnimeo, Pietro Delugas (September 2021)
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
MODULE globalmod
|
||||||
|
USE kinds, ONLY : dp
|
||||||
|
implicit none
|
||||||
|
!INTEGER, PARAMETER :: DP = selected_real_kind(14,200)
|
||||||
|
!
|
||||||
|
! method
|
||||||
|
character(len=50) :: method
|
||||||
|
!
|
||||||
|
! band indexes
|
||||||
|
integer :: Nb
|
||||||
|
!
|
||||||
|
! uniform grid of q-points in the IBZ (cart. coord. in units 2pi/alat)
|
||||||
|
integer :: Nq, Nlines
|
||||||
|
real(dp), allocatable :: q(:,:), eq(:,:)
|
||||||
|
!
|
||||||
|
! path of k-points in the IBZ
|
||||||
|
integer :: Nk
|
||||||
|
real(dp), allocatable :: k(:,:), ek(:,:), t(:)
|
||||||
|
!
|
||||||
|
! list of k-points for building the path in the BZ (cart. coord. in units 2pi/alat)
|
||||||
|
integer :: Nkl
|
||||||
|
integer, allocatable :: kln(:)
|
||||||
|
real(dp), allocatable :: kl(:,:)
|
||||||
|
!
|
||||||
|
! Crystal data
|
||||||
|
real(dp) :: at(3,3) ! real-space lattice translation vectors (cart. coord. in units of alat) a_i(:) = at(:,i)/alat
|
||||||
|
real(dp) :: bg(3,3) ! reciprocal-lattice (cart. coord. in units 2 pi/alat) b_i(:) = bg(:,i)/tpiba
|
||||||
|
integer :: Nsym ! number of symmetry operations of the point group of the crystal
|
||||||
|
real(dp), allocatable :: Op(:,:,:) ! 3x3 symmetry matrices: Op(Nsym,3,3) (cart. units)
|
||||||
|
real(dp), allocatable :: Op_tmp(:,:,:) ! this is just a buffer to convert Op in cartesian units.
|
||||||
|
! quite redundant here, but useful to use s_axis_to_cart without modifications
|
||||||
|
!
|
||||||
|
CONTAINS
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine read_input ()
|
||||||
|
use fouriermod, ONLY : NMax, allocate_fourier, NC, C, NUser, VecUser
|
||||||
|
use shepardmod, ONLY : PMetric, ScaleSphere
|
||||||
|
!
|
||||||
|
! read the input file and make all allocations
|
||||||
|
!
|
||||||
|
use qes_read_module, only: qes_read
|
||||||
|
use qes_types_module, only: band_structure_type, atomic_structure_type, symmetries_type, basis_set_type
|
||||||
|
use fox_dom
|
||||||
|
implicit none
|
||||||
|
integer :: ik, iq, ib, ikl, ivec
|
||||||
|
integer :: ilines
|
||||||
|
integer :: isym
|
||||||
|
character(len=50) :: string
|
||||||
|
!
|
||||||
|
type (node),pointer :: doc, pnode1, pnode2
|
||||||
|
type (band_structure_type) :: bandstr
|
||||||
|
type (atomic_structure_type) :: atstr
|
||||||
|
type (symmetries_type) :: symstr
|
||||||
|
type (basis_set_type) :: basisstr
|
||||||
|
!
|
||||||
|
! read data from the xml_file
|
||||||
|
doc => parsefile('pwscf.xml')
|
||||||
|
pnode1 => item(getElementsByTagname(doc, 'output'),0)
|
||||||
|
pnode2 => item(getElementsByTagname(pnode1, 'atomic_structure'),0)
|
||||||
|
call qes_read (pnode2, atstr)
|
||||||
|
pnode2 => item(getElementsByTagname(pnode1, 'symmetries'),0)
|
||||||
|
call qes_read (pnode2, symstr)
|
||||||
|
pnode2 => item(getElementsByTagname(pnode1, 'band_structure'),0)
|
||||||
|
call qes_read(pnode2, bandstr)
|
||||||
|
pnode2 => item(getElementsByTagname(pnode1, 'basis_set'),0)
|
||||||
|
call qes_read(pnode2, basisstr)
|
||||||
|
call destroy (doc)
|
||||||
|
nullify (doc, pnode1, pnode2)
|
||||||
|
|
||||||
|
read(*, *)
|
||||||
|
read(*, *) method
|
||||||
|
write(*,*) 'Interpolation method: ', method
|
||||||
|
if( TRIM(method).ne.'shepard'.and.TRIM(method).ne.'shepard-sphere'&
|
||||||
|
.and.TRIM(method).ne.'fourier'.and.TRIM(method).ne.'fourier-diff' ) then
|
||||||
|
write(*,*) 'ERROR: Wrong method ', method
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
! optionally add user-given star functions to the basis set
|
||||||
|
!
|
||||||
|
NUser = 0
|
||||||
|
read(*,*)
|
||||||
|
read(*,*) NUser
|
||||||
|
if(NUser.gt.0) then
|
||||||
|
write(*,*) NUser, ' user-given star functions found'
|
||||||
|
allocate( VecUser(3,NUser) )
|
||||||
|
VecUser = 0.0d0
|
||||||
|
do ivec = 1, NUser
|
||||||
|
read(*,*) VecUser(1:3,ivec)
|
||||||
|
write(*,'(3f12.6)') VecUser(1:3,ivec)
|
||||||
|
end do
|
||||||
|
elseif(NUser.eq.0) then
|
||||||
|
write(*,*) 'No user-given star functions provided'
|
||||||
|
else
|
||||||
|
write(*,*) 'ERROR: Wrong NUser'
|
||||||
|
write(*,*) ' Please provide non-negative NUser'
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
! read parameters for Shepard interpolation
|
||||||
|
!
|
||||||
|
read(*,*)
|
||||||
|
read(*,*) PMetric, ScaleSphere
|
||||||
|
write(*,*) 'PMetric: ', PMetric, 'ScaleSphere ', ScaleSphere
|
||||||
|
!
|
||||||
|
! read parameters for Fourier interpolation
|
||||||
|
!
|
||||||
|
read(*, *)
|
||||||
|
read(*, *) string, NMax
|
||||||
|
if( NMax.le.0) then
|
||||||
|
write(*,*) 'Wrong NMax: ', NMax
|
||||||
|
write(*,*) 'NMax must be greater than 0 '
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
read(*, *) string, NC
|
||||||
|
if( NC.le.0) then
|
||||||
|
write(*,*) 'Wrong NC: ', NC
|
||||||
|
write(*,*) 'NC must be greater than 0 '
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
Call allocate_fourier( )
|
||||||
|
read(*, *) string, C(1:NC)
|
||||||
|
write(*,*) NC, ' coefficients read for rho expansion: ', C(:)
|
||||||
|
!
|
||||||
|
! read system dependent quantities (crystal specifications uniform grid energies)
|
||||||
|
!
|
||||||
|
!
|
||||||
|
! read the list of Nkl special points
|
||||||
|
!
|
||||||
|
read(*, *)
|
||||||
|
read(*, *) Nkl
|
||||||
|
!
|
||||||
|
write(*,*) Nkl, ' special points read'
|
||||||
|
!
|
||||||
|
allocate( kl(3,Nkl), kln(Nkl ) )
|
||||||
|
!
|
||||||
|
do ikl = 1, Nkl
|
||||||
|
read(*, *) kl(:,ikl), kln(ikl)
|
||||||
|
write(*,'(3f12.6,I5)') kl(:,ikl), kln(ikl)
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
Nlines = Nkl - 1
|
||||||
|
Nk = sum(kln(1:Nlines)) + 1
|
||||||
|
write(*,*) 'Creating ', Nlines, ' lines connecting ', Nkl, ' special points with ', Nk, ' k-points'
|
||||||
|
!
|
||||||
|
allocate( k(3,Nk), t(Nk) )
|
||||||
|
!
|
||||||
|
Call build_kpath()
|
||||||
|
!
|
||||||
|
deallocate( kl, kln )
|
||||||
|
!
|
||||||
|
! read the uniform grid of q-points
|
||||||
|
!
|
||||||
|
Nq = bandstr%nks
|
||||||
|
Nb = bandstr%nbnd
|
||||||
|
!
|
||||||
|
write(*,*) Nq, ' points on the uniform grid, ', Nb, ' bands'
|
||||||
|
!write(*,'(A)') 'iq, q(iq, :), e(iq, :)'
|
||||||
|
!
|
||||||
|
allocate( q(3, Nq), eq(Nq, Nb), ek(Nk,Nb) )
|
||||||
|
!
|
||||||
|
do iq = 1, Nq
|
||||||
|
q(:,iq) = bandstr%ks_energies(iq)%k_point%k_point(:)
|
||||||
|
end do
|
||||||
|
do iq = 1, Nq
|
||||||
|
eq(iq,:) = bandstr%ks_energies(iq)%eigenvalues%vector(:)*27.211386245988
|
||||||
|
!write(*,'(I5,11f12.6)') iq, q(iq, :), eq(iq, :)
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
! read crystalline group specifications (direct and reciprocal vectors, symmetry operations)
|
||||||
|
!
|
||||||
|
at(1:3,1) = atstr%cell%a1 / atstr%alat
|
||||||
|
at(1:3,2) = atstr%cell%a2 / atstr%alat
|
||||||
|
at(1:3,3) = atstr%cell%a3 / atstr%alat
|
||||||
|
write(*,*) ' Crystal lattice vectors found '
|
||||||
|
write(*,*) 'Ra: ' , at(:,1)
|
||||||
|
write(*,*) 'Rb: ' , at(:,2)
|
||||||
|
write(*,*) 'Rc: ' , at(:,3)
|
||||||
|
bg(1:3,1) = basisstr%reciprocal_lattice%b1
|
||||||
|
bg(1:3,2) = basisstr%reciprocal_lattice%b2
|
||||||
|
bg(1:3,3) = basisstr%reciprocal_lattice%b3
|
||||||
|
write(*,*) ' Reciprocal lattice vectors found '
|
||||||
|
write(*,*) 'Ga: ' , bg(:,1)
|
||||||
|
write(*,*) 'Gb: ' , bg(:,2)
|
||||||
|
write(*,*) 'Gc: ' , bg(:,3)
|
||||||
|
Nsym = symstr%nsym
|
||||||
|
write(*,*) Nsym, ' symmetry operations found '
|
||||||
|
!
|
||||||
|
allocate( Op(3,3,Nsym), Op_tmp(3,3,Nsym) )
|
||||||
|
!
|
||||||
|
do isym = 1, Nsym
|
||||||
|
Op_tmp(:,:,isym) = reshape(symstr%symmetry(isym)%rotation%matrix,[3,3])
|
||||||
|
end do
|
||||||
|
Call s_axis_to_cart()
|
||||||
|
deallocate(Op_tmp)
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine read_input
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine build_kpath ()
|
||||||
|
!
|
||||||
|
! build the path of k-points connecting the Nkl special points
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
integer :: ik, iik, iline
|
||||||
|
real(dp) :: tk, dt, modt
|
||||||
|
!
|
||||||
|
k = 0.0d0
|
||||||
|
t = 0.0d0
|
||||||
|
! do the first point (t = 0)
|
||||||
|
ik = 1
|
||||||
|
t(ik) = 0.0d0
|
||||||
|
k(1, ik) = kl(1, 1)
|
||||||
|
k(2, ik) = kl(2, 1)
|
||||||
|
k(3, ik) = kl(3, 1)
|
||||||
|
! loop over lines
|
||||||
|
do iline = 1, Nlines
|
||||||
|
tk = 0.0d0
|
||||||
|
do iik = 1, kln(iline)
|
||||||
|
ik = ik + 1
|
||||||
|
modt = sqrt( (kl(1,iline+1)-kl(1,iline))**2 + (kl(2,iline+1)-kl(2,iline))**2 + (kl(3,iline+1)-kl(3,iline))**2 )
|
||||||
|
tk = tk + 1.0d0/dble(kln(iline))
|
||||||
|
dt = modt/dble(kln(iline))
|
||||||
|
t(ik) = t(ik-1) + dt
|
||||||
|
k(1,ik) = kl(1,iline) + (kl(1,iline+1) - kl(1,iline)) * tk
|
||||||
|
k(2,ik) = kl(2,iline) + (kl(2,iline+1) - kl(2,iline)) * tk
|
||||||
|
k(3,ik) = kl(3,iline) + (kl(3,iline+1) - kl(3,iline)) * tk
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine build_kpath
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine print_bands (label)
|
||||||
|
!
|
||||||
|
! print band structure
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
integer :: ik
|
||||||
|
character(len=*), intent(in) :: label
|
||||||
|
character(len=100) :: formt, filename
|
||||||
|
!
|
||||||
|
write(formt,'(A,I5,A)') '(', Nb+1 ,'f24.6)'
|
||||||
|
write(filename, '(A,A)') TRIM(label),'.dat'
|
||||||
|
!
|
||||||
|
write(*,*) 'writing band structure on ', filename
|
||||||
|
!
|
||||||
|
open(2, file=filename, status='unknown')
|
||||||
|
!
|
||||||
|
do ik = 1, Nk
|
||||||
|
write(2,formt) t(ik), ek(ik,:)
|
||||||
|
!write(*,formt) t(ik), ek(ik,:)
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
close(2)
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine print_bands
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine deallocate_global ()
|
||||||
|
implicit none
|
||||||
|
deallocate(q, eq, k, ek, t, Op)
|
||||||
|
end subroutine deallocate_global
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
!----------------------------------------------------------------------
|
||||||
|
SUBROUTINE s_axis_to_cart()
|
||||||
|
!----------------------------------------------------------------------
|
||||||
|
!! This routine transforms symmetry matrices expressed in the
|
||||||
|
!! basis of the crystal axis into rotations in cartesian axis.
|
||||||
|
!
|
||||||
|
IMPLICIT NONE
|
||||||
|
!
|
||||||
|
INTEGER :: isym
|
||||||
|
REAL(DP) :: sa(3,3), sb(3,3)
|
||||||
|
!
|
||||||
|
DO isym = 1,nsym
|
||||||
|
sa(:,:) = DBLE( Op_tmp(:,:,isym) )
|
||||||
|
sb = MATMUL( bg, sa )
|
||||||
|
Op(:,:,isym) = MATMUL( at, TRANSPOSE(sb) )
|
||||||
|
ENDDO
|
||||||
|
!
|
||||||
|
END SUBROUTINE s_axis_to_cart
|
||||||
|
END MODULE
|
||||||
|
|
@ -0,0 +1,89 @@
|
||||||
|
!
|
||||||
|
! Copyright (C) 2001-2022 Quantum ESPRESSO group
|
||||||
|
! This file is distributed under the terms of the
|
||||||
|
! GNU General Public License. See the file `License'
|
||||||
|
! in the root directory of the present distribution,
|
||||||
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
! Author: Ivan Carnimeo (September 2021)
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
!----------------------------------------------------------------------------
|
||||||
|
MODULE shepardmod
|
||||||
|
implicit none
|
||||||
|
save
|
||||||
|
!
|
||||||
|
integer, parameter :: dp = selected_real_kind(14,200)
|
||||||
|
!
|
||||||
|
integer :: PMetric ! metric for the (inverse) distance in the Shepard formula
|
||||||
|
!
|
||||||
|
real(dp) :: ScaleSphere ! scaling factor for the radius of the sphere for the modified method
|
||||||
|
!
|
||||||
|
CONTAINS
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
subroutine shepard(iwhat, Nb, Nq, q, eq, Nk, k, ek)
|
||||||
|
!
|
||||||
|
! compute the band structure with Shepard's interpolation
|
||||||
|
! iwhat = 1 ... basic Shepard method
|
||||||
|
! 2 ... modified method with the sphere radius
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: iwhat
|
||||||
|
integer, intent(in) :: Nq, Nk, Nb
|
||||||
|
real(dp), intent(in) :: q(3,Nq), k(3,Nk), eq(Nq,Nb)
|
||||||
|
real(dp), intent(out) :: ek(Nk,Nb)
|
||||||
|
! local variables
|
||||||
|
real(dp) :: w, d, dsum, esum, dthr, R
|
||||||
|
integer :: ib, iq, ik
|
||||||
|
!
|
||||||
|
if (iwhat.ne.1.and.iwhat.ne.2) then
|
||||||
|
write(*,*) 'wrong iwhat in shepard'
|
||||||
|
stop
|
||||||
|
else
|
||||||
|
write(*,'(A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||||
|
write(*,'(A)') 'Shepard interpolation method'
|
||||||
|
write(*,'(4(A,I5))') 'iwhat: ',iwhat, ' Nb: ',Nb, ' Nq: ',Nq, ' Nk: ',Nk
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
PMetric = 6 ! set a metric for distance
|
||||||
|
R = ScaleSphere * sqrt( (q(1,1)-q(1,2))**2 + (q(2,1)-q(2,2))**2 + (q(3,1)-q(3,2))**2 ) ! the radius of the search sphere is proportional to the uniform grid spacing
|
||||||
|
dthr = 0.00010d0**PMetric
|
||||||
|
!
|
||||||
|
ek = 0.0d0
|
||||||
|
!
|
||||||
|
do ib = 1, Nb
|
||||||
|
do ik = 1, Nk
|
||||||
|
!
|
||||||
|
dsum = 0.0d0
|
||||||
|
esum = 0.0d0
|
||||||
|
do iq = 1, Nq
|
||||||
|
d = abs(k(1,ik)-q(1,iq))**PMetric + abs(k(2,ik)-q(2,iq))**PMetric + abs(k(3,ik)-q(3,iq))**PMetric
|
||||||
|
if(d.gt.dthr) then
|
||||||
|
if(iwhat.eq.1) then
|
||||||
|
! basic Shepard method
|
||||||
|
w = 1.0d0/d
|
||||||
|
elseif(iwhat.eq.2) then
|
||||||
|
! search only inside the sphere R
|
||||||
|
w = (max(0.0d0, (R-d))/(R*d))**2
|
||||||
|
end if
|
||||||
|
dsum = dsum + w
|
||||||
|
esum = esum + w * eq(iq,ib)
|
||||||
|
else
|
||||||
|
ek(ik,ib) = eq(iq,ib)
|
||||||
|
write(*,*) ib, ik, iq, ' found', d
|
||||||
|
go to 10
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
ek(ik,ib) = esum / dsum
|
||||||
|
!
|
||||||
|
10 continue
|
||||||
|
!
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
end subroutine shepard
|
||||||
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
END MODULE
|
||||||
|
|
@ -594,6 +594,50 @@ SUBROUTINE MatCheck( Mat, n )
|
||||||
|
|
||||||
END SUBROUTINE MatCheck
|
END SUBROUTINE MatCheck
|
||||||
!
|
!
|
||||||
|
!--------------------------------------------------------------------
|
||||||
|
SUBROUTINE MatCheck_k(label, Mat, n, m )
|
||||||
|
!--------------------------------------------------------------------
|
||||||
|
! Compute different quantities of Mat
|
||||||
|
!
|
||||||
|
USE kinds, ONLY : dp
|
||||||
|
USE io_global, ONLY : stdout
|
||||||
|
!
|
||||||
|
IMPLICIT NONE
|
||||||
|
!
|
||||||
|
INTEGER, intent(in) :: n, m
|
||||||
|
!! the dimension of the matrix
|
||||||
|
COMPLEX(dp), intent(in) :: Mat(n,m)
|
||||||
|
!! the input matrix
|
||||||
|
CHARACTER(LEN=*) :: label
|
||||||
|
!
|
||||||
|
! ... local variables
|
||||||
|
!
|
||||||
|
INTEGER :: i, j
|
||||||
|
REAL(dp) :: tmp, MaxDiag, MaxOff, SumDiag, SumOff
|
||||||
|
!
|
||||||
|
MaxDiag = 0.0d0
|
||||||
|
MaxOff = 0.0d0
|
||||||
|
SumDiag = 0.0d0
|
||||||
|
SumOff = 0.0d0
|
||||||
|
tmp = 0.0d0
|
||||||
|
do i = 1, n
|
||||||
|
do j = 1, m
|
||||||
|
tmp = sqrt(dble(Mat(i,j)*conjg(Mat(i,j))))
|
||||||
|
if(i.eq.j) then
|
||||||
|
SumDiag = SumDiag + tmp
|
||||||
|
IF(tmp.gt.MaxDiag) MaxDiag = tmp
|
||||||
|
else
|
||||||
|
SumOff = SumOff + tmp
|
||||||
|
IF(tmp.gt.MaxOff) MaxOff = tmp
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
write(stdout,'(2A,2(A,I5))') 'Matrix ', TRIM(label), ' n: ', n, ' m: ', m
|
||||||
|
write(stdout,'(2(A,f12.6))') 'MaxAbsDiag =', MaxDiag, ' SumAbsDiag =', SumDiag
|
||||||
|
write(stdout,'(2(A,f12.6))') 'MaxAbsOff =', MaxOff, ' SumAbsOff =', SumOff
|
||||||
|
|
||||||
|
END SUBROUTINE MatCheck_k
|
||||||
|
!
|
||||||
!---------------------------------------------------------------------
|
!---------------------------------------------------------------------
|
||||||
SUBROUTINE PTSVD( Mat, m )
|
SUBROUTINE PTSVD( Mat, m )
|
||||||
!-------------------------------------------------------------------
|
!-------------------------------------------------------------------
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue