From 2659473de7d3ce83dbcd2276bc8a8fb788d23d85 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Sun, 20 Mar 2022 16:09:56 +0100 Subject: [PATCH] Integrate with AtomsBase (#558) --- Project.toml | 2 + docs/src/guide/input_output.jl | 5 +- examples/arbitrary_floattype.jl | 2 +- examples/ase.jl | 13 ++- examples/atomsbase.jl | 26 ++++++ examples/scf_callbacks.jl | 2 +- src/DFTK.jl | 7 +- src/Model.jl | 9 +- src/elements.jl | 9 +- src/external/ase.jl | 7 +- src/external/atomsbase.jl | 116 +++++++++++++++++++++++++ src/guess_density.jl | 11 ++- src/interpolation_transfer.jl | 5 +- src/pseudo/attach_psp.jl | 61 ++++++++++++++ src/pseudo/load_psp.jl | 24 ++---- src/standard_models.jl | 4 +- test/PspHgh.jl | 2 +- test/atomsbase.jl | 145 ++++++++++++++++++++++++++++++++ test/compute_bands.jl | 2 +- test/compute_density.jl | 2 +- test/energies_guess_density.jl | 2 +- test/list_psp.jl | 10 +++ test/load_psp.jl | 24 ------ test/runtests.jl | 3 +- test/scf_compare.jl | 2 +- test/silicon_lda.jl | 2 +- test/silicon_pbe.jl | 2 +- test/silicon_redHF.jl | 2 +- 28 files changed, 426 insertions(+), 75 deletions(-) create mode 100644 examples/atomsbase.jl create mode 100644 src/external/atomsbase.jl create mode 100644 src/pseudo/attach_psp.jl create mode 100644 test/atomsbase.jl create mode 100644 test/list_psp.jl delete mode 100644 test/load_psp.jl diff --git a/Project.toml b/Project.toml index 96f69b51a1..a365659374 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.4.7" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" +AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a" BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" Brillouin = "23470ee3-d0df-4052-8b1a-8cbd6363e7f0" Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" @@ -43,6 +44,7 @@ spglib_jll = "ac4a9f1e-bdb2-5204-990c-47c8b2f70d4e" [compat] AbstractFFTs = "1" +AtomsBase = "0.2.1" BlockArrays = "0.16.2" Brillouin = "0.5" Conda = "1" diff --git a/docs/src/guide/input_output.jl b/docs/src/guide/input_output.jl index 589ba1e137..948836b54c 100644 --- a/docs/src/guide/input_output.jl +++ b/docs/src/guide/input_output.jl @@ -40,7 +40,10 @@ magnetic_moments = load_magnetic_moments(qe_input) # not exposed inside the ASE datastructures. # See [Creating slabs with ASE](@ref) for more details. -atoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional="pbe")) for el in atoms] +atoms = map(atoms) do el + @assert el.symbol == :Fe + ElementPsp(:Fe, psp=load_psp("hgh/pbe/fe-q16.hgh")) +end; # Finally we run the calculation. diff --git a/examples/arbitrary_floattype.jl b/examples/arbitrary_floattype.jl index 29589ef527..3b73dc3ebf 100644 --- a/examples/arbitrary_floattype.jl +++ b/examples/arbitrary_floattype.jl @@ -29,7 +29,7 @@ using DFTK ## Setup silicon lattice a = 10.263141334305942 # lattice constant in Bohr lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]] -Si = ElementPsp(:Si, psp=load_psp(:Si, functional="lda")) +Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4")) atoms = [Si, Si] positions = [ones(3)/8, -ones(3)/8] diff --git a/examples/ase.jl b/examples/ase.jl index 866abae8e5..ddd4fd48b3 100644 --- a/examples/ase.jl +++ b/examples/ase.jl @@ -68,10 +68,17 @@ pyimport("ase.io").write("surface.png", surface * (3, 3, 1), using DFTK -atoms = load_atoms(surface) -atoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional="pbe")) for el in atoms] positions = load_positions(surface) -lattice = load_lattice(surface); +lattice = load_lattice(surface) +atoms = map(load_atoms(surface)) do el + if el.symbol == :Ga + ElementPsp(:Ga, psp=load_psp("hgh/pbe/ga-q3.hgh")) + elseif el.symbol == :As + ElementPsp(:As, psp=load_psp("hgh/pbe/as-q5.hgh")) + else + error("Unsupported element: $el") + end +end; # We model this surface with (quite large a) temperature of 0.01 Hartree # to ease convergence. Try lowering the SCF convergence tolerance (`tol`) diff --git a/examples/atomsbase.jl b/examples/atomsbase.jl new file mode 100644 index 0000000000..05aa5d720e --- /dev/null +++ b/examples/atomsbase.jl @@ -0,0 +1,26 @@ +using DFTK +using AtomsBase +using Unitful +using UnitfulAtomic +using PyCall + +# Construct system in the AtomsBase world +a = 10.26u"bohr" # Silicon lattice constant +lattice = a / 2 * [[0, 1, 1.], # Lattice as vector of vectors + [1, 0, 1.], + [1, 1, 0.]] +atoms = [:Si => ones(3)/8, :Si => -ones(3)/8] +system = periodic_system(atoms, lattice; fractional=true) + +# Use it inside DFTK: +# Attach pseudopotentials, construct appropriate model, +# discretise and solve. +system = attach_psp(system; family="hgh", functional="lda") + +model = model_LDA(system; temperature=1e-3) +basis = PlaneWaveBasis(model; Ecut=15, kgrid=[4, 4, 4]) +scfres = self_consistent_field(basis, tol=1e-8) + +# Get the DFTK model back into the AtomsBase world: +newsystem = atomic_system(model) +@show atomic_number(newsystem) diff --git a/examples/scf_callbacks.jl b/examples/scf_callbacks.jl index 1e356d9d39..2eecadad35 100644 --- a/examples/scf_callbacks.jl +++ b/examples/scf_callbacks.jl @@ -16,7 +16,7 @@ using DFTK using PyCall silicon = pyimport("ase.build").bulk("Si") -atoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional="lda")) +atoms = [ElementPsp(el.symbol, psp=load_psp("hgh/lda/si-q4.hgh")) for el in load_atoms(silicon)] positions = load_positions(silicon) lattice = load_lattice(silicon); diff --git a/src/DFTK.jl b/src/DFTK.jl index cb1d94e10b..d32b7a98c8 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -4,9 +4,9 @@ with plane-wave density-functional theory algorithms. """ module DFTK -using Printf -using Markdown using LinearAlgebra +using Markdown +using Printf using Requires using TimerOutputs using spglib_jll @@ -143,9 +143,11 @@ export guess_density export random_density export load_psp export list_psp +export attach_psp include("guess_density.jl") include("pseudo/load_psp.jl") include("pseudo/list_psp.jl") +include("pseudo/attach_psp.jl") export pymatgen_structure export ase_atoms @@ -154,6 +156,7 @@ export load_atoms export load_positions export load_magnetic_moments export run_wannier90 +include("external/atomsbase.jl") include("external/load_from_file.jl") include("external/ase.jl") include("external/pymatgen.jl") diff --git a/src/Model.jl b/src/Model.jl index d3b020f432..0ef6340154 100644 --- a/src/Model.jl +++ b/src/Model.jl @@ -168,9 +168,11 @@ function Model(lattice::AbstractMatrix{<:Quantity}, args...; kwargs...) Model(austrip.(lattice), args...; kwargs...) end -normalize_magnetic_moment(::Nothing) = Vec3{Float64}(zeros(3)) -normalize_magnetic_moment(mm::Number) = Vec3{Float64}(0, 0, mm) -normalize_magnetic_moment(mm::AbstractVector) = Vec3{Float64}(mm) + +normalize_magnetic_moment(::Nothing)::Vec3{Float64} = (0, 0, 0) +normalize_magnetic_moment(mm::Number)::Vec3{Float64} = (0, 0, mm) +normalize_magnetic_moment(mm::AbstractVector)::Vec3{Float64} = mm + """ :none if no element has a magnetic moment, else :collinear or :full @@ -211,7 +213,6 @@ function default_symmetries(lattice, atoms, positions, magnetic_moments, spin_po end - """ Maximal occupation of a state (2 for non-spin-polarized electrons, 1 otherwise). """ diff --git a/src/elements.jl b/src/elements.jl index 2a932e34f7..5a21adc05a 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -1,4 +1,5 @@ import PeriodicTable +using AtomsBase # Alias to avoid similarity of elements and Element in DFTK module namespace periodic_table = PeriodicTable.elements @@ -14,7 +15,7 @@ abstract type Element end charge_nuclear(::Element) = 0 """Chemical symbol corresponding to an element""" -atomic_symbol(::Element) = :X +AtomsBase.atomic_symbol(::Element) = :X # The preceeding functions are fallback implementations that should be altered as needed. """Return the total ionic charge of an atom type (nuclear charge - core electrons)""" @@ -46,7 +47,7 @@ struct ElementCoulomb <: Element end charge_ionic(el::ElementCoulomb) = el.Z charge_nuclear(el::ElementCoulomb) = el.Z -atomic_symbol(el::ElementCoulomb) = el.symbol +AtomsBase.atomic_symbol(el::ElementCoulomb) = el.symbol """ Element interacting with electrons via a bare Coulomb potential @@ -87,7 +88,7 @@ function ElementPsp(key; psp) end charge_ionic(el::ElementPsp) = el.psp.Zion charge_nuclear(el::ElementPsp) = el.Z -atomic_symbol(el::ElementPsp) = el.symbol +AtomsBase.atomic_symbol(el::ElementPsp) = el.symbol function local_potential_fourier(el::ElementPsp, q::T) where {T <: Real} q == 0 && return zero(T) # Compensating charge background @@ -108,7 +109,7 @@ struct ElementCohenBergstresser <: Element end charge_ionic(el::ElementCohenBergstresser) = 2 charge_nuclear(el::ElementCohenBergstresser) = el.Z -atomic_symbol(el::ElementCohenBergstresser) = el.symbol +AtomsBase.atomic_symbol(el::ElementCohenBergstresser) = el.symbol """ Element where the interaction with electrons is modelled diff --git a/src/external/ase.jl b/src/external/ase.jl index 31bbed4162..041079d15d 100644 --- a/src/external/ase.jl +++ b/src/external/ase.jl @@ -3,8 +3,11 @@ ase_cell(lattice) = pyimport("ase").cell.Cell(Array(lattice)' / austrip(1u"Å")) ase_cell(model::Model) = ase_cell(model.lattice) function ase_atoms(lattice_or_model, atoms, positions, magnetic_moments=[]) - # Collect Z magnetic moments - magmoms = [normalize_magnetic_moment(mom)[3] for mom in magnetic_moments] + if isempty(magnetic_moments) # Collect Z magnetic moments + magmoms = nothing + else + magmoms = [normalize_magnetic_moment(mom)[3] for mom in magnetic_moments] + end cell = ase_cell(lattice_or_model) symbols = string.(atomic_symbol.(atoms)) scaled_positions = reduce(hcat, positions)' diff --git a/src/external/atomsbase.jl b/src/external/atomsbase.jl new file mode 100644 index 0000000000..5493032cd2 --- /dev/null +++ b/src/external/atomsbase.jl @@ -0,0 +1,116 @@ +using AtomsBase +# Key functionality to integrate DFTK and AtomsBase + +function construct_system(lattice::AbstractMatrix, atoms::Vector, positions::Vector, magnetic_moments::Vector) + @assert length(atoms) == length(positions) + @assert isempty(magnetic_moments) || length(magnetic_moments) == length(atoms) + atomsbase_atoms = map(enumerate(atoms)) do (i, element) + kwargs = Dict{Symbol, Any}(:potential => element, ) + if !isempty(magnetic_moments) + kwargs[:magnetic_moment] = normalize_magnetic_moment(magnetic_moments[i]) + end + if element isa ElementPsp + kwargs[:pseudopotential] = element.psp.identifier + end + + position = lattice * positions[i] * u"bohr" + if atomic_symbol(element) == :X # dummy element ... should solve this upstream + Atom(:X, position; atomic_symbol=:X, atomic_number=0, atomic_mass=0u"u", kwargs...) + else + Atom(atomic_symbol(element), position; kwargs...) + end + end + periodic_system(atomsbase_atoms, collect(eachcol(lattice)) * u"bohr") +end + +""" + atomic_system(model::DFTK.Model, magnetic_moments=[]) + +Construct an AtomsBase atomic system from a DFTK model and associated magnetic moments. +""" +function AtomsBase.atomic_system(model::Model, magnetic_moments=[]) + construct_system(model.lattice, model.atoms, model.positions, magnetic_moments) +end + +""" + periodic_system(model::DFTK.Model, magnetic_moments=[]) + +Construct an AtomsBase atomic system from a DFTK model and associated magnetic moments. +""" +function AtomsBase.periodic_system(model::Model, magnetic_moments=[]) + atomic_system(model, magnetic_moments) +end + + +function parse_system(system::AbstractSystem{D}) where {D} + if !all(periodicity(system)) + error("DFTK only supports calculations with periodic boundary conditions.") + end + + # Parse abstract system and return data required to construct model + mtx = austrip.(reduce(hcat, bounding_box(system))) + T = eltype(mtx) + lattice = zeros(T, 3, 3) + lattice[1:D, 1:D] .= mtx + + # Cache for instantiated pseudopotentials. This is done to ensure that identical + # atoms are indistinguishable in memory, which is used in the Model constructor + # to deduce the atom_groups. + cached_pspelements = Dict{String, ElementPsp}() + atoms = map(system) do atom + if hasproperty(atom, :potential) + atom.potential + elseif hasproperty(atom, :pseudopotential) + get!(cached_pspelements, atom.pseudopotential) do + ElementPsp(atomic_symbol(atom); psp=load_psp(atom.pseudopotential)) + end + else + ElementCoulomb(atomic_symbol(atom)) + end + end + + positions = map(system) do atom + coordinate = zeros(T, 3) + coordinate[1:D] = lattice[1:D, 1:D] \ T.(austrip.(position(atom))) + Vec3{T}(coordinate) + end + + magnetic_moments = map(system) do atom + hasproperty(atom, :magnetic_moment) || return nothing + getproperty(atom, :magnetic_moment) + end + if all(m -> isnothing(m) || iszero(m) || isempty(m), magnetic_moments) + empty!(magnetic_moments) + else + magnetic_moments = normalize_magnetic_moment.(magnetic_moments) + end + + # TODO Use system to determine n_electrons + + (; lattice, atoms, positions, magnetic_moments) +end + + +function _call_with_system(f, system::AbstractSystem, args...; kwargs...) + @assert !(:magnetic_moments in keys(kwargs)) + parsed = parse_system(system) + f(parsed.lattice, parsed.atoms, parsed.positions, args...; + parsed.magnetic_moments, kwargs...) +end + + +""" + Model(system::AbstractSystem; kwargs...) + +AtomsBase-compatible Model constructor. Sets structural information (`atoms`, `positions`, +`lattice`, `n_electrons` etc.) from the passed `system`. +""" +Model(system::AbstractSystem; kwargs...) = _call_with_system(Model, system; kwargs...) + + +# Generate equivalent functions for AtomsBase +for fun in (:model_atomic, :model_DFT, :model_LDA, :model_PBE, :model_SCAN) + @eval function $fun(system::AbstractSystem, args...; kwargs...) + _call_with_system($fun, system, args...; kwargs...) + end +end diff --git a/src/guess_density.jl b/src/guess_density.jl index 22bc0d2a7c..9ebf700359 100644 --- a/src/guess_density.jl +++ b/src/guess_density.jl @@ -16,11 +16,12 @@ end @doc raw""" - guess_density(basis, magnetic_moments) + guess_density(basis, magnetic_moments=[]) + guess_density(basis, system) Build a superposition of atomic densities (SAD) guess density. -We take for the guess density a gaussian centered around the atom, of +We take for the guess density a Gaussian centered around the atom, of length specified by `atom_decay_length`, normalized to get the right number of electrons ```math \hat{ρ}(G) = Z \exp\left(-(2π \text{length} |G|)^2\right) @@ -31,7 +32,11 @@ The magnetic moments should be specified in units of ``μ_B``. function guess_density(basis::PlaneWaveBasis, magnetic_moments=[]) guess_density(basis, basis.model.atoms, basis.model.positions, magnetic_moments) end -@timing function guess_density(basis::PlaneWaveBasis, atoms, positions, magnetic_moments=[]) +function guess_density(basis::PlaneWaveBasis, system::AbstractSystem) + parsed = parse_system(system) + guess_density(basis, parsed.atoms, parsed.positions, parsed.magnetic_moments) +end +@timing function guess_density(basis::PlaneWaveBasis, atoms, positions, magnetic_moments) ρtot = _guess_total_density(basis, atoms, positions) if basis.model.n_spin_components == 1 ρspin = nothing diff --git a/src/interpolation_transfer.jl b/src/interpolation_transfer.jl index 3ccc83d7d0..5b371398b9 100644 --- a/src/interpolation_transfer.jl +++ b/src/interpolation_transfer.jl @@ -1,4 +1,5 @@ -using Interpolations +import Interpolations +import Interpolations: interpolate, extrapolate, scale, BSpline, Quadratic, OnCell using SparseArrays """ @@ -46,7 +47,7 @@ function interpolate_density(ρ_in::AbstractArray, grid_in, grid_out, lattice_in # interpolate ρ_in_supercell from grid grid_supercell to grid_out axes_in = (range(0, 1, length=grid_supercell[i]+1)[1:end-1] for i=1:3) - itp = interpolate(ρ_in_supercell, BSpline(Quadratic(Periodic(OnCell())))) + itp = interpolate(ρ_in_supercell, BSpline(Quadratic(Interpolations.Periodic(OnCell())))) sitp = scale(itp, axes_in...) ρ_interp = extrapolate(sitp, Periodic()) ρ_out = similar(ρ_in, grid_out) diff --git a/src/pseudo/attach_psp.jl b/src/pseudo/attach_psp.jl new file mode 100644 index 0000000000..59e9d60438 --- /dev/null +++ b/src/pseudo/attach_psp.jl @@ -0,0 +1,61 @@ +# Attach pseudopotentials to an atomic system + +""" + attach_psp(system::AbstractSystem, pspmap::AbstractDict) + +Return a new system with the `pseudopotential` property of all atoms set according +to the passed `pspmap`, which maps from the atomic symbol to a pseudopotential identifier. + +# Examples +Select pseudopotentials for all silicon and oxygen atoms in the system. +```julia-repl +julia> attach_psp(system, Dict(:Si => "hgh/lda/si-q4", :O => "hgh/lda/o-q6") +``` +""" +function attach_psp(system::AbstractSystem, pspmap::AbstractDict{Symbol,String}) + particles = map(system) do atom + symbol = atomic_symbol(atom) + + # Pseudo or explicit potential already set + if hasproperty(atom, :pseudopotential) || hasproperty(atom, :potential) + Atom(; atom) + elseif !(symbol in keys(pspmap)) + error("No pseudo identifier given for element $symbol.") + else + Atom(; atom, pseudopotential=pspmap[symbol]) + end + end + FlexibleSystem(system; particles) +end + + +function compute_pspmap(symbols::AbstractVector{Symbol}; core=:fullcore, kwargs...) + pspmap = map(unique(sort(symbols))) do symbol + list = list_psp(symbol; core, kwargs...) + if length(list) != 1 + error("Parameters passed do not uniquely identify a PSP file for element $symbol.") + end + symbol => list[1].identifier + end + Dict(pspmap...) +end +function compute_pspmap(system::AbstractSystem; kwargs...) + compute_pspmap(atomic_symbol(system); kwargs...) +end + + +""" + attach_psp(system::AbstractSystem; family=..., functional=..., core=...) + +For each atom look up a pseudopotential in the library using `list_psp`, which matches the +passed parameters and store its identifier in the `pseudopotential` property of all atoms. + +# Examples +Select HGH pseudopotentials for LDA XC functionals for all atoms in the system. +```julia-repl +julia> attach_psp(system; family="hgh", functional="lda") +``` +""" +function attach_psp(system::AbstractSystem; kwargs...) + attach_psp(system, compute_pspmap(system; kwargs...)) +end diff --git a/src/pseudo/load_psp.jl b/src/pseudo/load_psp.jl index ce91623f7a..742ba6a993 100644 --- a/src/pseudo/load_psp.jl +++ b/src/pseudo/load_psp.jl @@ -10,7 +10,7 @@ The file is searched in the directory `datadir_psp()` and by the `key`. If the `key` is a path to a valid file, the extension is used to determine the type of the pseudopotential file format and a respective class is returned. """ -function load_psp(key::AbstractString; identifier=nothing) +function load_psp(key::AbstractString) if startswith(key, "hgh/") || endswith(key, ".hgh") parser = parse_hgh_file extension = ".hgh" @@ -19,35 +19,23 @@ function load_psp(key::AbstractString; identifier=nothing) end Sys.iswindows() && (key = replace(key, "/" => "\\")) - if isfile(key) # Key is a file ... deduce identifier + if isfile(key) # Key is a file ... deduce identifier fullpath = key - if isnothing(identifier) - identifier = replace(key, "\\" => "/") - if startswith(identifier, datadir_psp()) - identifier = identifier[length(datadir_psp())+1:end] - end + identifier = replace(key, "\\" => "/") + if startswith(identifier, datadir_psp()) + identifier = identifier[length(datadir_psp())+1:end] end else # Not a file: Treat as identifier, add extension if needed - @assert isnothing(identifier) fullpath = joinpath(datadir_psp(), lowercase(key)) isfile(fullpath) || (fullpath = fullpath * extension) identifier = replace(lowercase(key), "\\" => "/") end if isfile(fullpath) - return parser(fullpath, identifier=identifier) + return parser(fullpath; identifier) else error("Could not find pseudopotential for identifier " * "'$identifier' in directory '$(datadir_psp())'") end end - - -function load_psp(element::Union{Symbol,Integer}; family="hgh", core=:fullcore, kwargs...) - list = list_psp(element; family=family, core=core, kwargs...) - if length(list) != 1 - error("Parameters passed to load_psp do not uniquely identify a PSP file.") - end - load_psp(list[1].path, identifier=list[1].identifier) -end diff --git a/src/standard_models.jl b/src/standard_models.jl index ef5ea4e245..3ad9790c82 100644 --- a/src/standard_models.jl +++ b/src/standard_models.jl @@ -1,4 +1,6 @@ -# Convenience functions to make standard models +# High-level convenience functions to make standard models +# Note: When adding a function here, also add a method taking an AbstractSystem +# to external/atomsbase.jl """ Convenience constructor, which builds a standard atomic (kinetic + atomic potential) model. diff --git a/test/PspHgh.jl b/test/PspHgh.jl index 47defba47f..780e13d963 100644 --- a/test/PspHgh.jl +++ b/test/PspHgh.jl @@ -80,7 +80,7 @@ end end @testset "Check qcut routines" begin - psp = load_psp(:Au, functional="pbe", family="hgh") + psp = load_psp("hgh/pbe/au-q11.hgh") ε = 1e-6 qcut = qcut_psp_local(Float64, psp) diff --git a/test/atomsbase.jl b/test/atomsbase.jl new file mode 100644 index 0000000000..9b4f83c3ed --- /dev/null +++ b/test/atomsbase.jl @@ -0,0 +1,145 @@ +using DFTK +using Unitful +using UnitfulAtomic +using AtomsBase +using Test +import DFTK: compute_pspmap + +@testset "DFTK -> AbstractSystem -> DFTK" begin + Si = ElementCoulomb(:Si) + C = ElementPsp(:C, psp=load_psp("hgh/pbe/c-q4.hgh")) + H = ElementPsp(:H, psp=load_psp("hgh/lda/h-q1.hgh")) + + lattice = randn(3, 3) + atoms = [Si, C, H, C] + positions = [rand(3) for _ in 1:4] + magnetic_moments = DFTK.normalize_magnetic_moment.(rand(4)) + + system = DFTK.construct_system(lattice, atoms, positions, magnetic_moments) + @test atomic_symbol(system) == [:Si, :C, :H, :C] + @test bounding_box(system) == collect(eachcol(lattice)) * u"bohr" + @test position(system) == [lattice * p * u"bohr" for p in positions] + + @test system[1].potential == Si + @test system[2].potential == C + @test system[3].potential == H + @test system[4].potential == C + @test !hasproperty(system[1], :pseudopotential) + @test system[2].pseudopotential == "hgh/pbe/c-q4.hgh" + @test system[3].pseudopotential == "hgh/lda/h-q1.hgh" + @test system[4].pseudopotential == "hgh/pbe/c-q4.hgh" + + parsed = DFTK.parse_system(system) + @test parsed.lattice == lattice + @test parsed.atoms == atoms + @test parsed.positions ≈ positions atol=1e-14 + @test parsed.magnetic_moments == magnetic_moments + + let system = attach_psp(system; family="hgh", functional="lda") + for i in 1:4 + @test system[i].potential == atoms[i] + @test system[i].magnetic_moment == magnetic_moments[i] + end + end + + for constructor in (Model, model_atomic, model_LDA, model_PBE, model_SCAN) + model = constructor(system) + @test model.spin_polarization == :collinear + newsys = atomic_system(model, magnetic_moments) + + @test atomic_symbol(system) == atomic_symbol(newsys) + @test bounding_box(system) == bounding_box(newsys) + @test boundary_conditions(system) == boundary_conditions(newsys) + @test maximum(maximum, position(system) - position(newsys)) < 1e-12u"bohr" + for (p, newp) in zip(system, newsys) + @test p.magnetic_moment == newp.magnetic_moment + @test p.potential == newp.potential + end + end +end + +@testset "AbstractSystem -> DFTK" begin + lattice = [12u"bohr" * rand(3) for _ in 1:3] + atoms = [:C => rand(3), :Si => rand(3), :H => rand(3), :C => rand(3)] + pos_units = last.(atoms) + pos_lattice = austrip.(reduce(hcat, lattice)) + system = periodic_system(atoms, lattice; fractional=true) + + let model = Model(system) + @test model.lattice == pos_lattice + @test model.positions ≈ pos_units atol=1e-14 + @test model.spin_polarization == :none + + @test length(model.atoms) == 4 + @test model.atoms[1] == ElementCoulomb(:C) + @test model.atoms[2] == ElementCoulomb(:Si) + @test model.atoms[3] == ElementCoulomb(:H) + @test model.atoms[4] == ElementCoulomb(:C) + end + + let system = attach_psp(system; family="hgh", functional="pbe") + @test length(system) == 4 + @test system[1].pseudopotential == "hgh/pbe/c-q4.hgh" + @test system[2].pseudopotential == "hgh/pbe/si-q4.hgh" + @test system[3].pseudopotential == "hgh/pbe/h-q1.hgh" + @test system[4].pseudopotential == "hgh/pbe/c-q4.hgh" + + parsed = DFTK.parse_system(system) + @test parsed.lattice == pos_lattice + @test parsed.positions ≈ pos_units atol=1e-14 + @test isempty(parsed.magnetic_moments) + + @test length(parsed.atoms) == 4 + @test parsed.atoms[1].psp.identifier == "hgh/pbe/c-q4.hgh" + @test parsed.atoms[2].psp.identifier == "hgh/pbe/si-q4.hgh" + @test parsed.atoms[3].psp.identifier == "hgh/pbe/h-q1.hgh" + @test parsed.atoms[4].psp.identifier == "hgh/pbe/c-q4.hgh" + end + + let system = attach_psp(system; family="hgh", functional="lda") + @test length(system) == 4 + @test system[1].pseudopotential == "hgh/lda/c-q4.hgh" + @test system[2].pseudopotential == "hgh/lda/si-q4.hgh" + @test system[3].pseudopotential == "hgh/lda/h-q1.hgh" + @test system[4].pseudopotential == "hgh/lda/c-q4.hgh" + + model = Model(system) + @test model.lattice == pos_lattice + @test model.positions ≈ pos_units atol=1e-14 + @test model.spin_polarization == :none + + @test length(model.atoms) == 4 + @test model.atoms[1].psp.identifier == "hgh/lda/c-q4.hgh" + @test model.atoms[2].psp.identifier == "hgh/lda/si-q4.hgh" + @test model.atoms[3].psp.identifier == "hgh/lda/h-q1.hgh" + @test model.atoms[4].psp.identifier == "hgh/lda/c-q4.hgh" + end +end + + +@testset "Check attach_psp routine selectively" begin + symbols = [:Cu, :Au, :Ni] + + let pspmap = compute_pspmap(symbols; functional="lda", family="hgh", core=:semicore) + @test pspmap[:Cu] == "hgh/lda/cu-q19.hgh" + @test pspmap[:Au] == "hgh/lda/au-q19.hgh" + @test pspmap[:Ni] == "hgh/lda/ni-q18.hgh" + end + + let pspmap = compute_pspmap(symbols; functional="lda", family="hgh", core=:fullcore) + @test pspmap[:Cu] == "hgh/lda/cu-q11.hgh" + @test pspmap[:Au] == "hgh/lda/au-q11.hgh" + @test pspmap[:Ni] == "hgh/lda/ni-q10.hgh" + end + + let pspmap = compute_pspmap([:Cu, :Au]; functional="pbe", family="hgh", core=:semicore) + @test pspmap[:Cu] == "hgh/pbe/cu-q19.hgh" + @test pspmap[:Au] == "hgh/pbe/au-q19.hgh" + end + + let pspmap = compute_pspmap(symbols; functional="pbe", family="hgh", core=:fullcore) + @test pspmap[:Cu] == "hgh/pbe/cu-q11.hgh" + @test pspmap[:Au] == "hgh/pbe/au-q11.hgh" + @test pspmap[:Ni] == "hgh/pbe/ni-q18.hgh" + end +end diff --git a/test/compute_bands.jl b/test/compute_bands.jl index 6fca135de0..18d9138ae5 100644 --- a/test/compute_bands.jl +++ b/test/compute_bands.jl @@ -122,7 +122,7 @@ end basis = PlaneWaveBasis(model, Ecut, testcase.kcoords, testcase.kweights) # Build Hamiltonian just from SAD guess - ρ0 = guess_density(basis, testcase.atoms, testcase.positions) + ρ0 = guess_density(basis) ham = Hamiltonian(basis; ρ=ρ0) # Check that plain diagonalization and compute_bands agree diff --git a/test/compute_density.jl b/test/compute_density.jl index 36d95d1281..e9a3582fdd 100644 --- a/test/compute_density.jl +++ b/test/compute_density.jl @@ -21,7 +21,7 @@ if mpi_nprocs() == 1 # not easy to distribute model = model_DFT(testcase.lattice, testcase.atoms, testcase.positions, :lda_xc_teter93; kwargs...) basis = PlaneWaveBasis(model, Ecut, kcoords, kweights; symmetries) - ham = Hamiltonian(basis; ρ=guess_density(basis, testcase.atoms, testcase.positions)) + ham = Hamiltonian(basis; ρ=guess_density(basis)) res = diagonalize_all_kblocks(lobpcg_hyper, ham, n_bands; tol) occ, εF = DFTK.compute_occupation(basis, res.λ) diff --git a/test/energies_guess_density.jl b/test/energies_guess_density.jl index 8f5d5d1b5b..647df56196 100644 --- a/test/energies_guess_density.jl +++ b/test/energies_guess_density.jl @@ -18,7 +18,7 @@ include("testcases.jl") [:lda_x, :lda_c_vwn]; symmetries=false) basis = PlaneWaveBasis(model; Ecut, kgrid, fft_size, kshift) - ρ0 = guess_density(basis, silicon.atoms, silicon.positions) + ρ0 = guess_density(basis) E, H = energy_hamiltonian(basis, nothing, nothing; ρ=ρ0) @test E["Hartree"] ≈ 0.3527293727197568 atol=5e-8 diff --git a/test/list_psp.jl b/test/list_psp.jl new file mode 100644 index 0000000000..9fdaadc59e --- /dev/null +++ b/test/list_psp.jl @@ -0,0 +1,10 @@ +using Test +using DFTK: load_psp, list_psp + +@testset "Check reading all HGH pseudos" begin + for record in list_psp() + psp = load_psp(record.identifier) + @test psp.identifier == record.identifier + @test psp.Zion == record.n_elec_valence + end +end diff --git a/test/load_psp.jl b/test/load_psp.jl deleted file mode 100644 index 7efd3c8e9d..0000000000 --- a/test/load_psp.jl +++ /dev/null @@ -1,24 +0,0 @@ -using Test -using DFTK: load_psp, list_psp - -@testset "Check reading all HGH pseudos" begin - for record in list_psp() - psp = load_psp(record.identifier) - @test psp.identifier == record.identifier - @test psp.Zion == record.n_elec_valence - end -end - -@testset "Check load_psp routine selectively" begin - psp_Cu = load_psp(:Cu, functional="lda", family="hgh", core=:semicore) - @test psp_Cu.identifier == "hgh/lda/cu-q19.hgh" - - psp_Au = load_psp(:Au, functional="pbe", family="hgh", core=:fullcore) - @test psp_Au.identifier == "hgh/pbe/au-q11.hgh" - - for element in (:Cu, :Ni, :Au) - sc = load_psp(element, functional="lda", family="hgh", core=:semicore) - fc = load_psp(element, functional="lda", family="hgh", core=:fullcore) - @test sc.Zion > fc.Zion - end -end diff --git a/test/runtests.jl b/test/runtests.jl index 8cb627ab34..362a7f70a9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -61,7 +61,8 @@ Random.seed!(0) include("PlaneWaveBasis.jl") include("Model.jl") include("interpolation.jl") - include("load_psp.jl") + include("atomsbase.jl") + include("list_psp.jl") include("PspHgh.jl") include("elements.jl") include("bzmesh.jl") diff --git a/test/scf_compare.jl b/test/scf_compare.jl index e945f9245c..fb22d7bd9a 100644 --- a/test/scf_compare.jl +++ b/test/scf_compare.jl @@ -38,7 +38,7 @@ include("testcases.jl") end # Run other SCFs with SAD guess - ρ0 = guess_density(basis, silicon.atoms, silicon.positions) + ρ0 = guess_density(basis) for solver in (scf_nlsolve_solver(), scf_damping_solver(1.2), scf_anderson_solver(), scf_CROP_solver()) @testset "Testing $solver" begin diff --git a/test/silicon_lda.jl b/test/silicon_lda.jl index 6b81467a3c..8d03b0e773 100644 --- a/test/silicon_lda.jl +++ b/test/silicon_lda.jl @@ -17,7 +17,7 @@ function run_silicon_lda(T ;Ecut=5, grid_size=15, spin_polarization=:none, kwarg ref_etot = -7.911817522631488 fft_size = fill(grid_size, 3) - Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.atnum, functional="lda", family="hgh")) + Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/lda/si-q4")) atoms = [Si, Si] if spin_polarization == :collinear diff --git a/test/silicon_pbe.jl b/test/silicon_pbe.jl index 703873fcd9..fdb16763f6 100644 --- a/test/silicon_pbe.jl +++ b/test/silicon_pbe.jl @@ -21,7 +21,7 @@ function run_silicon_pbe(T ;Ecut=5, grid_size=15, spin_polarization=:none, kwarg ref_etot = -7.854477356672080 fft_size = fill(grid_size, 3) - Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.atnum, functional="pbe", family="hgh")) + Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/pbe/si-q4")) atoms = [Si, Si] if spin_polarization == :collinear diff --git a/test/silicon_redHF.jl b/test/silicon_redHF.jl index cc504a7e65..fdece15d03 100644 --- a/test/silicon_redHF.jl +++ b/test/silicon_redHF.jl @@ -26,7 +26,7 @@ function run_silicon_redHF(T; Ecut=5, grid_size=15, spin_polarization=:none, kwa fft_size = fill(grid_size, 3) fft_size = DFTK.next_working_fft_size(T, fft_size) # ad-hoc fix for buggy generic FFTs - Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.atnum, functional="lda", family="hgh")) + Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/lda/si-q4")) atoms = [Si, Si] if spin_polarization == :collinear