Refactor atoms field in Model (JuliaMolSim#626)

docs/make.jl

@@ -14,8 +14,9 @@ SRCPATH   = joinpath(@__DIR__, "src")
 BUILDPATH = joinpath(@__DIR__, "build")
 ROOTPATH  = joinpath(@__DIR__, "..")
 CONTINUOUS_INTEGRATION = get(ENV, "CI", nothing) == "true"
-DFTKREV   = LibGit2.head(ROOTPATH)
-DFTKREPO  = "github.com/JuliaMolSim/DFTK.jl.git"
+DFTKREV    = LibGit2.head(ROOTPATH)
+DFTKBRANCH = try LibGit2.branch(LibGit2.GitRepo(ROOTPATH)) catch end
+DFTKREPO   = "github.com/JuliaMolSim/DFTK.jl.git"
 
 # Python and Julia dependencies needed for running the notebooks
 PYDEPS = ["ase", "pymatgen"]
@@ -126,7 +127,7 @@ makedocs(
 )
 
 # Dump files for managing dependencies in binder
-if CONTINUOUS_INTEGRATION
+if CONTINUOUS_INTEGRATION && DFTKBRANCH == "master"
     cd(BUILDPATH) do
         open("environment.yml", "w") do io
             print(io,

docs/src/advanced/conventions.md

@@ -60,8 +60,8 @@ auconvert(Å, 1.2)  # 1.2 Bohr in Ångström
     When setting up a calculation for DFTK
     one needs to ensure to convert to Bohr and atomic units.
     For some Python libraries (currently ASE and pymatgen)
-    DFTK directly ships conversion tools in form of the [`load_lattice`](@ref)
-    and [`load_atoms`](@ref) functions,
+    DFTK directly ships conversion tools in form of the `load_lattice`,
+    `load_positions` and `load_atoms` functions,
     which take care of such conversions. Examples which demonstrate this
     are [Creating slabs with ASE](@ref) and [Creating supercells with pymatgen](@ref).
 
@@ -80,8 +80,8 @@ reciprocal-space lattice vectors by ``B = 2\pi A^{-T}``.
     (notably Python and C) use row-major ordering.
     Care therefore needs to be taken to properly
     transpose the unit cell matrices ``A`` before using it with DFTK.
-    For the supported third-party packages `load_lattice`
-    and `load_atoms` again handle such conversion automatically.
+    For the supported third-party packages `load_lattice`,
+    `load_positions` and `load_atoms` again handle such conversion automatically.
 
 We use the convention that the unit cell in real space is
 ``[0, 1)^3`` in reduced coordinates and the unit cell in reciprocal

docs/src/advanced/data_structures.md

@@ -7,9 +7,10 @@ lattice = a / 2 * [[0 1 1.];
                    [1 0 1.];
                    [1 1 0.]]
 Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
-atoms = [Si => [ones(3)/8, -ones(3)/8]]
+atoms = [Si, Si]
+positions = [ones(3)/8, -ones(3)/8]
 
-model = model_LDA(lattice, atoms)
+model = model_LDA(lattice, atoms, positions)
 kgrid = [4, 4, 4]
 Ecut = 15
 basis = PlaneWaveBasis(model; Ecut, kgrid)
@@ -70,13 +71,13 @@ not limited to DFT. Convenience constructors are provided for common cases:
    [libxc](https://tddft.org/programs/libxc/functionals/) library,
    for example:
    ```
-   model_DFT(lattice, atoms, [:gga_x_pbe, :gga_c_pbe])
-   model_DFT(lattice, atoms, :lda_xc_teter93)
+   model_DFT(lattice, atoms, positions, [:gga_x_pbe, :gga_c_pbe])
+   model_DFT(lattice, atoms, positions, :lda_xc_teter93)
    ```
    where the latter is equivalent to `model_LDA`.
    Specifying no functional is the reduced Hartree-Fock model:
    ```
-   model_DFT(lattice, atoms, [])
+   model_DFT(lattice, atoms, positions, [])
    ```
 - `model_atomic`: A linear model, which contains no electron-electron interaction
   (neither Hartree nor XC term).

docs/src/advanced/symmetries.md

@@ -92,24 +92,25 @@ lattice = a / 2 * [[0 1 1.];
                    [1 0 1.];
                    [1 1 0.]]
 Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
-atoms = [Si => [ones(3)/8, -ones(3)/8]]
+atoms = [Si, Si]
+positions = [ones(3)/8, -ones(3)/8]
 Ecut = 5
 kgrid = [4, 4, 4]
 ```
 Let us demonstrate this in practice.
-We consider silicon, setup appropriately in the `lattice` and `atoms` objects
-as in [Tutorial](@ref) and to reach a fast execution, we take a small `Ecut` of `5`
+We consider silicon, setup appropriately in the `lattice`, `atoms` and `positions`
+objects as in [Tutorial](@ref) and to reach a fast execution, we take a small `Ecut` of `5`
 and a `[4, 4, 4]` Monkhorst-Pack grid.
 First we perform the DFT calculation disabling symmetry handling
 ```@example symmetries
-model_nosym = model_LDA(lattice, atoms; symmetries=false)
+model_nosym = model_LDA(lattice, atoms, positions; symmetries=false)
 basis_nosym = PlaneWaveBasis(model_nosym; Ecut, kgrid)
 scfres_nosym = @time self_consistent_field(basis_nosym, tol=1e-8)
 nothing  # hide
 ```
 and then redo it using symmetry (the default):
 ```@example symmetries
-model_sym = model_LDA(lattice, atoms)
+model_sym = model_LDA(lattice, atoms, positions)
 basis_sym = PlaneWaveBasis(model_sym; Ecut, kgrid)
 scfres_sym = @time self_consistent_field(basis_sym, tol=1e-8)
 nothing  # hide

docs/src/guide/input_output.jl

@@ -26,9 +26,10 @@
 
 using DFTK
 
-qe_input = "Fe_afm.pwi"
-atoms    = load_atoms(qe_input)
-lattice  = load_lattice(qe_input)
+qe_input  = "Fe_afm.pwi"
+atoms     = load_atoms(qe_input)
+positions = load_positions(qe_input)
+lattice   = load_lattice(qe_input)
 magnetic_moments = load_magnetic_moments(qe_input)
 
 # At this point a file of any format supported by ASE could be passed instead,
@@ -39,12 +40,11 @@ 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")) => position
-         for (el, position) in atoms];
+atoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional="pbe")) for el in atoms]
 
 # Finally we run the calculation.
 
-model = model_LDA(lattice, atoms, magnetic_moments=magnetic_moments, temperature=0.01)
+model = model_LDA(lattice, atoms, positions; magnetic_moments, temperature=0.01)
 basis = PlaneWaveBasis(model; Ecut=10, kgrid=(2, 2, 2))
 ρ0 = guess_density(basis, magnetic_moments)
 scfres = self_consistent_field(basis, tol=1e-4, ρ=ρ0, mixing=KerkerMixing());

docs/src/guide/parallelization.md

@@ -10,9 +10,10 @@ lattice = a / 2 * [[0 1 1.];
                    [1 0 1.];
                    [1 1 0.]]
 Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
-atoms = [Si => [ones(3)/8, -ones(3)/8]]
+atoms = [Si, Si]
+positions = [ones(3)/8, -ones(3)/8]
 
-model = model_LDA(lattice, atoms)
+model = model_LDA(lattice, atoms, positions)
 basis = PlaneWaveBasis(model; Ecut=5, kgrid=[2, 2, 2])
 
 DFTK.reset_timer!(DFTK.timer)

docs/src/guide/periodic_problems.jl

@@ -311,11 +311,12 @@ lattice = diagm([100., 0, 0])
 
 ## Place them at 20 and 80 in *fractional coordinates*,
 ## that is 0.2 and 0.8, since the lattice is 100 wide.
-nucleus = ElementGaussian()
-atoms = [nucleus => [[0.2, 0, 0], [0.8, 0, 0]]]
+nucleus   = ElementGaussian()
+atoms     = [nucleus, nucleus]
+positions = [[0.2, 0, 0], [0.8, 0, 0]]
 
 ## Assemble the model, discretize and build the Hamiltonian
-model = Model(lattice; atoms=atoms, terms=[Kinetic(), AtomicLocal()])
+model = Model(lattice; atoms, positions, terms=[Kinetic(), AtomicLocal()])
 basis = PlaneWaveBasis(model; Ecut=300, kgrid=(1, 1, 1));
 ham   = Hamiltonian(basis)
 

docs/src/guide/tutorial.jl

@@ -51,10 +51,11 @@ lattice = a / 2 * [[0 1 1.];  # Silicon lattice vectors
 Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
 
 ## Specify type and positions of atoms
-atoms = [Si => [ones(3)/8, -ones(3)/8]]
+atoms     = [Si, Si]
+positions = [ones(3)/8, -ones(3)/8]
 
 ## 2. Select model and basis
-model = model_LDA(lattice, atoms)
+model = model_LDA(lattice, atoms, positions)
 kgrid = [4, 4, 4]     # k-point grid (Regular Monkhorst-Pack grid)
 Ecut = 7              # kinetic energy cutoff
 ## Ecut = 190.5u"eV"  # Could also use eV or other energy-compatible units
@@ -96,9 +97,7 @@ plot(x, scfres.ρ[1, :, 1, 1], label="", xlabel="x", ylabel="ρ", marker=2)
 # for instance compute a band structure
 plot_bandstructure(scfres; kline_density=10)
 # or get the cartesian forces (in Hartree / Bohr)
-compute_forces_cart(scfres)[1]  # Select silicon forces
-# The `[1]` extracts the forces for the first kind of atoms,
-# i.e. `Si` (silicon) in the setup of the `atoms` list of step 1 above.
+compute_forces_cart(scfres)
 # As expected, they are almost zero in this highly symmetric configuration.
 
 # ## Where to go from here

examples/arbitrary_floattype.jl

@@ -30,10 +30,11 @@ using DFTK
 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"))
-atoms = [Si => [ones(3)/8, -ones(3)/8]]
+atoms = [Si, Si]
+positions = [ones(3)/8, -ones(3)/8]
 
 ## Cast to Float32, setup model and basis
-model = model_DFT(Array{Float32}(lattice), atoms, [:lda_x, :lda_c_vwn])
+model = model_DFT(Array{Float32}(lattice), atoms, positions, [:lda_x, :lda_c_vwn])
 basis = PlaneWaveBasis(model, Ecut=7, kgrid=[4, 4, 4])
 
 ## Run the SCF

examples/ase.jl

@@ -59,7 +59,7 @@ pyimport("ase.io").write("surface.png", surface * (3, 3, 1),
 #md # ```
 #nb # <img src="surface.png" width=500 height=500 />
 
-# Use the `load_atoms` and `load_lattice` functions
+# Use the `load_atoms`, `load_positions` and `load_lattice` functions
 # to convert to DFTK datastructures.
 # These two functions not only support importing ASE atoms into DFTK,
 # but a few more third-party datastructures as well.
@@ -69,15 +69,15 @@ 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")) => position
-         for (el, position) in atoms]
+atoms = [ElementPsp(el.symbol, psp=load_psp(el.symbol, functional="pbe")) for el in atoms]
+positions = load_positions(surface)
 lattice = load_lattice(surface);
 
 # We model this surface with (quite large a) temperature of 0.01 Hartree
 # to ease convergence. Try lowering the SCF convergence tolerance (`tol`)
 # or the `temperature` or try `mixing=KerkerMixing()`
 # to see the full challenge of this system.
-model = model_DFT(lattice, atoms, [:gga_x_pbe, :gga_c_pbe],
+model = model_DFT(lattice, atoms, positions, [:gga_x_pbe, :gga_c_pbe],
                   temperature=0.001, smearing=DFTK.Smearing.Gaussian())
 basis = PlaneWaveBasis(model; Ecut, kgrid)
 

examples/cohen_bergstresser.jl

@@ -12,11 +12,12 @@
 using DFTK
 
 Si = ElementCohenBergstresser(:Si)
+atoms = [Si, Si]
+positions = [ones(3)/8, -ones(3)/8]
 lattice = Si.lattice_constant / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
-atoms = [Si => [ones(3)/8, -ones(3)/8]];
 
-# Next we build the rather simple model and discretise it with moderate `Ecut`:
-model = Model(lattice; atoms=atoms, terms=[Kinetic(), AtomicLocal()])
+# Next we build the rather simple model and discretize it with moderate `Ecut`:
+model = Model(lattice; atoms, positions, terms=[Kinetic(), AtomicLocal()])
 basis = PlaneWaveBasis(model, Ecut=10.0, kgrid=(1, 1, 1));
 
 # We diagonalise at the Gamma point to find a Fermi level ...

examples/collinear_magnetism.jl

@@ -14,18 +14,18 @@ a = 5.42352  # Bohr
 lattice = a / 2 * [[-1  1  1];
                    [ 1 -1  1];
                    [ 1  1 -1]]
-Fe = ElementPsp(:Fe, psp=load_psp("hgh/lda/Fe-q8.hgh"))
-atoms = [Fe => [zeros(3)]];
+atoms     = [ElementPsp(:Fe, psp=load_psp("hgh/lda/Fe-q8.hgh"))]
+positions = [zeros(3)]
 
 # To get the ground-state energy we use an LDA model and rather moderate
 # discretisation parameters.
 
 kgrid = [3, 3, 3]  # k-point grid (Regular Monkhorst-Pack grid)
 Ecut = 15          # kinetic energy cutoff in Hartree
-model_nospin = model_LDA(lattice, atoms, temperature=0.01)
+model_nospin = model_LDA(lattice, atoms, positions, temperature=0.01)
 basis_nospin = PlaneWaveBasis(model_nospin; kgrid, Ecut)
 
-scfres_nospin = self_consistent_field(basis_nospin, tol=1e-6, mixing=KerkerDosMixing());
+scfres_nospin = self_consistent_field(basis_nospin; tol=1e-6, mixing=KerkerDosMixing());
 #-
 scfres_nospin.energies
 
@@ -44,7 +44,7 @@ scfres_nospin.energies
 # The structure (i.e. pair mapping and order) of the `magnetic_moments` array needs to agree
 # with the `atoms` array and `0` magnetic moments need to be specified as well.
 
-magnetic_moments = [Fe => [4, ]];
+magnetic_moments = [4];
 
 # !!! tip "Units of the magnetisation and magnetic moments in DFTK"
 #     Unlike all other quantities magnetisation and magnetic moments in DFTK
@@ -56,7 +56,7 @@ magnetic_moments = [Fe => [4, ]];
 # We repeat the calculation using the same model as before. DFTK now detects
 # the non-zero moment and switches to a collinear calculation.
 
-model = model_LDA(lattice, atoms, magnetic_moments=magnetic_moments, temperature=0.01)
+model = model_LDA(lattice, atoms, positions; magnetic_moments, temperature=0.01)
 basis = PlaneWaveBasis(model; Ecut, kgrid)
 ρ0 = guess_density(basis, magnetic_moments)
 scfres = self_consistent_field(basis, tol=1e-6; ρ=ρ0, mixing=KerkerDosMixing());

examples/compare_solvers.jl

@@ -8,27 +8,28 @@ lattice = a / 2 * [[0 1 1.];
                    [1 0 1.];
                    [1 1 0.]]
 Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
-atoms = [Si => [ones(3)/8, -ones(3)/8]]
+atoms     = [Si, Si]
+positions = [ones(3)/8, -ones(3)/8]
 
-model = model_LDA(lattice, atoms)
+model = model_LDA(lattice, atoms, positions)
 basis = PlaneWaveBasis(model; Ecut=5, kgrid=[3, 3, 3])
 tol = 1e-12
 
 ## SCF
-scfres_scf = self_consistent_field(basis, tol=tol,
+scfres_scf = self_consistent_field(basis; tol=tol,
                                    is_converged=DFTK.ScfConvergenceDensity(tol))
 
 ## Direct minimization
-scfres_dm = direct_minimization(basis, tol=tol)
+scfres_dm = direct_minimization(basis; tol=tol)
 
 ## Newton algorithm
 # start not too far from the solution to ensure convergence : we use here the
 # solution of a single iteration SCF
-scfres_start = self_consistent_field(basis, maxiter=1)
+scfres_start = self_consistent_field(basis; maxiter=1)
 # remove virtual orbitals
 ψ, _ = DFTK.select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation)
-scfres_newton = newton(basis, ψ, tol=tol)
+scfres_newton = newton(basis, ψ; tol=tol)
 
 println("|ρ_newton - ρ_scf| = ", norm(scfres_newton.ρ - scfres_scf.ρ))
-println("|ρ_newton - ρ_dm| = ", norm(scfres_newton.ρ - scfres_dm.ρ))
-println("|ρ_scf - ρ_dm| = ", norm(scfres_scf.ρ - scfres_dm.ρ))
+println("|ρ_newton - ρ_dm| =  ", norm(scfres_newton.ρ - scfres_dm.ρ))
+println("|ρ_scf - ρ_dm| =     ", norm(scfres_scf.ρ - scfres_dm.ρ))

examples/custom_potential.jl

@@ -34,36 +34,37 @@ lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];
 # ``0.5`` to break symmetry and get nonzero forces.
 x1 = 0.2
 x2 = 0.8
-nucleus = ElementGaussian(1.0, 0.5)
-atoms = [nucleus => [[x1, 0, 0], [x2, 0, 0]]];
+positions = [[x1, 0, 0], [x2, 0, 0]]
+gauss     = ElementGaussian(1.0, 0.5)
+atoms     = [gauss, gauss]
 
 # We setup a Gross-Pitaevskii model
 C = 1.0
 α = 2;
 n_electrons = 1  # Increase this for fun
 terms = [Kinetic(),
          AtomicLocal(),
-         LocalNonlinearity(ρ -> C * ρ^α),
-]
-model = Model(lattice; atoms=atoms, n_electrons=n_electrons, terms=terms,
+         LocalNonlinearity(ρ -> C * ρ^α)]
+model = Model(lattice; atoms, positions, n_electrons, terms,
               spin_polarization=:spinless);  # use "spinless electrons"
 
 # We discretize using a moderate Ecut and run a SCF algorithm to compute forces
-# afterwards. As there is no ionic charge associated to `nucleus` we have to specify
+# afterwards. As there is no ionic charge associated to `gauss` we have to specify
 # a starting density and we choose to start from a zero density.
 basis = PlaneWaveBasis(model; Ecut=500, kgrid=(1, 1, 1))
 ρ = zeros(eltype(basis), basis.fft_size..., 1)
-scfres = self_consistent_field(basis, tol=1e-8, ρ=ρ)
+scfres = self_consistent_field(basis; tol=1e-8, ρ=ρ)
 scfres.energies
+
 # Computing the forces can then be done as usual:
-hcat(compute_forces(scfres)...)
+compute_forces(scfres)
 
 # Extract the converged total local potential
 tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1
 
 # Extract other quantities before plotting them
-ρ = scfres.ρ[:, 1, 1, 1]  # converged density, first spin component
-ψ_fourier = scfres.ψ[1][:, 1];    # first k-point, all G components, first eigenvector
+ρ = scfres.ρ[:, 1, 1, 1]        # converged density, first spin component
+ψ_fourier = scfres.ψ[1][:, 1]   # first k-point, all G components, first eigenvector
 ψ = G_to_r(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]
 ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));
 

examples/custom_solvers.jl

@@ -8,10 +8,11 @@ lattice = a / 2 * [[0 1 1.];
                    [1 0 1.];
                    [1 1 0.]]
 Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
-atoms = [Si => [ones(3)/8, -ones(3)/8]]
+atoms = [Si, Si]
+positions =  [ones(3)/8, -ones(3)/8]
 
 ## We take very (very) crude parameters
-model = model_LDA(lattice, atoms)
+model = model_LDA(lattice, atoms, positions)
 basis = PlaneWaveBasis(model; Ecut=5, kgrid=[1, 1, 1]);
 
 # We define our custom fix-point solver: simply a damped fixed-point