Reactive scattering from a metal surface

Our implementation allows us to simulate vibrational (de-)excitation probability during reactive scattering events at metal surfaces for any diatomic molecule with a suitable model to describe energies and forces (and friction coefficients for MDEF simulations). Here, we investigate the reactive scattering of hydrogen on a Cu(111) metal surface as a prototypical example and we employ machine learning models for potential energy surface (MLIP) and electronic friction tensor (EFT).

To run this kind of simulation, a set of initial positions and velocities ($\mathbf{R}$ and $\mathbf{\dot{R}}$) with ro-vibrational quantum states $\nu$ and $j$ have to be generated (see EBK quantisation). With a specific ro-vibrational quantum state it is possible to compute different properties after molecular collision and energy transfer with the metal surface like the vibrational de-excitation probabilities discussed here.

To reproduce the state-to-state vibrational de-excitation probability results presented by [27] for this system, the same initial conditions can be generated with QuantisedDiatomic.generate_configurations setting the initial ro-vibrational quantum state to ($\nu=2, j=1$) as was explored in the original paper.

As shown earlier in the EBK documentation we are able to generate a semiclassically quantised distribution for a diatomic molecule on a collision course with a metal surface. In this example we follow the EBK example using machine learning potential to prepare our initial distribution and run our simulation.

Specifically, we have produced a set of initial conditions with different translational energy (translational_energy keyword) ranging from 0.04 to 0.4 eV, locating the hydrogen molecule 7 Å away from the metal surface (height keyword) with a normal incidence.

using NQCDynamics
using NQCDInterfASE
using Unitful
using NQCDynamics.InitialConditions: QuantisedDiatomic
using JLD2
using PythonCall

mace_calc = pyimport("mace.calculators")
io = pyimport("ase.io")
ase_units = pyimport("ase.units")

function mace_model(model_path, cur_atoms)
    calculator = mace_calc.MACECalculator(model_path=model_path, device="cpu", default_dtype="float32") # or "cuda"
    cur_atoms.set_calculator(calculator)
    model = ClassicalASEModel(cur_atoms)

    return model
end

model_path = "h2cu.model"
surface_structures_path = "surf_structures.xyz"

# settings
ν, J = 2, 1     # selected ro-vibrational quantum states  
nsamples = 10  # number of configurations      
Ek = 0.3u"eV"   # Translational energy [eV] ; range considered [0.0-0.4] eV
z = 7.0u"Å"     # Height [Å]  ; fixed at 7 Å
z = ang_to_au(z)
dof = 3

# Read surface structures (equilibrated at chosen temperature)
surface_ase = ase_io.read(surface_structures_path)
xatoms, xpositions, xcell = NQCDynamics.convert_from_ase_atoms(surface_ase[0])
set_periodicity!(xcell,[true,true,true])

# load PES model
model_pes = mace_model(model_path, atoms_path)
# set simulation
sim = Simulation(NQCDynamics.Atoms([:H,:H]), model_pes; cell=xcell)

configurations = QuantisedDiatomic.generate_configurations(sim, ν, J;
    samples=nsamples, 
    translational_energy=Ek, 
    height=z
)
v_h2 = first.(configurations)
r_h2 = last.(configurations)


# Re build the position and velocities for the whole system Cu56H2
n_atoms = length(atoms)
v = [zeros(dof,n_atoms) for i=1:length(r_h2)]
r = [zeros(dof,n_atoms) for i=1:length(r_h2)]

for i in 1:length(r_h2)
    surf_atoms, surf_positions, surf_cell = convert_from_ase_atoms(surface_ase[i])
    r[i][:,1:n_atoms-2] .= surf_positions
    r[i][:,n_atoms-1:n_atoms] .= r_h2[i]  

    v[i][:,1:n_atoms-2] .= austrip.(transpose(surface_ase[i].get_velocities().*ase_units.fs)*u"Å/fs")
    v[i][:,n_atoms-1:n_atoms] .= v_h2[i]
end


distribution = DynamicalDistribution(v, r, (dof,length(surface_ase[0])+2))
nothing # hide

In order to produce an unweighted distribution, the lateral and angular orientation are randomly selected within the unit cell. As an example of the spacial and orientation distribution generated with this module, a subset of data (300 configurations) is shown below. To run our production simulations, however, a set of 80,000 initial velocities and positions were used.

Data analysis and truncation function

Since we are interested in the dynamics only when the molecule is close to the surface, we can use a callback to terminate the simulation early to save us some time. This requires defining a function that returns true when we want the simulation to terminate. This means we can set our time span relatively long since we expect most simulations to terminate before reaching the time limit.

using Statistics: mean
using LinearAlgebra: norm

h2distance(p) = norm(p[:,ids_adsorbate[1]] .- p[:,ids_adsorbate[2]])

"Terminates simulation if returns `true`." 
mutable struct TrajectoryTerminator
    h2_indices
    ads_height_cutoff
    ads_dist_cutoff
    n_atoms_layer
end
function (t::TrajectoryTerminator)(u, t, integrator)::Bool
    R = get_positions(u)
    com_h2_z = minimum(R[3,t.h2_indices[1]:t.h2_indices[2]])
    top_surface_avg_z = mean(R[3,end-Int(t.n_atoms_layer)-1:end-2])
    zcom = au_to_ang(com_h2_z-top_surface_avg_z) # Convert vertical centre of mass to angstrom
    if zcom > t.ads_height_cutoff                         # Scattering event
        return true
    elseif au_to_ang(h2distance(R)) > t.ads_dist_cutoff   # Reactive event
        return true
    else
        return false
    end
end

termination_condition = TrajectoryTerminator(ids_adsorbate, 7.2, 2.25, 9)
terminate = DynamicsUtils.TerminatingCallback(termination_condition)
tspan = (0.0, 420.0u"fs")
nothing # hide

In this example, we consider the outcome a reactive event if the H-H bond length is larger than 2.5 Å in any point of during the trajectory and a scattering event if the molecule rebounds to a vertical distance from the metal surface greater than 8.1 Å.

MDEF with the LDFA

Now that we have set up the initial distribution and some of our simulation parameters, we can choose which form of friction we would like use. First, let's use the cube-based density implementation for LDFA friction provided by the FrictionProviders.jl. This takes a .cube file containing the electron density and will provide the friction during the dynamics. Here we initialize the MACE-based interatomic potential, together with the cube-based calculator.

# load cube EFT calculator
using FrictionProviders
density_model = CubeLDFA("path_to_cube.cube", cell)
model_eft = LDFAFriction(density_model, atoms; friction_atoms=ids_adsorbate)

Now we can pass all the variables defined so far to the Simulation and run multiple trajectories using run_dynamics.

model = CompositeFrictionModel(model_pes, model_eft)
sim = Simulation{MDEF}(atoms, model, cell=cell, temperature=300u"K")
ensemble = run_dynamics(sim, tspan, distribution; selection=1:nsamples,
    dt=0.1u"fs", output=OutputPosition, trajectories=nsamples, callback=terminate)

MDEF with machine-learned LDFA

Above, we used the cube-based LDFA interpretation of MDEF to perform the simulation. However, in this scheme, surface atoms have to stay fixed to match the cube densities. To include surface temperature effects machine learning models can be trained that predict densities at any surface configuration. Here, we run MDEF simulation employing two machine learning models, to predict adiabatic PES (based on MACE) and surface electron density (based on ACEpotentials.jl), utilizing FrictionProviders.jl, to run MDEF simulation.

using ACEpotentials
ace_model, ace_model_meta = ACEpotentials.load_model("ace_ldfa_model_path.json")
density_model = AdiabaticModels.ACEpotentialsModel(atoms, cell, ace_model) 
ace_density_model = AceLDFA(density_model)
model_eft = LDFAFriction(ace_density_model, atoms; friction_atoms=ids_adsorbate)

model = CompositeFrictionModel(model_pes, model_eft)

sim = Simulation{MDEF}(atoms, model, cell=cell, temperature=300u"K")
ensemble = run_dynamics(sim, tspan, distribution; selection=1:nsamples,
    dt=0.1u"fs", output=OutputPosition, trajectories=nsamples, callback=terminate)

MDEF with machine-learned TDPT (ODF)

Above, we used the LDFA interpretation of MDEF to perform the simulation. However, an alternative, TDPT (otherwise known as ODF) method can be used that provides full friction tensor. TDPT ML models can be incorporated in our simulation in a similar way as LDFA models, through FrictionProviders.jl. Here, we show how this can be done for ACEds (ACEfriction) models. Such models require the usage of JuLIP-type atoms. We can easily convert our ASE-type atoms into such format:

using ASE, JuLIP
atoms_ase_jl = ASE.ASEAtoms(atoms_ase)
atoms_julip = JuLIP.Atoms(atoms_ase_jl)

Having our atoms, ACEds (ACEfriction) EFT models can be then initialized and combined with previously loaded potential model.

using ACE
using ACEds.FrictionModels
using ACEds.FrictionModels: Gamma
aceds_model = ACEdsODF(read_dict(load_dict("aceds_odf_model_path.model")), Gamma, atoms_julip)
model_eft = ODFriction(ace_model; friction_atoms=ids_adsorbate)

model = CompositeFrictionModel(model_pes, model_eft)

sim = Simulation{MDEF}(atoms, model, cell=cell, temperature=300u"K")
ensemble = run_dynamics(sim, tspan, distribution; selection=1:nsamples,
    dt=0.1u"fs", output=OutputPosition, trajectories=nsamples, callback=terminate)

Visualisation

To show the effect of the truncation procedure, we have run 20 trajectories with and without the truncation function starting with an initial translation energy at 1.0 eV. For both figures, the total and kinetic energies are shown in the top panels along with the H-H distance and centre of mass z coordinate for each individual trajectory.