TDHF - fmauger1/QMol-grid GitHub Wiki
Time-dependent Hartree Fock (TDHF)
Time-dependent Hartree Fock corresponds to the time propagation of molecular orbitals under the Hartree-Fock mean-field Hamiltonian, which is obtained by running density-functional theory (DFT) computations with (i) an exact-exchange functional, (ii) no correlation functional, and (iii) no self-interaction correction (SIC) in the Hartree functional.
Examples
Because of the implicit (exact exchange) functional, TDHF require the use of general-purpose TDDFT propagators, which are implemented via extended phase-space symplectic split-operator schemes [Mauger 2025]; see the TDDFT documentation for a list of available schemes.
Dipole, dipole velocity, and dipole acceleration signals
We compare the high-order harmonic generation (HHG) spectrum computed from the dipole signals. First, we define the HF model and the ground-state configuration
% Molecular model
A1 = QMol_Va_softCoulomb('atom','A','charge',3,'position',-1.5);
A2 = QMol_Va_softCoulomb('atom','A','charge',3,'position', 1.5);
Vmol = QMol_DFT_Vext('atom',{A1,A2});
% HF model
HF = QMol_DFT_spinRes( ...
'xspan', -50:.1:50, ... small prime factors
'occupation', [2 2 2], ...
'externalPotential', Vmol, ...
'HartreePotential', QMol_DFT_Vh_conv, ...
'exchangeCorrelationPotential', QMol_DFT_Vx_XX_conv);
% DFT ground state
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(HF);
Then, we define the external field -- ramped up for 2 laser cycles followed by a constant-amplitude plateau -- and TDHF 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_TDDFT_abs_CAP('shape','sin^1/8','length',10,'amplitude',.5);
% Propagator
dt = 5e-2; % Propagation time step
TDHF = QMol_TDDFT_ESSO_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.
'saveDipoleOrbitalIndex', 'all', ...
'absorbingBoundary', ABC );
TDHF.propagate(HF);
Finally, we plot the HHG spectrum computed from the dipole and the contributions from the various HF orbitals
% Spectral analysis
W = zeros(size(TDHF.outDipole.time));
ind = TDHF.outDipole.time > T;
W(ind) = sin(pi*(0:sum(ind)-1)/(sum(ind)-1)).^2; % Hann window
% Total signal
F = fourierTool.fftGrid(TDHF.outDipole.time);
D = fft(W.*TDHF.outDipole.total).*F.^2;
% Orbital-resolved signals
O = NaN(size(TDHF.outDipole.orbital_x));
for k = 1:numel(HF.occupation)
O(k,:) = fft(TDHF.outDipole.orbital_x(k,:).*W).*F.^2;
end
% Plot results
ind = F > 0; % Positive frequencies
F = F(ind)/w;
D = abs(D(ind)).^2;
O = abs(O(:,ind)).^2;
figure
subplot(2,1,1), set(gca,'yscale','log'), hold on
plot(F,D,'-' ,'LineWidth',3,'DisplayName','dipole')
hold off
%xlim(TDDFT.outEnergyDFT.time([1 end]))
xlim([0 120]), xlabel('harmonic order')
ylim([1e-6 1]), ylabel('intensity (arb. units)')
title('total harmonic signal')
legend show
subplot(2,1,2), set(gca,'yscale','log'), hold on
plot(F,O(3,:),'-' ,'LineWidth',3,'DisplayName','HFO #3')
plot(F,O(2,:),'--','LineWidth',2,'DisplayName','HFO #2')
plot(F,O(1,:),':' ,'LineWidth',2,'DisplayName','HFO #1')
hold off
%xlim(TDDFT.outEnergyDFT.time([1 end]))
xlim([0 120]), xlabel('harmonic order')
ylim([1e-6 1]), ylabel('intensity (arb. units)')
title('orbital-resolved harmonic signal')
legend show
HF and orbital energies
Here we check the conservation of the HF energy in a spin-polarized molecular model where an electron is suddenly removed form a localized-core orbital and subjected to an external driving field in the dipole approximation and length gauge. First, we define the HF model and compute the (neutral) ground-state configuration
% Molecular model
A1 = QMol_Va_softCoulomb('atom','X','charge',4,'position',-4);
A2 = QMol_Va_softCoulomb('atom','A','charge',2,'position', 0);
A3 = QMol_Va_softCoulomb('atom','A','charge',2,'position', 3);
Vmol = QMol_DFT_Vext('atom',{A1,A2,A3});
% HF model
HF = QMol_DFT_spinPol( ...
'xspan', -50:.1:50, ... small prime factors
'occupation', {[1 1 1 1],[1 1 1 1]}, ...
'externalPotential', Vmol, ...
'HartreePotential', QMol_DFT_Vh_conv, ...
'exchangeCorrelationPotential', QMol_DFT_Vx_XX_conv);
% DFT ground state
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(HF);
Next, we remove the core orbital from the down-spin channel and define the external driving field
% Remove core orbital
HF.orbital.set('orbitalDown',HF.KSO.orbitalDown(:,2:end));
HF.set('occupation',{[1 1 1 1],[1 1 1]});
% External potential-vector field
E0 = convertLaser.intensity2amplitude(1e13);
w = convertLaser.wavelength2frequency(800e-9);
EField = @(t) E0*sin(w*t);
Finally we propagate the TDHF dynamics using the Forest-Ruth symplectic split-operator scheme
% TDHF propagator
TDHF = QMol_TDDFT_ESSO_4FR( ...
'time', 0:10:100, ...
'timeStep', 2e-2, ...
'electricField', EField, ...
'externalFieldGauge', 'length', ...
'saveExternalField', true, ...
'saveEnergyDFT', true, ...
'saveEnergyDFTtime', 1, ... save DFT energy every 1 a.u.
'absorbingBoundary', QMol_TDDFT_abs_CAP);
TDHF.propagate(HF);
With this, we can plot the variation of the total HF energy and its various components.
% Plot results
figure
subplot(2,1,1), hold on
plot(TDHF.outEnergyDFT.time,TDHF.outEnergyDFT.total-TDHF.outEnergyDFT.total(1),...
'-','LineWidth',2,'DisplayName','total');
hold off
xlim(TDHF.outEnergyDFT.time([1 end]))
xlabel('time (a.u.)')
ylabel('error (a.u.)')
title('total HF energy')
subplot(2,1,2), hold on
plot(TDHF.outEnergyDFT.time,sum(TDHF.outEnergyDFT.kinetic,1)-sum(TDHF.outEnergyDFT.kinetic(:,1)),...
'-','LineWidth',2,'DisplayName','kinetic')
plot(TDHF.outEnergyDFT.time,sum(TDHF.outEnergyDFT.external,1)-sum(TDHF.outEnergyDFT.external(:,1)),...
'-','LineWidth',2,'DisplayName','external')
plot(TDHF.outEnergyDFT.time,sum(TDHF.outEnergyDFT.externalField,1)-sum(TDHF.outEnergyDFT.externalField(:,1)),...
'-','LineWidth',2,'DisplayName','external field')
plot(TDHF.outEnergyDFT.time,TDHF.outEnergyDFT.Hartree-TDHF.outEnergyDFT.Hartree(1),...
'-','LineWidth',2,'DisplayName','Hartree')
plot(TDHF.outEnergyDFT.time,TDHF.outEnergyDFT.exchangeCorrelation-TDHF.outEnergyDFT.exchangeCorrelation(1),...
'-','LineWidth',2,'DisplayName','exchange correlation')
plot(TDHF.outEnergyDFT.time,TDHF.outEnergyDFT.autonomization-TDHF.outEnergyDFT.autonomization(1),...
'-','LineWidth',2,'DisplayName','autonomization')
hold off
xlim(TDHF.outEnergyDFT.time([1 end]))
xlabel('time (a.u.)')
ylabel('variation (a.u.)')
title('HF energy components')
legend show
The sudden change in the HF energy error is 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 spin-restricted molecular model driven by an external field in the dipole approximation and velocity gauge. First, we define the HF model and compute the (neutral) ground-state configuration
% Molecular model
A1 = QMol_Va_softCoulomb('atom','X','charge',4,'position',-4);
A2 = QMol_Va_softCoulomb('atom','A','charge',2,'position', 0);
A3 = QMol_Va_softCoulomb('atom','A','charge',2,'position', 3);
Vmol = QMol_DFT_Vext('atom',{A1,A2,A3});
% HF model
HF = QMol_DFT_spinRes( ...
'xspan', -50:.1:50, ... small prime factors
'occupation', [2 2 2 2], ...
'externalPotential', Vmol, ...
'HartreePotential', QMol_DFT_Vh_conv, ...
'exchangeCorrelationPotential', QMol_DFT_Vx_XX_conv);
% DFT ground state
SCF = QMol_DFT_SCF_Anderson('mixing',.4);
SCF.solveSCF(HF);
Next, we define the external potential-vector field and propagate the TDHF 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;
% TDDFT propagator
ABC = QMol_TDDFT_abs_CAP('length',40);
TDHF = QMol_TDDFT_ESSO_4FR( ...
'time', -T:2*pi/w:T, ...
'timeStep', 5e-2, ...
'potentialVector', AField, ...
'externalFieldGauge', 'velocity', ...
'saveExternalField', true, ...
'saveIonization', true, ...
'saveIonizationOrbitalIndex', 'all', ...
'saveIonizationTime', 2, ... save ionization yield every 2 a.u.
'absorbingBoundary', ABC);
TDHF.propagate(HF);
Finally, we plot the total and HF-orbital resolved ionization rates
% Plot results
T = .5*(TDHF.outIonization.time(2:end)+TDHF.outIonization.time(1:end-1));
figure
subplot(2,1,1), hold on
plot(TDHF.outIonization.time,TDHF.outIonization.potentialVector,...
'-','LineWidth',2,'DisplayName','A field');
hold off
xlim(TDHF.outIonization.time([1 end]))
xlabel('time (a.u.)')
ylabel('field (a.u.)')
title('potential vector')
subplot(2,1,2), hold on
plot(T,(TDHF.outIonization.total(2:end)-TDHF.outIonization.total(1:end-1))/(T(2)-T(1)),...
'-','LineWidth',2,'DisplayName','total')
plot(T,(TDHF.outIonization.orbital(1,2:end)-TDHF.outIonization.orbital(1,1:end-1))/(T(2)-T(1)),...
':','LineWidth',2,'DisplayName','HFO #1')
plot(T,(TDHF.outIonization.orbital(2,2:end)-TDHF.outIonization.orbital(2,1:end-1))/(T(2)-T(1)),...
':','LineWidth',2,'DisplayName','HFO #2')
plot(T,(TDHF.outIonization.orbital(3,2:end)-TDHF.outIonization.orbital(3,1:end-1))/(T(2)-T(1)),...
':','LineWidth',2,'DisplayName','HFO #3')
plot(T,(TDHF.outIonization.orbital(4,2:end)-TDHF.outIonization.orbital(4,1:end-1))/(T(2)-T(1)),...
':','LineWidth',2,'DisplayName','FHO #4')
hold off
xlim(TDHF.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 of the one-body density to store the spin density, down-sampled and over a reduced domain, to illustrate how this feature can be used as an alternative to storing the entire one-body density. We also use an installable output function of the HF orbitals to follow the current density in the TDHF dynamics. Here we consider a spin-polarized molecular model in which we introduce a localized ionization. First, we define the HF model and compute the (neutral) ground-state configuration
% Molecular model
A1 = QMol_Va_softCoulomb('atom','X','charge',2,'position',-3);
A2 = QMol_Va_softCoulomb('atom','A','charge',2,'position', 0);
A3 = QMol_Va_softCoulomb('atom','A','charge',2,'position', 3);
Vmol = QMol_DFT_Vext('atom',{A1,A2,A3});
% DFT model
HF = QMol_DFT_spinPol( ...
'xspan', -50:.1:50, ... small prime factors
'occupation', {[1 1 1],[1 1 1]}, ...
'externalPotential', Vmol, ...
'HartreePotential', QMol_DFT_Vh_conv, ...
'exchangeCorrelationPotential', QMol_DFT_Vx_XX_conv);
% DFT ground state
SCF = QMol_DFT_SCF_Anderson;
SCF.solveSCF(HF);
Next, we alter the electronic structure to emulate the localized ionization (by mixing the highest two orbitals in the down-spin channel), define the output function handles and propagate the TDHF dynamics using the Forest-Ruth symplectic split-operator scheme
% Remove core orbital
HF.orbital.set('orbitalDown',[HF.KSO.KSOdw(:,1) (HF.KSO.KSOdw(:,2)+HF.KSO.KSOdw(:,3))/sqrt(2)]);
HF.set('occupation',{[1 1 1],[1 1]});
% Installable output functions
ind = HF.disc.x >= -15 & HF.disc.x <= 15 & mod(1:numel(HF.disc.x),2) == 1; % 1/2 grid points between -15 and 15
funRho = @(rho,t) rho.densityUp(ind) - rho.densityDown(ind);
funKSO = @(KSO,t) real(-1i*conj(KSO.KSOup(:,1)).*ifft(HF.disc.D.*fft(KSO.KSOup(:,1)))) * HF.occ{1}(1) + ...
real(-1i*conj(KSO.KSOup(:,2)).*ifft(HF.disc.D.*fft(KSO.KSOup(:,2)))) * HF.occ{1}(2) + ...
real(-1i*conj(KSO.KSOup(:,3)).*ifft(HF.disc.D.*fft(KSO.KSOup(:,3)))) * HF.occ{1}(3) + ...
real(-1i*conj(KSO.KSOdw(:,1)).*ifft(HF.disc.D.*fft(KSO.KSOdw(:,1)))) * HF.occ{2}(1) + ...
real(-1i*conj(KSO.KSOdw(:,2)).*ifft(HF.disc.D.*fft(KSO.KSOdw(:,2)))) * HF.occ{2}(2);
% TDDFT propagator
ABC = QMol_TDDFT_abs_CAP('length',20);
TDHF = QMol_TDDFT_ESSO_4FR( ...
'time', 0:10:100, ...
'timeStep', 2e-2, ...
'saveOutputFunctionDensity', funRho, ...
'saveOutputFunctionDensityTime', 1, ... save spin density every 1 a.u.
'saveOutputFunctionOrbital', funKSO, ...
'saveOutputFunctionOrbitalTime', 1, ... save current every 2 a.u.
'absorbingBoundary', ABC);
TDHF.propagate(HF);
With this, we can plot the results of the installable output functions
% Plot results
figure
subplot(2,1,1), hold on
imagesc(TDHF.outOutputFunctionDensity.time,HF.x(ind),TDHF.outOutputFunctionDensity.result)
hold off
xlim(TDHF.outOutputFunctionDensity.time([1 end]))
ylim([-10 10])
xlabel('time (a.u.)')
ylabel('position (a.u.)')
title('spin density')
colorbar vert
subplot(2,1,2), hold on
imagesc(TDHF.outOutputFunctionOrbital.time,HF.x,TDHF.outOutputFunctionOrbital.result)
hold off
xlim(TDHF.outOutputFunctionOrbital.time([1 end]))
ylim([-10 10])
xlabel('time (a.u.)')
ylabel('position (a.u.)')
title('current density')
colorbar vert
References
[Mauger 2025] Extended symplectic split operator scheme: F. Mauger and C. Chandre, "Extended phase-space symplectic integration for electronic dynamics," in preparation (2025).
Notes
The results displayed in this documentation page were generated using version 01.23 of the QMol-grid package.