Calculating potential energy surfaces, forces and non-adiabatic couplings with the Newns-Anderson Hamiltonian

Note

This tutorial focuses using NQCModels.jl to build a Newns-Anderson model of a molecule chemisorped to a surface. In the process we will also see how one can use NQCCalculators.jl to evaulate and update important properties of the model.

The Newns-Anderson Hamiltonian

The Newns-Anderson model is a lattice Hamiltonian that parameterises a small set of impurity states" coupled to a large number ofbath” states. Within the context of surface science the impurity" states are taken to represent the electronic eigenstates of a molecular adsorbate and thebath” states the band structure that characterises the material on which it is adsorbed. In second-quantised form of the model Hamiltonian reads \[ \hat{H} = \sum_{a}{\varepsilon_{a}\hat{n}_{a}} + \sum_{k}{\varepsilon_{k}\hat{n}_{k}} + \sum_{a,k}{V_{ak}(\hat{c}_{a}^\dagger\hat{c}_{k} + \hat{c}_{k}^\dagger\hat{c}_{a})}. \] Here, the subscript \(a\) runs over the adsorbate states and \(k\) the bath states. \(\varepsilon_a\) and \(\varepsilon_k\) are the energies of their respective states which the double sum over \(V_{ak}\) encodes the couplings between every given pair of surface and adsorbate states.

The model implemented in NQCModels.jl has a single adsorbate energy state coupled to a bath of energy states clustered around the fermi-level representing a either a metallic or semiconducting system (depending on the chosen discretisation scheme). The model is also taken to be in the ``wideband limit” where the \(k\) dependence of the couplings is dropped and the density of states (DOS) is approximated via variable sparsity of sampling in the bath state discritisation (Shenvi2009). In after these simplifications the Hamiltonian is able to be written down in a simplified matrix form, \[ \begin{bmatrix} \varepsilon_{a} & V_1 & V_2 & \cdots & V_k \\ V_1 & \varepsilon_1 & 0 & \cdots & 0 \\ V_2 & 0 & \varepsilon_2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ V_k & 0 & 0 & \cdots & \varepsilon_k \end{bmatrix} \]

Dependencies

The following packages are all that is required to complete this tutorial. If you don’t want to add them to your default Julia environment simply type ] activate 'NQCD Newns-Anderson tutorial' into the REPL before adding the packages.

using Unitful, UnitfulAtomic
using DelimitedFiles
using LinearAlgebra
using Combinatorics
using NQCModels
using NQCCalculators
using CairoMakie
using ColorSchemes, Colors

Building a Newns-Anderson model

Newns-Anderson model in NQCModels is given its other common name, AndersonHolstein() and requires a QuantumModel as input that paramterises both the quantum states of the adsorbate far from the surface and its coupling strength as a function of its distance from the surface. Alongside this impurity_model input, AndersonHolstein() also requries the input of a WideBandBathDiscretisation object which is used to construct the bath energies and is generated automatically by a collection of discretisation functions.

The first thing to do when building a Newns-Anderson model with NQCModels is to choose what kind of impurity model and bath discretisation you would like. Some sensible defaults are the ErpenbeckThoss() model and a ShenviGaussLegendre() discretisation.

impurity_model = NQCModels.QuantumModels.ErpenbeckThoss(;Γ=0.01) #Gamma controls the coupling strength with the surface and can be properly parameterised with DFT, but 0.01 is a sensible default for the purpose of testing

bandmin = austrip(-270 * u"eV") #austrip converts unitful quantities into atomic units, the default units for NQCD
bandmax = austrip(270 * u"eV")
nstates = 100 #this controls the number of bath states, after including the additional single adsorbate state the AndersonHolstein Hamiltonain becomes a (nstates + 1, nstates + 1) matrix

bath = ShenviGaussLegendre(nstates, bandmin, bandmax)

NAmodel = AndersonHolstein(impurity_model, bath)

This code creates an instance of the Newns-Anderson model, stored as a julia struct. This struct contains all of the information you’ve input along with the number of electrons (calculated automatically by the discretisation method), the fermi level (defaults to 0.0 but can be set via the kwarg fermi_level) and an additional quantity saved as couplings_rescale that we will cover later on in the tutorial. For now we shall see what we can achieve with the NQCModels and the NAmodel struct.

NQCCalculators and Adiabatic Potential Energy Surfaces

The most obvious first step is to construct the matrix representation of the Newns-Anderson model that we saw above this can be using the NQCModels function potential(). This function takes the model we have already built along with a position argument that encodes the distance between the adsorbate and the surface. This position is wrapped in hcat() as in general NQCD functions take position as a matrix.

Hamiltonian_mat = NQCModels.potential(NAmodel, hcat(1.0))

This matrix should have the same structure as the matrix shown above, a structure characteristic of the Anderson model in the diabatic basis. Finding the adiabatic eigenstates of the system at the given height can be easily achieved using the LinearAlgebra.eigen() function. Of more interest in dynamics is the values and derivatives of these eigenstates at range of heights, the potential energy surfaces (PESs) of the system.

Before plotting the potential energy surfaces it would be useful to briefly discuss the applications of NQCModels to dynamics simulations and in doing so mention how it interfaces with the suite of dynamics software available in NQCDynamics.jl. The interface between the dynamics algorithms and the library of Hamiltonians and analytic potentials stored in NQCModels is handled by a third package called NQCCalculators.jl, this package uses the model to build a new type of struct called a cache. This cache contains the current instances of all the values that can be calculated from the model and can be automatically updated at each time-step using the update_cache!() function. To create a cache one uses the Create_Cache() function which takes as input the model, the number of atoms subject to the model and a type T that should almost always be set to Float64

NAcache = Create_Cache(NAmodel, 1, Float64)
NAcache.potential
update_cache!(NAcache, hcat(1.0))
NAcache.potential

Building Adiabatic Potential Energy Surfaces

We are now ready to generate a plot of the adiabatic potnetial energy surfaces of our model, for both ground and excited states. We have already shown how to construct a model of an adsorbate on a metal surface, so for this section we will instead plot the adiabatic eigenstates of a molecule chemisorbed to a semiconductor. For this we employ the LuHertlMaurer() model, an analytic potential fit to DFT data for H/Ge(111) and the GapGaussLegendre() discretisation which uses a modification of the GaussLegendre algorithm to create a set of bathstates for a surface with bandgap.

impurity_model = NQCModels.QuantumModels.LuHertlMaurer()
bandmin = austrip(-25 * u"eV")
bandmax = austrip(25 * u"eV")
gapwidth = austrip(0.49 * u"eV")
nstates = 150

bath = GapGaussLegendre(nstates, bandmin, bandmax, gapwidth)

NAmodel = AndersonHolstein(impurity_model, bath; couplings_rescale = 2.5)
NAcache = Create_Cache(NAmodel, 1, Float64)

Now we have built the model and constructed a cache for it, the steps required to plot the adiabatic excited states are the following:

Generate a ground state electronic occupation vector, a fock state where the first \(N_{\text{el}}\) elements are 1 and the rest of it are 0.
Generate a large number of combinatoric permutations of excited states with a fixed number of electrons \(N_{\text{el}}\) and within some given active space in order to obtain a random collection of low energy excited states.
Diagonalise the Newns-Anderson Hamiltonian \(H_{\text{NA}}\) at a given height \(x\) and sum the eigenvalues \(\varepsilon_i\) of occupied states to obtain the total many-body energy of a given excited state.
Looping over each of the generated excited states and a variety of positions \(x\), repeat step 3 to produce the excited state potential energy surfaces.

These steps are all performed by the adiabatic_surfaces() function below which will output a sorted list of PESs in the form of a Vector{<:Vector} (a vector of vectors).

"""
    adiabatic_surfaces computes the potential energy surfaces, the nonadiabatic couplings between them and the forces due to their derivatives. 
    The surfaces are sorted according to their energy farthest from the surface.

    Output:
        energies_raw : sorted adiabatic PES energies
        U_0s : the state independent, external potential at each point
        configurations : configuration vectors sorted according to the adiabatic PES energy at the last grid point
        couplings: the non-adiabatic coupling matrix between the single-electron adiabatic states at each evaluated position
        GS_forces: the forces due to the groundstate PES at each evaluated position

"""
function adiabatic_surfaces(NAcache, x_ang::AbstractArray{Float64})

    ## The number of bath states
    nstates = NQCModels.nstates(NAcache)

    ## The number of electrons 
    nelectrons = NQCModels.nelectrons(NAcache) #NQCModels.nelectrons(model) -> Int(nstates/2)

    x = austrip.(x_ang * u"Å")
    matrix_x = hcat.(x)
    U_0s= zero(x)
    ∂U_0s= zero(x)
    eigs = [zero(NAcache.eigen.values) for _ in CartesianIndices(x)]
    
    #NAcache.nonadiabatic_coupling is stored as a Matrix{<:Matrix}, here the exterior matrix is 1x1 so below we only access the internal structure to store the values
    couplings = [zero(NAcache.nonadiabatic_coupling[1]) for _ in CartesianIndices(x)]
    derivative_eigen = [zero(NAcache.eigen.values) for _ in CartesianIndices(x)]

    for (i, r) in enumerate(matrix_x)
        update_cache!(NAcache, r)
        U_0s[i] = ustrip(auconvert(u"eV", NQCModels.state_independent_potential(NAcache.model, r)))
        ∂U_0s[i] = ustrip(auconvert(u"eV/Å", NQCModels.state_independent_derivative(NAcache.model, r)[1]))
        eigs[i] .= NAcache.eigen.values # length of eigs is 200(grid points); eigs[1] has 101 eigenstates
        couplings[i] .= NAcache.nonadiabatic_coupling[1]
        derivative_eigen[i] = eigvals(NAcache.derivative[1]) #-sum(eigen(NAcache.derivative[1]).values[1:nelectrons])
        #println(i)
    end

    # groundstate configuration
    config = vcat(ones(nelectrons), zeros(nstates-nelectrons))

    perms = multiset_permutations(config, nstates)
    configurations = []
    forces = []
    energies = []
    sortpoints = []
    """
        This setting makes sure that each PES has the same second quantization configuration
    """
    # loop over all the configurational permutations from groundstate to higher excited states
    for config in Iterators.take(perms, 20000)
        E = [sum(e[config .!= 0]) for e in eigs] # sum the eigenvalues of the occupied states at each grid point
        F = [-sum(d[config .!= 0]) for d in derivative_eigen] # calculate forces via summing across the eigenvalue derivatives of occupied eigenstates at each grid point
        push!(energies, ustrip.(auconvert.(u"eV", E)))
        push!(forces, ustrip.(auconvert.(u"eV/Å", F)))
        push!(sortpoints, E[end]) # record the adiabatic PES energy at the last grid point
        push!(configurations, config)
    end
    perm = sortperm(sortpoints)
    energies_raw = energies[perm] # sort the energies according to the adiabatic PES energy at the last grid point
    forces_raw = forces[perm]
    configurations = configurations[perm] # sort the configuration according to the adiabatic PES energy at the last grid point
    return energies_raw, U_0s, configurations, couplings, forces_raw, ∂U_0s
end

** Plotting the Adiabatic PESs **

Now that a function has been written to find the ground and excited state adiabatic potential energy surfaces, it’s time to set up some functions to visualise the surfaces we can calculate.

function plot_surfaces!(NAcache, ax, x_ang::AbstractArray{Float64}, gapwidth, groundstate_align_zero::Bool, reduce_PES::Bool)
    energies_raw, U_0s = adiabatic_surfaces(NAcache, x_ang)[1:2]

    if groundstate_align_zero
        groundstate_energy_end = energies_raw[1][end] + U_0s[end]
        energies_raw = broadcast(.-, energies_raw,  groundstate_energy_end) #broadcast to every PESs
    end

    if reduce_PES
        plotted_indices = vcat(collect(1:2),collect(3:12:20),collect(21:50:2000))
    else
        plotted_indices = collect(1:2000)
    end
    # extract only the indexed PESs
    energies = energies_raw[plotted_indices] # energies_raw[1] groundstate 

    @info "The depth of the ground state well is $(minimum(energies[1] .+ U_0s) - (energies[1] .+ U_0s)[end]) eV"

    # first excited state - ground state
    f_minus_g = energies[2] .- energies[1]
    gap_max = maximum(max.(f_minus_g))
    gap_min = minimum(min.(f_minus_g))

    @info "The maximum gap of ground state and first excited state is $gap_max eV"
    @info "The minimum gap of ground state and first excited state is $gap_min eV"

    color_range = []
    ## Plot PES ##
    for e in energies
        c = (e[end] - energies[1][end]) / (energies[end][end] - energies[1][end])
        push!(color_range, c)
    end
    colors = [ColorSchemes.get(colormap, 1-c) for c in color_range]

    excited_PES_energy_5A = []
    for (i, e) in Iterators.reverse(enumerate(energies))
        y = e .+ U_0s # NA Hamiltonian = adiabatic eigenvalues + U_0
        lines!(ax, x_ang, y, color=colors[i], linewidth= i==1 ? 4 : 1)

        if y[end] < 1.92 && y[end] > 0.4
            push!(excited_PES_energy_5A, y[end])
        end
    end

    ground_end = energies[1][end] + U_0s[end]

    DFT = readdlm("resources/Newns-Anderson-tutorial/DFT_HGe111_restatom.txt")[:,3:4]
    DFT[:,2] = DFT[:,2] .- DFT[end,2]
    aglinment =  ground_end - DFT[end,2] # adiabatic ground state - DFT ground state
    y_DFT = DFT[:,2] .+ aglinment # align the DFT ground state to the adiabatic ground states

    lines!(ax, DFT[:,1], y_DFT, label="DFT PES",linestyle=:dash, color=:blue, linewidth=3)

    ## bandgap##
    hlines!(ax, [ground_end + ustrip(auconvert(u"eV", gapwidth))], color=:black, linestyle=:dash, label="Band Gap", linewidth=3)
    hlines!(ax, [ground_end], color=:black, linestyle=:dash,linewidth=3)
end

function CustomAxis(f; kwargs...)
        my_ax = Axis(f; kwargs...)
        extra_ax = Axis(f; xaxisposition=:top, yaxisposition=:right, kwargs...)
        
        linkaxes!(my_ax, extra_ax)
        hidespines!(extra_ax)
        hidedecorations!(extra_ax, ticks=false, minorticks=false)

        return my_ax
    end

function plot_adiabatic_PES(NAcache, positions::AbstractVector, gapwidth::Float64; groundstate_align_zero::Bool=true, reduce_PES::Bool=false)
    ## set plotting theme ##
    colorscheme = ColorScheme(parse.(Colorant, ["#045275", "#089099", "#7CCBA2", "#FCDE9C", "#F0746E", "#DC3977", "#7C1D6F"]))
    colormap = ColorScheme(parse.(Colorant, ["#FCE1A4","#FABF7B","#F08F6E","#E05C5C","#D12959","#AB1866","#6E005F"]))

    ## create plotting figure ##
    ylimitslow = -2
    ylimitsup = 4

    fig = Figure(figure_padding=(5, 5, 5, 5), size=(480, 720/MathConstants.golden),fontsize=17)

    # PES axis
    ax1 = CustomAxis(fig[1,1]; xlabel="x / Å", limits=(1, 6, ylimitslow, ylimitsup), ylabel = "PES: Energy /eV", xgridvisible = false, ygridvisible = false)
    ax1.title = ""

    # plot DFT groundstate and adiabatic surfaces
    plot_surfaces!(NAcache, ax1, positions, gapwidth, groundstate_align_zero, reduce_PES)

    # label
    Legend(fig[1,1], ax1, tellwidth=false, tellheight=false, valign=:top, halign=:right, margin=(5, 5, 5, 5), orientation=:horizontal)
    return fig
end

Now we have constructed the model and written the functions needed to perform the calculations and plotting of the potential energy surfaces, all that remains is to set the range of positions over which the calculation is run and call the functions.

x_ang = range(0.8, 6, length=200) #range of distances in Å for the plot to be calculated over
plot_adiabtic_PES(NAcache, x_ang, gapwidth)

Calling these functions together should return this plot:

Adiabatic potential surfaces of \(H_{\text{NA}}\) configurated with the Lu-Hertl-Maurer model and Gapped the Gauss Legendre discretisation with a 0.49 eV bandgap

Plotting the Forces and Non-Adiabatic Couplings

This is no the only thing we can do with the adiabatic_surfaces!() function, however. Using update_cache!() also gives the user direct access to the Hamiltonian derivative which can be used to calculate forces and non-adiabatic couplings, values which are vital to propagating non-adiabatic dynamics. To show how these values change over the trajectory we can write a simple plotting script to visualise the net forces due to the potential energy surfaces and the nonadiabtic couplings between the ground and first excited states.

function plot_derivatives_and_couplings(NAcache, positions)
    ## Create Figure ##
    fig = Figure(figure_padding=(5, 5, 5, 5), size=(480, 720/MathConstants.golden),fontsize=17)

    # generate axis
    ylimitslow = 0.0
    ylimitsup = 3.0
    ax1 = CustomAxis(fig[1,1]; limits=(0.8, 6.0, -3.0, 3.0), ylabel = "Forces (eV/Å)", xgridvisible = false, ygridvisible = false)
    ax2 = CustomAxis(fig[2,1]; xlabel="x (Å)", limits=(0.8, 6.0, 0.0, 0.04), ylabel = "Couplings (Å⁻¹)", xgridvisible = false, ygridvisible = false)
    ax1.title = "Forces due to PESs"
    ax2.title = "Non-adiabaitc Couplings"

    ## Calculate Forces and Couplings##
    _, _, _, couplings, forces_raw, ∂U_0s = adiabatic_surfaces(NAcache, positions)
    plotted_indices = vcat(collect(1:2),collect(3:12:20),collect(21:50:2000))
    forces = forces_raw[plotted_indices]

    ## Plot Forces ##
    color_range = []
    for f in forces
        c = (f[end] - forces[1][1])
        push!(color_range, c)
    end
    colors = [ColorSchemes.get(colormap, 1-c) for c in color_range]

    for (i, f) in Iterators.reverse(enumerate(forces))
        f = f .+ ∂U_0s
        lines!(ax1, positions, f, color=colors[i], linewidth = 0.2)#i==1 ? 4 : 1)
        println(length(forces))
    end

    ## Plot Couplings ##
    c = []
    for (x,_) in enumerate(couplings)
        push!(c, couplings[x][1,2])
    end
    c = ustrip.(auconvert.(u"Å^-1", c))
    println(c)
    lines!(ax2, positions, c, color=colors[18], linewidth = 1.0)
    

    return fig
end

plot_derivatives_and_couplings(NAcache, x_ang)

These plotting routines should produce the following plot:

Non-adiabatic couplings and Forces due to potential energy surfaces of \(H_{\text{NA}}\) configured with the Lu-Hertl-Maurer model and Gapped the Gauss Legendre discretisation with a 0.49 eV bandgap

Note

Reference:

Philip W. Anderson. Localized magnetic states in metals. Physical Review 124.1 (1961): 41.
Newns, D. M. Self-Consistent Model of Hydrogen Chemisorption. Physical Review 178, no. 3 (1969): 1123–35.
Xuexun Lu, Nils Hertl, and Reinhard J. Maurer. A Haldane–Anderson Hamiltonian Model for Hyperthermal Hydrogen Scattering from a Semiconductor Surface. arXiv:2508.13360 (2025).
Xuexun Lu, Matthew Larkin and Nils Hertl. System-bath models. NQCD Documentation.