|
| 1 | +@testsetup module Regression |
| 2 | +using DFTK |
| 3 | +using Unitful |
| 4 | +using UnitfulAtomic |
| 5 | +using AtomsBase |
| 6 | +using ..TestCases: magnesium |
| 7 | + |
| 8 | +high_symmetry = let |
| 9 | + a = 4.474 |
| 10 | + lattice = [[0, a, a], [a, 0, a], [a, a, 0]]u"bohr" |
| 11 | + x = 6.711 |
| 12 | + y = 2.237 |
| 13 | + atoms = [ |
| 14 | + Atom(:Cu, [0, 0, 0]u"bohr", magnetic_moment=0), |
| 15 | + Atom(:O, [x, y, x]u"bohr", magnetic_moment=0), |
| 16 | + Atom(:O, [x, y, y]u"bohr", magnetic_moment=0), |
| 17 | + ] |
| 18 | + system = periodic_system(atoms, lattice) |
| 19 | + merge(DFTK.parse_system(system), (; temperature=0.03, Ecut=10, kgrid=[4,4,4], |
| 20 | + n_electrons=45)) |
| 21 | +end |
| 22 | +high_kpoints = merge(magnesium, (; kgrid=[13,13,13], Ecut=10)) |
| 23 | +high_Ecut = merge(magnesium, (; kgrid=[4,4,4], Ecut=60)) |
| 24 | +testcases = (; high_symmetry, high_kpoints, high_Ecut) |
| 25 | +end |
| 26 | + |
| 27 | +@testitem "Hamiltonian application" tags=[:regression] setup=[TestCases, Regression] begin |
| 28 | + using DFTK |
| 29 | + using LinearAlgebra |
| 30 | + using BenchmarkTools |
| 31 | + using .Main: SUITE |
| 32 | + |
| 33 | + for testcase in Regression.testcases |
| 34 | + model = Model(testcase.lattice, testcase.atoms, testcase.positions; |
| 35 | + testcase.temperature, terms=[Kinetic()]) |
| 36 | + basis = PlaneWaveBasis(model; testcase.Ecut, testcase.kgrid) |
| 37 | + |
| 38 | + n_electrons = testcase.n_electrons |
| 39 | + n_bands = div(n_electrons, 2, RoundUp) |
| 40 | + ψ = [Matrix(qr(randn(ComplexF64, length(G_vectors(basis, kpt)), n_bands)).Q) |
| 41 | + for kpt in basis.kpoints] |
| 42 | + filled_occ = DFTK.filled_occupation(model) |
| 43 | + occupation = [filled_occ * rand(n_bands) for _ = 1:length(basis.kpoints)] |
| 44 | + occ_scaling = n_electrons / sum(sum(occupation)) |
| 45 | + occupation = [occ * occ_scaling for occ in occupation] |
| 46 | + |
| 47 | + (; ham) = energy_hamiltonian(basis, ψ, occupation) |
| 48 | + |
| 49 | + SUITE["ham"] = BenchmarkGroup() |
| 50 | + SUITE["ham"] = @benchmarkable for ik = 1:length($(basis.kpoints)) |
| 51 | + $(ham.blocks)[ik]*$ψ[ik] |
| 52 | + end |
| 53 | + end |
| 54 | +end |
| 55 | + |
| 56 | +@testitem "Single SCF step" tags=[:regression] setup=[TestCases, Regression] begin |
| 57 | + using DFTK |
| 58 | + using BenchmarkTools |
| 59 | + using .Main: SUITE |
| 60 | + |
| 61 | + for testcase in Regression.testcases |
| 62 | + model = model_LDA(testcase.lattice, testcase.atoms, testcase.positions; |
| 63 | + testcase.temperature) |
| 64 | + basis = PlaneWaveBasis(model; testcase.Ecut, testcase.kgrid) |
| 65 | + SUITE["scf"] = BenchmarkGroup() |
| 66 | + SUITE["scf"] = @benchmarkable self_consistent_field($basis; tol=1e5) |
| 67 | + end |
| 68 | +end |
| 69 | + |
| 70 | +@testitem "Density + symmetrization" tags=[:regression] setup=[TestCases, Regression] begin |
| 71 | + using DFTK |
| 72 | + using BenchmarkTools |
| 73 | + using .Main: SUITE |
| 74 | + |
| 75 | + for testcase in Regression.testcases |
| 76 | + model = model_LDA(testcase.lattice, testcase.atoms, testcase.positions; |
| 77 | + testcase.temperature) |
| 78 | + basis = PlaneWaveBasis(model; testcase.Ecut, testcase.kgrid) |
| 79 | + scfres = self_consistent_field(basis; tol=10) |
| 80 | + |
| 81 | + ψ, occupation = DFTK.select_occupied_orbitals(basis, scfres.ψ, scfres.occupation; |
| 82 | + threshold=1e-6) |
| 83 | + |
| 84 | + SUITE["density", "ρ", "sym"] = BenchmarkGroup() |
| 85 | + SUITE["density", "ρ"] = @benchmarkable compute_density($basis, $ψ, $occupation) |
| 86 | + SUITE["density", "sym"] = @benchmarkable DFTK.symmetrize_ρ($basis, $(scfres.ρ)) |
| 87 | + end |
| 88 | +end |
| 89 | + |
| 90 | +@testitem "Basis construction" tags=[:regression] setup=[TestCases, Regression] begin |
| 91 | + using DFTK |
| 92 | + using BenchmarkTools |
| 93 | + using .Main: SUITE |
| 94 | + |
| 95 | + for testcase in Regression.testcases |
| 96 | + model = model_LDA(testcase.lattice, testcase.atoms, testcase.positions; |
| 97 | + testcase.temperature) |
| 98 | + SUITE["basis"] = BenchmarkGroup() |
| 99 | + SUITE["basis"] = @benchmarkable PlaneWaveBasis($model; Ecut=$(testcase.Ecut), |
| 100 | + kgrid=$(testcase.kgrid)) |
| 101 | + end |
| 102 | +end |
| 103 | + |
| 104 | +@testitem "Sternheimer" tags=[:regression] setup=[TestCases, Regression] begin |
| 105 | + using DFTK |
| 106 | + using BenchmarkTools |
| 107 | + using .Main: SUITE |
| 108 | + |
| 109 | + for testcase in Regression.testcases |
| 110 | + model = model_LDA(testcase.lattice, testcase.atoms, testcase.positions; |
| 111 | + testcase.temperature) |
| 112 | + basis = PlaneWaveBasis(model; testcase.Ecut, testcase.kgrid) |
| 113 | + scfres = self_consistent_field(basis; tol=10) |
| 114 | + |
| 115 | + rhs = DFTK.compute_projected_gradient(basis, scfres.ψ, scfres.occupation) |
| 116 | + SUITE["sternheimer"] = BenchmarkGroup() |
| 117 | + SUITE["sternheimer"] = @benchmarkable DFTK.solve_ΩplusK_split($scfres, $rhs) |
| 118 | + end |
| 119 | +end |
| 120 | + |
| 121 | +@testitem "Response with AD" tags=[:regression] setup=[TestCases, Regression] begin |
| 122 | + using DFTK |
| 123 | + using BenchmarkTools |
| 124 | + using LinearAlgebra |
| 125 | + using ForwardDiff |
| 126 | + using .Main: SUITE |
| 127 | + |
| 128 | + function make_basis(ε::T; a=10., Ecut=30) where {T} |
| 129 | + lattice=T(a) * I(3) # lattice is a cube of ``a`` Bohrs |
| 130 | + # Helium at the center of the box |
| 131 | + atoms = [ElementPsp(:He; psp=load_psp("hgh/lda/He-q2"))] |
| 132 | + positions = [[1/2, 1/2, 1/2]] |
| 133 | + |
| 134 | + model = model_DFT(lattice, atoms, positions, [:lda_x, :lda_c_vwn]; |
| 135 | + extra_terms=[ExternalFromReal(r -> -ε * (r[1] - a/2))], |
| 136 | + symmetries=false) |
| 137 | + PlaneWaveBasis(model; Ecut, kgrid=[1, 1, 1]) # No k-point sampling on isolated system |
| 138 | + end |
| 139 | + |
| 140 | + # dipole moment of a given density (assuming the current geometry) |
| 141 | + function dipole(basis, ρ) |
| 142 | + @assert isdiag(basis.model.lattice) |
| 143 | + a = basis.model.lattice[1, 1] |
| 144 | + rr = [a * (r[1] - 1/2) for r in r_vectors(basis)] |
| 145 | + sum(rr .* ρ) * basis.dvol |
| 146 | + end |
| 147 | + |
| 148 | + # Function to compute the dipole for a given field strength |
| 149 | + function compute_dipole(ε; tol=1e-8, kwargs...) |
| 150 | + scfres = self_consistent_field(make_basis(ε; kwargs...); tol) |
| 151 | + dipole(scfres.basis, scfres.ρ) |
| 152 | + end |
| 153 | + |
| 154 | + SUITE["response", "ad"] = BenchmarkGroup() |
| 155 | + SUITE["response", "ad"] = @benchmarkable ForwardDiff.derivative($compute_dipole, 0.0) |
| 156 | +end |
0 commit comments