InitialConditions

InitialConditions

Functions and types for generating initial conditions for simulations.

run_advancedhmc_sampling(sim::Simulation, r, n_samples;
    target_acceptance=0.5, kwargs...)

Perform Hamiltonian Monte Carlo sampling for the simulation sim using AdvancedHMC.jl.

run_advancedmh_sampling(sim, r, steps, σ; movement_ratio=nothing, movement_ratio_internal=nothing, kwargs...)

Sample the configuration space for the simulation sim starting from r.

Total number of steps is given by steps and σ is the dictionary of step sizes for each species.

movement_ratio denotes the fraction of system moved each step. internal_ratio works as for movement_ratio but for the internal modes of the ring polymer. For movement_ratio = 0, every degree of freedom is moved at each step, if movement_ratio = 1, then nothing will happen.

If neither arguments are defined, default behaviour is to move one atom (and one ring polymer normal mode) per step on average.

Further kwargs are passed to AdvancedMH.sample to allow for extra functionality.

QuantisedDiatomic

This module exports two user facing functions:

• generate_configurations Creates a set of velocities and positions for diatomic molecule with specified vibrational ν and rotational J quantum numbers.

• quantise_diatomic Obtains vibrational ν and rotational J quantum numbers for a diatomic molecule with a given set of velocities and positions.

The central concept of this module is the EBK procedure which is nicely detailed here: [4]

Inspired by VENUS96: [28]

Randomly rotate each column of two 3*N matrix, same rotation for all columns.

calculate_diatomic_energy(model::AdiabaticModel, bond_length::Real;
    height=10, normal_vector=[0, 0, 1])

Returns potential energy of diatomic with bond_length at height from surface.

Orients molecule parallel to the surface at the specified height, the surface is assumed to intersect the origin. This requires that the model implicitly provides the surface, or works fine without one.

classical_rotation_energy(J::Union{Int, CartesianIndex}, config::Any, ind1::Union{Int, CartesianIndex}, ind2::Union{Int, CartesianIndex}, sim::Simulation)

Classical rotation energy of a rigid diatomic rotor in Hartree

classical_translational_energy(config::Any, ind1::Union{Int, CartesianIndex}, ind2::Union{Int, CartesianIndex}, sim::Simulation)

Returns the classical translational energy in Hartree

combine_slab_and_molecule(atom_indices, molecule, slab)

Revert the transformation separate_slab_and_molecule

Randomly orient molecule in space for a given bond length and radial momentum

find_total_energy(V::EffectivePotential, ν)

Returns the energy associated with the specified quantum numbers

generate_configurations(sim, ν, J; samples=1000, height=10, normal_vector=[0, 0, 1],
    translational_energy=0, direction=[0, 0, -1], position=[0, 0, height])

Generate positions and momenta for given quantum numbers

translational_energy, direction and position specify the kinetic energy in a specific direction with the molecule placed with centre of mass at position.

Keyword arguments height and normal_vector become relevant if the potential requires specific placement of the molecule. These allow the molecule to be placed at a distance height in the direction normal_vector when performing potential evaluations.

harmonic_vibration_energy(ν::Union{Int, CartesianIndex}, k::Float, ind1::Union{Int, CartesianIndex}=1, ind2::Union{Int, CartesianIndex}=2, sim::Simulation)

Vibrational energy of a harmonic oscillator with the force constant k and vibrational level ν.

quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve::BindingCurve;
show_timer=false, reset_timer=false,
height=10, normal_vector=[0, 0, 1], atom_indices=[1,2], max_translation=1) 
)

Quantise the vibrational and rotational degrees of freedom for the specified positions and velocities using the BindingCurve specified. A binding curve will be automatically generated if you do not supply one.

If the potential can be evaluated for the diatomic only, independent of position, supplying a Simulation for just the diatomic will speed up evaluation.

When evaluating the potential, the molecule is moved to height in direction normal_vector. If the potential is independent of centre of mass position, this has no effect. Otherwise, be sure to modify these parameters to give the intended behaviour.

If a Simulation with a PeriodicCell is supplied, periodic copies of the diatomic atoms will be used if positions are close to cell boundaries. Set max_translation to the radius of surrounding unit cells to search. (e.g. 1 if positions are already wrapped around cell boundaries)

quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; 
bond_lengths=0.5:0.01:5.0, 
height=10.0, 
surface_normal=[0, 0, 1.0], 
atom_indices=[1, 2], 
max_translation=1,
show_timer=false, reset_timer=false
)

Quantise the vibrational and rotational degrees of freedom for the specified positions and velocities.

If the potential can be evaluated for the diatomic only, independent of position, supplying a Simulation for just the diatomic will speed up evaluation.

When evaluating the potential, the molecule is moved to height in direction normal_vector. If the potential is independent of centre of mass position, this has no effect. Otherwise, be sure to modify these parameters to give the intended behaviour.

If a Simulation with a PeriodicCell is supplied, periodic copies of the diatomic atoms will be used if positions are close to cell boundaries. Set max_translation to the radius of surrounding unit cells to search. (e.g. 1 if positions are already wrapped around cell boundaries)

quantise_diatomic(sim::Simulation, v::Vector{Matrix}, r::Vector{Matrix};
bond_lengths=0.5:0.01:5.0, 
height=10.0, 
surface_normal=[0, 0, 1.0], 
atom_indices=[1, 2], 
show_timer=false, reset_timer=false
)

Quantise the vibrational and rotational degrees of freedom of multiple atomic configurations given as a vector of velocity matrices and a vector of position matrices.

If the potential can be evaluated for the diatomic only, independent of position, supplying a Simulation for just the diatomic will speed up evaluation.

When evaluating the potential, the molecule is moved to height in direction normal_vector. If the potential is independent of centre of mass position, this has no effect. Otherwise, be sure to modify these parameters to give the intended behaviour.

If a Simulation with a PeriodicCell is supplied, periodic copies of the diatomic atoms will be used if positions are close to cell boundaries. Set max_translation to the radius of surrounding unit cells to search. (e.g. 1 if positions are already wrapped around cell boundaries)

Specify show_timer=true for performance timings of the EBK quantisation process and reset_timer=true to see timings for each individual quantisation.

Pick a random bond length and corresponding radial momentum that matches the radial probability distribution.

separate_slab_and_molecule(atom_indices, r)

Get the coordinates of the molecule and slab separately.

sqrt_avoid_negatives(x)

Same as sqrt but returns zero(x) if x is negative. Used here just in case the endpoints are a little off.

MetropolisHastings

Sampling of the initial conditions using the Metropolis-Hastings Markov chain Monte Carlo method.

Included within is the ability to sample the canonical distribution for adiabatic classical and ring polymer systems.

Usage involves creating an instance of an AbstractSystem{MonteCarlo} and calling run_monte_carlo_sampling.

Parameters for Monte carlo simulations.

Container for storing simulation quantities

acceptance_probability(system::AbstractSystem{MonteCarlo})

Return the Metropolis-Hastings acceptance probability.

apply_random_perturbation!(system::AbstractSystem{MonteCarlo}, R::AbstractMatrix, atom::Integer)

Randomly perturb the xyz coordinates of a single atom.

assess_proposal!(system::AbstractSystem{MonteCarlo}, Rᵢ, Rₚ, output, i)

Update the energy, check for acceptance, and update the output.

propose_centroid_move!(system::RingPolymerSystem{MonteCarlo}, Rᵢ::Array{T, 3}, Rₚ::Array{T, 3}) where {T<:AbstractFloat}

Propose a move for the ring polymer centroid for one atom.

propose_move!(system::System{MonteCarlo}, Rᵢ::Matrix{T}, Rₚ::Matrix{T}) where {T<:AbstractFloat}

Propose simple cartesian move for a single atom.

propose_normal_mode_move!(system::RingPolymerSystem{MonteCarlo}, Rᵢ::Array{T, 3}, Rₚ::Array{T, 3}) where {T}

Propose a move for a single normal mode for a single atom.

run_main_loop!(system::System{MonteCarlo}, Rᵢ::Matrix, Rₚ::Matrix, output::MonteCarloOutput)

Main loop for classical systems.

Main loop for ring polymer systems.

run_monte_carlo_sampling(sim::AbstractSimulation, monte::MonteCarloParameters, R0)

Perform Monte Carlo sampling for the system defined by the sim and monte parameters.

From the initial positions specified R0 the system will be explored using the Metropolis-Hastings algorithm.

write_output!(output::MonteCarloOutput, Rᵢ::AbstractArray, energy::AbstractFloat, i::Integer)

Store the current configuration and associated energy.