/
radial_orbitals.jl
61 lines (48 loc) · 1.96 KB
/
radial_orbitals.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
# TODO: Replace by FuncArray et al.
const RadialOrbitalArray{T,N,B<:Basis} = Applied{<:Any, typeof(*), <:Tuple{BasisOrRestricted{B},AbstractArray{T,N}}}
const AdjointRadialOrbitalArray{T,N,B<:Basis} = Applied{<:Any, typeof(*), <:Tuple{Adjoint{<:Any,<:AbstractArray{T,N}},AdjointBasisOrRestricted{B}}}
function Base.similar(A::RadialOrbitalArray)
R,v = A.args
applied(*, R, similar(v))
end
"""
RadialOrbital
A radial orbital is represented using basis coupled with a vector of
expansion coefficients with respect to that basis; the basis
implemented as an `AbstractQuasimatrix`.
"""
const RadialOrbital{T,B} = RadialOrbitalArray{T,1,B}
const AdjointRadialOrbital{T,B} = AdjointRadialOrbitalArray{T,1,B}
Base.axes(v::RadialOrbital, p::Int) = p == 1 ? axes(v)[1] : 1
"""
RadialOrbitals
A collection of radial orbitals is instead represented using a
matrix of such expansion coefficients, where each matrix column
corresponds to a single radial orbital.
"""
const RadialOrbitals{T,B<:Basis} = RadialOrbitalArray{T,2,B}
const AdjointRadialOrbitals{T,B<:Basis} = AdjointRadialOrbitalArray{T,2,B}
const RadialOperator{B<:Basis,M<:AbstractMatrix} =
Applied{<:Any,typeof(*),<:Tuple{BasisOrRestricted{B},M,AdjointBasisOrRestricted{B}}}
function LinearAlgebra.adjoint(op::RadialOperator)
Ac,M,B = op.args
applied(*, B', M', parent(Ac))
end
function Base.show(io::IO, ro::RadialOperator{B,M}) where {B,M}
write(io, "RadialOperator(")
show(io, M)
write(io, ")")
end
Base.iszero(op::RadialOperator) = iszero(op.args[2])
Base.iszero(op::MulQuasiArray{<:Any,2,<:Tuple{
<:BasisOrRestricted,<:AbstractMatrix,<:AdjointBasisOrRestricted}}) =
iszero(op.applied)
matrix(o::RadialOperator) = o.args[2]
function Base.zero(o::RadialOperator)
A,B,C = o.args
applied(*, A, Diagonal(zeros(eltype(B),size(B,1))), C)
end
function Base.:(+)(a::RadialOperator{B}, b::RadialOperator{B}) where B
R = first(a.args)
applied(*, R, matrix(a) + matrix(b), R')
end