main
1module GExt
2""" Perform extrapolation of density matrices using their geometric structure with
3 least-squares method.
4
5 nbas: number of basis functions
6 nocc: number of occupied orbitals (half the number of electrons)
7 nmat: number of elements used to perform the extrapolation plus one (the first element is
8 supposed to contain the point we use as a reference, and we do not need its
9 associated positions)
10 nqm: number of QM atoms
11"""
12
13using LinearAlgebra, DelimitedFiles, Printf
14# To make sure the guesses we write are precise enough
15Base.show(io::IO, f::Float64) = @printf(io, "%24.16e", f)
16
17export gext
18
19function gext(; ε=1e-5)
20 """ Main function which does the extrapolation """
21 cmat = get_cmat("cmat.dat")
22 smat = get_smat("smat.dat")
23 hpos = get_hpos("hpos.dat")
24
25 # Make sure that we have the correct number of elements
26 @assert size(smat, 2) == size(cmat, 3) == size(hpos, 3)
27 nbas, nocc, nmat = size(cmat)
28
29 # Perform orthonormalisation
30 C = Array{Float64}(undef, nbas, nocc, nmat)
31 for i in 1:nmat
32 C[:, :, i] = sqrt(tril_to_sym(smat[:, i]))*cmat[:, :, i]
33 end
34
35 # The reference point
36 γ_ref = C[:, :, 1]
37
38 # Project the points on the tangent space at γ_ref
39 Γs = logGrs(γ_ref, C)
40
41 # Create a matrix of previous descriptors
42 poly = [get_coulomb(hpos[:, :, imat]) for imat in 1:nmat-1]
43 P = [poly[i][j] for i in eachindex(poly), j in eachindex(poly[1])]
44
45 # Descriptor of the current point
46 poly_test = get_coulomb(hpos[:, :, nmat])
47
48 # Optional stablisation
49 ds = hcat(poly_test, zeros(nmat-1))
50 Ps = hcat(P, ε*Matrix(I, nmat-1, nmat-1))
51
52 # Least-squares on the Coulomb descriptors
53 α = ds'*pinv(Ps)
54
55 # Extrapolation on the tangent space
56 Γguess = sum(α[i]*Γs[:, :, i+1] for i in 1:nmat-1)
57 # Retract the point back to the Grassmannian
58 Cguess = expGr(γ_ref, Γguess)
59
60 # Compute the idempotent guess
61 pguess = 2*Cguess*Cguess'
62 lguess = vec_tril(pguess)
63
64 # Write the guess on file
65 writedlm("guess", lguess)
66end
67
68function get_coulomb(hpos)
69 """ Return the Coulomb descriptor, which is the lower triangular part of the symmetric
70 Coulomb matrix """
71 charges = hpos[1, :]
72 positions = hpos[2:4, :]
73
74 nqm = size(positions, 2)
75 G = Matrix{Float64}(undef, nqm, nqm)
76 for i in 1:nqm
77 for j in 1:nqm
78 if i == j
79 G[i, i] = 0.5*charges[i]^(2.4)
80 else
81 G[i, j] = charges[i]*charges[j]/norm(positions[:, i] - positions[:, j])
82 end
83 end
84 end
85
86 return vec_tril(G)
87end
88
89#############################################
90#### Exponential and logarithm functions ####
91#############################################
92
93function expGr(γ_ref, Γ)
94 """ Perform the exponential of the tangent vector Γ at base point γ_ref """
95 d = svd(Γ)
96 Y = γ_ref * d.V * cos(Diagonal(d.S)) + (d.U) * sin(Diagonal(d.S))
97 # Optional, to make sure that rounding errors do not prevent Y from being orthogonal
98 return Matrix(qr(Y).Q)
99end
100
101function logGr(γ_ref, γ)
102 """ Perform the logarithm of the point γ at base point γ_ref """
103 Z = γ'*γ_ref
104 W = Z \ (γ'-Z*γ_ref')
105 d = svd(W, full = false)
106 A = (d.V)*atan(Diagonal(svd(W, full = false).S))*(d.U)'
107 return A
108end
109
110function logGrs(γ_ref, γs)
111 """ Perform the logarithms on a list of points """
112 nbas, nocc, nmat = size(cmat)
113
114 Γs = Array{Float64}(undef, nbas, nocc, nmat)
115 for iγ in 1:nmat
116 Γs[:, :, iγ] = logGr(γ_ref, γs[:, :, iγ])
117 end
118
119 return Γs
120end
121
122#############################################
123#### Helper functions ####
124#############################################
125
126function tril_to_sym(v)
127 """ Convert a vector representing the lower triangular part of a symmetric matrix to
128 a symmetric matrix """
129 n = Int(floor(sqrt(2*length(v))))
130 M = zeros(n, n)
131
132 k = 0
133 for i in 1:n
134 for j in 1:i
135 k += 1
136 M[i, j] = v[k]
137 M[j, i] = v[k]
138 end
139 end
140
141 return M
142end
143
144function vec_tril(M)
145 """ Extract the the lower triangular part of a square matrix """
146 m, n = size(M)
147 @assert m == n
148 l = n*(n+1) ÷ 2
149 v = zeros(l)
150
151 k = 0
152 for i in 1:n
153 for j in 1:i
154 v[k + j] = M[i, j]
155 end
156 k += i
157 end
158
159 return v
160end
161
162#############################################
163#### Data processing functions ####
164#############################################
165
166function get_cmat(file)
167 """ Placeholder for a function that returns an array cmat of size nbas×nocc×nmat
168 containing molecular orbitals matrices """
169end
170
171function get_smat(file)
172 """ Placeholder for a function that returns an array smat of size nbas×(nbas+1)/2×nmat
173 containing lower diagonal of overlap matrices """
174end
175
176function get_hpos(file)
177 """ Placeholder for a function that returns an array smat of size 4×nqm×nmat
178 containing charges and positions of the QM atoms for the previous positions as well as
179 the current for which we want the extrapolation """
180end
181
182end