TDCI - fmauger1/QMol-grid GitHub Wiki
Time-dependent configuration interaction (TDCI)
The QMol-grid package can simulate the time evolution of TDCI
$$ i\partial_t \psi (x;t)=\hat{\mathcal{H}} \psi (x;t), ~~~~~~ (1) $$
where $\hat{\mathcal{H}}$ is the TDCI Hamiltonian operator, potentially including external driving fields. We refer to the corresponding documentation pages for details about how to properly define configuration interaction (CI) models and available parameters.
Note: TDCI is a recently added feature to the QMol-grid package. Expect some adjustments as we get more experience running TDCI for various systems.
TDCI components
For clarity and streamlined future development of TDCI capabilities in the QMol-grid package, we sort TDCI component (classes) into different topical groups.
External driving field
External driving fields provide the explicit time dependence in the Hamiltonian operator
QMol_extField_dipoleis a container class for external driving fields in the dipole approximation
Absorbing boundaries
Absorbing boundaries emulate the exit from the simulation domain of outgoing wave packets from the simulation domain
QMol_TDCI_abs_CAPenables user-defined absorbing boundaries using a complex absorbing potential (CAP)QMol_TDCI_abs_CAP_spacebuilds the CAP based on the spatial overlap between of orbital basis in the CI model and spatial mask for the absorbing boundariesQMol_TDCI_abs_maskenables user-defined mask-type absorbing boundariesQMol_TDCI_abs_mask_spacebuilds the mask based on the spatial overlap between of orbital basis in the CI model and spatial mask for the absorbing boundaries
Mask-type absorbing boundaries are applied at the end of each propagation step and are not scaled with the time step while CAPs are integrated as a complex component to the potential and automatically scaled with the propagation time step.
Memory profiler
The choice of the CI model together with the set of observables to save, and how often, during the TDCI propagation can greatly affect the memory resources required to carry out simulation
QMol_CI_profilerprovides estimates of the memory requirements for a given TDCI simulation
See the CI documentation for tips and examples as to how to optimize simulations.
TDCI propagators
Currently, the QMol-grid package only support symplectic split-operator TDCI propagators
Symplectic split-operator propagators
QMol_TDCI_SSO_2Simplements the 2nd order Strang splitting (a.k.a. Verlet) schemeQMol_TDCI_SSO_4BMimplements the 4th order Blanes and Moan optimized schemeQMol_TDCI_SSO_4FRimplements the 4th order Forest-Ruth symplectic schemeQMol_TDCI_SSO_6BMimplements the 6th order Blanes and Moan optimized scheme
Symplectic split-operator TDCI propagators are build from the Hamiltonian description of the dynamics, similarly to the discussion in [Mauger 2024] and with the two components in the split being the CI and dipole-coupling matrices.
Propagation results
All TDCI propagators in the QMol-Grid package take a CI model object as input and evolve it over the user-specified time interval and action of external driving field(s). During the propagation, they also support on-the-fly computation and storage of the following observables:
- Dipole, dipole velocity, and dipole acceleration: Each can individually be enabled and define its own sampling time and provide the signal contribution from (selected) individual wave functions. The dipole velocity and acceleration are self-consistently computed using analytical formula (i.e., not from taking the time derivative of the dipole).
- CI and wave function energies: Tracking the CI energy can be used as a proxy for evaluating the accuracy of the numerical simulations (since the energy should be conserved). In addition to the total energy, activating the CI-energy feature keeps track of the variations in field-free Hamiltonian (kinetic + potential), external-field contribution, and autonomization [Mauger 2024] components. The instantantous wave function energy is defined as $
\left\langle \psi (t)\left|\hat{\mathcal{H}} \right|\psi (t)\right\rangle$ . Both the CI and wave-function energies are directly computed from the CI-object getEnergy method -- see the corresponding documentation page for details. - Ionization signal: Keep track of how much density is absorbed on the boundaries in each of the wave-function channels - warning: this is a mostly untested feature.
- Wave functions: Saving the wave functions is generally discouraged as they generally lead to very large memory requirements; saving the Schrödinger-equation object into separate MATLAB files or using the installable output functions are recommended, where possible -- see below.
Additionally, one can
- Add the external-field information to each of the requested output results structures.
- Save the intermediate-time CI-model object into a separate MATLAB file. This is the best option if the intermediate wave functions are required for post-processing.
- Save the results of an installable output function: An installable output function enables on-the-fly custom analysis and saving of the TDCI dynamics. The output function must return a scalar, vector, or array of arbitrary size and shape but must be consistent throughout the propagation. The output results are appended in an array with one extra dimension compared to the output, and the last index labels the output time - warning: this is a mostly untested feature.
- Generate a restart file from which the propagation can be restarted in case it gets stopped before the end of the propagation.
All output can individually be enabled and define their own sampling times, which need not be uniformly distributed in time. The respective results are stored in output structure properties in the TDSE-propagator object; each of these structures keep a record of the precise time at which the results have been sampled (as these may be slightly different from the requested input -- see the respective propagator documentation pages for details).
Examples
TDCI propagators provides general support for on-the-fly computation and storage of various commonly used observables and output results. We give sample scripts for some of these.
Dipole, dipole velocity, and dipole acceleration signals
We compare the absorption cross section computed from the dipole, dipole-velocity, and dipole-acceleration signals. First, we define the molecular model and DFT calculation to build the orbital basis
% Molecular model
x = -50:.1:50;
Va = QMol_Va_softCoulomb('name','atom','charge',4);
% DFT calculation of the orbital basis
Vext = QMol_DFT_Vext('atom',{Va});
Vh = QMol_DFT_Vh_conv;
Vxc = {QMol_DFT_Vx_LDA_soft,QMol_DFT_Vc_LDA_soft};
DFT = QMol_DFT_spinRes( ...
'xspan', x, ...
'occupation', [2 2 0 0 0], ...
'externalPotential', Vext, ...
'HartreePotential', Vh, ...
'exchangeCorrelationPotential', Vxc, ...
'selfInteractionCorrection', 'ADSIC' );
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(DFT);
yielding
(...)
=== Orbital energies =====================================================
Orbital Occ. (elec.) Energy (-eV) Error(a.u.)
------- ------------ ------------ -----------
1 2.00 42.404 9.284e-13
2 2.00 15.014 1.460e-12
3 0.00 6.515 2.418e-12
4 0.00 3.489 2.792e-12
5 0.00 2.030 8.338e-13
------------------------------------------------------------------------
=== DFT-component energies ===============================================
Component Energy (a.u.) Energy (eV)
----------- ------------- -------------
Kinetic 1.378 37.490
External -11.652 -317.061
Hartree 3.070 83.529
Exch.-corr. -0.211 -5.735
----------- ------------- -------------
Total -7.415 -201.778
----------------------------------------------
All the occupied and virtual orbitals have converged and can thus be used to build the CI model, which we do next (at the CISD level)
% Set the CI model
CI = QMol_CI_conv(DFT);
CI.set('type','CISD');
CI.setConfigurationBasis;
CI.computeCImatrix;
CI.sparsify;
yielding
(... [building the configuration basis] ...)
##########################################################################
=== Set the configuration-state basis ====================================
* Initialization DONE
* Configuration-state basis
Type = CISD
Total spin = 0
Reference(s) = -1 -2 1 2
Active = 3 4 5
-5 -4 -3
55 configuration-state basis elements
V-01.23.002 (06/10/2025) G. Visentin & F. Mauger
##########################################################################
(... [calculating the CI matrix] ...)
##########################################################################
=== CI-matrix calculation ================================================
Calculation progress (%)
Calculation step 0 20 40 60 80 100
---------------------- ----------------------------------------------
Core integrals ||||||||||||||||||||||||||||||||||||||||||||||
Dipole coupling elements ||||||||||||||||||||||||||||||||||||||||||||||
Two-electron integrals ||||||||||||||||||||||||||||||||||||||||||||||
CI matrix ||||||||||||||||||||||||||||||||||||||||||||||
------------------------------------------------------------------------
CI-matrix calculation complete
##########################################################################
Note that we sparsify the CI matrix, which speeds up TDCI calculations. Next, we identify bright excited states (dipole coupled to the CI ground state) and grond initial state for the cross-section calculation
% Identify bright and dark excited states
CI.computeGroundState(10);
E_CI = CI.getEnergy; % ground- and excited-state energies
DX = sum(CI.waveFunction(:,1).*CI.DXmatrix*CI.waveFunction,1); % dipole coupling strength
% Ground state
CI.computeGroundState(1);
Then, we define the external field with a $\cos^4$ envelope, total duration of 3 fs, and 12 eV central frequency, and the TDCI propagator using the Forest-Ruth symplectic split-operator scheme
% External field
E_0 = convertLaser.intensity2amplitude(1e11);
T = convertUnit.second2au(3e-15);
t_0 = .6*T;
w_0 = convertUnit.ev2au(12);
EF = @(t) E_0 * cos(pi/T*min(max(t-t_0,-.5*T),.5*T)).^4 ...
.* cos(w_0*(t-t_0));
t = 0:100:1000;
% Absorbing boundaries
ABC = QMol_TDCI_abs_CAP_space('shape','sin^2','distance',20,'length',10,'amplitude',1);
% Propagator
TDCI = QMol_TDCI_SSO_4FR( ...
'time', t, ...
'timeStep', 2e-2, ...
'electricField', EF, ...
'saveExternalField', true, ...
'saveDipole', true, ... <x> position expectation value
'saveDipoleTime', 0.1, ...
'saveDipoleWaveFunctionIndex', 'all', ...
'saveDipoleVelocity', true, ... <p> velocity expectation value
'saveDipoleVelocityTime', 'dipole', ... same sampling as for the dipole signal
'saveDipoleVelocityWaveFunctionIndex', 'all', ...
'saveDipoleAcceleration', true, ... <a> acceleration expectation value
'saveDipoleAccelerationTime', 'dipole', ...
'saveDipoleAccelerationWaveFunctionIndex','all',...
'absorbingBoundary', ABC); % absorbing boundaries
TDCI.propagate(CI); % Run TDCI calculation
Finally, we plot the cross section computed from the dipole, dipole-velocity, and dipole-acceleration signals
% Recover data
t = TDCI.outDipole.time; % time sampling
E = TDCI.outDipole.electricField; % electric field
X = TDCI.outDipole.waveFunction_x; % <x>
P = TDCI.outDipoleVelocity.waveFunction_x; % <p>
A = TDCI.outDipoleAcceleration.waveFunction_x; % <a>
% Fourier transforms
W = cos(.5*pi*(t-t(1))/(t(end)-t(1))).^2; % window
N = max(numel(t),1e5); % increase spectral resolution
freq = fourierTool.fftGrid(t,N); % frequency sampling
FE = fft(W.*E,N);
FX = fft(W.*(X-X(1)),N); % remove permanent dipole
FP = fft(W.*(P-P(1)),N);
FA = fft(W.*(A-A(1)),N);
% Absorption cross sections
Dw = convertUnit.ev2au(5); % cross-section bandwidth
ind = (freq >= w_0-.5*Dw) & (freq <= w_0+.5*Dw); % selected frequencies
freq = freq(ind);
s_X = imag(FX(ind)https://github.com/fmauger1/QMol-grid/wiki/FE(ind)).*freq;
s_P = -real(FP(ind)https://github.com/fmauger1/QMol-grid/wiki/FE(ind));
s_A = -imag(FA(ind)https://github.com/fmauger1/QMol-grid/wiki/FE(ind))https://github.com/fmauger1/QMol-grid/wiki/freq;
% Plot absorption cross section
freq = convertUnit.au2ev(freq);
figure
set(gca,'box','on','FontSize',12,'LineWidth',2)
hold on
plot(freq,s_X,'-' ,'LineWidth',2,'DisplayName','\sigma_{<x>}')
plot(freq,s_P,'--','LineWidth',2,'DisplayName','\sigma_{<p>}')
plot(freq,s_A,':' ,'LineWidth',2,'DisplayName','\sigma_{<a>}')
xline(convertUnit.au2ev(E_CI(abs(DX)>1e-10)-E_CI(1)),'k-','LineWidth',2,'DisplayName','Th.')
xlabel('frequency (eV)')
xlim(freq([1 end]))
ylabel('abs. cross section (arb. units)')
legend show
CI and wave function energies
Here we check the conservation of the CI energy for a driven model starting from a non-stationarry initial wave function. First, we model and use DFT to calculate the molecular orbital basis for the CISD expansion
% Molecular model
X_C = [-8.95, -6.65, -3.75, -1.45, 1.45, 3.75, 6.65, 8.95];
Vext = QMol_DFT_Vext('atom',{ ...
QMol_Va_softCoulomb('name','C1','position',X_C(1)), ...
QMol_Va_softCoulomb('name','C2','position',X_C(2)), ...
QMol_Va_softCoulomb('name','C3','position',X_C(3)), ...
QMol_Va_softCoulomb('name','C4','position',X_C(4)), ...
QMol_Va_softCoulomb('name','C5','position',X_C(5)), ...
QMol_Va_softCoulomb('name','C6','position',X_C(6)), ...
QMol_Va_softCoulomb('name','C7','position',X_C(7)), ...
QMol_Va_softCoulomb('name','C8','position',X_C(8)) ...
});
occ = [2 2 2 2 0 0 0 0];
% DFT calculation of the molecular orbital basis
Vh = QMol_DFT_Vh_conv;
Vxc = {QMol_DFT_Vx_LDA_soft,QMol_DFT_Vc_LDA_soft};
DFT = QMol_DFT_spinRes(...
'xspan', -100:.2:100, ...
'occupation', occ, ...
'externalPotential', Vext, ...
'HartreePotential', Vh, ...
'exchangeCorrelationPotential', Vxc, ...
'selfInteractionCorrection', 'ADSIC' );
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(DFT);
% CI model
CI = QMol_CI_conv(DFT,'type','CISD');
CI.setConfigurationBasis;
CI.computeCImatrix;
CI.sparsify;
For the initial condition, we consider the single Slater determinant formed by the HOMO $\to$ LUMO excitation in the up spin channel
% Initial condition
psi_0 = all(CI.configurationBasis == [-1 -2 -3 -4 1 2 3 5],2);
CI.set('waveFunction',psi_0+0); % Convert logical to double
CI.initialize;
Next, we define the external driving field
% External field
E0 = convertLaser.intensity2amplitude(1e13);
w = convertLaser.wavelength2frequency(800e-9);
EField = @(t) E0*sin(w*t);
Finally we propagate the TDCI dynamics using the Forest-Ruth symplectic split-operator scheme
% TDCI propagation
TDCI = QMol_TDCI_SSO_4FR( ...
'time', 0:10:100, ...
'timeStep', 1e-2, ...
'electricField', EField, ...
'externalFieldGauge', 'length', ...
'saveExternalField', true, ...
'saveEnergyCI', true, ...
'saveEnergyCItime', 1);
TDCI.propagate(CI);
With this, we can plot the variation of the total CI energy and its various components.
% Plot results
figure
subplot(2,1,1), set(gca,'box','on','FontSize',12,'LineWidth',2)
hold on
plot(TDCI.outEnergyCI.time,TDCI.outEnergyCI.total-TDCI.outEnergyCI.total(1),'LineWidth',2,'DisplayName','total')
hold off, legend show
xlim(TDCI.outEnergyCI.time([1 end]))
ylabel('error (a.u.)')
title('total CI energy')
subplot(2,1,2), set(gca,'box','on','FontSize',12,'LineWidth',2)
hold on
plot(TDCI.outEnergyCI.time,TDCI.outEnergyCI.Hamiltonian-TDCI.outEnergyCI.Hamiltonian(1),'LineWidth',2,'DisplayName','kinetic + potential')
plot(TDCI.outEnergyCI.time,TDCI.outEnergyCI.externalField,'-','LineWidth',2,'DisplayName','external field')
plot(TDCI.outEnergyCI.time,TDCI.outEnergyCI.autonomization,'-','LineWidth',2,'DisplayName','autonomization')
hold off, legend show
xlim(TDCI.outEnergyCI.time([1 end]))
ylabel('variation (a.u.)')
title('CI energy components')
Note: to keep the example as simple as possible, we did not include any absorbing boundary conditions. It is generally advisable to have some when considering externally driven systems, to avoid spurious reflections of the wave function due to the finite dimension of the CI model.
References
[Mauger 2024] F. Mauger, C. Chandre, M.B. Gaarde, K. Lopata, and K.J. Schafer, "Hamiltonian formulation and symplectic split-operator schemes for time-dependent density-functional-theory equations of electron dynamics in molecules," Communications in Nonlinear Science and Numerical Simulation 129, 107685 (2024).
[Visentin 2025] G. Visentin and F. Mauger, "Configuration-interaction calculations with density-functional theory molecular orbitals for modeling valence- and core-excited states in molecules," arXiv:2509.08245 (2025).
Notes
The results displayed in this documentation page were generated using version 01.23 of the QMol-grid package.
- TDCI features were introduced in version 01.23