%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=16)
class Trekant():
def __init__(self,points,col='cornflowerblue'):
self.points = points.copy() # .copy() as otherwise accessible from outside
self._color = col
self.triangle = None
self.cross = None
def __str__(self):
return 'Trekant {} {}'.format(self.points,self.color)
@property
def cm(self):
return np.mean(self.points,axis=0)
def draw(self,ax):
cm = self.cm
if self.triangle is None:
self.triangle = ax.fill(self.points[:,0],self.points[:,1],color=self.color)[0]
self.cross = ax.plot(cm[0],cm[1],marker='x',
markeredgecolor='k',
markersize=20,
markeredgewidth=3)[0]
else:
self.triangle.update({'xy': self.points})
self.cross.set_data(cm[0],cm[1])
def move(self, r):
self.points += r
self.draw()
@property
def color(self):
return self._color
@color.setter
def color(self,col):
self._color = col
if self.triangle is not None:
self.triangle.set(facecolor=self.color)
fig, ax = plt.subplots(figsize=(6,6))
ax.set_aspect('equal', adjustable='box')
ax.grid()
ax.set_xlim([0,8])
ax.set_ylim([0,5])
ax.set_xticks(range(9))
P1 = np.array([2,1])
P2 = np.array([4,3])
P3 = np.array([1,4])
points = np.array([P1,P2,P3])
t = Trekant(points)
t.draw(ax)
class Walker():
def __init__(self):
pass
def __str__(self):
return ##
def distance(self):
pass
def move(self):
theta = np.random.
delta_r = np.array
self.r +=
def plot_all_walkers():
pass
walkers = [Walker() for _ in range(Nwalkers)]
plot_all_walkers(ax,walkers)
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=16)
class Walker():
def __init__(self,col='c'):
self.r = np.array([0.,0.])
self.col = col
def __repr__(self): # ikke nødvendig, men giver pænt output i Jupyter Notebook
return 'Walker({})'.format(self.col)
def __str__(self):
return 'x={}, y={}, r={}'.format(self.r[0],self.r[1],self.distance())
def distance(self):
return np.linalg.norm(self.r)
def move(self):
theta = np.random.rand() * 2 * np.pi
delta_r = np.array([np.cos(theta),np.sin(theta)])
self.r += delta_r
def plot_all_walkers(ax,walkers_to_plot):
artist = ax.scatter([w.r[0] for w in walkers_to_plot],
[w.r[1] for w in walkers_to_plot],
color=[w.col for w in walkers_to_plot])
return artist
fig, axs = plt.subplots(1,2,figsize=(8,4))
for ax in axs:
ax.set_aspect('equal', adjustable='box')
ax.grid()
ax.set_xlim([-25,25])
ax.set_ylim([-25,25])
# create the walkers in a list for convenience
Nwalkers = 100
walkers = [Walker() for _ in range(Nwalkers)]
# let one of them have a different color
walkers[-1].col = 'orange'
plot_all_walkers(axs[0],walkers)
Ntimesteps = 100
average_distances = []
for _ in range(Ntimesteps):
for w in walkers:
w.move()
distances = [w.distance() for w in walkers]
average_distances.append(np.mean(distances))
plot_all_walkers(axs[1],walkers)
fig.tight_layout()
fig.savefig('uge4_walkers_1.png')
fig, axs = plt.subplots(1,2,figsize=(10,4))
for ax in axs:
ax.grid()
ts = list(range(Ntimesteps))
axs[0].plot(ts,average_distances)
axs[1].plot(np.log10(ts[1:]),np.log10(average_distances[1:]))
axs[0].set_xlabel('$t$')
axs[1].set_xlabel('$\log(t)$')
axs[0].set_ylabel('$d$')
axs[1].set_ylabel('$\log(d)$')
fig.tight_layout()
fig.savefig('uge4_walkers_2.png')
# start over for animation
fig, ax = plt.subplots(figsize=(4,4))
ax.set_aspect('equal', adjustable='box')
ax.grid()
ax.set_xlim([-25,25])
ax.set_ylim([-25,25])
# new walkers
walkers = [Walker() for _ in range(Nwalkers)]
# let one of them have a different color
walkers[-1].col = 'orange'
artist = plot_all_walkers(ax,walkers)
def update(i):
for w in walkers:
w.move()
artist.set_offsets([w.r for w in walkers])
return [artist]
from matplotlib import animation
plt.rc('animation', html='jshtml')
anim = animation.FuncAnimation(fig,
update,
frames=100,
interval=10,
blit=True)
anim
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font', size=16)
class WaterMolecule():
def __init__(self,r,theta=None):
self.r = r.copy()
if theta is None:
theta = np.random.rand() * 2 * np.pi
self.theta = theta
self.artists = None
self.ax = None
@property
def theta(self):
return self._theta
@theta.setter
def theta(self,new_theta):
self._theta = new_theta
self.calculate_rHs()
def __str__(self):
return 'x={}, y={}, theta={}'.format(self.r[0],self.r[1],self.theta)
def calculate_rHs(self):
x = self.r[0] + np.cos(self.theta)
y = self.r[1] + np.sin(self.theta)
self.r_H1 = np.array([x,y])
theta0 = 109 / 180 * np.pi
x = self.r[0] + np.cos(self.theta + theta0)
y = self.r[1] + np.sin(self.theta + theta0)
self.r_H2 = np.array([x,y])
def rattle(self):
self.theta += (np.random.rand() - 0.5) * 0.5
def draw(self,ax=None):
# create if first time with this Axis
if ax is not None and self.ax is not ax:
self.ax = ax
artist1 = ax.plot([],[],marker='o',
markerfacecolor='red',
markeredgecolor='k',
markersize=40)[0]
artist2 = ax.plot([],[],marker='o',
markerfacecolor='w',
markeredgecolor='k',
markersize=25)[0]
artist3 = ax.plot([],[],marker='o',
markerfacecolor='w',
markeredgecolor='k',
markersize=25)[0]
self.artists = [artist1,artist2,artist3]
self.artists[0].set_data(self.r[0],self.r[1])
self.artists[1].set_data(self.r_H1[0],self.r_H1[1])
self.artists[2].set_data(self.r_H2[0],self.r_H2[1])
fig, axs = plt.subplots(1,2,figsize=(12,6))
# create the molecules
Natoms_x, Natoms_y = 3,3
Natoms = Natoms_x * Natoms_y
lsize = 3.5
xs,ys = np.meshgrid(np.linspace(-lsize,lsize,Natoms_x),
np.linspace(-lsize,lsize,Natoms_y))
molecules = [WaterMolecule([x,y]) for x,y in zip(xs.flatten(),ys.flatten())]
# plot the molecules
for molecule in molecules:
molecule.draw(axs[0])
axs[0].set_aspect('equal', adjustable='box')
axs[0].grid()
axs[0].set_xlim([-5,5])
axs[0].set_ylim([-5,5])
from matplotlib import animation
plt.rc('animation', html='jshtml')
Uplot = axs[1].plot([],[])[0]
axs[1].set_xlim([0,100])
axs[1].set_ylim([-5.85,-5.6])
def calculate_total_potential_energy(molecules):
W_OO = 0
W_HO = -0.2
W_HH = 0.1
u = 0
for m1 in molecules:
r1_O = m1.r
r1_H1 = m1.r_H1
r1_H2 = m1.r_H2
for m2 in [m for m in molecules if m != m1]:
r2_O = m2.r
r2_H1 = m2.r_H1
r2_H2 = m2.r_H2
for r1,r2,W in [(r1_O,r2_O,W_OO),
(r1_O,r2_H1,W_HO),
(r1_O,r2_H2,W_HO),
(r1_H1,r2_O,W_HO),
(r1_H1,r2_H1,W_HH),
(r1_H1,r2_H2,W_HH),
(r1_H2,r2_O,W_HO),
(r1_H2,r2_H1,W_HH),
(r1_H2,r2_H2,W_HH),
]:
u += W / np.linalg.norm(np.array(r1)-np.array(r2))
return u
Us = []
def update(i):
if i < 50:
delta_U0 = 10
else:
delta_U0 = .000000001
U = calculate_total_potential_energy(molecules)
Us.append(U)
if 5 < i < 95:
for _ in range(2):
# try to move each atom once and plot again
for molecule in molecules:
theta_before = molecule.theta
molecule.rattle()
U_after = calculate_total_potential_energy(molecules)
delta_U = U_after - U
if delta_U < 0 or np.random.rand() < np.exp(-delta_U/delta_U0):
# accept
U = U_after
molecule.draw()
else:
# undo
molecule.theta = theta_before
Uplot.set_data(range(len(Us)),Us)
return np.array([m.artists for m in molecules]).flatten().tolist() + [Uplot]
anim = animation.FuncAnimation(fig,
update,
frames=100,
interval=10,
blit=True)
anim
Choose which booklet to go to: