Transverse field Ising model
import numpy as np
def solve_tfim_mz(L, J, h):
# 1. Define Pauli matrices
sz = np.array([[1, 0], [0, -1]], dtype=complex)
sx = np.array([[0, 1], [1, 0]], dtype=complex)
id2 = np.eye(2)
def get_op(op, i, L):
"""Constructs an operator acting on site i in a chain of length L."""
res = 1
for j in range(L):
res = np.kron(res, op if i == j else id2)
return res
# 2. Construct Hamiltonian: H = J * sum(sz_i sz_{i+1}) + h * sum(sx_i)
dim = 2**L
H = np.zeros((dim, dim), dtype=complex)
# Interaction term
for i in range(L - 1):
H += J * (get_op(sz, i, L) @ get_op(sz, i+1, L))
# Transverse field term
for i in range(L):
H += h * get_op(sx, i, L)
# 3. Exact Diagonalization to find ground state
evals, evecs = np.linalg.eigh(H)
psi0 = evecs[:, 0]
# 4. Calculate sqrt(<Mz^2>) / L
Mz_op = np.zeros((dim, dim), dtype=complex)
for i in range(L):
Mz_op += get_op(sz, i, L)
Mz_sq_exp = np.vdot(psi0, (Mz_op @ Mz_op) @ psi0).real
return np.sqrt(Mz_sq_exp) / L
# Parameters
L, J, h = 10, -4, 1.0
result = solve_tfim_mz(L, J, h)
print(f"Numerical mz for L={L}: {result:.5f}")
Transverse field Ising model
program tfim_bitmask
implicit none
integer, parameter :: L = 12
integer, parameter :: dim = 2**L
real(8), parameter :: Jex = -4.0d0
real(8), parameter :: h = 1.0d0
complex(8), dimension(dim, dim) :: Ham
complex(8), dimension(dim) :: psi0
real(8), dimension(dim) :: w
complex(8), dimension(dim*2) :: work
real(8), dimension(dim*3) :: rwork
integer :: i, j, site, info, lwork
integer :: state, flipped_state, bit1, bit2
real(8) :: mz_sq_val, total_sz
Ham = (0.0d0, 0.0d0)
lwork = dim * 2
! 1. Construct Hamiltonian using Bit-Manipulation
do state = 0, dim - 1
! --- Transverse Field Term (h * sum sigma_x) ---
! sigma_x flips the bit at the specified site
do site = 0, L - 1
flipped_state = ieor(state, 2**site)
Ham(flipped_state + 1, state + 1) = Ham(flipped_state + 1, state + 1) + h
end do
! --- Ising Interaction Term (J * sum sigma_z_i * sigma_z_{i+1}) ---
do site = 0, L - 2
! btest returns .true. (Up) or .false. (Down)
if (btest(state, site) .eqv. btest(state, site + 1)) then
! Spins are same (Up-Up or Down-Down)
Ham(state + 1, state + 1) = Ham(state + 1, state + 1) + Jex
else
! Spins are different (Up-Down or Down-Up)
Ham(state + 1, state + 1) = Ham(state + 1, state + 1) - Jex
end if
end do
end do
! 2. Diagonalize using LAPACK (zheev)
print *, "Diagonalizing matrix of size:", dim
call zheev('V', 'U', dim, Ham, dim, w, work, lwork, rwork, info)
if (info /= 0) then
print *, "Error in diagonalization"
stop
end if
! 3. Calculate RMS Magnetization: sqrt(<(sum sz)^2>) / L
! The ground state vector is the first column of H after zheev
psi0 = Ham(:, 1)
mz_sq_val = 0.0d0
do state = 0, dim - 1
total_sz = 0.0d0
do site = 0, L - 1
if (btest(state, site)) then
total_sz = total_sz + 1.0d0 ! .true. is Up (+1)
else
total_sz = total_sz - 1.0d0 ! .false. is Down (-1)
end if
end do
mz_sq_val = mz_sq_val + (abs(psi0(state + 1))**2) * (total_sz**2)
end do
print '(A, I2, A, F10.6)', "Result for L=", L, " is mz = ", sqrt(mz_sq_val) / L
end program tfim_bitmask
Next project:
If a ray is reflected inside a rectangular box (whose inside is mirror) and if the ray falls on a side and get reflected, it will come back to the starting point after finite number of reflections if the slipe of y=mx+c ray , m=rational number. If m is irrational, then the ray will never come back to the original starting point and the whole box will light up i.e., it will fill the box inside.
I shall do it by using c/c++ or Python or Mathematica and show the result here in future.
Animation To show the path of Moon around the Sun. Its not helical as I thought on pen and paper. Its a simply connected single loop, not a perfect circle, not an ellipse, but a closed curve.
I think if irrational numbers are involved in periodicity, it will never come to the same spot again.
I am running all my python code in various methods because i am using Windows, Ubuntu and MacOS
In Ubuntu and MacOS its by creating a virtual environment
In Windows, I am doing it by IDLE.