Merge branch 'master' into inexact_krylov

Author: mfherbst

benchmark/Project.toml

@@ -2,7 +2,9 @@
 AtomsIO = "1692102d-eeb4-4df9-807b-c9517f998d44"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DFTK = "acf6eb54-70d9-11e9-0013-234b7a5f5337"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LibGit2 = "76f85450-5226-5b5a-8eaa-529ad045b433"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"

benchmark/benchmarks.jl

@@ -13,3 +13,4 @@ const SUITE = BenchmarkGroup()
 @include_case SUITE "silicon"
 @include_case SUITE "aluminium_rattled"
 @include_case SUITE "SrVO3"
+@include_case SUITE "aluminium12"

benchmark/cases/Al12.extxyz

@@ -0,0 +1,14 @@
+12
+pbc=[T, T, T] Lattice="12.15000000 0.00000000 0.00000000 0.00000000 4.05000000 0.00000000 0.00000000 0.00000000 4.05000000" Properties=species:S:1:pos:R:3:velocities:R:3:Z:I:1
+Al         0.00000000       0.00000000       0.00000000         0.00000000       0.00000000       0.00000000         13
+Al         0.00000000       2.02500000       2.02500000         0.00000000       0.00000000       0.00000000         13
+Al         2.02500000       0.00000000       2.02500000         0.00000000       0.00000000       0.00000000         13
+Al         2.02500000       2.02500000       0.00000000         0.00000000       0.00000000       0.00000000         13
+Al         4.05000000       0.00000000       0.00000000         0.00000000       0.00000000       0.00000000         13
+Al         4.05000000       2.02500000       2.02500000         0.00000000       0.00000000       0.00000000         13
+Al         6.07500000       0.00000000       2.02500000         0.00000000       0.00000000       0.00000000         13
+Al         6.07500000       2.02500000       0.00000000         0.00000000       0.00000000       0.00000000         13
+Al         8.10000000       0.00000000       0.00000000         0.00000000       0.00000000       0.00000000         13
+Al         8.10000000       2.02500000       2.02500000         0.00000000       0.00000000       0.00000000         13
+Al        10.12500000       0.00000000       2.02500000         0.00000000       0.00000000       0.00000000         13
+Al        10.12500000       2.02500000       0.00000000         0.00000000       0.00000000       0.00000000         13

benchmark/cases/aluminium12.jl

@@ -0,0 +1,26 @@
+# This is a 12-atom Al supercell without any rattling.
+
+module aluminium12
+include("common.jl")
+
+const SUITE = BenchmarkGroup(["metal", "PBE", "fast"])
+function setup_basis(; symmetries=true, kgrid=(3, 7, 7))
+    system = load_system(joinpath(@__DIR__, "Al12.extxyz"))
+    pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.stringent.upf")
+    model  = model_DFT(system; pseudopotentials, functionals=PBE(), symmetries,
+                               temperature=1e-3, smearing=Smearing.Gaussian())
+    PlaneWaveBasis(model; Ecut=30, kgrid)
+end
+
+#=
+basis   = setup_basis()
+scfres  = setup_dummy_scfres(basis)
+=#
+
+# Response benchmark
+scfres_nosym = setup_dummy_scfres(setup_basis(; symmetries=false, kgrid=(1, 2, 2));
+                                  mixing=KerkerMixing())
+perturb = setup_atomic_perturbation(scfres_nosym)
+SUITE["response"] = @benchmarkable bm_response($scfres_nosym, $perturb) evals=1 samples=1
+
+end

benchmark/cases/aluminium_rattled.jl

@@ -14,8 +14,8 @@ function setup_basis()
     PlaneWaveBasis(model; Ecut=64, kgrid=(1, 1, 1))
 end
 
-basis  = setup_basis()
-scfres = setup_dummy_scfres(basis)
+basis   = setup_basis()
+scfres  = setup_dummy_scfres(basis)
 
 # Assemble benchmarks
 SUITE["setup"] = @benchmarkable setup_basis() evals=1 samples=3

benchmark/cases/common.jl

@@ -1,14 +1,46 @@
 using AtomsIO
 using BenchmarkTools
 using DFTK
+using ForwardDiff
+using LinearAlgebra
 using PseudoPotentialData
 using Random
 
-function setup_dummy_scfres(basis)
+function setup_dummy_scfres(basis; kwargs...)
     Random.seed!(42)
-    self_consistent_field(basis; maxiter=1, callback=identity)
+    self_consistent_field(basis; maxiter=1, callback=identity, kwargs...)
 end
 
+function setup_atomic_perturbation(scfres; n_displaced=4)
+    # Displace up to n_displaced atoms into a random direction.
+    #
+    basis = scfres.basis
+    model = scfres.basis.model
+    @assert length(model.symmetries) == 1  # Symmetries disabled in model
+
+    Random.seed!(1234)
+    n_displaced = min(n_displaced, length(model.positions))
+    displaced_atoms = shuffle(1:length(model.positions))[1:n_displaced]
+    R = map(1:length(model.positions)) do i
+        if i in displaced_atoms
+            normalize!(-ones(3) + 2 * rand(3))
+        else
+            zeros(3)
+        end
+    end
+
+    function displacement(ε::T) where {T}
+        newpositions = basis.model.positions + ε * R
+        modelV = Model(Matrix{T}(model.lattice), model.atoms, newpositions;
+                       terms=[DFTK.AtomicLocal(), DFTK.AtomicNonlocal()], symmetries=false)
+        basisV = PlaneWaveBasis(modelV; basis.Ecut, basis.kgrid)
+        DFTK.total_local_potential(Hamiltonian(basisV))
+    end
+    δV  = ForwardDiff.derivative(displacement, 0.0)
+    δHψ = -DFTK.multiply_ψ_by_blochwave(scfres.basis, scfres.ψ, δV, zeros(3))
+end
+
+
 #
 # Benchmarks
 #
@@ -21,10 +53,15 @@ end
 
 bm_compute_forces(scfres) = compute_forces_cart(scfres)
 
+function bm_response(scfres, δHψ; tol=1e-6, kwargs...)
+    Random.seed!(101)
+    DFTK.solve_ΩplusK_split(scfres, δHψ; tol, callback=identity, kwargs...)
+end
+
 # TODO other benchmarks to include:
 #  - Hamiltonian application
 #  - Diagonalisation
-#  - Response
+#  - Elastic constant computation
 
 #
 # Benchmark collections

benchmark/cases/silicon.jl

@@ -4,7 +4,7 @@ module silicon
 include("common.jl")
 
 const SUITE = BenchmarkGroup(["insulator", "LDA", "fast"])
-function setup_basis()
+function setup_basis(; symmetries=true, kgrid=(8, 8, 8))
     a = 10.26  # Silicon lattice constant in Bohr
     lattice = a / 2 * [[0 1 1.];
                        [1 0 1.];
@@ -13,8 +13,8 @@ function setup_basis()
     atoms     = [Si, Si]
     positions = [ones(3)/8, -ones(3)/8]
 
-    model = model_DFT(lattice, atoms, positions; functionals=LDA())
-    PlaneWaveBasis(model; Ecut=35, kgrid=(8, 8, 8))
+    model = model_DFT(lattice, atoms, positions; functionals=LDA(), symmetries)
+    PlaneWaveBasis(model; Ecut=35, kgrid)
 end
 
 basis  = setup_basis()
@@ -23,4 +23,10 @@ scfres = setup_dummy_scfres(basis)
 # Assemble benchmarks
 SUITE["setup"] = @benchmarkable setup_basis() evals=1 samples=3
 add_default_benchmarks!(SUITE, basis, scfres)
+
+# Response benchmark
+scfres_nosym = setup_dummy_scfres(setup_basis(; symmetries=false, kgrid=(4, 4, 4)))
+perturb = setup_atomic_perturbation(scfres_nosym)
+SUITE["response"] = @benchmarkable bm_response($scfres_nosym, $perturb; tol=1e-8) evals=1 samples=1
+
 end