API

Settings for the dynamics.

Dynamics(dt :: T = 1.0 #femtosecodns, totaltime :: T # In femtosecodns, algorithm :: String = "generalized-verlet", kongliuratio :: T = 1.0, verbose :: Bool = true, evolvecorrelators :: Bool = true, correlated :: Bool = true, seed :: Int64, N :: Int64,

# Save the data each
save_each :: Int64 = 1,
save_filename :: String,
save_correlators :: Bool)

The settings for the simulation. dt and total_time are in femtoseconds

• dt is the time step [in femtoseconds]
• total_time is the total simulation time [in femtoseconds]
• algorithm is the integration algorithm to use. Either "generalized-verlet" or "semi-implicit-verlet"
• kong_liu_ratio is the ratio exploits the importance sampling.
• verbose is a flag to print out information during the simulation
• evolve_correlators is a flag to evolve the correlators <RR>, <PP> and <RP>
• correlated if true, the correlated approach for computing ensemble averages is used
• seed is the seed for the random number generator
• N is the number of stochastic configurations
• save_filename is the name of the file where to save the data
• save_correlators is a flag to save the correlators information
• save_each is the number of steps between each save of the data

Electric Field.

Base.@kwdef mutable struct ElectricField{T <: AbstractFloat} 
    fun :: Function #Time in fs, unit 
    Zeff :: Matrix{T} 
    edir :: Vector{T} #Must have unit norm
    eps :: Matrix{T}
end

This structure contains the information about the external IR electric field.

• fun is the function that describes the electric field as a function of time
• Zeff is the effective charge matrix
• edir is the direction of the electric field
• eps is the dielectric constant matrix

Ensemble average information stocastic displacements and forces have index i, j, coordinate i, configuration j

Ensemble(rho0 :: WignerDistribution{T}, positions :: Matrix{T}, # Positions are multiplied by the square root of the masses energies :: Vector{T}, forces :: Matrix{T}, # Forces are divided by the squareroot of masses stress :: Matrix{T}, # Stress in eV/A^3 sschaenergies :: Vector{T}, sschaforces :: Matrix{T}, nconfigs :: Int32, weights :: Vector{T}, temperature :: T, correlated :: Bool, y0 :: Matrix{T}, #unitcell :: Matrix{T})

The WignerDistribution.

It can be initialized using the method $init_from_dyn$. Alternatively, one can initialize a generic one with the field

WignerDistribution(n_atoms; type = Float64, n_dims=3)

The structure contains the following fields:

Base.@kwdef mutable struct WignerDistribution{T<: AbstractFloat}
    R_av    :: Vector{T}
    P_av    :: Vector{T}
    masses  :: Vector{T}
    n_atoms :: Int
    n_modes :: Int
    RR_corr :: Matrix{T}
    PP_corr :: Matrix{T}
    RP_corr :: Matrix{T}
    alpha :: Matrix{T}
    beta :: Matrix{T}
    gamma   :: Matrix{T}
    
    settings :: GeneralSettings
    λs :: Vector{T}
    λs_vect :: Matrix{T}
    evolve_correlators :: Bool
    cell :: Matrix{T}
    atoms :: Vector{String}
end

Note that all the variable here are with a tilde (mass rescaled) So that we can use linear-algebra on them quickly.

calculate_ensemble!(ensemble::Ensemble{T}, crystal) where {T<:AbstractFloat}

Evaluate energies, forces, and stresses for each configuration in an Ensemble using a given calculator ( assigned to crystal). Results are written back into the ensemble in place.

Arguments

• ensemble::Ensemble{T}: ensemble of stochastic configurations. Updated in place.
• crystal: external calculator object (e.g. ASE/PyCall wrapper) providing:
  • crystal.positions (settable Cartesian coordinates)
  • get_potential_energy()
  • get_forces() (returns N_atoms × 3 array)
  • get_stress() (returns stress tensor

Units

• Energies are converted to Rydberg (* CONV_RY).
• Forces are flattened, divided by √masses, and converted with CONV_RY/CONV_BOHR.
• Stress is assumed to be in eV/ų and stored directly.

Returns

• The updated ensemble (modification is in place, return value is mainly for chaining).

classicevolution!(Rs::Vector{T}, Ps::Vector{T}, dt::T, clfor, part) where {T <: AbstractFloat}

Classic dynamics.

eulerstep!(wigner:: WignerDistribution{T}, dt :: T, avgfor :: Vector{T}, d2V_dr2 :: Matrix{T}) where {T <: AbstractFloat}

Explicit Euler method. Implemented for testing and didactic use. MUST BE AVOIDED, it is unconditionally unstable.

function fixedstep!(wigner :: WignerDistribution{T}, dt :: T, avgfor :: Vector{T}, d2V_dr2 :: Matrix{T},part ) where {T <: AbstractFloat}

Evolves the quantum state freezing the nuclear degrees of freedom. This corresponds to a classic evolution in which however the forces are computed as ensemble averages using the initial quantum distribution.

semiimpliciteulerstep!(wigner:: WignerDistribution{T}, dt :: T, avgfor :: Vector{T}, d2V_dr2 :: Matrix{T}) where {T <: AbstractFloat}

Experimental.

generalized_verlet_step!(rho :: WignerDistribution{T}, dt :: T, avg_for :: Vector{T}, d2V_dr2 :: Matrix{T}, bc0 :: Vector{T} ,part ) where {T <: AbstractFloat}

Evolves the means and second–order correlators of a WignerDistribution by one generalized-Verlet step. This is algorithm has O(dt^3) error in time step. The algorithnm is illustrated in the work from Libbi et al., npj Computat. Mater. 11, 102 (2025), and consists of the explicit evolution of the variables ⟨R⟩, ⟨P⟩ and ⟨RR⟩, and of an implicit evolution for the variables ⟨PP⟩ and ⟨RP⟩. The latter is performed iteratively.

Arguments

• rho::WignerDistribution{T}: state (means & correlators), updated in place.
• dt::T: time step.
• avg_for::Vector{T}: current average force ⟨f⟩.
• d2V_dr2::Matrix{T}: Hessian (∂²V). Assumed symmetric in practice.
• part::Int: stage selector (1 = predictor, 2 = corrector).

Given the following variables name ⟨R⟩ = Rav, ⟨P⟩ = Pav, ⟨f⟩ = avgfor, ⟨RR⟩ = RRcorr, ⟨PP⟩ = PPcorr, ⟨RP⟩ = RPcorr, ⟨d²V⟩ = d2V_dr2

• part == 1 (predictor):

  • Means (half-kick)
    • ⟨R⟩ ← ⟨R⟩ + dt⟨P⟩ + ½dt²*⟨f⟩
    • ⟨P⟩ ← ⟨P⟩ + ½dt⟨f⟩
  • Correlators
    • let K = ⟨RR⟩ * ⟨d²V⟩
    • ⟨RR⟩ ← ⟨RR⟩ + dt(⟨RP⟩ + ⟨RP⟩ᵀ) − ½dt²(K + Kᵀ) + dt²⟨PP⟩
    • save ⟨RP⟩₀ = ⟨RP⟩
    • ⟨RP⟩ ← ⟨RP⟩ + ½dt(⟨PP⟩ − K)
    • ⟨PP⟩ ← ⟨PP⟩ − ½dt(⟨d²V⟩⟨RP⟩₀ + (⟨d²V⟩⟨RP⟩₀)ᵀ)

• part == 2 (corrector / finishing kick):

  • Means
    • ⟨P⟩ ← ⟨P⟩ + ½dt⟨f⟩
  • Correlators (fixed-point refinement of ⟨PP⟩ and ⟨RP⟩)
    • pre-update: ⟨RP⟩ ← ⟨RP⟩ − ½dt(⟨RR⟩*⟨d²V⟩)
    • iterate a few times (typically 4):
      • ⟨RP⟩₁ ← ⟨RP⟩ + ½dt⟨PP⟩₁
      • ⟨PP⟩₁ ← ⟨PP⟩ − ½dt(⟨d²V⟩⟨RP⟩₁ + (⟨d²V⟩⟨RP⟩₁)ᵀ)
    • write back final ⟨PP⟩, ⟨RP⟩

generate_ensemble!(N, ensemble:: Ensemble{T}, wigner_distribution :: WignerDistribution{T}) where {T <: AbstractFloat}

Generate or update an ensemble of stochastic configurations for nuclear dynamics, based on a Wigner distribution.

Arguments

• N::Int: Number of configurations to generate (must be even).
• ensemble::Ensemble{T}: Target ensemble structure. Arrays are either updated in place if their size already matches N, or reallocated otherwise.
• wigner_distribution::WignerDistribution{T}: Wigner distribution that encodes the covariance matrices, eigenvalues, and eigenvectors used to generate the ensemble.

Behavior

• If ensemble.correlated == true, correlated configurations are generated from ensemble.y0 using the square root of the correlator matrix.
• Otherwise, random Gaussian vectors are drawn (using randn) and scaled by √λ or 1/√λ depending on whether ρ.evolve_correlators is set.
• In MPI runs, random numbers are generated on rank 0 and broadcast to all ranks so that every processor works with the same ensemble.

Returns

• The updated ensemble (modifications are in place).

get_average_energy(ensemble :: Ensemble{T}) where {T <: AbstractFloat}

Computes the ensemble average of the potential energy.

get_average_forces(ensemble :: Ensemble{T}) where {T <: AbstractFloat}

Computes the ensemble average of the forces.

get_average_stress(ensemble :: Ensemble{T}, wigner :: WignerDistribution{T}) where {T <: AbstractFloat}

Computes the ensemble average of the stress. The result is returned in Ry/Bohr^3

function getaverages!(avgfor :: Vector{T}, d2vdr2 :: Matrix{T}, ensemble :: Ensemble{T}, wignerdistribution :: WignerDistribution{T}) where {T <: AbstractFloat}

Evaluate the average force and d2v_dr2 from the ensemble and the wigner distribution. The weights on the ensemble should have been already updated.

get_classic_forces(wigner_distribution:: WignerDistribution{T}, crystal) where {T <: AbstractFloat}

Compute the forces at the centroids, i.e., the classical forces acting on the crystal.

get_kong_liu(ensemble :: Ensemble{T}) where {T <: AbstractFloat}

Compute the Kong-Liu ration

get_random_y(N, N_modes, settings :: Dynamics{T}) where {T <: AbstractFloat}

Generates a random matrix of numbers sampled from a normal distribution. Each row corresponds to a phonon mode, and contains N random configurations (or N/2 if evenodd = true). For parallel execution, the matrix ymui is broadcast so that all processors use the same random variables.

get_total_energy(ensemble :: Ensemble{T}, wigner_distribution :: WignerDistribution{T}) where {T <: AbstractFloat}

Computes the total quantum energy, as a sum of the quantum potential energy (avgene), quantum kynetic energy (Krho), and classic kynetic energy (K_point)

function get_λs(RR_corr :: Matrix{T}) where {T <: AbstractFloat}

Get the eigenvalues of <RR>.

init_ensemble_from_python(py_ensemble, settings::Dynamics{T}) where {T<:AbstractFloat}

Construct a Julia Ensemble object from a Python-side ensemble exported via PyCall.

This routine converts positions, forces, stresses, and energies from the Python py_ensemble (typically created by CellConstructor / SSCHA Python codes) into the internal representation used by QuantumGaussianDynamics.

Arguments

• py_ensemble: A Python object with fields:
  • current_dyn – Python dynamical matrix object
  • T0 – temperature (in K)
  • N – number of configurations
  • xats – atomic random displacements, shape (N, Nat, 3)
  • energies – raw energies for each configuration
  • sscha_energies – SSCHA energies
  • forces – atomic forces, shape (N, Nat, 3)
  • sscha_forces – SSCHA forces, shape (N, Nat, 3)
  • stresses – stress tensors, shape (N, 3, 3)

Conversions performed

• Positions: (Nat, 3, N) → (3*Nat, N). Converted to Bohr and multiplied by √masses.
• Energies: copied directly from py_ensemble.energies and sscha_energies.
• Forces: (Nat, 3, N) → (3*Nat, N). Converted to Ry/Bohr and divided by √masses.
• Stresses: (3,3) tensors packed into 6×N Voigt form: 1=xx, 2=yy, 3=zz, 4=yz, 5=xz, 6=xy.
• Weights: initialized to 1.0 for all configurations.
• Random seed: if settings.seed != 0, reseeds RNG for reproducibility.
• y0: initial random Gaussian vectors, generated with get_random_y. If settings.correlated==false, filled with zeros of the correct shape.

Returns

• Ensemble{T}: A new ensemble populated with data from Python, ready for use in QuantumGaussianDynamics simulations.

init_from_dyn(dyn, TEMPERATURE::T, settings::Dynamics{T}) where {T<:AbstractFloat}

Initialize a WignerDistribution object starting from a dynamical matrix.

This routine constructs the equilibrium nuclear quantum state (in terms of positions, momenta, and correlation matrices) associated with the phonon modes of the provided dynamical matrix at a given temperature.

Arguments

• dyn: A dynamical matrix object from SSCHA calculations
• TEMPERATURE::T: The target temperature.
• settings::Dynamics{T}: Simulation settings.

Details

• Builds a supercell from dyn and computes the number of modes and atoms.
• Diagonalizes the supercell dynamical matrix to obtain phonon frequencies and eigenvectors (polarizations).
• Computes Gaussian width parameters (alpha, beta) and correlators (RR_corr, PP_corr, RP_corr) at the specified temperature.
• Initializes average positions (R_av) and momenta (P_av), rescaled by atomic masses.
• Removes translational acoustic modes from the correlators/eigenvalues.

Returns

• rho::WignerDistribution: An initialized Wigner distribution representing the quantum nuclear state corresponding to the given dynamical matrix and temperature.

Notes

• Frequencies are assumed to be in Rydberg units.
• Cell vectors and coordinates are converted to Bohr units.
• If settings.evolve_correlators == false, eigen-decomposition of alpha is used; otherwise, the RR_corr matrix is diagonalized (deprecated).

integrate!(wigner::WignerDistribution{T},
           ensemble::Ensemble{T},
           settings::Dynamics{T},
           crystal,
           efield::ElectricField{T}) where {T<:AbstractFloat}

Time-integrate the TDSCHA equations of motion. This mutates wigner (⟨R⟩, ⟨P⟩, correlators) and ensemble (weights, possibly regenerated configurations), and writes diagnostics to disk at the cadence settings.save_each.

Integration scheme

• Uses the algorithm specified by settings.algorithm: "euler" | "semi-implicit-euler" | "semi-implicit-verlet" | "fixed" | "generalized-verlet" | "none".
• At each step:
  1. Build average force and Hessian‐like term from the current ensemble: get_averages!.
  2. Add external field forces get_external_forces(t_fs, efield, wigner).
  3. Advance quantum variables with the selected integrator.
  4. Advance classical centroids (Rs, Ps) via classic_evolution!.
  5. Update eigenstructure of RR and reweight the ensemble update_weights!.
  6. If the Kong–Liu ratio falls below settings.kong_liu_ratio, regenerate the ensemble and recompute forces.

Files written (rank 0 only)

All prefixed by settings.save_filename * "settings.dt-settings.total_time-settings.N": All in Rydberg units, NOT scaled by masses.

• .pos: ⟨R⟩ and ⟨P⟩ plus total energy.
• .rho: RR correlator.
• .for: average force and “classical” force, and d²V/dR².
• .cl: classical centroids and classical energy.
• .ext: external field force.
• .str: average Voigt stress.
• .bet: PP correlator.
• .gam: RP correlator.

Returns

nothing. The function mutates wigner, ensemble, and writes output files.

Errors

Throws ArgumentError if settings.algorithm is unknown.

Example

```julia w = WignerDistribution{Float64}(; ... ) # prepared from a phonon dyn ens = Ensemble{Float64}(; ... ) # configurations & weights ef = ElectricField{Float64}(; fun=t->exp(-(t-50)^2/200), Zeff=Z, edir=û, eps=ϵ) integrate!(w, ens, Dynamics(dt=1.0, totaltime=200.0, algorithm="generalized-verlet", kongliuratio=0.8, verbose=true, correlated=true, seed=1, N=size(ens.positions,2), saveeach=50, save_filename="run"), crystal, ef)

pulse(t, A, w, t0, sig)

Gaussian wavepacket pulse

Remove acoustic sum rule from eigenvalue and eigenvectors

semiimpliciteulerstep!(wigner:: WignerDistribution{T}, dt :: T, avgfor :: Vector{T}, d2V_dr2 :: Matrix{T}) where {T <: AbstractFloat}

Semi-implicit Euler method. Implemented for testing and didactic use. Not efficient, use semi-implicit Verlet instead.

function semi_implicit_verlet_step!(rho:: WignerDistribution{T}, dt :: T, avg_for :: Vector{T}, d2V_dr2 :: Matrix{T}, part ) where {T <: AbstractFloat}

Evolves the quantum means and correlators of a WignerDistribution by one semi-implicit Verlet integration step. This algorithms has error O(dt^2) on the integration of the correlators.

This is a two-stage scheme:

• part == 1 updates the positions (R_av) and momenta (P_av) using the current force avg_for.
• part == 2 completes the step with the updated force, and evolves the second-order correlators (RR_corr, RP_corr, PP_corr).

Arguments

• ρ::WignerDistribution{T}: system state (means and correlators), mutated in place.
• dt::T: time step.
• avg_for::Vector{T}: average force vector ⟨F⟩.
• d2V_dr2::Matrix{T}: Hessian (second derivatives of the potential).
• part::Int: integration stage:

Method

Given the following variables name ⟨R⟩ = Rav, ⟨P⟩ = Pav, ⟨f⟩ = avgfor, ⟨RR⟩ = RRcorr, ⟨PP⟩ = PPcorr, ⟨RP⟩ = RPcorr, ⟨d²V⟩ = d2V_dr2

• In part == 1:
  • ⟨R⟩ ← ⟨R⟩ + dt·⟨P⟩ + ½·dt²·⟨f⟩
  • ⟨P⟩ ← ⟨P⟩ + ½·dt·⟨f⟩
• In part == 2:
  • ⟨P⟩ ← ⟨P⟩ + ½·dt·⟨f⟩
  • ⟨RP⟩ ← ⟨RP⟩ + dt·⟨PP⟩ − dt·⟨RR⟩·⟨d²V⟩
  • ⟨PP⟩ ← ⟨PP⟩ − dt·(⟨d²V⟩·⟨RP⟩ + (⟨d²V⟩·⟨RP⟩)ᵀ)
  • ⟨RR⟩ ← ⟨RR⟩ + dt·(⟨RP⟩ + ⟨RP⟩ᵀ)

sin_field(t, A, w)

Sinusoidal external field.

Update the weights of the ensemble according with the new wigner distribution