For spectrocopy EIG, added the possibility of using Arnoldi iteration techniques
implemented in ARPACK for calculating a subset of eigenvalues rather than all of
them or just the largest one.
This is invoked by setting algorithm='arnoldi' in the constructor MSSPEC.
calc = MSSPEC(spectroscopy='EIG', algorithm='arnoldi')
This version finds NEV eigenvalues with the largest magnitudes. The value of
NEV is set in the code to be 2% of the total number. The algorithm also sets
the number of basis vectors, NCV, to be 2*NEV. Options for allowing the user
to set these values as well as other convergence criteria should be added in
later versions.
Files src/msspec/calculator.py and src/msspec/parameters.py have been modified
to accept algorithm='arnoldi'.
A new directory, src/msspec/spec/fortran/eig/ar, has been created and contains
various Fortran files implementing the Arnoldi iteration. Some files are written
in Fortran90 free-format and use some features of the Fortran 2003 standard.
Compilation requires the ARPACK library to be installed and accessible using -larpack
Compilation rules and options have been introduced/modified where necessary in the
files src/options.mk, src/msspec/spec/fortran/Makefile and
src/msspec/spec/fortran/eig_ar.mk
This version has been tested for the CaO example.
This commit is contained in:
kmdunseath
parent
f5a2b9779c
commit
deb350a520
|
|
@ -95,6 +95,7 @@ from msspec.parameters import TMatrixParameters
|
||||||
from msspec.phagen.fortran.libphagen import main as do_phagen
|
from msspec.phagen.fortran.libphagen import main as do_phagen
|
||||||
from msspec.spec.fortran import _eig_mi
|
from msspec.spec.fortran import _eig_mi
|
||||||
from msspec.spec.fortran import _eig_pw
|
from msspec.spec.fortran import _eig_pw
|
||||||
|
from msspec.spec.fortran import _eig_ar
|
||||||
from msspec.spec.fortran import _phd_mi_noso_nosp_nosym
|
from msspec.spec.fortran import _phd_mi_noso_nosp_nosym
|
||||||
from msspec.spec.fortran import _phd_se_noso_nosp_nosym
|
from msspec.spec.fortran import _phd_se_noso_nosp_nosym
|
||||||
from msspec.spec.fortran import _phd_ce_noso_nosp_nosym
|
from msspec.spec.fortran import _phd_ce_noso_nosp_nosym
|
||||||
|
|
@ -418,6 +419,8 @@ class _MSCALCULATOR(Calculator):
|
||||||
do_spec = _eig_mi.run
|
do_spec = _eig_mi.run
|
||||||
elif self.global_parameters.algorithm == 'power':
|
elif self.global_parameters.algorithm == 'power':
|
||||||
do_spec = _eig_pw.run
|
do_spec = _eig_pw.run
|
||||||
|
elif self.global_parameters.algorithm == 'arnoldi':
|
||||||
|
do_spec = _eig_ar.run
|
||||||
else:
|
else:
|
||||||
LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not "
|
LOGGER.error("\'{}\' spectroscopy with \'{}\' algorithm is not "
|
||||||
"an allowed combination.".format(self.global_parameters.spectroscopy,
|
"an allowed combination.".format(self.global_parameters.spectroscopy,
|
||||||
|
|
@ -986,8 +989,8 @@ class _EIG(_MSCALCULATOR):
|
||||||
_MSCALCULATOR.__init__(self, spectroscopy='EIG', algorithm=algorithm,
|
_MSCALCULATOR.__init__(self, spectroscopy='EIG', algorithm=algorithm,
|
||||||
polarization=polarization, dichroism=dichroism,
|
polarization=polarization, dichroism=dichroism,
|
||||||
spinpol=spinpol, folder=folder, txt=txt)
|
spinpol=spinpol, folder=folder, txt=txt)
|
||||||
if algorithm not in ('inversion', 'power'):
|
if algorithm not in ('inversion', 'power', 'arnoldi'):
|
||||||
LOGGER.error("Only the 'inversion' or the 'power' algorithms "
|
LOGGER.error("Only the 'inversion', 'power' or 'arnoldi' algorithms "
|
||||||
"are supported in EIG spectroscopy mode")
|
"are supported in EIG spectroscopy mode")
|
||||||
exit(1)
|
exit(1)
|
||||||
self.iodata = iodata.Data('EIG Simulation')
|
self.iodata = iodata.Data('EIG Simulation')
|
||||||
|
|
|
||||||
|
|
@ -770,7 +770,8 @@ class GlobalParameters(BaseParameters):
|
||||||
Parameter('algorithm', types=str, allowed_values=('expansion',
|
Parameter('algorithm', types=str, allowed_values=('expansion',
|
||||||
'inversion',
|
'inversion',
|
||||||
'correlation',
|
'correlation',
|
||||||
'power'),
|
'power',
|
||||||
|
'arnoldi'),
|
||||||
default='expansion', doc="""
|
default='expansion', doc="""
|
||||||
You can choose the algorithm used for the computation of the scattering path operator.
|
You can choose the algorithm used for the computation of the scattering path operator.
|
||||||
|
|
||||||
|
|
@ -778,6 +779,7 @@ class GlobalParameters(BaseParameters):
|
||||||
- '**expansion**', for the Rehr-Albers series expansion
|
- '**expansion**', for the Rehr-Albers series expansion
|
||||||
- '**correlation**', for the correlation-expansion algorithm
|
- '**correlation**', for the correlation-expansion algorithm
|
||||||
- '**power**', for the power method approximation scheme (only for spectroscopy='EIG')
|
- '**power**', for the power method approximation scheme (only for spectroscopy='EIG')
|
||||||
|
- '**arnoldi**', for computing multiple eigenvalues using Arnoldi iteration (only for spectroscopy='EIG')
|
||||||
|
|
||||||
The series expansion algorithm is well suited for high energy since the number of terms
|
The series expansion algorithm is well suited for high energy since the number of terms
|
||||||
required decreases as the energy increases. The exact solution is obtained by the matrix inversion
|
required decreases as the energy increases. The exact solution is obtained by the matrix inversion
|
||||||
|
|
|
||||||
|
|
@ -1,6 +1,6 @@
|
||||||
.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw comp_curve clean
|
.PHONY: all phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve clean
|
||||||
|
|
||||||
all: phd_se phd_mi phd_ce eig_mi eig_pw comp_curve
|
all: phd_se phd_mi phd_ce eig_mi eig_pw eig_ar comp_curve
|
||||||
|
|
||||||
phd_se:
|
phd_se:
|
||||||
@+$(MAKE) -f phd_se_noso_nosp_nosym.mk all
|
@+$(MAKE) -f phd_se_noso_nosp_nosym.mk all
|
||||||
|
|
@ -17,6 +17,9 @@ eig_mi:
|
||||||
eig_pw:
|
eig_pw:
|
||||||
@+$(MAKE) -f eig_pw.mk all
|
@+$(MAKE) -f eig_pw.mk all
|
||||||
|
|
||||||
|
eig_ar:
|
||||||
|
@+$(MAKE) -f eig_ar.mk all
|
||||||
|
|
||||||
comp_curve:
|
comp_curve:
|
||||||
@+$(MAKE) -f comp_curve.mk all
|
@+$(MAKE) -f comp_curve.mk all
|
||||||
|
|
||||||
|
|
@ -26,4 +29,5 @@ clean::
|
||||||
@+$(MAKE) -f phd_ce_noso_nosp_nosym.mk $@
|
@+$(MAKE) -f phd_ce_noso_nosp_nosym.mk $@
|
||||||
@+$(MAKE) -f eig_mi.mk $@
|
@+$(MAKE) -f eig_mi.mk $@
|
||||||
@+$(MAKE) -f eig_pw.mk $@
|
@+$(MAKE) -f eig_pw.mk $@
|
||||||
|
@+$(MAKE) -f eig_ar.mk $@
|
||||||
@+$(MAKE) -f comp_curve.mk $@
|
@+$(MAKE) -f comp_curve.mk $@
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,249 @@
|
||||||
|
!==============================================================================!
|
||||||
|
module arnoldi_mod
|
||||||
|
!==============================================================================!
|
||||||
|
implicit none
|
||||||
|
private
|
||||||
|
!
|
||||||
|
integer, parameter :: sp = kind(1.0)
|
||||||
|
integer, parameter :: dp = kind(1.0d0)
|
||||||
|
integer, parameter :: zp = kind((0.0d0,0.0d0))
|
||||||
|
!
|
||||||
|
complex(zp), parameter :: one = complex(1.0d0, 0.0d0)
|
||||||
|
complex(zp), parameter :: zero = complex(0.0d0, 0.0d0)
|
||||||
|
!
|
||||||
|
! Public procedures
|
||||||
|
!
|
||||||
|
public :: arnoldi_iteration
|
||||||
|
!
|
||||||
|
! Private data
|
||||||
|
!
|
||||||
|
! ARPACK's debug common block
|
||||||
|
!
|
||||||
|
integer :: logfil = 6, ndigit = -3, mcaupd = 1
|
||||||
|
integer :: mgetv0 = 0, msaupd = 0, msaup2 = 0, msaitr = 0, mseigt = 0, &
|
||||||
|
msapps = 0, msgets = 0, mseupd = 0, mnaupd = 0, mnaup2 = 0, &
|
||||||
|
mnaitr = 0, mneigh = 0, mnapps = 0, mngets = 0, mneupd = 0, &
|
||||||
|
mcaup2 = 0, mcaitr = 0, mceigh = 0, mcapps = 0, mcgets = 0, &
|
||||||
|
mceupd = 0
|
||||||
|
!
|
||||||
|
common /debug/ logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt, &
|
||||||
|
msapps, msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh, &
|
||||||
|
mnapps, mngets, mneupd, mcaupd, mcaup2, mcaitr, mceigh, &
|
||||||
|
mcapps, mcgets, mceupd
|
||||||
|
!==============================================================================!
|
||||||
|
contains
|
||||||
|
!==============================================================================!
|
||||||
|
! Public procedures
|
||||||
|
!==============================================================================!
|
||||||
|
subroutine arnoldi_iteration (ndim, a, nev, d)
|
||||||
|
!
|
||||||
|
! Use ARPACK routines to find a few eigenvalues (lambda) and corresponding
|
||||||
|
! eigenvectors (x) for the standard eigenvalue problem:
|
||||||
|
!
|
||||||
|
! A*x = lambda*x
|
||||||
|
!
|
||||||
|
! where A is a general NDIM by NDIM complex matrix
|
||||||
|
!
|
||||||
|
! This subroutine is based on an ARPACK test program adapted by Logan Boudet
|
||||||
|
! as part of his M1 project in Rennes in 2022.
|
||||||
|
!
|
||||||
|
integer, intent(in) :: ndim
|
||||||
|
complex(zp), intent(in) :: a(:,:)
|
||||||
|
integer, intent(inout) :: nev ! number of eigenvalues required
|
||||||
|
complex(zp), intent(out) :: d(:) ! vector of required eigenvalues
|
||||||
|
!
|
||||||
|
! Local data
|
||||||
|
!
|
||||||
|
character(1), parameter :: bmat = "I" ! standarg eigenvalue problem
|
||||||
|
character(2), parameter :: which = "LM" ! find NEV eigenvalues of
|
||||||
|
! ! largest magnitude
|
||||||
|
!
|
||||||
|
integer :: ido, ierr, info, j, lworkl, nconv, ncv
|
||||||
|
integer :: iparam(11)
|
||||||
|
integer :: ipntr(14)
|
||||||
|
real(dp) :: tol
|
||||||
|
complex(zp) :: sigma
|
||||||
|
logical :: rvec
|
||||||
|
!
|
||||||
|
real(dp), allocatable :: rd(:,:)
|
||||||
|
real(dp), allocatable :: rwork(:)
|
||||||
|
complex(zp), allocatable :: ax(:)
|
||||||
|
complex(zp), allocatable :: resid(:)
|
||||||
|
complex(zp), allocatable :: v(:,:)
|
||||||
|
complex(zp), allocatable :: workd(:)
|
||||||
|
complex(zp), allocatable :: workev(:)
|
||||||
|
complex(zp), allocatable :: workl(:)
|
||||||
|
logical, allocatable :: select(:)
|
||||||
|
!
|
||||||
|
! External BLAS/LAPACK functions used
|
||||||
|
!
|
||||||
|
real(dp), external :: dznrm2, dlapy2
|
||||||
|
!
|
||||||
|
!
|
||||||
|
write(6,*) "----------------- BEGIN OF ARNOLDI ITERATION -----------------"
|
||||||
|
!
|
||||||
|
! NCV is is the largest number of basis vectors that will be used in the
|
||||||
|
! Implicitly Restarted Arnoldi Process. Work per major iteration is
|
||||||
|
! proportional to NDIM*NCV*NCV.
|
||||||
|
!
|
||||||
|
! Note: we must have NCV >= NEV + 2, and preferably NCV >= 2*NEV
|
||||||
|
!
|
||||||
|
! ncv = max(ceiling(1.125*nev + 15), nev+3)
|
||||||
|
ncv = 2*nev
|
||||||
|
!
|
||||||
|
iparam(11) = 0
|
||||||
|
ipntr(14) = 0
|
||||||
|
!
|
||||||
|
! stopping criteria; machine precision is used if tol <= 0
|
||||||
|
!
|
||||||
|
tol = 0.0_dp
|
||||||
|
!
|
||||||
|
! Algorithm mode
|
||||||
|
!
|
||||||
|
iparam(1) = 1 ! exact shift strategy
|
||||||
|
iparam(3) = 300 ! maximum number of iterations
|
||||||
|
iparam(7) = 1 ! use mode1 of ZNAUPD
|
||||||
|
!
|
||||||
|
lworkl = ncv*(3*ncv + 5)
|
||||||
|
allocate(rwork(ncv))
|
||||||
|
allocate(resid(ndim), v(ndim,ncv))
|
||||||
|
allocate(workd(3*ndim), workev(2*ncv), workl(lworkl))
|
||||||
|
allocate(select(ncv))
|
||||||
|
!
|
||||||
|
! IDO is the reverse communication parameter used to determine action to be taken
|
||||||
|
! on return from ZNAUPD. Its initial value must be 0.
|
||||||
|
!
|
||||||
|
ido = 0
|
||||||
|
!
|
||||||
|
! On entry, INFO == 0 instructs ZNAUPD to use a random starting vector.
|
||||||
|
! To specify a particular starting vector, set INFO to a non-zero value.
|
||||||
|
! The required startng vector should then be supplied in RESID.
|
||||||
|
!
|
||||||
|
info = 0
|
||||||
|
!
|
||||||
|
! Main reverse communication loop
|
||||||
|
!
|
||||||
|
do
|
||||||
|
call znaupd(ido, bmat, ndim, which, nev, tol, resid, ncv, v, ndim, &
|
||||||
|
iparam, ipntr, workd, workl, lworkl, rwork, info)
|
||||||
|
|
||||||
|
if (abs(ido) /= 1) exit
|
||||||
|
!
|
||||||
|
! Matrix-vector multiplication y = A*x
|
||||||
|
! Initial vector x is stored starting at workd(ipntr(1))
|
||||||
|
! The result y should be stored in workd(ipntr(2))
|
||||||
|
!
|
||||||
|
call matvec(a, ndim, ndim, workd(ipntr(1)), workd(ipntr(2)))
|
||||||
|
!
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
if (info < 0) then
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "Error: znaupd returned info = ", info
|
||||||
|
write(6,*) "Check the documentation of znaupd for more information"
|
||||||
|
write(6,*)
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
! No fatal errors, post-process with ZNEUPD to extract computed eigenvalues.
|
||||||
|
! Eigenvectors may be also computed by setting RVEC = .TRUE.)
|
||||||
|
!
|
||||||
|
if (info == 1) then
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "Maximum number of iterations reached."
|
||||||
|
write(6,*)
|
||||||
|
else if (info == 3) then
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "No shifts could be applied during implicit Arnoldi update"
|
||||||
|
write(6,*) "Try increasing NCV."
|
||||||
|
write(6,*)
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
rvec = .false.
|
||||||
|
!
|
||||||
|
call zneupd(rvec, 'A', select, d, v, ndim, sigma, workev, bmat, ndim, &
|
||||||
|
which, nev, tol, resid, ncv, v, ndim, iparam, ipntr, workd, &
|
||||||
|
workl, lworkl, rwork, ierr)
|
||||||
|
!
|
||||||
|
if (ierr /= 0) then
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "Error: zneupd returned ierr = ", ierr
|
||||||
|
write(6,*) "Check the documentation of zneupd for more information"
|
||||||
|
write(6,*)
|
||||||
|
stop
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
! Eigenvalues are returned in the one dimensional array D and if RVEC == .TRUE.
|
||||||
|
! the corresponding eigenvectors are returned in the first NCONV == IPARAM(5)
|
||||||
|
! columns of the two dimensional array V
|
||||||
|
!
|
||||||
|
nconv = iparam(5)
|
||||||
|
!
|
||||||
|
if (rvec) then
|
||||||
|
!
|
||||||
|
! Compute the residual norm || A*x - lambda*x || for the NCONV accurately
|
||||||
|
! computed eigenvalues and eigenvectors
|
||||||
|
!
|
||||||
|
allocate(rd(ncv,3), ax(ndim))
|
||||||
|
!
|
||||||
|
do j = 1, nconv
|
||||||
|
call matvec(a, ndim, ndim, v(:,j), ax)
|
||||||
|
call zaxpy(ndim, -d(j), v(:,j), 1, ax, 1)
|
||||||
|
rd(j,1) = real(d(j), dp)
|
||||||
|
rd(j,2) = aimag(d(j))
|
||||||
|
rd(j,3) = dznrm2(ndim, ax, 1) / dlapy2(rd(j,1), rd(j,2))
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
call dmout(6, nconv, 3, rd, 2*nev, -6, &
|
||||||
|
"Ritz values (Real, Imag) and relative residuals")
|
||||||
|
!
|
||||||
|
deallocate(rd, ax)
|
||||||
|
!
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "SUMMARY"
|
||||||
|
write(6,*) "======="
|
||||||
|
write(6,*)
|
||||||
|
write(6,*) "Size of the matrix is ", ndim
|
||||||
|
write(6,*) "The number of Ritz values requested is ", nev
|
||||||
|
write(6,*) "The number of Arnoldi vectors generated (NCV) is ", ncv
|
||||||
|
write(6,*) "Portion of the spectrum: ", which
|
||||||
|
write(6,*) "The number of converged Ritz values is ", nconv
|
||||||
|
write(6,*) "The number of implicit Arnoldi update iterations taken is ", iparam(3)
|
||||||
|
write(6,*) "The number of OP*x is ", iparam(9)
|
||||||
|
write(6,*) "The convergence criterion is ", tol
|
||||||
|
write(6,*)
|
||||||
|
!
|
||||||
|
nev = nconv
|
||||||
|
!
|
||||||
|
write(6,*) "------------------ END OF ARNOLDI ITERATION ------------------"
|
||||||
|
!
|
||||||
|
deallocate(rwork)
|
||||||
|
deallocate(resid, v)
|
||||||
|
deallocate(workd, workev, workl)
|
||||||
|
deallocate(select)
|
||||||
|
!
|
||||||
|
return
|
||||||
|
end subroutine arnoldi_iteration
|
||||||
|
!==============================================================================!
|
||||||
|
! Private procedures
|
||||||
|
!==============================================================================!
|
||||||
|
subroutine matvec (a, n, lda, x, y)
|
||||||
|
!
|
||||||
|
! Compute the matrix-vector product a*x, storing the result in y
|
||||||
|
!
|
||||||
|
complex(zp), intent(in) :: a(:,:)
|
||||||
|
integer, intent(in) :: n
|
||||||
|
integer, intent(in) :: lda
|
||||||
|
complex(zp), intent(in) :: x(*)
|
||||||
|
complex(zp), intent(out) :: y(*)
|
||||||
|
!
|
||||||
|
!
|
||||||
|
call zgemv('n', n, n, one, a, lda, x, 1, zero, y, 1)
|
||||||
|
!
|
||||||
|
return
|
||||||
|
end subroutine matvec
|
||||||
|
!==============================================================================!
|
||||||
|
end module arnoldi_mod
|
||||||
|
!==============================================================================!
|
||||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,291 @@
|
||||||
|
!==============================================================================!
|
||||||
|
subroutine eig_mat_ar (je, e_kin)
|
||||||
|
!==============================================================================!
|
||||||
|
! This subroutine stores the G_o T kernel matrix and computes a subset of the
|
||||||
|
! the eigenvalues with largest magnitude using Arnoldi iteration methods from
|
||||||
|
! ARPACK
|
||||||
|
!
|
||||||
|
use dim_mod, only: n_gaunt, nl_m
|
||||||
|
use coor_mod, only: natyp, ncorr, n_prot, sym_at
|
||||||
|
use outfiles_mod, only: outfile2
|
||||||
|
use outunits_mod, only: iuo1, iuo2
|
||||||
|
use trans_mod, only: lmax, vk, tl
|
||||||
|
!
|
||||||
|
use arnoldi_mod, only: arnoldi_iteration
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
!
|
||||||
|
integer, parameter :: sp = kind(1.0)
|
||||||
|
integer, parameter :: dp = kind(1.0d0)
|
||||||
|
integer, parameter :: cp = kind((0.0,0.0))
|
||||||
|
integer, parameter :: zp = kind((0.0d0,0.0d0))
|
||||||
|
!
|
||||||
|
! Subroutine arguments
|
||||||
|
!
|
||||||
|
integer, intent(in) :: je
|
||||||
|
real(sp), intent(in) :: e_kin
|
||||||
|
!
|
||||||
|
! Local data
|
||||||
|
!
|
||||||
|
integer, parameter :: ibess = 3
|
||||||
|
integer, parameter :: nprint = 10
|
||||||
|
real(dp), parameter :: pi = acos(-1.0_dp)
|
||||||
|
real(dp), parameter :: small = 0.0001_dp
|
||||||
|
complex(zp), parameter :: ic = complex(0.0d0,1.0d0)
|
||||||
|
complex(zp), parameter :: zeroc = complex(0.0d0,0.0d0)
|
||||||
|
complex(zp), parameter :: four_pi_i = 4.0_dp*pi*ic
|
||||||
|
!
|
||||||
|
integer :: jatl, jlin, jtyp, jnum, lj, lmj, mj, nbtypj
|
||||||
|
integer :: katl, klin, ktyp, knum, lk, lmk, mk, nbtypk
|
||||||
|
integer :: j, jp, l, lin, l_max, l_min, m
|
||||||
|
integer :: ndim, n_dot, n_eig, nev, nfin, nltwo, npr, n_xmax
|
||||||
|
real(sp) :: eig, xj, yj, zj, xmax_l
|
||||||
|
real(dp) :: attkj, xkj, ykj, zkj, rkj, zdkj, krkj
|
||||||
|
complex(zp) :: expkj, sum_l, tlk
|
||||||
|
!
|
||||||
|
integer, save :: iout2, iout3
|
||||||
|
!
|
||||||
|
real(sp), allocatable :: w1(:), w2(:)
|
||||||
|
real(dp), allocatable :: gnt(:)
|
||||||
|
complex(zp), allocatable :: hl1(:), sm(:,:), ylm(:,:), w(:)
|
||||||
|
character(:), allocatable :: outfile, path
|
||||||
|
!
|
||||||
|
!
|
||||||
|
if (je == 1) then
|
||||||
|
!
|
||||||
|
! Name of second output file where eigenvalues will be written
|
||||||
|
!
|
||||||
|
n_dot = index(outfile2, '.')
|
||||||
|
outfile = outfile2(1:n_dot)//'egv'
|
||||||
|
path = outfile2(1:n_dot)//'pth'
|
||||||
|
open(newunit=iout2, file=outfile, status='unknown')
|
||||||
|
open(newunit=iout3, file=path, status='unknown')
|
||||||
|
!
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
! Construction of the multiple scattering kernel matrix G_o T.
|
||||||
|
! Elements are stored using a linear index LINJ representing (J,LJ)
|
||||||
|
!
|
||||||
|
! First compute Go T array dimension
|
||||||
|
!
|
||||||
|
jlin = 0
|
||||||
|
do jtyp = 1, n_prot
|
||||||
|
nbtypj = natyp(jtyp)
|
||||||
|
lmj = lmax(jtyp,je)
|
||||||
|
do jnum = 1, nbtypj
|
||||||
|
do lj = 0, lmj
|
||||||
|
do mj = -lj, lj
|
||||||
|
jlin = jlin + 1
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
ndim = jlin
|
||||||
|
write(6,*) "GoT matrix sm has dimension ", ndim
|
||||||
|
!
|
||||||
|
allocate(sm(ndim,ndim))
|
||||||
|
sm = zeroc
|
||||||
|
!
|
||||||
|
nltwo = 2*nl_m
|
||||||
|
allocate(ylm(0:nltwo, -nltwo:nltwo))
|
||||||
|
ylm = zeroc
|
||||||
|
!
|
||||||
|
allocate(hl1(0:nltwo))
|
||||||
|
hl1 = zeroc
|
||||||
|
!
|
||||||
|
allocate(gnt(0:n_gaunt))
|
||||||
|
gnt = 0.0_dp
|
||||||
|
!
|
||||||
|
jlin = 0
|
||||||
|
do jtyp = 1, n_prot
|
||||||
|
nbtypj = natyp(jtyp)
|
||||||
|
lmj = lmax(jtyp,je)
|
||||||
|
do jnum = 1, nbtypj
|
||||||
|
jatl = ncorr(jnum,jtyp)
|
||||||
|
xj = sym_at(1,jatl)
|
||||||
|
yj = sym_at(2,jatl)
|
||||||
|
zj = sym_at(3,jatl)
|
||||||
|
do lj = 0, lmj
|
||||||
|
do mj = -lj, lj
|
||||||
|
jlin = jlin + 1
|
||||||
|
!
|
||||||
|
klin=0
|
||||||
|
do ktyp = 1, n_prot
|
||||||
|
nbtypk = natyp(ktyp)
|
||||||
|
lmk = lmax(ktyp,je)
|
||||||
|
do knum = 1, nbtypk
|
||||||
|
katl = ncorr(knum,ktyp)
|
||||||
|
!
|
||||||
|
if (katl /= jatl) then
|
||||||
|
xkj = real(sym_at(1,katl) - xj, dp)
|
||||||
|
ykj = real(sym_at(2,katl) - yj, dp)
|
||||||
|
zkj = real(sym_at(3,katl) - zj, dp)
|
||||||
|
rkj = sqrt(xkj*xkj + ykj*ykj + zkj*zkj)
|
||||||
|
krkj = real(vk(je), dp)*rkj
|
||||||
|
attkj = exp(-aimag(cmplx(vk(je)))*rkj)
|
||||||
|
expkj = (xkj + ic*ykj)/rkj
|
||||||
|
zdkj = zkj/rkj
|
||||||
|
call sph_har2(2*nl_m, zdkj, expkj, ylm, lmj+lmk)
|
||||||
|
call besphe2(lmj+lmk+1, ibess, krkj, hl1)
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
do lk = 0,lmk
|
||||||
|
l_min = abs(lk-lj)
|
||||||
|
l_max = lk + lj
|
||||||
|
tlk = cmplx(tl(lk,1,ktyp,je))
|
||||||
|
do mk = -lk, lk
|
||||||
|
klin = klin + 1
|
||||||
|
!
|
||||||
|
sm(klin,jlin) = zeroc
|
||||||
|
sum_l = zeroc
|
||||||
|
if (katl /= jatl) then
|
||||||
|
call gaunt2(lk, mk, lj, mj, gnt)
|
||||||
|
do l = l_min, l_max, 2
|
||||||
|
m = mj - mk
|
||||||
|
if (abs(m) <= l) then
|
||||||
|
sum_l = sum_l + (ic**l)*hl1(l)*ylm(l,m)*gnt(l)
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
sum_l = sum_l*attkj*four_pi_i
|
||||||
|
else
|
||||||
|
sum_l = zeroc
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
sm(klin,jlin) = tlk*sum_l
|
||||||
|
!
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
deallocate(ylm, hl1, gnt)
|
||||||
|
!
|
||||||
|
! Compute subset of eigenvalues of SM using ARPACK
|
||||||
|
!
|
||||||
|
! NEV is the number of eigenvalues required, set to 2% of NDIM
|
||||||
|
!
|
||||||
|
nev = ceiling(0.02*ndim)
|
||||||
|
allocate(w(nev))
|
||||||
|
!
|
||||||
|
call arnoldi_iteration(ndim, sm, nev, w)
|
||||||
|
!
|
||||||
|
deallocate(sm)
|
||||||
|
!
|
||||||
|
! Save results to filestream for easy access from python
|
||||||
|
!
|
||||||
|
call save_eigenvalues(w, nev, e_kin)
|
||||||
|
!
|
||||||
|
! Write results to OUTFILE on unit IOUT2
|
||||||
|
!
|
||||||
|
write(iout2,75)
|
||||||
|
write(iout2,110)
|
||||||
|
write(iout2,80) e_kin
|
||||||
|
write(iout2,110)
|
||||||
|
write(iout2,75)
|
||||||
|
write(iout2,105)
|
||||||
|
write(iout2,75)
|
||||||
|
!
|
||||||
|
allocate(w1(nev), w2(nev))
|
||||||
|
!
|
||||||
|
n_eig = 0
|
||||||
|
xmax_l = 0.0
|
||||||
|
n_xmax = 0
|
||||||
|
do lin = 1, nev
|
||||||
|
eig = real(abs(w(lin)))
|
||||||
|
write(iout2,100) real(w(lin)), aimag(w(lin)), eig
|
||||||
|
if ((eig-xmax_l) > 0.0001) n_xmax = lin
|
||||||
|
xmax_l = max(xmax_l, eig)
|
||||||
|
w1(lin) = eig
|
||||||
|
if (eig > 1.000) then
|
||||||
|
n_eig = n_eig + 1
|
||||||
|
end if
|
||||||
|
end do
|
||||||
|
!
|
||||||
|
write(iout2,75)
|
||||||
|
write(iout2,85) xmax_l
|
||||||
|
write(iout2,90) n_xmax
|
||||||
|
write(iout2,95) w(n_xmax)
|
||||||
|
write(iout2,75)
|
||||||
|
!
|
||||||
|
! Summarize results in main output file
|
||||||
|
!
|
||||||
|
call ordre(nev, w1, nfin, w2)
|
||||||
|
!
|
||||||
|
write(iuo1,5)
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,15) w2(1)
|
||||||
|
write(iuo1,20) w2(nfin)
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,10)
|
||||||
|
!
|
||||||
|
if (n_eig >= 1) then
|
||||||
|
if (n_eig == 1) then
|
||||||
|
write(iuo1,25) ndim
|
||||||
|
else
|
||||||
|
write(iuo1,30) n_eig, ndim
|
||||||
|
end if
|
||||||
|
end if
|
||||||
|
!
|
||||||
|
write(iuo1,65) n_xmax
|
||||||
|
write(iuo1,70) w(n_xmax)
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iout3,100) real(w(n_xmax)), aimag(w(n_xmax))
|
||||||
|
!
|
||||||
|
npr = nprint/5
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,35) 5*npr
|
||||||
|
write(iuo1,10)
|
||||||
|
do jp = 0, npr-1
|
||||||
|
j = 5*jp
|
||||||
|
write(iuo1,40) w2(j+1), w2(j+2), w2(j+3), w2(j+4), w2(j+5)
|
||||||
|
enddo
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,45) w2(1)
|
||||||
|
write(iuo2,*) e_kin, w2(1)
|
||||||
|
if (n_eig == 0) then
|
||||||
|
write(iuo1,50)
|
||||||
|
else
|
||||||
|
write(iuo1,55)
|
||||||
|
end if
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,10)
|
||||||
|
write(iuo1,60)
|
||||||
|
!
|
||||||
|
deallocate(w, w1, w2)
|
||||||
|
!
|
||||||
|
return
|
||||||
|
!
|
||||||
|
5 format(/,11X,'----------------- EIGENVALUE ANALYSIS ','-----------------')
|
||||||
|
10 format(11X,'-',54X,'-')
|
||||||
|
15 format(11X,'-',14X,'MAXIMUM MODULUS : ',F9.6,13X,'-')
|
||||||
|
20 format(11X,'-',14X,'MINIMUM MODULUS : ',F9.6,13X,'-')
|
||||||
|
25 format(11X,'-',6X,'1 EIGENVALUE IS > 1 OF A TOTAL OF ',I8,6X,'-')
|
||||||
|
30 format(11X,'-',4X,I5,' EIGENVALUES ARE > 1 OF A TOTAL OF ',I8,2X,'-')
|
||||||
|
35 format(11X,'-',11X,'THE ',I3,' LARGEST EIGENVALUES ARE :',11X,'-')
|
||||||
|
40 format(11X,'-',6X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,2X,F7.4,5X,'-')
|
||||||
|
45 format(11X,'-',4X,'SPECTRAL RADIUS OF THE KERNEL MATRIX : ',F8.5,3X,'-')
|
||||||
|
50 format(11X,'-',5X,'---> THE MULTIPLE SCATTERING SERIES ','CONVERGES',4X,'-')
|
||||||
|
55 format(11X,'-',10X,'---> NO CONVERGENCE OF THE MULTIPLE',9X,'-',/,11X,'-', &
|
||||||
|
18X,'SCATTERING SERIES',19X,'-')
|
||||||
|
60 format(11X,'----------------------------------------','----------------',/)
|
||||||
|
65 format(11X,'-',5X,' LABEL OF LARGEST EIGENVALUE : ',I5,8X,'-')
|
||||||
|
70 format(11X,'-',5X,' LARGEST EIGENVALUE : ','(',F6.3,',',F6.3,')',8X,'-')
|
||||||
|
75 format(' ')
|
||||||
|
80 format(' KINETIC ENERGY : ',F7.2,' eV')
|
||||||
|
85 format(' LARGEST MODULUS OF EIGENVALUE : ',F6.3)
|
||||||
|
90 format(' LABEL OF LARGEST EIGENVALUE : ',I5)
|
||||||
|
95 format(' LARGEST EIGENVALUE : (',F6.3,',',F6.3,')')
|
||||||
|
100 format(5X,F9.5,2X,F9.5,2X,F9.5)
|
||||||
|
105 format(7X,'EIGENVALUES :',3X,'MODULUS :')
|
||||||
|
110 format(2X,'-------------------------------')
|
||||||
|
!==============================================================================!
|
||||||
|
end subroutine eig_mat_ar
|
||||||
|
!==============================================================================!
|
||||||
|
|
@ -0,0 +1,103 @@
|
||||||
|
C
|
||||||
|
C=======================================================================
|
||||||
|
C
|
||||||
|
SUBROUTINE EIGDIF_AR
|
||||||
|
C
|
||||||
|
C This subroutine computes some of the eigenvalues of the
|
||||||
|
C multiple scattering matrix using Arnoldi iteration as
|
||||||
|
C implemented in ARPACK.
|
||||||
|
C
|
||||||
|
C Last modified : 21 June 2023
|
||||||
|
C
|
||||||
|
C INCLUDE 'spec.inc'
|
||||||
|
USE DIM_MOD
|
||||||
|
USE CONVTYP_MOD
|
||||||
|
USE COOR_MOD, NTCLU => NATCLU, NTP => NATYP
|
||||||
|
USE DEBWAL_MOD
|
||||||
|
USE EIGEN_MOD, NE => NE_EIG, E0 => E0_EIG, EFIN => EFIN_EIG
|
||||||
|
USE OUTFILES_MOD
|
||||||
|
USE OUTUNITS_MOD
|
||||||
|
USE RESEAU_MOD
|
||||||
|
USE TESTS_MOD
|
||||||
|
USE TRANS_MOD
|
||||||
|
USE VALIN_MOD, E1 => E0, PHLUM => PHILUM
|
||||||
|
C
|
||||||
|
COMPLEX IC,ONEC
|
||||||
|
COMPLEX TLT(0:NT_M,4,NATM,NE_M)
|
||||||
|
C
|
||||||
|
C
|
||||||
|
DATA CONV /0.512314/
|
||||||
|
C
|
||||||
|
IC=(0.,1.)
|
||||||
|
ONEC=(1.,0.)
|
||||||
|
C
|
||||||
|
OPEN(UNIT=IUO2, FILE=OUTFILE2, STATUS='UNKNOWN')
|
||||||
|
C
|
||||||
|
C Loop over the energies
|
||||||
|
C
|
||||||
|
DO JE=1,NE
|
||||||
|
IF(NE.GT.1) THEN
|
||||||
|
ECIN=E0+FLOAT(JE-1)*(EFIN-E0)/FLOAT(NE-1)
|
||||||
|
ELSEIF(NE.EQ.1) THEN
|
||||||
|
ECIN=E0
|
||||||
|
ENDIF
|
||||||
|
CALL LPM(ECIN,XLPM,*6)
|
||||||
|
XLPM1=XLPM/A
|
||||||
|
IF(IPRINT.GT.0) WRITE(IUO1,56) A,XLPM1
|
||||||
|
IF(ITL.EQ.0) THEN
|
||||||
|
VK(JE)=SQRT(ECIN+VINT)*CONV*A*ONEC
|
||||||
|
VK2(JE)=CABS(VK(JE)*VK(JE))
|
||||||
|
ENDIF
|
||||||
|
GAMMA=1./(2.*XLPM1)
|
||||||
|
IF(IPOTC.EQ.0) THEN
|
||||||
|
VK(JE)=VK(JE)+IC*GAMMA
|
||||||
|
ENDIF
|
||||||
|
IF(I_MFP.EQ.0) THEN
|
||||||
|
VK(JE)=CMPLX(REAL(VK(JE)))
|
||||||
|
VK2(JE)=CABS(VK(JE)*VK(JE))
|
||||||
|
ENDIF
|
||||||
|
IF(I_VIB.EQ.1) THEN
|
||||||
|
IF(IDCM.GE.1) WRITE(IUO1,22)
|
||||||
|
DO JAT=1,N_PROT
|
||||||
|
IF(IDCM.EQ.0) THEN
|
||||||
|
XK2UJ2=VK2(JE)*UJ2(JAT)
|
||||||
|
ELSE
|
||||||
|
XK2UJ2=VK2(JE)*UJ_SQ(JAT)
|
||||||
|
WRITE(IUO1,23) JAT,UJ_SQ(JAT)*A*A
|
||||||
|
ENDIF
|
||||||
|
CALL DWSPH(JAT,JE,XK2UJ2,TLT,I_VIB)
|
||||||
|
DO LAT=0,LMAX(JAT,JE)
|
||||||
|
TL(LAT,1,JAT,JE)=TLT(LAT,1,JAT,JE)
|
||||||
|
ENDDO
|
||||||
|
ENDDO
|
||||||
|
ENDIF
|
||||||
|
C
|
||||||
|
C Eigenvalue calculation
|
||||||
|
C
|
||||||
|
ckmd IF(I_PWM.EQ.0) THEN
|
||||||
|
ckmd CALL EIG_MAT_MS(JE,ECIN)
|
||||||
|
ckmd ELSE
|
||||||
|
ckmd CALL SPEC_RAD_POWER(JE,ECIN)
|
||||||
|
ckmd ENDIF
|
||||||
|
C
|
||||||
|
call eig_mat_ar(je, ecin)
|
||||||
|
C
|
||||||
|
C
|
||||||
|
C End of the loop on the energy
|
||||||
|
C
|
||||||
|
ENDDO
|
||||||
|
GOTO 7
|
||||||
|
C
|
||||||
|
6 WRITE(IUO1,55)
|
||||||
|
C
|
||||||
|
22 FORMAT(16X,'INTERNAL CALCULATION OF MEAN SQUARE DISPLACEMENTS',/,
|
||||||
|
&25X,' BY DEBYE UNCORRELATED MODEL:',/)
|
||||||
|
23 FORMAT(21X,'ATOM TYPE ',I5,' MSD = ',F8.6,' ANG**2')
|
||||||
|
55 FORMAT(///,12X,' <<<<<<<<<< THIS VALUE OF ILPM IS NOT',
|
||||||
|
&'AVAILABLE >>>>>>>>>>')
|
||||||
|
56 FORMAT(4X,'LATTICE PARAMETER A = ',F6.3,' ANGSTROEMS',4X,
|
||||||
|
*'MEAN FREE PATH = ',F6.3,' * A',//)
|
||||||
|
C
|
||||||
|
7 RETURN
|
||||||
|
C
|
||||||
|
END
|
||||||
|
|
@ -0,0 +1,22 @@
|
||||||
|
SUBROUTINE RUN(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_,
|
||||||
|
& NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_,
|
||||||
|
& NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_,
|
||||||
|
& N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
|
||||||
|
|
||||||
|
USE DIM_MOD
|
||||||
|
IMPLICIT INTEGER (A-Z)
|
||||||
|
CF2PY INTEGER, INTENT(IN,COPY) :: NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_
|
||||||
|
CF2PY INTEGER, INTENT(IN,COPY) :: NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_
|
||||||
|
CF2PY INTEGER, INTENT(IN,COPY) :: NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_
|
||||||
|
CF2PY INTEGER, INTENT(IN,COPY) :: N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_
|
||||||
|
|
||||||
|
CALL ALLOCATION(NATP_M_, NATCLU_M_, NAT_EQ_M_, N_CL_L_M_,
|
||||||
|
& NE_M_, NL_M_, LI_M_, NEMET_M_, NO_ST_M_, NDIF_M_, NSO_M_,
|
||||||
|
& NTEMP_M_, NODES_EX_M_, NSPIN_M_, NTH_M_, NPH_M_, NDIM_M_,
|
||||||
|
& N_TILT_M_, N_ORD_M_, NPATH_M_, NGR_M_)
|
||||||
|
|
||||||
|
CALL DO_MAIN()
|
||||||
|
CALL CLOSE_ALL_FILES()
|
||||||
|
|
||||||
|
END SUBROUTINE
|
||||||
|
|
||||||
|
|
@ -0,0 +1,14 @@
|
||||||
|
memalloc_src := memalloc/dim_mod.f memalloc/modules.f memalloc/allocation.f
|
||||||
|
cluster_gen_src := $(wildcard cluster_gen/*.f)
|
||||||
|
common_sub_src := $(wildcard common_sub/*.f)
|
||||||
|
renormalization_src := $(wildcard renormalization/*.f)
|
||||||
|
#eig_common_src := $(wildcard eig/common/*.f)
|
||||||
|
eig_common_src := $(filter-out eig/common/lapack_eig.f, $(wildcard eig/common/*.f))
|
||||||
|
eig_ar_src := $(wildcard eig/ar/*.f)
|
||||||
|
eig_ar_src_f90 := $(wildcard eig/ar/*.f90)
|
||||||
|
|
||||||
|
SRCS = $(memalloc_src) $(cluster_gen_src) $(common_sub_src) $(renormalization_src) $(eig_common_src) $(eig_ar_src_f90) $(eig_ar_src)
|
||||||
|
MAIN_F = eig/ar/main.f
|
||||||
|
SO = _eig_ar.so
|
||||||
|
|
||||||
|
include ../../../options.mk
|
||||||
|
|
@ -31,7 +31,7 @@ IFORT_FFLAGS_DBG =
|
||||||
################################################################################
|
################################################################################
|
||||||
# F2PY CONFIGURATION #
|
# F2PY CONFIGURATION #
|
||||||
################################################################################
|
################################################################################
|
||||||
F2PYFLAGS = --opt=-O2 -llapack
|
F2PYFLAGS = --opt=-O2 -llapack -larpack
|
||||||
F2PYFLAGS_DBG = --debug-capi --debug
|
F2PYFLAGS_DBG = --debug-capi --debug
|
||||||
################################################################################
|
################################################################################
|
||||||
|
|
||||||
|
|
@ -103,7 +103,7 @@ endif
|
||||||
|
|
||||||
FFLAGS = $($(PREFIX)_FFLAGS$(SUFFIX))
|
FFLAGS = $($(PREFIX)_FFLAGS$(SUFFIX))
|
||||||
|
|
||||||
OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS))))
|
OBJS = $(addprefix $(BUILDDIR)/, $(patsubst %.f90,%.o, $(patsubst %.f,%.o, $(filter-out $(MAIN_F), $(SRCS)))))
|
||||||
|
|
||||||
.PHONY: clean obj all info
|
.PHONY: clean obj all info
|
||||||
|
|
||||||
|
|
@ -141,6 +141,12 @@ $(BUILDDIR)/%.o: %.f
|
||||||
$(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION)
|
$(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION)
|
||||||
|
|
||||||
|
|
||||||
|
$(BUILDDIR)/%.o: %.f90
|
||||||
|
@echo "Compiling $@..."
|
||||||
|
mkdir -p $(basename $@)
|
||||||
|
$(FC) $(FFLAGS) -J $(BUILDDIR) -I $(BUILDDIR) -fPIC -o $@ -c $^ $(OUPUT_REDIRECTION)
|
||||||
|
|
||||||
|
|
||||||
$(SO): $(OBJS) $(MAIN_F)
|
$(SO): $(OBJS) $(MAIN_F)
|
||||||
@echo "building Python binding $@..."
|
@echo "building Python binding $@..."
|
||||||
mkdir -p $(BUILDDIR)
|
mkdir -p $(BUILDDIR)
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue