Skip to content
Snippets Groups Projects
Commit ad0ff6bc authored by ziskaj00's avatar ziskaj00
Browse files

Update sheet6.py

parent ab481ff8
Branches
No related tags found
No related merge requests found
...@@ -5,6 +5,7 @@ import matplotlib.pyplot as plt ...@@ -5,6 +5,7 @@ import matplotlib.pyplot as plt
import numpy as np import numpy as np
from numpy.typing import ArrayLike from numpy.typing import ArrayLike
from scipy.constants import Boltzmann as k_B from scipy.constants import Boltzmann as k_B
from tqdm import tqdm
# Problem 6.1 # Problem 6.1
# a) # a)
...@@ -53,7 +54,7 @@ def verlet_thermostat( ...@@ -53,7 +54,7 @@ def verlet_thermostat(
f, U[0] = force(r[0], L) f, U[0] = force(r[0], L)
for i in range(1, N): for i in tqdm(range(1, N)):
# first half step # first half step
p[i] = p[i-1] + 1/2 * f * dt p[i] = p[i-1] + 1/2 * f * dt
r[i] = r[i-1] + p[i] / m * dt r[i] = r[i-1] + p[i] / m * dt
...@@ -94,7 +95,7 @@ def verlet( ...@@ -94,7 +95,7 @@ def verlet(
f, U[0] = force(r[0], L) f, U[0] = force(r[0], L)
for i in range(1, N): for i in tqdm(range(1, N)):
# first half step # first half step
p[i] = p[i-1] + 1/2 * f * dt p[i] = p[i-1] + 1/2 * f * dt
r[i] = r[i-1] + p[i] / m * dt r[i] = r[i-1] + p[i] / m * dt
...@@ -179,7 +180,7 @@ T = 1 / (2 * m0) * np.sum(p**2, axis=(1, 2)) ...@@ -179,7 +180,7 @@ T = 1 / (2 * m0) * np.sum(p**2, axis=(1, 2))
H = U + T H = U + T
dH = H - np.mean(H) dH = H - np.mean(H)
c_V = np.mean(dH**2) / (N * k_BT**2) c_V = np.mean(dH**2) / (N * k_BT**2)
dc_V = bp.std(dH**2) / (N * k_BT**2) dc_V = np.std(dH**2) / (N * k_BT**2)
print(f"c_V = {c_V} ± {dc_V} k_B") print(f"c_V = {c_V} ± {dc_V} k_B")
# c) # c)
...@@ -210,5 +211,5 @@ T = 1 / (2 * m0) * np.sum(p**2, axis=(1, 2)) ...@@ -210,5 +211,5 @@ T = 1 / (2 * m0) * np.sum(p**2, axis=(1, 2))
H = U + T H = U + T
dH = H - np.mean(H) dH = H - np.mean(H)
c_V = np.mean(dH**2) / (N * k_BT**2) c_V = np.mean(dH**2) / (N * k_BT**2)
dc_V = bp.std(dH**2) / (N * k_BT**2) dc_V = np.std(dH**2) / (N * k_BT**2)
print(f"c_V = {c_V} ± {dc_V} k_B") print(f"c_V = {c_V} ± {dc_V} k_B")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment