在这种情况下势只是一个乘法算符即可以直接指数化的对角矩阵。演化的步骤是知道t时刻的波函数将波函数傅立叶变换到动量空间乘以动量算符然后逆变换到坐标空间乘以势能项接着变换到动量空间乘以动量算符最后变换回来即可得到t时刻的波函数 。clear all,clc; % 空间网格 % 空间总长度2457.6a.u.,步长0.15a.u. xmin -1228.8; xmax 1228.8; Nx 16384; dx (xmax - xmin) / Nx; x linspace(xmin, xmax, Nx); % 时间网格 % 总时间1600a.u.,步长0.1a.u. T 800; % 总时间 (a.u.) dt 0.1; % 时间步长 (a.u.) Nt round(T / dt); % 时间步数 t (0:Nt-1) * dt; % 动量空间网格 % 为了分裂算符法中动量算符直接与波函数作用对动量空间作fftshift dp 2*pi / (xmax - xmin); p (-Nx/2:Nx/2-1) * dp; p fftshift(p); % 调整为自然序 % 原子参数软库仑势 a 1; % 原子序数 b 1; % 软化参数 Vatom -a./ sqrt(x.^2 b); % 激光脉冲 % E E0·A(t)·cos(wtø),取峰值强度1e14W/cm2,中心波长800nm包络为高斯型 a0 5.29177210903e-11; % 玻尔半径 (m) epsilon 8.854187817e-12; % 真空介电常数(F/m) Eh 4.3597447222071e-18; % 哈特里能量 (J) c 299792458; % 光速(m/s) hbar 1.054571817e-34; % 约化普朗克常数 (J·s) au_time hbar/Eh; % 原子时间单位 (s) % 激光参数 peak_intensity 1e14; % 峰值强度 (W/cm²) FWHM 10e-15; % FWHM wavelength 800e-9; FWHM FWHM / au_time; lambda wavelength / a0; c0 c/2.1876912633/1e6; w0 2*pi*c0 / lambda; % 中心频率(a.u.) E0 sqrt(2*peak_intensity*1e4/epsilon/c)/5.1422e11; %电场振幅(a.u.) tau FWHM / (2*sqrt(2*log(2))); E E0 * exp(-(t-T/2).^2 / (2 * tau^2)*4*log(2)) .* cos(w0 * (t-T/2)); %电场形式 % 吸收边界参数 % 左右各50a.u.范围内快速衰减 x0 50; absorb ones(Nx, 1); idx_left find(x x(1) x0, 1, last); idx_right find(x x(end) - x0, 1, first); absorb(1:idx_left) exp(-(x(1) x0 - x(1:idx_left)).^2 / (2*x0^2)); absorb(idx_right:end) exp(-(x(idx_right:end) - (x(end) - x0)).^2 / (2*x0^2)); %% 虚时间演化求基态波函数 max_steps 400; % 演化时间 dt 0.1; sigma 10; % 初始化波函数 psi exp(-x.^2/(2*sigma^2)); norm_x sqrt(sum(abs(psi).^2*dx)); psi (psi/norm_x); % 归一化 % 虚时间演化 for i 1:max_steps % 动能演化 psi_p fft(psi); psi_p psi_p .* exp(-dt*p.^2/4); psi ifft(psi_p); % 势能演化 psi psi .* exp(-dt*Vatom); % 动能演化 psi_p fft(psi); psi_p psi_p .* exp(-dt*p.^2/4); psi ifft(psi_p); % 归一化波函数 norm_psi sqrt(sum(abs(psi).^2)*dx); psi psi / norm_psi; end %% 计算基态能量 % 取实部并归一化 psi real(psi); norm_psi sqrt(sum(psi.^2)*dx); psi psi / norm_psi; % 势能期望值 V_expect sum(psi.^2 .* Vatom) * dx; % 动能期望值 psi_p fft(psi); d2_psi ifft(-p.^2 .* psi_p); T_expect -0.5 * real(sum(conj(psi).*d2_psi)) * dx; % 总能量 E_ground T_expect V_expect; fprintf(基态能量 %.4f\n, E_ground); %% % 分裂算符时间演化 % 预计算偶极矩 dipole zeros(Nt, 1); for n 1:Nt V_total Vatom - x * E(n); % 总势能 % 画图势场变化和波函数实部变化 % subplot(2,1,1) % plot(x,V_total); % xlim([-100,100]);ylim([-3,1]) % pause(0.01); % subplot(2,1,2) % plot(x,real(psi)); % xlim([-10,10]);ylim([-2,1]); pause(0.01); % 时间演化 psi_p fft(psi); psi_p psi_p .* exp(-1i * dt/2 * p.^2 / 2); psi ifft(psi_p); psi psi .* exp(-1i * dt/2 * V_total); psi_p fft(psi); psi_p psi_p .* exp(-1i * dt/2 * p.^2 / 2); psi ifft(psi_p); % 吸收边界 psi psi .* absorb; % 偶极矩 x ∫ψ* x ψ dx dipole(n) trapz(x, conj(psi) .* x .* psi); end % 加窗函数减少频谱泄漏 window hann(Nt); dipole dipole .* window; % 傅里叶变换 N_fft 2^nextpow2(2*Nt); dw fft(dipole, N_fft); w 2*pi*(0:N_fft-1) / (N_fft * dt); % 中心角频率 (a.u.) % 转换为谐波级次 harmonic_order w / w0; % 结果可视化 % 激光电场 figure(Name, Laser Field); plot(t, E, b, LineWidth, 1.5); xlabel(Time (a.u.)); ylabel(Electric Field (a.u.)); title(Laser Field); grid on; % 偶极矩 figure(Name, Dipole); plot(t, real(dipole), r, LineWidth, 1.5); xlabel(Time (a.u.)); ylabel(Dipole (a.u.)); title(Dipole); grid on; % 高次谐波谱 figure(Name, HHG Spectrum); max_harmonic 100; % 显示的最大谐波级次 final_order find(harmonic_order max_harmonic, 1, last); semilogy(harmonic_order(1:final_order), abs(dw(1:final_order)).^2, k, LineWidth, 1.5); xlabel(Harmonic Order (w/w0)); ylabel(Intensity (a.u.)); title(High-Harmonic Generation Spectrum(I1e15W/cm2)); xlim([0 max_harmonic]); grid on;