|
| 1 | +""" |
| 2 | + EDIIS |
| 3 | +""" |
| 4 | +mutable struct EDIIS <: Accelerator |
| 5 | + state::Tuple{DiisState, DiisState} |
| 6 | + sync_spins::Bool |
| 7 | + conditioning_threshold::Float64 |
| 8 | + coefficient_threshold::Float64 |
| 9 | + |
| 10 | + function EDIIS(problem::ScfProblem; n_diis_size = 5, sync_spins = true, conditioning_threshold = 1e-14, coefficient_threshold = 1e-6, kwargs...) |
| 11 | + stateα = DiisState(n_diis_size) |
| 12 | + stateβ = DiisState(n_diis_size) |
| 13 | + new((stateα, stateβ), sync_spins, conditioning_threshold, coefficient_threshold) |
| 14 | + end |
| 15 | +end |
| 16 | + |
| 17 | +""" |
| 18 | + Calculates coefficients |
| 19 | +""" |
| 20 | +function diis_solve_coefficients(::EDIIS, B::AbstractArray; energybuffer::AbstractArray, kwargs...) |
| 21 | + m = size(B,1) |
| 22 | + E = map(energies -> energies["total"], energybuffer) |
| 23 | + |
| 24 | + function f(x) |
| 25 | + c = x.^2/sum(x.^2) |
| 26 | + E'*c - 1/2 * c'*B*c |
| 27 | + end |
| 28 | + |
| 29 | + guess = ones(m) / m |
| 30 | + guess[1] = 1 |
| 31 | + |
| 32 | + options = Optim.Options(x_tol = 10e-5) |
| 33 | + res = optimize(f, guess, BFGS(), options; inplace = false, autodiff = :forward) |
| 34 | + x = Optim.minimizer(res) |
| 35 | + c = x.^2/sum(x.^2) |
| 36 | + |
| 37 | + # If number of iterations in optimization is zero, reuse old matrix |
| 38 | + if Optim.iterations(res) == 0 |
| 39 | + c = zeros(m) |
| 40 | + c[1] = 1 |
| 41 | + end |
| 42 | + |
| 43 | + return c, 0 |
| 44 | +end |
| 45 | + |
| 46 | +""" |
| 47 | + Linear System Matrix for the EDIIS accelerator. |
| 48 | +""" |
| 49 | +function diis_build_matrix(::EDIIS, state::DiisState) |
| 50 | + @assert state.n_diis_size > 0 |
| 51 | + @assert state.iterate.length > 0 |
| 52 | + |
| 53 | + # B has dimension m <= state.n_diis_size, since in the |
| 54 | + # beginning of the iteration we do not have the full number of |
| 55 | + # previous iterates yet. |
| 56 | + |
| 57 | + m = min(state.n_diis_size, length(state.iterate)) |
| 58 | + |
| 59 | + B = zeros(m,m) |
| 60 | + |
| 61 | + # Since the Matrix B is symmetric, we only have to calculate |
| 62 | + # the upper triagonal and can use the Julia object 'Symmetric' |
| 63 | + # to fill the lower triagonal of the matrix. This also allows |
| 64 | + # Julia to use optimized algorithems for symmetric matrices. |
| 65 | + # |
| 66 | + # The matrix B is filled in the following form: |
| 67 | + # |
| 68 | + # B[1,1] B[1,2] B[1,3] … B[1,m] |
| 69 | + # 0 B[2,2] B[2,3] … B[2,m] |
| 70 | + # 0 0 B[3,3] … B[3,m] |
| 71 | + # ⋮ ⋮ ⋮ … ⋮ |
| 72 | + # 0 0 0 … B[m,m] |
| 73 | + # |
| 74 | + # Note, that the values B[1,1] … B[1,m-1] will become the values |
| 75 | + # of B[2,2] … B[2,m] in the next iteration. This means, that we effectively |
| 76 | + # only have to calculate the first row and can reuse all consecutive ones. |
| 77 | + # |
| 78 | + # We would like to store the rows in such a way, that the storage variable |
| 79 | + # can be mapped to a row immediately. Since the values are shifted to the |
| 80 | + # right in each iteration and the last element is discarded it makes sense |
| 81 | + # to use a Circular Buffer to hold each row. |
| 82 | + # |
| 83 | + # Since the last row is also discarded after the new row is calculated we |
| 84 | + # use a Circular Buffer again to store the rows. |
| 85 | + |
| 86 | + # Definition of matrix elements |
| 87 | + b(i,j) = tr((state.iterate[i] - state.iterate[j]) * (state.density[i] - state.density[j])) |
| 88 | + |
| 89 | + # Fill the first row with newly calculated values and cache them |
| 90 | + # in a newly created Circular Buffer |
| 91 | + newValues = CircularBuffer{Any}(state.n_diis_size) |
| 92 | + map(j -> push!(newValues, b(1,j)), 1:m) |
| 93 | + fill!(newValues, 0) |
| 94 | + |
| 95 | + # Push newly calculated row to the row buffer. We use the iterationstate to |
| 96 | + # store these. |
| 97 | + pushfirst!(state.iterationstate, newValues) |
| 98 | + |
| 99 | + # Now fill all rows with cached values, |
| 100 | + # push a '0' on each buffer to prepare it for the next iteration |
| 101 | + for i in 1:m |
| 102 | + B[i,1:m] = state.iterationstate[i][1:m] |
| 103 | + |
| 104 | + # Since we want to use this buffer as the 2nd row of A in the next |
| 105 | + # iteration we need the following layout of the buffer |
| 106 | + # 0 A[1,1] A[1,2] … A[1,m-1] |
| 107 | + # so we need to push a 0 to the beginning |
| 108 | + pushfirst!(state.iterationstate[i], 0) |
| 109 | + end |
| 110 | + |
| 111 | + return Symmetric(B) |
| 112 | +end |
| 113 | + |
| 114 | +function needs_error(::EDIIS) |
| 115 | + false |
| 116 | +end |
| 117 | + |
| 118 | +function needs_density(::EDIIS) |
| 119 | + true |
| 120 | +end |
| 121 | + |
| 122 | +""" |
| 123 | + When synchronizing spins the resulting DIIS matrices have to be added |
| 124 | + together but the constraint must be kept as is. |
| 125 | +""" |
| 126 | +function merge_diis_matrices_spin_blocks(::EDIIS, B1::AbstractArray, B2::AbstractArray) |
| 127 | + B1 + B2 |
| 128 | +end |
0 commit comments