\(\text{Li}_3 \text{ Model}\)

[1]:

import pyvcham
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

2026-02-03 21:33:43.516665: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: SSE4.1 SSE4.2, in other operations, rebuild TensorFlow with the appropriate compiler flags.
2026-02-03 21:33:47,391 [INFO] root: Logging configuration successfully loaded.

[2]:

"""
This is an example of a Linear Vibronic Coupling (LVC) calculation using the vcham package.
The molecular system studied is Li3, which is a triatomic molecule consisting of three lithium atoms in D3h symmetry.
The GS is doubly degenerate at this symmetry. We will observe E x e Jahn-Teller (JT) effect in this system.

We will build the database from the abinitio data. Both the database and the displacement
data will be saved in a list of numpy arrays. This is because it might be that the
range of the displacement data is not the same for all coordinates.

"""

[2]:

'\nThis is an example of a Linear Vibronic Coupling (LVC) calculation using the vcham package.\nThe molecular system studied is Li3, which is a triatomic molecule consisting of three lithium atoms in D3h symmetry.\nThe GS is doubly degenerate at this symmetry. We will observe E x e Jahn-Teller (JT) effect in this system.\n\nWe will build the database from the abinitio data. Both the database and the displacement\ndata will be saved in a list of numpy arrays. This is because it might be that the\nrange of the displacement data is not the same for all coordinates.\n\n'

[3]:

# Define the number of states and normal modes
number_states = 2
number_normal_modes = 3

all_q = []
all_data = []

# Define the normal modes to process
normal_modes_to_process = [1, 1, 3]
# This is done because the first two modes are degenerate (1, 1) and the third mode is different (3).

# Loop over selected normal modes
for normal_mode in normal_modes_to_process:
    mode_data = []

    # For each state (root 1 for ground state, root 2 for valence-excited in this case)
    for root in range(1, 3):
        # Construct the filename and load the data
        filename = f'li3_abinitio/gs{root}_v{normal_mode}.dat'
        data = np.genfromtxt(filename)
        # Assume data columns: first column = displacement (x), second column = energy (y)
        displacements = data[:, 0] / 10  # Scale displacement
        energies = data[:, 1]

        # Create a DataFrame for clarity (optional)
        df = pd.DataFrame({'disp': displacements, 'gs': energies})
        # Shift the energies and convert from atomic units to eV
        adjusted_energies = (df["gs"].to_numpy() + 22.35342994) * pyvcham.constants.AU_TO_EV
        mode_data.append(adjusted_energies)

    # Convert the data for this mode to a numpy array and store the displacement vector
    all_data.append(np.array(mode_data))
    all_q.append(df["disp"].to_numpy())

# The dimensions of the all_q and all_data lists are:
# all_q: (number_normal_modes, number of displacements)
# all_data: (number_normal_modes, number_states, number of displacements)

[4]:

for mode_idx in range(3):
    plt.figure(figsize=(8, 5))

    # Plot data for each state
    for state_idx in range(number_states):
        plt.scatter(
            all_q[mode_idx],
            all_data[mode_idx][state_idx],
            s=0.5,
            label=f"State {state_idx}"
        )

    plt.xlabel(f"$Q_{{{mode_idx+1}}}$")
    plt.ylabel("Energy [a.u]")
    plt.title(f"Normal Mode {mode_idx} Ab initio Data")
    plt.legend(fontsize=8)
    plt.show()

[5]:

# Define vibrational frequencies in cm^-1 and convert to eV
vib_freq = np.array([389.55, 389.55, 308.52]) * pyvcham.constants.CM1_TO_AU * pyvcham.constants.AU_TO_EV

# Define Jahn-Teller effects
jt_effects = [
    {
        'mode': 0,                 # Normal mode index where the JT effect is active
        'state_pairs': [(0, 1)],   # State pair(s) coupled by the JT effect
        'types': ['Exe'],          # JT effect type(s)
        'active': True             # This mode is actively optimized
    },
    {
        'mode': 1,                 # Another normal mode where the same JT effect applies
        'state_pairs': [(0, 1)],   # Same state pairs as in mode 0
        'types': ['Exe'],          # Same JT type
        'active': False,           # This mode is inactive; parameters will be copied from mode 0
        'source': 0                # Indicates that parameters from mode 0 should be reused
    }
]

# Define diabatic function types for each mode based on the abinitio data.
# Antimorse means that the asymptotic behavior is at negative displacement.
diab_f_each_mode = ["ho", "ho", "antimorse"]
# For each mode, replicate the chosen diabatic function for all states.
diabatic_functions = [[mode_type] * number_states for mode_type in diab_f_each_mode]

# Create the VCSystem object
system = pyvcham.VCSystem(
    vc_type="linear",
    units="eV",
    number_normal_modes=3,  # 3N-5 vibrational modes
    number_states=2,        # Degenerate ground state (2 states total)
    coupling_with_gs=True,  # Couplings with respect to the ground state non neglected
    # symmetry_point_group="Cs",
    # symmetry_states=["Ap", "App"],
    # symmetry_modes=["App", "Ap", "Ap"],
    symmetry_point_group="C2v",
    symmetry_states=["A1", "B2"],
    symmetry_modes=["A1", "B2", "A1"],
    vib_freq=vib_freq,
    diab_funct=diabatic_functions,
    displacement_vector=all_q,
    database_abinitio=all_data,
    jt_effects=jt_effects
)
# Note: Energy shifts are nonzero because of the calculation issues in the abinitio data (the distances are not optimized for the method). Li3 is
# has not D3h symmetry at the equilibrium distance.

2026-02-03 21:33:47,916 [INFO] pyvcham.symm_vcham: Jahn-Teller effect detected, applying 2 effects
2026-02-03 21:33:47,918 [INFO] pyvcham.vcham_system: Computed vertical energy shifts: [0.         0.01619812]

[6]:

system.symmetry_matrix

[6]:

array([[[0., 2.],
        [2., 0.]],

       [[2., 0.],
        [0., 2.]],

       [[1., 0.],
        [0., 1.]]])

[7]:

import importlib
importlib.reload(pyvcham.lvc)

[7]:

<module 'pyvcham.lvc' from '/Users/emilioerc/Desktop/pyvcham/src/pyvcham/lvc.py'>

[8]:

# Initialize the LVC Hamiltonian for each normal mode
# and optimize the parameters

for mode in range(system.number_normal_modes):
    model = pyvcham.LVCHam(normal_mode=mode, VCSystem=system, nepochs=10000)
    model.initialize_params()
    model.initialize_loss_function()
    model.optimize()

2026-02-03 21:33:47,954 [INFO] pyvcham.lvc:

----- Initializing LVC Hamiltonian Builder -----
2026-02-03 21:33:47,955 [INFO] pyvcham.lvc: Normal mode: 0
2026-02-03 21:33:47,963 [INFO] pyvcham.lvc: JT off-diagonal pairs: [[0, 1]]
2026-02-03 21:33:47,968 [INFO] pyvcham.lvc: Total parameters to optimize: 2
2026-02-03 21:33:47,969 [INFO] pyvcham.lvc: Optimizing mode 0...
<unknown>:3: SyntaxWarning: invalid escape sequence '\c'
2026-02-03 21:33:48,788 [INFO] pyvcham.lvc: Step 0, Loss: 8.535259
2026-02-03 21:33:48,855 [INFO] pyvcham.lvc: Step 100, Loss: 6.840737
2026-02-03 21:33:48,906 [INFO] pyvcham.lvc: Step 200, Loss: 5.148804
2026-02-03 21:33:48,957 [INFO] pyvcham.lvc: Step 300, Loss: 3.461690
2026-02-03 21:33:49,008 [INFO] pyvcham.lvc: Step 400, Loss: 1.787459
2026-02-03 21:33:49,061 [INFO] pyvcham.lvc: Step 500, Loss: 0.202611
2026-02-03 21:33:49,114 [INFO] pyvcham.lvc: Step 600, Loss: 0.000159
2026-02-03 21:33:49,164 [INFO] pyvcham.lvc: Step 700, Loss: 0.000159
2026-02-03 21:33:49,183 [INFO] pyvcham.lvc: Early stopping triggered at step 735 (no improvement for 150 steps). Best loss: 0.000159

2026-02-03 21:33:49,580 [INFO] pyvcham.vcham_system: Stored JT parameters for mode 0: {'on': array([-0.02821373], dtype=float32), 'off': array([0.04355564], dtype=float32)}
2026-02-03 21:33:49,582 [INFO] pyvcham.lvc: Saved JT parameters for mode 0.
2026-02-03 21:33:49,583 [INFO] pyvcham.lvc: Optimization finished – Final reported loss: 0.000159
2026-02-03 21:33:49,583 [INFO] pyvcham.lvc: Optimization completed in 1.21 seconds.
2026-02-03 21:33:49,586 [INFO] pyvcham.lvc:

----- Initializing LVC Hamiltonian Builder -----
2026-02-03 21:33:49,586 [INFO] pyvcham.lvc: Normal mode: 1
2026-02-03 21:33:49,589 [INFO] pyvcham.lvc: JT off-diagonal pairs: []
2026-02-03 21:33:49,595 [INFO] pyvcham.lvc: Total parameters to optimize: 2
2026-02-03 21:33:49,595 [INFO] pyvcham.lvc: Skipping optimization: mode inactive or no parameters.

2026-02-03 21:33:49,845 [INFO] pyvcham.lvc:

----- Initializing LVC Hamiltonian Builder -----
2026-02-03 21:33:49,846 [INFO] pyvcham.lvc: Normal mode: 2
2026-02-03 21:33:49,848 [INFO] pyvcham.lvc: JT off-diagonal pairs: []
2026-02-03 21:33:49,849 [INFO] pyvcham.lvc: n_var_list (non-JT): [3, 3]
2026-02-03 21:33:49,851 [INFO] pyvcham.lvc: Total parameters to optimize: 6
2026-02-03 21:33:49,852 [INFO] pyvcham.lvc: Optimizing mode 2...
2026-02-03 21:33:50,255 [INFO] pyvcham.lvc: Step 0, Loss: 6044.231445
2026-02-03 21:33:50,319 [INFO] pyvcham.lvc: Step 100, Loss: 904.519592
2026-02-03 21:33:50,373 [INFO] pyvcham.lvc: Step 200, Loss: 136.350250
2026-02-03 21:33:50,423 [INFO] pyvcham.lvc: Step 300, Loss: 18.566757
2026-02-03 21:33:50,478 [INFO] pyvcham.lvc: Step 400, Loss: 0.944561
2026-02-03 21:33:50,531 [INFO] pyvcham.lvc: Step 500, Loss: 0.008205
2026-02-03 21:33:50,589 [INFO] pyvcham.lvc: Step 600, Loss: 0.007505
2026-02-03 21:33:50,646 [INFO] pyvcham.lvc: Step 700, Loss: 0.006866
2026-02-03 21:33:50,697 [INFO] pyvcham.lvc: Step 800, Loss: 0.006284
2026-02-03 21:33:50,748 [INFO] pyvcham.lvc: Step 900, Loss: 0.005768
2026-02-03 21:33:50,798 [INFO] pyvcham.lvc: Step 1000, Loss: 0.005311
2026-02-03 21:33:50,850 [INFO] pyvcham.lvc: Step 1100, Loss: 0.004907
2026-02-03 21:33:50,900 [INFO] pyvcham.lvc: Step 1200, Loss: 0.004549
2026-02-03 21:33:50,949 [INFO] pyvcham.lvc: Step 1300, Loss: 0.004229
2026-02-03 21:33:51,001 [INFO] pyvcham.lvc: Step 1400, Loss: 0.003941
2026-02-03 21:33:51,055 [INFO] pyvcham.lvc: Step 1500, Loss: 0.003681
2026-02-03 21:33:51,113 [INFO] pyvcham.lvc: Step 1600, Loss: 0.003442
2026-02-03 21:33:51,167 [INFO] pyvcham.lvc: Step 1700, Loss: 0.003222
2026-02-03 21:33:51,220 [INFO] pyvcham.lvc: Step 1800, Loss: 0.002993
2026-02-03 21:33:51,271 [INFO] pyvcham.lvc: Step 1900, Loss: 0.002708
2026-02-03 21:33:51,323 [INFO] pyvcham.lvc: Step 2000, Loss: 0.002467
2026-02-03 21:33:51,377 [INFO] pyvcham.lvc: Step 2100, Loss: 0.002264
2026-02-03 21:33:51,427 [INFO] pyvcham.lvc: Step 2200, Loss: 0.002090
2026-02-03 21:33:51,478 [INFO] pyvcham.lvc: Step 2300, Loss: 0.001938
2026-02-03 21:33:51,531 [INFO] pyvcham.lvc: Step 2400, Loss: 0.001802
2026-02-03 21:33:51,591 [INFO] pyvcham.lvc: Step 2500, Loss: 0.001678
2026-02-03 21:33:51,650 [INFO] pyvcham.lvc: Step 2600, Loss: 0.001563
2026-02-03 21:33:51,708 [INFO] pyvcham.lvc: Step 2700, Loss: 0.001454
2026-02-03 21:33:51,760 [INFO] pyvcham.lvc: Step 2800, Loss: 0.001351
2026-02-03 21:33:51,811 [INFO] pyvcham.lvc: Step 2900, Loss: 0.001253
2026-02-03 21:33:51,864 [INFO] pyvcham.lvc: Step 3000, Loss: 0.001159
2026-02-03 21:33:51,914 [INFO] pyvcham.lvc: Step 3100, Loss: 0.001069
2026-02-03 21:33:51,967 [INFO] pyvcham.lvc: Step 3200, Loss: 0.000983
2026-02-03 21:33:52,018 [INFO] pyvcham.lvc: Step 3300, Loss: 0.000901
2026-02-03 21:33:52,072 [INFO] pyvcham.lvc: Step 3400, Loss: 0.000823
2026-02-03 21:33:52,132 [INFO] pyvcham.lvc: Step 3500, Loss: 0.000750
2026-02-03 21:33:52,199 [INFO] pyvcham.lvc: Step 3600, Loss: 0.000682
2026-02-03 21:33:52,259 [INFO] pyvcham.lvc: Step 3700, Loss: 0.000619
2026-02-03 21:33:52,320 [INFO] pyvcham.lvc: Step 3800, Loss: 0.000561
2026-02-03 21:33:52,381 [INFO] pyvcham.lvc: Step 3900, Loss: 0.000509
2026-02-03 21:33:52,437 [INFO] pyvcham.lvc: Step 4000, Loss: 0.000462
2026-02-03 21:33:52,491 [INFO] pyvcham.lvc: Step 4100, Loss: 0.000420
2026-02-03 21:33:52,547 [INFO] pyvcham.lvc: Step 4200, Loss: 0.000383
2026-02-03 21:33:52,598 [INFO] pyvcham.lvc: Step 4300, Loss: 0.000352
2026-02-03 21:33:52,701 [INFO] pyvcham.lvc: Step 4400, Loss: 0.000325
2026-02-03 21:33:52,753 [INFO] pyvcham.lvc: Step 4500, Loss: 0.000304
2026-02-03 21:33:52,802 [INFO] pyvcham.lvc: Step 4600, Loss: 0.000286
2026-02-03 21:33:52,853 [INFO] pyvcham.lvc: Step 4700, Loss: 0.000272
2026-02-03 21:33:52,906 [INFO] pyvcham.lvc: Step 4800, Loss: 0.000261
2026-02-03 21:33:52,958 [INFO] pyvcham.lvc: Step 4900, Loss: 0.000253
2026-02-03 21:33:53,011 [INFO] pyvcham.lvc: Step 5000, Loss: 0.000247
2026-02-03 21:33:53,063 [INFO] pyvcham.lvc: Step 5100, Loss: 0.000243
2026-02-03 21:33:53,110 [INFO] pyvcham.lvc: Step 5200, Loss: 0.000240
2026-02-03 21:33:53,158 [INFO] pyvcham.lvc: Step 5300, Loss: 0.000238
2026-02-03 21:33:53,237 [INFO] pyvcham.lvc: Step 5400, Loss: 0.000237
2026-02-03 21:33:53,282 [INFO] pyvcham.lvc: Step 5500, Loss: 0.000236
2026-02-03 21:33:53,328 [INFO] pyvcham.lvc: Step 5600, Loss: 0.000236
2026-02-03 21:33:53,379 [INFO] pyvcham.lvc: Step 5700, Loss: 0.000236
2026-02-03 21:33:53,428 [INFO] pyvcham.lvc: Step 5800, Loss: 0.000236
2026-02-03 21:33:53,477 [INFO] pyvcham.lvc: Step 5900, Loss: 0.000244
2026-02-03 21:33:53,499 [INFO] pyvcham.lvc: Early stopping triggered at step 5943 (no improvement for 150 steps). Best loss: 0.000236

2026-02-03 21:33:53,787 [INFO] pyvcham.lvc: Optimization finished – Final reported loss: 0.000236
2026-02-03 21:33:53,788 [INFO] pyvcham.lvc: Optimization completed in 3.65 seconds.

[9]:

# Define general system information for output
general_data = {
    "molecule": "Li3",
    "calculation_info": "Ground state and first excited state",
    "method": "CASSCF(3,7)",
    "basis": "cc-pvtz",
    "software": "OpenMolcas",
    "additional_info": "-"
}
# Add the reference geometry to the system
ref_geom = "li3_abinitio/geo_v1_00.xyz"
system.add_geometry(ref_geom)

# Define the JSON filename for saving the VCSystem data
filename_json = "results/lvc_li3.json"

# Convert the VCSystem to JSON and save it
pyvcham.utils.VCSystem_to_json(system, general_data, filename_json, rewrite=True)

# Define the output operator filename
output_op = "results/lvc_li3.op"

# Generate the MCTDH operator file from the JSON data
pyvcham.utils.json_to_mctdh(infile=filename_json, outfile=output_op)

2026-02-03 21:34:00,266 [INFO] pyvcham.vcham_system: Reference geometry updated.
2026-02-03 21:34:00,268 [INFO] pyvcham.utils: Warning: Overwriting existing file results/lvc_li3.json
2026-02-03 21:34:00,269 [INFO] pyvcham.utils: Data successfully saved to results/lvc_li3.json
2026-02-03 21:34:00,271 [INFO] pyvcham.utils: No interactions found.
2026-02-03 21:34:00,271 [INFO] pyvcham.utils: MCTDH operator file successfully written to: results/lvc_li3.op