Code documentation

VCHam Object

class pyvcham.vcham_system.VCSystem(elements=[], reference_geometry=None, vc_type='linear', units='eV', number_normal_modes=None, number_states=None, displacement_vector=[], database_abinitio=[], coupling_with_gs=True, symmetry_point_group=None, totally_sym_irrep=None, symmetry_modes=[], symmetry_states=[], vib_freq=None, energy_shift=None, vcham=[], diab_funct=[], symmetry_matrix=[], jt_effects=[], dipole_matrix=[], shared_states=None)

Bases: object

Initialize the VCSystem object.

Parameters:

• elements (List[str])

• reference_geometry (ndarray | None)

• vc_type (str)

• units (str)

• number_normal_modes (int | None)

• number_states (int | None)

• displacement_vector (List[ndarray | Any])

• database_abinitio (List[ndarray | Any])

• coupling_with_gs (bool)

• symmetry_point_group (str | None)

• totally_sym_irrep (str | None)

• symmetry_modes (List[str])

• symmetry_states (List[str])

• vib_freq (ndarray | None)

• energy_shift (ndarray | None)

• vcham (List[ndarray])

• diab_funct (List[str])

• symmetry_matrix (List[ndarray])

• jt_effects (List[Dict[str, Any]])

• dipole_matrix (List[ndarray])

• shared_states (Dict[int, List[Tuple[int, int]]])

elements

List of element symbols in the system.

Type:

list of str, optional

reference_geometry

Reference geometry of the system, given as a NumPy array of shape (n_atoms, 3).

Type:

np.ndarray, optional

vc_type

Type of the vibrational coupling system, e.g., “linear”, “quadratic”, etc.

Type:

str, default “linear”

units

Units for the energy and frequency values, e.g., “eV”, “Hartree”, etc.

Type:

str, default “eV”

number_normal_modes

Number of normal modes in the system.

Type:

int, optional

number_states

Number of electronic states considered in the system.

Type:

int, optional

displacement_vector

List of displacement vectors for each normal mode, where each vector is a NumPy array.

Type:

list of np.ndarray, optional

database_abinitio

List of ab initio data arrays, where each array corresponds to a state and contains energy values for different displacements. It is a list of NumPy arrays, each of shape (n_states, n_displacements). Overall shape is (number_normal_modes, number_states, n_displacements).

Type:

list of np.ndarray, optional

coupling_with_gs

Whether the system is coupled with the ground state (GS). If not provided, it will be initialized as True. It is useful for highly excited states.

Type:

bool, optional

symmetry_point_group

The point group symmetry of the system, e.g., “C2v”, “Cs”, etc.

Type:

str, optional

totally_sym_irrep

The totally symmetric irreducible representation of the system. It is automatically computed from the symmetry point group.

Type:

str, optional

symmetry_modes

List of irreps corresponding to the normal modes.

Type:

list of str, optional

symmetry_states

List of irreps corresponding to the electronic states.

Type:

list of str, optional

vib_freq

Vibrational frequencies for the normal modes, given as a NumPy array.

Type:

np.ndarray, optional

energy_shift

Energy shifts for the states, given as a NumPy array. It is computed from the ab initio database.

Type:

np.ndarray, optional

vcham

List of vibrational coupling Hamiltonian matrices, where each matrix is a NumPy array. It is used to store the Hamiltonian for each normal mode.

Type:

list of np.ndarray, optional

diab_funct

List of diabatic function names for each normal mode.

Type:

list of str, optional

symmetry_matrix

List of symmetry matrices, where each matrix is a NumPy array. It is used to store the symmetry operations for the system. It is computed from the symmetry variables.

Type:

list of np.ndarray, optional

jt_effects

List of Jahn-Teller (JT) effects, where each effect is a dictionary containing: - ‘mode’: int, the index of the normal mode. - ‘state_pairs’: list, pairs of states coupled by the JT effect. - ‘types’: list, JT effect types for each state pair. Optionally, if a JT effect is inactive (‘active’: False), it must include: - ‘source’: int, the mode index to copy parameters from.

Type:

list of dict, optional

dipole_matrix

List of dipole matrices, where each matrix is a NumPy array. It is used to store the dipole moments for the system. The dimensions of the dipole matrix should be (n_states, n_states, 3).

Type:

list of np.ndarray, optional

add_dipole_matrix(dipole_matrix)

Add or update the dipole matrix.

Parameters:

dipole_matrix (np.ndarray) – The dipole matrix to be added or updated. It should be a 3D numpy array where each slice corresponds to a state and each vector corresponds to a mode. The dipole matrix should have dimensions (n_states, n_states, 3). The last dimension represents the x, y, z components of the dipole moment.

Return type:

None

add_geometry(geometry)

Add or update the reference geometry. This method processes the provided geometry, ensures it is centered at the center of mass (COM), and updates the system’s reference geometry.

Parameters:

geometry (list of tuple or str) – The geometry data to be added or updated. It can be provided in two formats: 1. A list of tuples, where each tuple contains an element symbol and its coordinates: - e.g., [(“H”, [0.0, 0.0, 0.0]), (“O”, [0.0, 0.0, 1.0]), …]. 2. A string path to an XYZ file containing the geometry.

Return type:

None

Linear vibronic coupling Hamiltonian module

class pyvcham.lvc.LVCHam(normal_mode, VCSystem, funct_guess=None, nepochs=3000, optimization=True)

Bases: object

Builder and optimizer for a Linear Vibronic Coupling (LVC) Hamiltonian.

This class extracts parameters from a VCSystem object, prepares both diagonal (on-diagonal) and off-diagonal terms, and performs gradient-based optimization.

Parameters:

• normal_mode (int)

• VCSystem (Any)

• funct_guess (Any | None)

• nepochs (int)

• optimization (bool)

initialize_loss_function(fn='huber', **kwargs)

Set up the loss function for optimization.

Parameters:

fn (str)

Return type:

None

initialize_params(lambda_guess=0.1, jt_guess=0.01, kappa_guess=0.1)

Initialize TensorFlow variables for optimization based on symmetry and Jahn-Teller (JT) effects.

Parameters:

• lambda_guess (float) – Initial guess for off-diagonal lambda parameters. Defaults to 0.1.

• jt_guess (float) – Initial guess for JT off-diagonal parameters. Defaults to 0.01.

• kappa_guess (float) – Initial guess for on-diagonal kappa parameters. Defaults to 0.1.

Return type:

None

optimize()

Perform gradient-based optimization with a smart early-stopping mechanism.

Return type:

None

Coupling Functions

This module provides coupling functions for the vibronic coupling Hamiltonian.

Currently, the module supports linear coupling, with a design that allows extension to additional coupling types.

New coupling functions can be added by defining a new function and updating the COUPLING_TYPES, coupling_funct, n_var, and initial_guesses dictionaries accordingly.

Functions

• linear_coupling(): Computes the linear coupling function.

Variables

• COUPLING_TYPES: List of supported coupling types.

• coupling_funct: Dictionary mapping coupling types to their implementation functions.

• n_var: Dictionary specifying the number of parameters per coupling function.

• initial_guesses: Dictionary of default initial parameter guesses for each coupling type.

pyvcham.couplingfunct.linear_coupling(q, k1)

Compute the linear coupling function.

The linear coupling is defined as:

\[f(q) = k_1 \cdot q\]

where \(q\) is the displacement along a normal mode, and \(k_1\) is the linear coupling constant.

Parameters:

• q (tf.Tensor) – Tensor of displacements, typically along a normal mode, with dtype tf.float32.

• k1 (tf.Tensor or float) – Coupling constant as a tensor or scalar. If a float is provided, it is converted to a tf.Tensor with dtype tf.float32.

Returns:

Tensor of computed linear coupling values with dtype tf.float32.

Return type:

tf.Tensor

Diabatic Functions

This module contains the diabatic functions.

New diabatic functions can be added by defining a new function and updating the potential_functions, n_var, initial_guesses, and kappa_compatible dictionaries accordingly.

Available Potential Types

• "ho": Harmonic oscillator

• "morse": General Morse potential

• "antimorse": Anti-Morse potential (uses general_morse with negative alpha)

• "quartic": General quartic potential

Functions

• harmonic_oscillator()

• general_quartic_potential()

• general_morse()

Variables

• potential_functions

• n_var

• initial_guesses

• kappa_compatible

pyvcham.diabfunct.general_morse(q, params, gs=False)

Compute the general Morse potential.

The Morse potential models anharmonic bond stretching in diatomic molecules, widely used in quantum chemistry for potential energy surfaces.

The potential is:

\[V(q) = D_e \left(1 - e^{-\alpha (q - q_0)}\right)^2 - D_e \left(1 - e^{-\alpha q_0}\right)^2\]

where: - \(D_e\) is the dissociation energy. - \(\alpha\) is the range parameter (well width). - \(q_0\) is the equilibrium bond distance. - The offset ensures \(V(0) = 0\).

Parameters:

• q (tf.Tensor) – Tensor of displacements (bond lengths), broadcast-compatible shape, dtype tf.float32.

• params (list) – List with three elements: - De: Dissociation energy, scalar or tensor. - alpha: Range parameter, scalar or tensor. - q0: Equilibrium bond distance, scalar or tensor. All broadcast-compatible with q.

• gs (bool)

Returns:

Tensor of potential energy values, matching q’s shape.

Return type:

tf.Tensor

pyvcham.diabfunct.general_quartic_potential(q, omega, params)

Compute the general quartic potential for a diabatic state.

This function computes a quartic potential with quadratic and quartic terms, ideal for modeling anharmonic vibrational effects, such as bond stretching, in quantum chemistry.

The potential is defined as:

\[V(q) = \frac{1}{2} \omega^2 q^2 + \frac{1}{2} k_2 q^2 + \frac{1}{24} k_3 q^4\]

where: - \(q\) is the displacement along a normal mode. - \(\omega\) is the mode frequency. - \(k_2\) is an additional quadratic coefficient. - \(k_3\) is the quartic coefficient introducing anharmonicity.

Parameters:

• q (tf.Tensor) – Tensor of displacements, broadcast-compatible shape, dtype tf.float32.

• omega (tf.Tensor) – Tensor of mode frequency, scalar or broadcast-compatible with q, dtype tf.float32.

• params (list) – List with two elements: - k2: Quadratic coefficient. - k3: Quartic coefficient. Both scalar or broadcast-compatible with q.

Returns:

Tensor of potential energy values, matching q’s shape.

Return type:

tf.Tensor

pyvcham.diabfunct.harmonic_oscillator(q, omega, params)

Compute the harmonic oscillator potential.

The potential is given by:

\[V(q) = \frac{1}{2} \omega q^2 + \frac{1}{2} k_1 q^2\]

where: - \(q\) is the displacement along a normal mode. - \(\omega\) is the frequency of the mode. - \(k_1\) is an additional quadratic coupling constant.

Parameters:

• q (tf.Tensor) – Tensor of displacements along a normal mode, with shape suitable for broadcasting. Must be dtype tf.float32.

• omega (tf.Tensor) – Tensor representing the mode frequency, scalar or broadcast-compatible with q. Must be dtype tf.float32.

• params (list) – List with one element, k1, the quadratic constant, scalar or broadcast-compatible with q.

Returns:

Tensor of diabatic potential energy values, matching the shape of q.

Return type:

tf.Tensor

Constants

This module contains physical constants and conversion factors.

pyvcham.constants.H

Planck’s constant (J·s)

pyvcham.constants.C

Speed of light in vacuum (m/s)

pyvcham.constants.AU_TO_J

Hartree to Joule conversion factor (J)

pyvcham.constants.CM1_TO_AU

Conversion factor from wavenumbers (cm^-1) to atomic units

pyvcham.constants.AU_TO_EV

Conversion factor from atomic units to electron volts (eV)

Utilities

Utility functions for PyVCHAM package. This module provides various utility functions for handling molecular geometries, rotating geometries, encoding data to JSON, and processing VCHAM parameters.

class pyvcham.utils.DipoleInteraction(molecules, pairs_interacting=None)

Bases: object

A class representing dipole interactions between pairs of molecules.

Parameters:

• molecules (List[Molecule])

• pairs_interacting (List[Tuple[int, int]] | None)

molecules

List of Molecule objects involved in the dipole interactions.

Type:

List[‘Molecule’]

pairs_interacting

List of tuples specifying which pairs of molecules (by their indices in the molecules list) interact via dipole interactions. If None, all possible unique pairs are used.

Type:

Optional[List[Tuple[int, int]]]

class pyvcham.utils.Molecule(molecule_idx, interacting_states=None, CM=None, rot_angles=None)

Bases: object

Represents a molecule in the VCHAM system.

Parameters:

• molecule_idx (int)

• interacting_states (List[int] | None)

• CM (ndarray)

• rot_angles (Tuple[float, float, float])

molecule_idx

Index of the molecule in the system.

Type:

int

interacting_states

List of interacting states for the molecules.

Type:

Optional[List[int]]

CM

Center of mass of the molecule in Bohr units.

Type:

np.ndarray

rot_angles

Euler angles for the molecule’s rotation (degrees).

Type:

Tuple[float, float, float]

pyvcham.utils.VCSystem_to_json(VCSystem, general_data={}, output_name='vcham_data.json', rewrite=False)

Convert a VCSystem object to JSON and save to a file.

Parameters:

• VCSystem (Any)

• general_data (Dict[str, Any])

• output_name (str)

• rewrite (bool)

Return type:

None

pyvcham.utils.eulerMatrix(xhi, theta, phi)

Generate a rotation matrix for Euler angles in ZYZ convention. This function constructs a 3x3 rotation matrix based on the provided Euler angles.

Parameters:

• xhi (float) – Euler angle in radians for the first rotation about the Z-axis.

• theta (float) – Euler angle in radians for the second rotation about the Y-axis.

• phi (float) – Euler angle in radians for the third rotation about the Z-axis.

Returns:

A 3x3 rotation matrix corresponding to the Euler angles in ZYZ convention.

Return type:

np.ndarray

pyvcham.utils.json_to_mctdh(infile, outfile)

Reads a JSON file and writes a single MCTDH operator file.

Parameters:

• infile (str) – JSON file path.

• outfile (str) – Output MCTDH operator file path.

Return type:

None

pyvcham.utils.mctdh_to_json(infile, outfile, template_json=None)

Reads an MCTDH operator (.op) file and converts it into a JSON file (.json).

If a template_json file is provided (or if a .json file with the same base name as infile exists in the same directory), non-numerical metadata (such as symmetries, point groups, and elements) are preserved while all parameter values (frequencies, energy shifts, diabatic curve parameters, kappa, and lambda) are updated from the .op file. Otherwise, a valid standalone JSON dictionary is reconstructed.

Parameters:

• infile (str) – Input MCTDH operator (.op) file path.

• outfile (str) – Output JSON (.json) file path.

• template_json (Optional[str], optional) – Optional path to a reference JSON file to update parameters in-place.

Return type:

None

pyvcham.utils.merge_jsons(infiles, outfile, molecules=None, interactions=None, rewrite=False)

Reads multiple JSON files and writes a single JSON file file by concatenating them with proper indexing offsets and possible coupling terms coming from the functions specified in the functions list.

Parameters:

• infiles (Union[str, List[str]]) – One or more JSON file paths.

• outfile (str) – Output MCTDH operator file path.

• molecules (List[Molecule])

• interactions (Any | None)

• rewrite (bool)

Return type:

None

pyvcham.utils.op_to_json(infile, outfile, template_json=None)

Reads an MCTDH operator (.op) file and converts it into a JSON file (.json).

If a template_json file is provided (or if a .json file with the same base name as infile exists in the same directory), non-numerical metadata (such as symmetries, point groups, and elements) are preserved while all parameter values (frequencies, energy shifts, diabatic curve parameters, kappa, and lambda) are updated from the .op file. Otherwise, a valid standalone JSON dictionary is reconstructed.

Parameters:

• infile (str) – Input MCTDH operator (.op) file path.

• outfile (str) – Output JSON (.json) file path.

• template_json (Optional[str], optional) – Optional path to a reference JSON file to update parameters in-place.

Return type:

None

pyvcham.utils.plot_surfaces_from_json(json_file, mode, surface_type='diabatic', points=100, range_q=(-5.0, 5.0), molecule_idx=0, ax=None)

Plots the diabatic or adiabatic potential energy surfaces along a specified normal mode using the parameters stored in a PyVCHAM JSON file.

Parameters:

• json_file (str) – Path to the JSON file generated by PyVCHAM.

• mode (int) – Index of the normal mode to plot.

• surface_type (str, optional) – Type of surface to plot: “diabatic” or “adiabatic”. Default is “diabatic”.

• points (int, optional) – Number of points to evaluate the surfaces along the normal mode. Default is 100.

• range_q (Tuple[float, float], optional) – The (min, max) range of the dimensionless normal mode displacement to plot. Default is (-5.0, 5.0).

• molecule_idx (int, optional) – Index of the molecule if the JSON contains multiple molecules. Default is 0.

• ax (matplotlib.axes.Axes, optional) – An optional Axes object to plot on. If None, a new figure and axes are created.

Returns:

The axes object containing the plot.

Return type:

matplotlib.axes.Axes

pyvcham.utils.rotate_dipole_matrix(dipole_matrix, xhi, theta, phi)

Rotate a dipole moment matrix using Euler angles in ZYZ convention.

pyvcham.utils.rotate_geometry(geometry, xhi, theta, phi)

Rotate a 3D geometry using Euler angles in ZYZ convention.