Quantum mechanics in the browser

[timqian]

使用 javascript 写一个求解 Kohn–Sham equations 的库

[timqian]

step1:

convert this piece of code to javascript. 用 Davidson 算法 求 H 原子的几个本征值和本征向量，与文献做对比。

Things to learn

• 求解矩阵本征值的方法(wikipedia)
• Davidson method 求最小的 n 个本征向量
• 求解高维离散空间中哈密顿量的方法(wikipedia)

[timqian]

repo: https://github.com/timqian/DFT.js

[timqian]

matlab code solving atom in a box (1 dimension)

g = 100;
X = linspace(-1, 1, g);         % one dimensiton space lattice
h = X(2) - X(1);                % latice spacing
e = ones(g,1);               
L = spdiags([e -2*e e], -1:1, g, g) / h^2; % 1D finite difference Laplacian
H = -0.5 * L;                             % Hamiltonian of the atom
[PSIs,Es] = eigs(H);                       % eigenvalues of the atom

level = 2;         % energy level
disp(['Energy for the atom' num2str(Es(level,level)*27.21, 5) 'eV']); %display result
scatter(X, PSIs(:,level).^2);

%scatter([1,2,3,4,5,6],[-Es(1,1),-Es(2,2),-Es(3,3),-Es(4,4),-Es(5,5),-Es(6,6)])

http://octave-online.net/

[timqian]

matlab code solving H atom

g = 30; g3 = g^3;
p = linspace(-4, 4, g);         % one dimensiton space lattice
[X, Y, Z] = meshgrid(p, p, p);  % three dimension space lattice
h = p(2) - p(1);                % latice spacing
X = X(:); Y = Y(:); Z = Z(:);   % all elements of arraty as a single column
R = sqrt(X.^2 + Y.^2 + Z.^2);   % distance from the center
Vext = -1 ./ R;                 % potential energy
e = ones(g,1);               
L = spdiags([e -2*e e], -1:1, g, g) / h^2; % 1D finite difference Laplacian
I = speye(g);
L3 = kron(kron(L,I), I) + kron(kron(I, L), I) + kron(kron(I, I), L);  % extend Laplacian to 3 D
H = -0.5 * L3 + spdiags(Vext, 0, g3, g3);  % Hamiltonian of H atom
[PSIs,Es] = eigs(H);                      % Smallest eigenvalue of H
degree = 5
disp(['Total energy for Hydrogen atom' num2str(Es(degree,degree)*27.21, 5) 'eV']); %display result
%PSI_3 = reshape(PSI, [g,g,g]);
scatter3(X,Y,Z,PSIs(:,degree).^2 *13000);