TDSE - fmauger1/QMol-grid GitHub Wiki
Time-dependent Schrödinger equation (TDSE)
The QMol-grid package can simulate the time evolution of the TDSE
$$ i\partial _t \psi (x;t) = \hat{\mathcal{H}} \psi (x;t),~~~~~~(1) $$
where $\hat{\mathcal{H}}$ is the TDSE Hamiltonian operator, potentially including external driving fields. We refer to the corresponding documentation pages for details about how to properly define Schrödinger-equation (SE) models and available parameters.
TDSE components
For clarity and streamlined future development of TDSE capabilities in the QMol-grid package, we sort TDSE 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 near the edges of the simulation domain, therefore avoiding (reducing) reflections on the boundaries
QMol_TDSE_abs_CAPenables absorbing boundaries using a complex absorbing potential (CAP) at the edges of the domainQMol_TDSE_abs_maskenables absorbing boundaries using a mask function at the edges of the domain
Memory profiler
The choice of the Schrödinger-equation model together with the set of observables to save, and how often, during the TDSE propagation can greatly affect the memory resources required to carry out simulation
QMol_SE_profilerprovides estimates of the memory requirements for a given TDSE simulation
See the Schrödinger-equation documentation for tips and examples as to how to optimize simulations.
TDSE propagators
The QMol-grid package supports various types of TDSE propagators
General-purpose propagators
General-purpose propagators are the most versatile schemes with few, if any, restrictions on the Schrödinger-equation model but require more resources
- No general-purpose propagators have been implemented yet
Propagators for basis-set models
For Schrödinger-equation models that use a basis-set discretization scheme
- No basis-set propagators have been implemented yet
Symplectic split-operator propagators
For Schrödinger-equation models discretized on a grid
QMol_TDSE_SSO_2Simplements the 2nd order Strang splitting (a.k.a. Verlet) schemeQMol_TDSE_SSO_4BMimplements the 4th order Blanes and Moan optimized schemeQMol_TDSE_SSO_4FRimplements the 4th order Forest-Ruth symplectic schemeQMol_TDSE_SSO_6BMimplements the 6th order Blanes and Moan optimized scheme
Symplectic split-operator TDSE propagators are build from the Hamiltonian description of the dynamics, similarly to the discussion in [Mauger 2024].
Propagation results
All TDSE propagators in the QMol-Grid package take a Schrödinger-equation 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).
- Schrödinger-equation and wave function energies: Tracking the Schrödinger-equation 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 Schrödinger-equation-energy feature keeps track of the variations in the kinetic, potential, and autonomization [Mauger 2024] components. The instantaneous wave function energy is defined as $
\left\langle \psi (t)\left|\hat{\mathcal{H}} \right|\psi (t)\right\rangle$ . Both the Schrödinger-equation and wave-function energies are directly computed from the Schrödinger-equation-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.
- Wave functions and projection of the wave functions: Each can individually be enabled and define its own sampling time. 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. Alternatively, saving the projection of the wave functions onto a specific basis generally has more tractable memory requirements, to the cost of possibly losing some details of the wave functions and the time to compute the projections.
Additionally, one can
- Add the external-field information to each of the requested output results structures.
- Save the intermediate-time Schrödinger-equation-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 TDSE 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.
- 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
TDSE 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 high-order harmonic generation (HHG) spectrum computed from the dipole, dipole-velocity, and dipole-acceleration signals. First, we define the Schrödinger-equation model and the ground-state configuration
% Molecular model
Va_1 = QMol_Va_softCoulomb('atom','A','charge',1,'position',-2.5,'softeningParameter',1.5);
Va_2 = QMol_Va_softCoulomb('atom','A','charge',1,'position', 2.5,'softeningParameter',1.5);
% Schrödinger-equation model
V = QMol_SE_V('atom',{Va_1,Va_2});
SE = QMol_SE( ...
'xspan', -50:.1:50, ... small prime factors
'numberWaveFunction', 2, ...
'potential', V);
% Ground state
GSS = QMol_SE_eigs;
GSS.computeGroundState(SE);
Then, we define the external field -- ramped up for 2 laser cycles followed by a constant-amplitude plateau -- and TDSE propagator using the Forest-Ruth symplectic split-operator scheme
% External field
E0 = convertLaser.intensity2amplitude(2e13);
w = convertLaser.wavelength2frequency(2000e-9);
T = 2*2*pi/w;
EField = @(t) E0*sin(w*t).*sin(.5*pi*min(t,T)/T).^2;
t = 0:2*pi/w:T+5*2*pi/w;
% Absorbing boundaries
ABC = QMol_TDSE_abs_CAP('shape','sin^1/8','length',10,'amplitude',.5);
% Propagator
dt = 5e-2; % Propagation time step
TDSE = QMol_TDSE_SSO_4FR( ...
'time', t, ...
'timeStep', dt, ...
'electricField', EField, ...
'saveExternalField', true, ...
'saveDipole', true, ...
'saveDipoleTime', -floor(.5/dt), ... equally-spaced time sampling close to 0.5 a.u.
'saveDipoleWaveFunctionIndex', 'all', ...
'saveDipoleVelocity', true, ...
'saveDipoleVelocityTime', 'dipole', ... same sampling as for the dipole signal
'saveDipoleVelocityWaveFunctionIndex', 'all', ...
'saveDipoleAcceleration', true, ...
'saveDipoleAccelerationTime', 'dipole', ... same sampling as for the dipole signal
'saveDipoleAccelerationWaveFunctionIndex', 'all', ...
'absorbingBoundary', ABC );
TDSE.propagate(SE);
Finally, we plot the HHG spectra computed from the dipole, dipole-velocity, and dipole-acceleration signals from each of the wave functions
% Spectral analysis
W = zeros(size(TDSE.outDipole.time));
ind = TDSE.outDipole.time > T;
W(ind) = sin(pi*(0:sum(ind)-1)/(sum(ind)-1)).^2; % Hann window
F = fourierTool.fftGrid(TDSE.outDipole.time);
D = NaN(size(TDSE.outDipole.waveFunction_x));
V = NaN(size(TDSE.outDipoleVelocity.waveFunction_x));
A = NaN(size(TDSE.outDipoleAcceleration.waveFunction_x));
for k = 1:SE.numberWaveFunction
D(k,:) = fft(TDSE.outDipole.waveFunction_x(k,:).*W).*F.^2;
V(k,:) = fft(TDSE.outDipoleVelocity.waveFunction_x(k,:).*W).*F;
A(k,:) = fft(TDSE.outDipoleAcceleration.waveFunction_x(k,:).*W);
end
% Plot results
ind = F > 0; % Positive frequencies
F = F(ind)/w;
D = abs(D(:,ind)).^2;
V = abs(V(:,ind)).^2;
A = abs(A(:,ind)).^2;
figure
subplot(2,1,1), set(gca,'yscale','log'), hold on
plot(F,D(1,:),'-' ,'LineWidth',3,'DisplayName','dipole')
plot(F,V(1,:),'--' ,'LineWidth',3,'DisplayName','velocity')
plot(F,A(1,:),':' ,'LineWidth',3,'DisplayName','acceleration')
hold off
xlim([0 120]), xlabel('harmonic order')
ylim([1e-6 1]*1e-6), ylabel('intensity (arb. units)')
title('wave function 1')
legend show
subplot(2,1,2), set(gca,'yscale','log'), hold on
plot(F,D(2,:),'--','LineWidth',2,'DisplayName','dipole')
plot(F,V(2,:),'--' ,'LineWidth',3,'DisplayName','velocity')
plot(F,A(2,:),':' ,'LineWidth',3,'DisplayName','acceleration')
hold off
xlim([0 120]), xlabel('harmonic order')
ylim([1e-6 1]*1e-4), ylabel('intensity (arb. units)')
title('wave function 2')
legend show
Schrödinger-equation and wave function energies
Here we check the conservation of the Schrödinger-equation energy in where one of the molecular centers is suddenly shifted and subjected to an external driving field in the dipole approximation and length gauge. First, we define the Schrödinger-equation model and compute the (neutral) ground-state configuration
% Molecular model
Va_1 = QMol_Va_softCoulomb('atom','X','charge',1.2,'position',-4);
Va_2 = QMol_Va_softCoulomb('atom','A','charge',.8,'position', 0);
Va_3 = QMol_Va_softCoulomb('atom','A','charge',.8,'position', 3);
% Schrödinger-equation model
V = QMol_SE_V('atom',{Va_1,Va_2,Va_3});
SE = QMol_SE( ...
'xspan', -50:.1:50, ... small prime factors
'numberWaveFunction', 3, ...
'potential', V);
% Ground state
GSS = QMol_SE_eigs;
GSS.computeGroundState(SE);
Next, we move the "X" atomic center and define the external driving field
% Move the X center
Va_1.set('position',-3.5);
% External potential-vector field
E0 = convertLaser.intensity2amplitude(1e13);
w = convertLaser.wavelength2frequency(800e-9);
EField = @(t) E0*sin(w*t);
Finally we propagate the TDSE dynamics using the Forest-Ruth symplectic split-operator scheme
% TDSE propagator
TDSE = QMol_TDSE_SSO_4FR( ...
'time', 0:10:100, ...
'timeStep', 2e-2, ...
'electricField', EField, ...
'externalFieldGauge', 'length', ...
'saveExternalField', true, ...
'saveEnergySE', true, ...
'saveEnergySETime', 1, ... save SE energy every 1 a.u.
'absorbingBoundary', QMol_TDSE_abs_CAP);
TDSE.propagate(SE);
With this, we can plot the variation of the total Schrödinger-equation energy and its various components.
% Plot results
figure
subplot(2,1,1), hold on
plot(TDSE.outEnergySE.time,TDSE.outEnergySE.total-TDSE.outEnergySE.total(1),...
'-','LineWidth',2,'DisplayName','total');
hold off
xlim(TDSE.outEnergySE.time([1 end]))
xlabel('time (a.u.)')
ylabel('error (a.u.)')
title('total Schrödinger-equation energy')
subplot(2,1,2), hold on
plot(TDSE.outEnergySE.time,sum(TDSE.outEnergySE.kinetic,1)-sum(TDSE.outEnergySE.kinetic(:,1)),...
'-','LineWidth',2,'DisplayName','kinetic')
plot(TDSE.outEnergySE.time,sum(TDSE.outEnergySE.potential,1)-sum(TDSE.outEnergySE.potential(:,1)),...
'-','LineWidth',2,'DisplayName','potential')
plot(TDSE.outEnergySE.time,sum(TDSE.outEnergySE.externalField,1)-sum(TDSE.outEnergySE.externalField(:,1)),...
'-','LineWidth',2,'DisplayName','external field')
plot(TDSE.outEnergySE.time,TDSE.outEnergySE.autonomization-TDSE.outEnergySE.autonomization(1),...
'-','LineWidth',2,'DisplayName','autonomization')
hold off
xlim(TDSE.outEnergySE.time([1 end]))
xlabel('time (a.u.)')
ylabel('variation (a.u.)')
title('Schrödinger-equation energy components')
legend show
The sudden change in the Schrödinger-equation energy error is probably caused by small portions of the wave function being absorbed on the edges of the simulation domain.
Ionization signal
Here we track the ionization signal of a molecular model driven by an external field in the dipole approximation and velocity gauge. First, we define the Schrödinger-equation model and compute the (neutral) ground-state configuration
% Molecular model
Va_1 = QMol_Va_softCoulomb('atom','X','charge',.7,'position',-4);
Va_2 = QMol_Va_softCoulomb('atom','A','charge',.6,'position', 0);
Va_3 = QMol_Va_softCoulomb('atom','A','charge',.6,'position', 3);
% Schrödinger-equation model
V = QMol_SE_V('atom',{Va_1,Va_2,Va_3});
SE = QMol_SE( ...
'xspan', -50:.1:50, ... small prime factors
'numberWaveFunction', 3, ...
'potential', V);
% Ground state
GSS = QMol_SE_eigs;
GSS.computeGroundState(SE);
Next, we define the driving potential-vector field and propagate the TDSE dynamics using the Forest-Ruth symplectic split-operator scheme
% External potential-vector field
E0 = convertLaser.intensity2amplitude(1e14);
w = convertLaser.wavelength2frequency(800e-9);
T = 5*2*pi/w;
AField = @(t) E0/w*sin(w*t).*cos(.5*pi*t/T).^2;
% TDSE propagator
ABC = QMol_TDSE_abs_CAP('length',40);
TDSE = QMol_TDSE_SSO_4FR( ...
'time', -T:2*pi/w:T, ...
'timeStep', 5e-2, ...
'potentialVector', AField, ...
'externalFieldGauge', 'velocity', ...
'saveExternalField', true, ...
'saveIonization', true, ...
'saveIonizationWaveFunctionIndex', 'all', ...
'saveIonizationTime', 2, ... save ionization yield every 2 a.u.
'absorbingBoundary', ABC);
TDSE.propagate(SE);
Finally, we plot the ionization rates for each of the wave functions
% Plot results
T = .5*(TDSE.outIonization.time(2:end)+TDSE.outIonization.time(1:end-1));
figure
subplot(2,1,1), hold on
plot(TDSE.outIonization.time,TDSE.outIonization.potentialVector,...
'-','LineWidth',2,'DisplayName','A field');
hold off
xlim(TDSE.outIonization.time([1 end]))
xlabel('time (a.u.)')
ylabel('field (a.u.)')
title('potential vector')
subplot(2,1,2), hold on
plot(T,(TDSE.outIonization.waveFunction(1,2:end)-TDSE.outIonization.waveFunction(1,1:end-1))/(T(2)-T(1)),...
':','LineWidth',2,'DisplayName','wfcn #1')
plot(T,(TDSE.outIonization.waveFunction(2,2:end)-TDSE.outIonization.waveFunction(2,1:end-1))/(T(2)-T(1)),...
':','LineWidth',2,'DisplayName','wfcn #2')
plot(T,(TDSE.outIonization.waveFunction(3,2:end)-TDSE.outIonization.waveFunction(3,1:end-1))/(T(2)-T(1)),...
':','LineWidth',2,'DisplayName','wfcn #3')
hold off
xlim(TDSE.outIonization.time([1 end]))
xlabel('time (a.u.)')
ylabel('variation (a.u.)')
title('Ionization rate')
legend show
Installable output functions
We use an installable output function to store the current density in each of the wave functions, down-sampled and over a reduced domain, to illustrate how this feature can be used as an alternative to storing the entire wave functions. First, we define the Schrödinger-equation model and compute the ground-state configuration
% Molecular model
Va_1 = QMol_Va_softCoulomb('atom','X','charge',1.2,'position',-4);
Va_2 = QMol_Va_softCoulomb('atom','A','charge',.8,'position', 0);
Va_3 = QMol_Va_softCoulomb('atom','A','charge',.8,'position', 3);
% Schrödinger-equation model
V = QMol_SE_V('atom',{Va_1,Va_2,Va_3});
SE = QMol_SE( ...
'xspan', -50:.1:50, ... small prime factors
'numberWaveFunction', 3, ...
'potential', V);
% Ground state
GSS = QMol_SE_eigs;
GSS.computeGroundState(SE);
Next, we move the "X" atomic center, define the output function handle and propagate the TDSE dynamics using the Forest-Ruth symplectic split-operator scheme
% Move the X center
Va_1.set('position',-3.5);
% Installable output functions
ind = SE.disc.x >= -15 & SE.disc.x <= 15 & mod(1:numel(SE.disc.x),2) == 1; % 1/2 grid points between -15 and 15
cFcn = @(wfcn) real(-1i*conj(wfcn).*ifft(SE.disc.D.*fft(wfcn)));
fcn = @(wfcn,t) [subsref(cFcn(wfcn.wfcn(:,1)),struct('type','()','subs',{{ind}})), ...
subsref(cFcn(wfcn.wfcn(:,2)),struct('type','()','subs',{{ind}}))];
% fcn = @(wfcn,t) [real(-1i*conj(wfcn.wfcn(:,1)).*ifft(SE.disc.D.*fft(wfcn.wfcn(:,1)))), ...
% real(-1i*conj(wfcn.wfcn(:,2)).*ifft(SE.disc.D.*fft(wfcn.wfcn(:,2))))];
% TDSE propagator
ABC = QMol_TDSE_abs_CAP('length',20);
TDSE = QMol_TDSE_SSO_4FR( ...
'time', 0:10:100, ...
'timeStep', 2e-2, ...
'saveOutputFunction', fcn, ...
'saveOutputFunctionTime', 1, ... save current every 1 a.u.
'absorbingBoundary', ABC);
TDSE.propagate(SE);
With this, we can plot the results of the installable output functions
% Plot results
figure
subplot(2,1,1), hold on
imagesc(TDSE.outOutputFunction.time,SE.x(ind),...
reshape(TDSE.outOutputFunction.result(:,1,:),sum(ind),numel(TDSE.outOutputFunction.time)))
hold off
xlim(TDSE.outOutputFunction.time([1 end]))
ylim([-10 10])
xlabel('time (a.u.)')
ylabel('position (a.u.)')
title('wave function 1 current')
colorbar vert
subplot(2,1,2), hold on
imagesc(TDSE.outOutputFunction.time,SE.x(ind),...
reshape(TDSE.outOutputFunction.result(:,2,:),sum(ind),numel(TDSE.outOutputFunction.time)))
hold off
xlim(TDSE.outOutputFunction.time([1 end]))
ylim([-10 10])
xlabel('time (a.u.)')
ylabel('position (a.u.)')
title('wave function 2 current')
colorbar vert
Notes for developers
The QMol_TDSE abstract class defines the common interface and run-time documentation for TDSE propagators -- see its documentation page for details on how to implement new TDSE propagators.
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).
Notes
The results displayed in this documentation page were generated using version 01.21 of the QMol-grid package.
- TDSE features were introduced in version 01.20