Possible flag-based interface to reduce duplication

Project.toml

@@ -6,12 +6,14 @@ version = "0.2.0"
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
+CubicSplines = "9c784101-8907-5a6d-9be6-98f00873c89b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 EzXML = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PeriodicTable = "7b2266bf-644c-5ea3-82d8-af4bbd25a884"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
@@ -22,10 +24,12 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 Bessels = "0.2"
+CubicSplines = "0.2"
 DocStringExtensions = "0.9"
 EzXML = "1"
 Interpolations = "0.14"
 OffsetArrays = "1"
+OrderedCollections = "1"
 PeriodicTable = "1"
 Polynomials = "3"
 PrettyTables = "2"

docs/src/api.md

@@ -13,24 +13,52 @@ PsPFile
 UpfFile
 Psp8File
 HghFile
+
 AbstractPsP
+
 NumericPsP{T}
 NormConservingPsP{T}
 UltrasoftPsP{T}
 ProjectorAugmentedWavePsP{T}
+
 AnalyticalPsP
 HghPsP{T}
+
+AbstractPsPQuantity
+
+AbstractDensity
+ValenceChargeDensity
+CoreChargeDensity
+
+AbstractProjector
+NumericProjector
+NumericState
+
+PsPCoupling
+BetaCoupling
+AugmentationCoupling
+
+AbstractPotential
+NumericLocalPotential
+AugmentationFunction
+
+EvaluationSpace
+RealSpace
+FourierSpace
 ```
 
 ## Functions
 
 ```@docs
 load_psp_file
 load_psp
+load_family_psp_files
+load_family_psps
 list_families
-load_family
+list_family_psps
+show_family_summary
 show_family_periodic_table
-show_family_list
+show_family_table
 
 identifier
 element
@@ -41,43 +69,21 @@ relativistic_treatment
 formalism
 valence_charge
 atomic_charge
-n_projector_radials
-n_projector_angulars
-n_chi_function_radials
-n_chi_function_angulars
 is_norm_conserving
 is_ultrasoft
 is_paw
 has_spin_orbit
-has_core_density
-has_valence_density
-has_chi_functions
-local_potential_cutoff_radius
-projector_cutoff_radius
-chi_function_cutoff_radius
-valence_charge_density_cutoff_radius
-core_charge_density_cutoff_radius
-pseudo_cutoff_radius
-projector_radial_indices
-chi_function_radial_indices
-projector_coupling
-pseudo_energy_correction
-augmentation_coupling
-local_potential_real
-local_potential_fourier
-projector_real
-projector_fourier
-chi_function_real
-chi_function_fourier
-valence_charge_density_real
-valence_charge_density_fourier
-core_charge_density_real
-core_charge_density_fourier
-augmentation_real
-augmentation_fourier
+n_radials
+n_angulars
+has_quantity
+get_quantity
+cutoff_radius
+psp_quantity_evaluator
+psp_energy_correction
 
 PseudoPotentialIO.hankel_transform
 PseudoPotentialIO.build_interpolator_real
+PseudoPotentialIO.trapezoid
 PseudoPotentialIO.simpson
 PseudoPotentialIO.fast_sphericalbesselj
 ```

docs/src/tutorial.jl

@@ -1,42 +1,53 @@
 # # Tutorial
 #
-# Here, we'll take a look at the basic usage of PseudoPotentialIO.jl: searching for, loading, and investigating the contents of pseudopotentials.
+# Here, we'll take a look at the basic usage of PseudoPotentialIO.jl: searching for,
+# loading, and investigating the contents of pseudopotentials.
 
 using PseudoPotentialIO
 using CairoMakie
 using Colors
 
 # ## 1. Searching for pseudopotentials
-# PseudoPotentialIO provides easy access to a variety of pseudopotential families with zero setup using Julia [Artifacts](https://docs.julialang.org/en/v1/stdlib/Artifacts/), [LazyArtifacts.jl](https://github.com/JuliaPackaging/LazyArtifacts.jl), and [PseudoLibrary](https://github.com/JuliaMolSim/PseudoLibrary).
+# PseudoPotentialIO provides easy access to a variety of pseudopotential families with zero
+# setup using Julia [Artifacts](https://docs.julialang.org/en/v1/stdlib/Artifacts/),
+# [LazyArtifacts.jl](https://github.com/JuliaPackaging/LazyArtifacts.jl), and
+# [PseudoLibrary](https://github.com/JuliaMolSim/PseudoLibrary).
 # In order to list the available pre-bundled families, use `list_families`
 
-list_families(with_info=true)
+println(list_families())
 
-# Because families are downloaded lazily, detailed information on some of the families is missing (they need to be downloaded first).
+# Because families are downloaded lazily, detailed information on some of the families is
+# missing (they need to be downloaded first).
 # In order to download a pseudopotential family, you can load it using `load_family`.
 #
 # !!! note "Loading pseudopotential families"
 #     `load_family` can also load all the pseudopotentials in a local directory!
 #
-# You can then see an overview of which elements the family supports with `show_family_periodic_table`.
+# You can then see a summary of the family with`show_family_summary` and which elements it
+# supports with `show_family_periodic_table`.
 
-family = load_family("hgh_lda_upf");
+family = load_family_psp_files("hgh_lda_upf");
 show_family_periodic_table(family)
 
 # For more detailed information, use `show_family_list`.
 # You can restrict the output by providing a list of elements that you're interested in.
 
-show_family_list(family)  # Show all the pseudopotentials
-show_family_list(family; elements=["Ba", "Ti", "O"])  # Only show the pseudos for Ba, Ti, and O
+show_family_table(family)  # Show all the pseudopotentials
+show_family_table(family; elements=["Ba", "Ti", "O"])  # Only show the pseudos for Ba, Ti, and O
 
 # ## 2. Loading pseudopotential files
-# To load an individual pseudopotential _file_, use `load_psp_file`, specifying the family name or directory and the filename of the pseudopotential
+# To load an individual pseudopotential _file_, use `load_psp_file`, specifying the family
+# name or directory and the filename of the pseudopotential
 
 Ba_psp_file = load_psp_file("hgh_lda_upf", "Ba.pz-sp-hgh.UPF")
 
-# PseudoPotentialIO distinguishes between pseudopotential _files_ and the pseudopotentials themselves.
-# Structures like `HghFile`, `UpfFile`, and `Psp8File` correspond to pseudopotential file formats and make the quantities that these files contain directly available, with no unit conversion or processing.
-# For example, we can take a look at the author field from the header in the barium UPF pseudopotential we just loaded
+# PseudoPotentialIO distinguishes between pseudopotential _files_ and the pseudopotentials
+# themselves.
+# Structures like `HghFile`, `UpfFile`, and `Psp8File` correspond to pseudopotential file
+# formats and make the quantities that these files contain directly available, with no unit
+# conversion or processing.
+# For example, we can take a look at the author field from the header in the barium UPF
+# pseudopotential we just loaded
 
 Ba_psp_file.header.author
 
@@ -45,39 +56,55 @@ Ba_psp_file.header.author
 propertynames(Ba_psp_file)
 
 # ## 3. Loading pseudopotentials
-# Once you've decided that you would like to use a given pseudopotential for a calculation, either convert its `PsPFile` struct a corresponding pseudopotential structure
+# Once you've decided that you would like to use a given pseudopotential for a calculation,
+# either convert its `PsPFile` struct a corresponding pseudopotential structure
 
 Ba_psp_from_File = load_psp(Ba_psp_file)
 
 # , or load the pseudopotential structure directly from the file using `load_psp`
 
 Ba_psp_from_disk = load_psp("hgh_lda_upf", "Ba.pz-sp-hgh.UPF")
 
-# We can confirm that these pseudopotentials are identical
-Ba_psp_from_File == Ba_psp_from_disk
-
-# This procedure has processed the contents of the UPF file (on disk) or the `UpfFile` struct into a common and consistent data representation that PseudoPotentialIO uses for calculations. We can see that the contents have changed by looking at the property names of our new `UpfPsP` struct
+# This procedure has processed the contents of the UPF file (on disk) or the `UpfFile`
+# struct into a common and consistent data representation that PseudoPotentialIO uses for
+# calculations. We can see that the contents have changed by looking at the property names
+# of our new `UpfPsP` struct
 
 propertynames(Ba_psp_from_disk)
 
 # ## 4. Inspecting pseudopotential quantities
-# One thing that we might want to do with a processed pseudopotential is to visualize some of the quantities it contains.
+# One thing that we might want to do with a processed pseudopotential is to visualize some
+# of the quantities it contains.
 # Let's plot the Kleinman-Bylander projectors from a PseudoDojo barium pseudopotential
 
 Ba_psp = load_psp("pd_nc_sr_pbesol_standard_0.4.1_upf", "Ba.upf");
 let
     linestyles = [:solid, :dash, :dot]
     colors = Colors.JULIA_LOGO_COLORS
     fig = Figure(); ax = Axis(fig[1,1], xlabel="r [a₀]", ylabel="β(r)")
-    for l in angular_momenta(Ba_psp)                  # Iterate over each angular momentum 0:lmax
+
+    ## Iterate over each angular momentum 0:lmax
+    for l in angular_momenta(Ba_psp)
         color = colors[l+1]
-        for n in projector_radial_indices(Ba_psp, l)  # Iterate over each projector at l 1:nmax
+        ## Iterate over each projector at l 1:nmax
+        for n in 1:n_radials(NumericProjector(), Ba_psp, l)
             linestyle = linestyles[n]
-            r²βln = Ba_psp.β[l][n]                    # Projector multiplied by r²
-            i_rc_ln = lastindex(r²βln)                # Cutoff radius index
-            βln = r²βln ./ Ba_psp.r[1:i_rc_ln].^2     # Remove the r² prefactor
-            lines!(ax, Ba_psp.r[1:i_rc_ln], βln, label="|β[$l][$n]⟩",
-                   linestyle=linestyle, color=color)
+
+            ## Get a callable (interpolator or function) which lets us evaluate the
+            ## pseudopotential quantity in real- or Fourier-space
+            r²βln = psp_quantity_evaluator(
+                RealSpace(),      # Select real- or Fourier-space using [Real,Fourier]Space()
+                NumericProjector(),  # Select the quantity, here a non-local β-projector
+                Ba_psp, l, n)
+
+            ## Real-space evaluators for numeric pseudos fail outside their cutoff radius,
+            ## which we can find using the `cutoff_radius` function
+            rgrid = 0:0.01:cutoff_radius(NumericProjector(), Ba_psp, l, n)
+
+            ## Remove the r² prefactor
+            βln = r²βln.(rgrid) ./ rgrid.^2
+
+            lines!(ax, rgrid, βln, label="|β[$l][$n]⟩", linestyle=linestyle, color=color)
         end
     end
     axislegend()