Trajectory Methods¶

The classes documented here are responsible for running molecular dynamics simulations with different algorithms.

Bases: object

Class to propagate a single FSSH trajectory.

This class implements the Fewest Switches Surface Hopping (FSSH) algorithm for nonadiabatic molecular dynamics simulations. It handles the propagation of both nuclear and electronic degrees of freedom, including surface hopping events between electronic states.

Initialize the SurfaceHoppingMD class.

Parameters:

model (Any) – Model object defining problem
x0 (ArrayLike) – Initial position
v0 (ArrayLike) – Initial velocity
rho0 (ArrayLike, int, or str) – Initial density matrix or state index. If an integer, populates a single state. If a matrix, populates a density matrix (requires state0 for active state).
tracer (Any, optional) – spawn from TraceManager to collect results, by default “default”
queue (Any, optional) – Trajectory queue, by default None
strict_option_check (bool, optional) – Whether to strictly check options, by default True
propagator (str or dict, optional) – The propagator to use for nuclear motion. Can be a string (e.g., ‘vv’, ‘fssh’) or a dictionary with more options. Default is ‘vv’.
last_velocity (array-like, optional) – The velocity from the previous step, used for restarts. Default is zeros.
bounds (tuple or list, optional) – Tuple or list of (lower, upper) bounds for the simulation box. Used to determine if the trajectory is inside a region. Default is None.
duration (dict, optional) – Dictionary controlling simulation duration (overrides max_steps, max_time, etc.). Default is auto-generated.
dt (float, optional) – Time step for nuclear propagation (in atomic units). Default is fs_to_au.
t0 (float, optional) – Initial time. Default is 0.0.
previous_steps (int, optional) – Number of previous steps (for restarts). Default is 0.
trace_every (int, optional) – Interval (in steps) at which to record trajectory data. Default is 1.
max_steps (int, optional) – Maximum number of simulation steps. Default is 1000000.
max_time (float, optional) – Maximum simulation time. Default is 1e25.
seed_sequence (int or numpy.random.SeedSequence, optional) – Seed or SeedSequence for random number generation. Default is None.
outcome_type (str, optional) – Type of outcome to record (e.g., ‘state’). Default is ‘state’.
electronics (object, optional) – Initial electronic state object. Default is None.
electronic_integration (str, optional) – Method for integrating electronic equations (‘exp’ or ‘linear-rk4’). Default is ‘exp’.
max_electronic_dt (float, optional) – Maximum time step for electronic integration. Default is 0.1.
starting_electronic_intervals (int, optional) – Initial number of intervals for electronic integration. Default is 4.
weight (float, optional) – Statistical weight of the trajectory. Default is 1.0.
restarting (bool, optional) – Whether this is a restarted trajectory. Default is False.
hopping_probability (str, optional) – Method for computing hopping probability (‘tully’ or ‘poisson’). Default is ‘tully’.
zeta_list (list, optional) – List of pre-determined random numbers for hopping decisions. Default is [].
state0 (int, optional) – Initial electronic state (used if rho0 is a matrix). Required if rho0 is not scalar.
hopping_method (str, optional) – Hopping method: ‘cumulative’, ‘cumulative_integrated’, or ‘instantaneous’. Default is ‘cumulative’.

recognized_options = ['propagator', 'last_velocity', 'bounds', 'dt', 't0', 'previous_steps', 'duration', 'max_steps', 'max_time', 'seed_sequence', 'outcome_type', 'trace_every', 'electronics', 'electronic_integration', 'max_electronic_dt', 'starting_electronic_intervals', 'weight', 'restarting', 'hopping_probability', 'zeta_list', 'state0', 'hopping_method']¶

Initialize the SurfaceHoppingMD class.

Parameters:

model (Any) – Model object defining problem
x0 (ArrayLike) – Initial position
v0 (ArrayLike) – Initial velocity
rho0 (ArrayLike, int, or str) – Initial density matrix or state index. If an integer, populates a single state. If a matrix, populates a density matrix (requires state0 for active state).
tracer (Any, optional) – spawn from TraceManager to collect results, by default “default”
queue (Any, optional) – Trajectory queue, by default None
strict_option_check (bool, optional) – Whether to strictly check options, by default True
propagator (str or dict, optional) – The propagator to use for nuclear motion. Can be a string (e.g., ‘vv’, ‘fssh’) or a dictionary with more options. Default is ‘vv’.
last_velocity (array-like, optional) – The velocity from the previous step, used for restarts. Default is zeros.
bounds (tuple or list, optional) – Tuple or list of (lower, upper) bounds for the simulation box. Used to determine if the trajectory is inside a region. Default is None.
duration (dict, optional) – Dictionary controlling simulation duration (overrides max_steps, max_time, etc.). Default is auto-generated.
dt (float, optional) – Time step for nuclear propagation (in atomic units). Default is fs_to_au.
t0 (float, optional) – Initial time. Default is 0.0.
previous_steps (int, optional) – Number of previous steps (for restarts). Default is 0.
trace_every (int, optional) – Interval (in steps) at which to record trajectory data. Default is 1.
max_steps (int, optional) – Maximum number of simulation steps. Default is 1000000.
max_time (float, optional) – Maximum simulation time. Default is 1e25.
seed_sequence (int or numpy.random.SeedSequence, optional) – Seed or SeedSequence for random number generation. Default is None.
outcome_type (str, optional) – Type of outcome to record (e.g., ‘state’). Default is ‘state’.
electronics (object, optional) – Initial electronic state object. Default is None.
electronic_integration (str, optional) – Method for integrating electronic equations (‘exp’ or ‘linear-rk4’). Default is ‘exp’.
max_electronic_dt (float, optional) – Maximum time step for electronic integration. Default is 0.1.
starting_electronic_intervals (int, optional) – Initial number of intervals for electronic integration. Default is 4.
weight (float, optional) – Statistical weight of the trajectory. Default is 1.0.
restarting (bool, optional) – Whether this is a restarted trajectory. Default is False.
hopping_probability (str, optional) – Method for computing hopping probability (‘tully’ or ‘poisson’). Default is ‘tully’.
zeta_list (list, optional) – List of pre-determined random numbers for hopping decisions. Default is [].
state0 (int, optional) – Initial electronic state (used if rho0 is a matrix). Required if rho0 is not scalar.
hopping_method (str, optional) – Hopping method: ‘cumulative’, ‘cumulative_integrated’, or ‘instantaneous’. Default is ‘cumulative’.

classmethod restart(model, log, **options) → SurfaceHoppingMD¶

Restart a simulation from a previous trajectory log.

Parameters:

model (Any) – Model object defining problem
log (Any) – Previous trajectory log
**options (Any) – Additional options for the simulation

Returns:

New instance initialized from the log data

Return type:

SurfaceHoppingMD

update_weight(weight: float) → None¶

Update weight held by trajectory and by trace.

Parameters:: weight (float) – New weight value

clone() → SurfaceHoppingMD¶

Clone existing trajectory for spawning.

Returns:: Copy of current object
Return type:: SurfaceHoppingMD

random() → float64¶

Get random number for hopping decisions.

Returns:: Uniform random number between 0 and 1
Return type:: np.float64

currently_interacting() → bool¶

Determine whether trajectory is currently inside an interaction region.

Returns:: True if trajectory is inside interaction region, False otherwise
Return type:: bool

duration_initialize(options: Dict[str, Any]) → None¶

Initialize variables related to continue_simulating.

Parameters:: options (Dict[str, Any]) – Dictionary with options for simulation duration

continue_simulating() → bool¶

Decide whether trajectory should keep running.

Returns:: True if trajectory should keep running, False if it should finish
Return type:: bool

trace(force: bool = False) → None¶

Add results from current time point to tracing function.

Only adds snapshot if nsteps%trace_every == 0, unless force=True.

Parameters:: force (bool, optional) – Force snapshot regardless of trace_every, by default False

snapshot() → Dict[str, Any]¶

Collect data from run for logging.

Returns:: Dictionary with all data from current time step
Return type:: Dict[str, Any]

kinetic_energy() → float64¶

Calculate kinetic energy.

Returns:: Kinetic energy
Return type:: np.float64

potential_energy(electronics: ElectronicModel_ = None) → floating¶

Calculate potential energy.

Parameters:: electronics (ElectronicModel, optional) – electronic states from current step, by default None
Returns:: Potential energy
Return type:: np.floating

total_energy(electronics: ElectronicModel_ = None) → floating¶

Calculate total energy (kinetic + potential).

Parameters:: electronics (ElectronicModel, optional) – Electronic states from current step, by default None
Returns:: Total energy
Return type:: np.floating

Calculate nonadiabatic coupling matrix.

Parameters:

electronics (ElectronicModel, optional) – electronic states from current step, by default None
velocity (ArrayLike, optional) – Velocity used to compute NAC, by default None

Returns:

NAC matrix

Return type:

ArrayLike

Calculate kinetic energy along given momentum mode.

Parameters:: direction (ArrayLike) – Numpy array defining direction
Returns:: Kinetic energy along specified direction
Return type:: np.float64

draw_new_zeta() → float¶

Get a new zeta value for hopping.

First checks the input list of zeta values in self.zeta_list. If no value is found in zeta_list, then a random number is pulled.

Returns:: Zeta value for hopping
Return type:: float

Determine if a hop with given rescale direction and energy change is allowed.

Parameters:

direction (ArrayLike) – Momentum unit vector
dE (float) – Change in energy such that Enew = Eold + dE

Returns:

Whether the hop is allowed

Return type:

bool

direction_of_rescale(source: int, target: int, electronics: ElectronicModel_ = None) → ndarray¶

Return direction in which to rescale momentum.

Parameters:

source (int) – Active state before hop
target (int) – Active state after hop
electronics (ElectronicModel, optional) – Electronic model information (used to pull derivative coupling), by default None

Returns:

Unit vector pointing in direction of rescale

Return type:

np.ndarray

Update velocity by rescaling the momentum in the specified direction and amount.

Parameters:

direction (ArrayLike) – The direction of the momentum to rescale
reduction (np.floating) – How much kinetic energy should be damped

Compute the Hamiltonian used to propagate the electronic wavefunction.

Parameters:

last_electronics (ElectronicModel) – Electronic states at previous time step
this_electronics (ElectronicModel) – Electronic states at current time step
velo (ArrayLike, optional) – Velocity at midpoint between current and previous time steps, by default None

Returns:

Nonadiabatic coupling Hamiltonian at midpoint between current and previous time steps

Return type:

np.ndarray

propagate_electronics(last_electronics: ElectronicModel_, this_electronics: ElectronicModel_, dt: floating) → None¶

Propagate density matrix from t to t+dt.

The propagation assumes the electronic energies and couplings are static throughout. This will only be true for fairly small time steps.

Parameters:

last_electronics (ElectronicModel) – Electronic states at t
this_electronics (ElectronicModel) – Electronic states at t+dt
dt (np.floating) – Time step

surface_hopping(last_electronics: ElectronicModel_, this_electronics: ElectronicModel_)¶

Compute probability of hopping, generate random number, and perform hops.

Parameters:

last_electronics (ElectronicModel) – Electronic states at previous time step
this_electronics (ElectronicModel) – Electronic states at current time step

hopper(gkndt: ndarray) → List[Dict[str, float]]¶

Determine whether and where to hop.

Parameters:: gkndt (np.ndarray) – Array of individual hopping probabilities
Returns:: List of dictionaries with (target_state, weight) pairs
Return type:: List[Dict[str, float]]

hop_update(hop_from, hop_to)¶

Handle any extra operations that need to occur after a hop.

Parameters:

hop_from (int) – State before hop
hop_to (int) – State after hop

hop_to_it(hop_targets: List[Dict[str, float | int]], electronics: ElectronicModel_ = None) → None¶

Hop from the current active state to the given state, including rescaling the momentum.

Parameters:

hop_targets (List[Dict[str, Union[float, int]]]) – List of (target, weight) pairs
electronics (ElectronicModel, optional) – Electronic states for current step, by default None

simulate() → Trace¶

Run the surface hopping molecular dynamics simulation.

Returns:: Trace of trajectory
Return type:: Trace

class Ehrenfest(*args: Any, **kwargs: Any)¶

Bases: SurfaceHoppingMD

Ehrenfest dynamics.

This class implements Ehrenfest dynamics, where the electronic degrees of freedom are treated quantum mechanically and the nuclear degrees of freedom are treated classically. The force on the nuclei is computed as the expectation value of the force operator over the electronic density matrix.

Initialize Ehrenfest dynamics.

Parameters:

*args (Any) – Positional arguments passed to SurfaceHoppingMD
**kwargs (Any) – Keyword arguments passed to SurfaceHoppingMD

__init__(*args: Any, **kwargs: Any)¶

Initialize Ehrenfest dynamics.

Parameters:

*args (Any) – Positional arguments passed to SurfaceHoppingMD
**kwargs (Any) – Keyword arguments passed to SurfaceHoppingMD

potential_energy(electronics: ElectronicModel_ = None) → floating¶

Calculate Ehrenfest potential energy.

The potential energy is computed as the trace of the product of the density matrix and the Hamiltonian.

Parameters:: electronics (ElectronicModel, optional) – ElectronicStates from current step, by default None
Returns:: Potential energy
Return type:: np.floating

surface_hopping(last_electronics: ElectronicModel_, this_electronics: ElectronicModel_)¶

Handle surface hopping.

In Ehrenfest dynamics, surface hopping is not performed as the electronic states are treated quantum mechanically.

Parameters:

last_electronics (ElectronicModel) – Electronic states from previous step
this_electronics (ElectronicModel) – Electronic states from current step

class AugmentedFSSH(*args: Any, **options: Any)¶

Bases: SurfaceHoppingMD

Augmented-FSSH (A-FSSH) dynamics, by Subotnik and coworkers.

Initial implementation based on original paper:: Subotnik, Shenvi JCP 134, 024105 (2011); doi: 10.1063/1.3506779

Initialize the SurfaceHoppingMD class.

Parameters:

model (Any) – Model object defining problem
x0 (ArrayLike) – Initial position
v0 (ArrayLike) – Initial velocity
rho0 (ArrayLike, int, or str) – Initial density matrix or state index. If an integer, populates a single state. If a matrix, populates a density matrix (requires state0 for active state).
tracer (Any, optional) – spawn from TraceManager to collect results, by default “default”
queue (Any, optional) – Trajectory queue, by default None
strict_option_check (bool, optional) – Whether to strictly check options, by default True
propagator (str or dict, optional) – The propagator to use for nuclear motion. Can be a string (e.g., ‘vv’, ‘fssh’) or a dictionary with more options. Default is ‘vv’.
last_velocity (array-like, optional) – The velocity from the previous step, used for restarts. Default is zeros.
bounds (tuple or list, optional) – Tuple or list of (lower, upper) bounds for the simulation box. Used to determine if the trajectory is inside a region. Default is None.
duration (dict, optional) – Dictionary controlling simulation duration (overrides max_steps, max_time, etc.). Default is auto-generated.
dt (float, optional) – Time step for nuclear propagation (in atomic units). Default is fs_to_au.
t0 (float, optional) – Initial time. Default is 0.0.
previous_steps (int, optional) – Number of previous steps (for restarts). Default is 0.
trace_every (int, optional) – Interval (in steps) at which to record trajectory data. Default is 1.
max_steps (int, optional) – Maximum number of simulation steps. Default is 1000000.
max_time (float, optional) – Maximum simulation time. Default is 1e25.
seed_sequence (int or numpy.random.SeedSequence, optional) – Seed or SeedSequence for random number generation. Default is None.
outcome_type (str, optional) – Type of outcome to record (e.g., ‘state’). Default is ‘state’.
electronics (object, optional) – Initial electronic state object. Default is None.
electronic_integration (str, optional) – Method for integrating electronic equations (‘exp’ or ‘linear-rk4’). Default is ‘exp’.
max_electronic_dt (float, optional) – Maximum time step for electronic integration. Default is 0.1.
starting_electronic_intervals (int, optional) – Initial number of intervals for electronic integration. Default is 4.
weight (float, optional) – Statistical weight of the trajectory. Default is 1.0.
restarting (bool, optional) – Whether this is a restarted trajectory. Default is False.
hopping_probability (str, optional) – Method for computing hopping probability (‘tully’ or ‘poisson’). Default is ‘tully’.
zeta_list (list, optional) – List of pre-determined random numbers for hopping decisions. Default is [].
state0 (int, optional) – Initial electronic state (used if rho0 is a matrix). Required if rho0 is not scalar.
hopping_method (str, optional) – Hopping method: ‘cumulative’, ‘cumulative_integrated’, or ‘instantaneous’. Default is ‘cumulative’.

__init__(*args: Any, **options: Any)¶

Initialize the SurfaceHoppingMD class.

Parameters:

model (Any) – Model object defining problem
x0 (ArrayLike) – Initial position
v0 (ArrayLike) – Initial velocity
rho0 (ArrayLike, int, or str) – Initial density matrix or state index. If an integer, populates a single state. If a matrix, populates a density matrix (requires state0 for active state).
tracer (Any, optional) – spawn from TraceManager to collect results, by default “default”
queue (Any, optional) – Trajectory queue, by default None
strict_option_check (bool, optional) – Whether to strictly check options, by default True
propagator (str or dict, optional) – The propagator to use for nuclear motion. Can be a string (e.g., ‘vv’, ‘fssh’) or a dictionary with more options. Default is ‘vv’.
last_velocity (array-like, optional) – The velocity from the previous step, used for restarts. Default is zeros.
bounds (tuple or list, optional) – Tuple or list of (lower, upper) bounds for the simulation box. Used to determine if the trajectory is inside a region. Default is None.
duration (dict, optional) – Dictionary controlling simulation duration (overrides max_steps, max_time, etc.). Default is auto-generated.
dt (float, optional) – Time step for nuclear propagation (in atomic units). Default is fs_to_au.
t0 (float, optional) – Initial time. Default is 0.0.
previous_steps (int, optional) – Number of previous steps (for restarts). Default is 0.
trace_every (int, optional) – Interval (in steps) at which to record trajectory data. Default is 1.
max_steps (int, optional) – Maximum number of simulation steps. Default is 1000000.
max_time (float, optional) – Maximum simulation time. Default is 1e25.
seed_sequence (int or numpy.random.SeedSequence, optional) – Seed or SeedSequence for random number generation. Default is None.
outcome_type (str, optional) – Type of outcome to record (e.g., ‘state’). Default is ‘state’.
electronics (object, optional) – Initial electronic state object. Default is None.
electronic_integration (str, optional) – Method for integrating electronic equations (‘exp’ or ‘linear-rk4’). Default is ‘exp’.
max_electronic_dt (float, optional) – Maximum time step for electronic integration. Default is 0.1.
starting_electronic_intervals (int, optional) – Initial number of intervals for electronic integration. Default is 4.
weight (float, optional) – Statistical weight of the trajectory. Default is 1.0.
restarting (bool, optional) – Whether this is a restarted trajectory. Default is False.
hopping_probability (str, optional) – Method for computing hopping probability (‘tully’ or ‘poisson’). Default is ‘tully’.
zeta_list (list, optional) – List of pre-determined random numbers for hopping decisions. Default is [].
state0 (int, optional) – Initial electronic state (used if rho0 is a matrix). Required if rho0 is not scalar.
hopping_method (str, optional) – Hopping method: ‘cumulative’, ‘cumulative_integrated’, or ‘instantaneous’. Default is ‘cumulative’.

compute_delF(this_electronics)¶

Compute the difference in forces between states.

Parameters:: this_electronics (ElectronicModel) – Current electronic state information.
Returns:: Matrix of force differences between states.
Return type:: numpy.ndarray

advance_delR(last_electronics, this_electronics)¶

Propagate delR using Eq. (29) from Subotnik 2011 JCP.

Parameters:

last_electronics (ElectronicModel) – Electronic state information from previous step.
this_electronics (ElectronicModel) – Electronic state information from current step.

advance_delP(last_electronics, this_electronics)¶

Propagate delP using Eq. (31) from Subotnik JCP 2011.

Parameters:

last_electronics (ElectronicModel) – Electronic state information from previous step.
this_electronics (ElectronicModel) – Electronic state information from current step.

direction_of_rescale(source: int, target: int, electronics: ElectronicModel_ = None) → ndarray¶

Return direction in which to rescale momentum.

In Subotnik JCP 2011, they suggest to use the difference between the momenta on delP.

Parameters:

source (int) – Active state before hop.
target (int) – Active state after hop.
electronics (ElectronicModel, optional) – Electronic model information (ignored), by default None.

Returns:

Unit vector pointing in direction of rescale.

Return type:

numpy.ndarray

gamma_collapse(electronics: ElectronicModel_ = None) → ndarray¶

Compute probability of collapse to each electronic state.

Uses Eq. (55) in Subotnik JCP 2011. This formula has some major problems and is tweaked or abandoned in future papers.

Parameters:: electronics (ElectronicModel, optional) – Electronic model information (forces), by default None.
Returns:: Array of collapse probabilities for each state.
Return type:: numpy.ndarray

surface_hopping(last_electronics: ElectronicModel_, this_electronics: ElectronicModel_) → None¶

Specialized version of surface_hopping that handles collapsing.

Parameters:

last_electronics (ElectronicModel) – ElectronicStates from previous step.
this_electronics (ElectronicModel) – ElectronicStates from current step.

hop_update(hop_from, hop_to)¶

Shift delR and delP after hops.

Parameters:

hop_from (int) – State index before hop.
hop_to (int) – State index after hop.

class EvenSamplingTrajectory(*args, **options)¶

Bases: SurfaceHoppingMD

Trajectory surface hopping with even sampling of phase space.

This class extends SurfaceHoppingMD to implement even sampling of phase space during trajectory surface hopping simulations.

Parameters:

*args (tuple) – Arguments passed to SurfaceHoppingMD constructor
**options (dict) – Options passed to SurfaceHoppingMD constructor
spawn_stack (SpawnStack, optional) – Stack for managing trajectory spawning
quadrature (dict, optional) – Quadrature parameters for sampling
mcsamples (int, optional) – Number of Monte Carlo samples

Initialize the EvenSamplingTrajectory.

Parameters:

*args (tuple) – Arguments passed to SurfaceHoppingMD constructor
**options (dict) – Options passed to SurfaceHoppingMD constructor
spawn_stack (SpawnStack, optional) – Stack for managing trajectory spawning
quadrature (dict, optional) – Quadrature parameters for sampling
mcsamples (int, optional) – Number of Monte Carlo samples

recognized_options = ['propagator', 'last_velocity', 'bounds', 'dt', 't0', 'previous_steps', 'duration', 'max_steps', 'max_time', 'seed_sequence', 'outcome_type', 'trace_every', 'electronics', 'electronic_integration', 'max_electronic_dt', 'starting_electronic_intervals', 'weight', 'restarting', 'hopping_probability', 'zeta_list', 'state0', 'hopping_method', 'spawn_stack', 'quadrature', 'mcsamples']¶

__init__(*args, **options)¶

Initialize the EvenSamplingTrajectory.

Parameters:

*args (tuple) – Arguments passed to SurfaceHoppingMD constructor
**options (dict) – Options passed to SurfaceHoppingMD constructor
spawn_stack (SpawnStack, optional) – Stack for managing trajectory spawning
quadrature (dict, optional) – Quadrature parameters for sampling
mcsamples (int, optional) – Number of Monte Carlo samples

clone(spawn_stack: Any | None = None) → EvenSamplingTrajectory¶

Create a clone of the current trajectory.

Parameters:: spawn_stack (Any, optional) – New spawn stack to use for the clone, by default None
Returns:: Cloned trajectory instance.
Return type:: EvenSamplingTrajectory

Determine whether and where to hop based on probabilities.

Parameters:: probs (ArrayLike) – Array of individual hopping probabilities.
Returns:: List of dictionaries containing target state and hop weight information.
Return type:: List[Dict[str, Union[int, float]]]

hop_to_it(hop_to: List[Dict[str, Any]], electronics: ElectronicModel_ = None) → None¶

Handle hopping by spawning new trajectories.

This method spawns new trajectories instead of enacting hops directly. It must: - Copy current trajectory - Initiate hops on the copied trajectories - Make no changes to current trajectory - Set next threshold for spawning

Parameters:

hop_to (List[Dict[str, Any]]) – List of states and associated weights on which to hop.
electronics (ElectronicModel, optional) – Model class for electronic structure calculations, by default None

Trajectory Methods¶

Table of Contents

Previous topic

Next topic

This Page