DensityMatrix method

dm.jl

#!/usr/bin/env julia

using LinearAlgebra;
using Printf;

import Random;
Random.seed!(1234);

N  = 1_000_000;
δt = 0.1;
ħ  = 0.658212;

# initial density matrix
ρ0 = diagm([1.0 + 0im, 1, 0]);
# ρ0[1, 2] = 1.0;
# ρ0[2, 1] = 1.0;

# hamiltonian initialize
# 1.0      0.325977 0.549051
# 0.325977 2.0      0.218587
# 0.549051 0.218587 3.0
H   = diagm([1.0, 2, 3]);
H[1,2] = rand();
H[1,3] = rand();
H[2,3] = rand();
H[2,1] = H[1,2];
H[3,1] = H[1,3];
H[3,2] = H[2,3];

# liouville equation
# dρ/dt = -i/ħ [H,ρ]
# U = exp(-iHt/ħ)
# ρ(t) = U(t) ρ(0) U'(t)

ρt = zeros(Float64, 3, N);
ρijt = zeros(Float64, 3, 3, N);

for i in 1:N
    t = i * δt;
    U = exp(-im*t/ħ*H);
    ρ = U * ρ0 * U';

    ρijt[:, :, i] = abs2.(ρ);
    ρt[:,i] =  abs2.(diag(ρ));
end

@printf "min(ρ11) = %9.6f , max(ρ11) = %9.6f\n" minimum(ρt[1,:]) maximum(ρt[1,:])
@printf "min(ρ22) = %9.6f , max(ρ22) = %9.6f\n" minimum(ρt[2,:]) maximum(ρt[2,:])
@printf "min(ρ33) = %9.6f , max(ρ33) = %9.6f\n" minimum(ρt[3,:]) maximum(ρt[3,:])

ρmax = zeros(Float64, 3, 3);
ρmin = zeros(Float64, 3, 3);
for j in 1:3
    for i in 1:3
        ρmax[i,j] = maximum(ρijt[i,j,:])
        ρmin[i,j] = minimum(ρijt[i,j,:])
    end
end

@info "" ρmax ρmin

wf.jl

#!/usr/bin/env julia

using LinearAlgebra;
using Printf;

import Random;
Random.seed!(1234);

N  = 1_000_000;
δt = 0.1;
ħ  = 0.658212;

# initial state
ψ0 = [1.0 + 0im, 1, 0];

# hamiltonian initialize
# 1.0      0.325977 0.549051
# 0.325977 2.0      0.218587
# 0.549051 0.218587 3.0
H   = diagm([1.0, 2, 3]);
H[1,2] = rand();
H[1,3] = rand();
H[2,3] = rand();
H[2,1] = H[1,2];
H[3,1] = H[1,3];
H[3,2] = H[2,3];

# TDSE
# ψ' = exp(-im Ht/ħ) ψ

ρt = zeros(Float64, 3, N);

for i in 1:N
    t = i * δt;
    ρt[:,i] = abs2.(exp(-im*H*t/ħ) * ψ0);
end

@printf "min(ρ11) = %9.6f , max(ρ11) = %9.6f\n" minimum(ρt[1,:]) maximum(ρt[1,:])
@printf "min(ρ22) = %9.6f , max(ρ22) = %9.6f\n" minimum(ρt[2,:]) maximum(ρt[2,:])
@printf "min(ρ33) = %9.6f , max(ρ33) = %9.6f\n" minimum(ρt[3,:]) maximum(ρt[3,:])