In [3]:
%%HTML
<link rel="stylesheet" type="text/css" href="custom.css">
In [2]:
%matplotlib inline
from matplotlib import pyplot as plt
import pickle
from ase import Atoms
from ase.cluster import Cluster, Icosahedron
from ase import units as un
import ase.io
from IPython.core.display import Image
import numpy as np
from numpy import array, sum, std, mean
import scipy.interpolate as inter
import gzip
Naszym celem jest zbadanie zachowania warstwowych nanocząstek Fe-Pt w wysokich temperaturach. Mamy nadzieję wyznaczyć granicę stabilności termicznej takich nanocząstek. Poniższa animacja pokazuje funkcję kalorymetryczną ($E_P - E_K$) nanocząstki jako funkcję temperatury, razem z odcinkami animacji przedstawiającymi wygląd nanocząstki w danej temperaturze.
Aby móc prowadzić rachunki musimy najpierw zdefiniować badaną strukturę. W naszym przypadku będzie to ikosaedryczna (dwudziestościenna) cząstka złożona z 55 atomów (1Pt+12Fe+42Pt). Szczęśliwie, biblioteka którą się posługujemy (ASE - Atomistic Simulation Environment) zawiera funkcję budującą dwudziestościenne cząstki.
In [3]:
ico=Icosahedron("Fe", 3); ico.center(vacuum=2.0)
ase.io.write('fe-ico.pov', ico, show_unit_cell=2, rotation='115y,15x', run_povray=True)
Image(filename='fe-ico.png')
Out[3]:
In [4]:
def LC_ico(a1,a2,ln):
'''
Stwórz strukturę warstwową rzędu ln
o naprzemiennych warstwach atomów a1, a2
'''
ico=Icosahedron(a1,ln)
icon=ico.get_atomic_numbers()
for l in range(ln,0,-1):
if (ln-l)%2 :
a=a1
else :
a=a2
il=Icosahedron(a,l)
an=il.get_atomic_numbers()
icon[:len(an)]=an
ico.set_atomic_numbers(icon)
return ico
Upewnijmy się że nasza procedura działa poprawnie rysując stworzoną komórkę w przekroju. W tym celu wykorzystamy możliwość wyboru elementów listy przy pomocy tablicy logicznej (fancy indexing). Taka konstrukcja jest specyficzna dla języka Python.
In [5]:
lc=LC_ico("Fe","Pt",3); lc.center(vacuum=2.0) # Warstwowa nanocząstka
cut=lc.get_center_of_mass()[2] # Płaszczyzna cięcia (CM)
ase.io.write('cebula.pov',
lc[[a.index for a in lc if a.position[2]>=cut]], # Selekcja atomów
show_unit_cell=2, rotation='135y,15x',
canvas_height=300,run_povray=True)
Image(filename='cebula.png')
Out[5]:
Równanie ruchu atomu o indeksie $i$:
$$ \frac{d^2{\mathbf x}_i}{{dt}^2}=\frac{{\mathbf F}_i}{m_i} $$rozwiązywane jest klasycznymi metodami numerycznego całkowania równan ruchu (np. algorytm Verleta), siły ${\mathbf F}_i$ pochodzą z rachunku ab-initio (siły Hellmanna-Feynmana).
Dla naszej struktury przeprowadzamy serię obliczeń w których nadajemy atomom różną średnią energię kinetyczną - czyli temperaturę:
$$ T=\frac{2\langle E_k \rangle}{3k_B}. $$W rezultacie otrzymujemy trajektorie układu, które poddajemy dalszej analizie.
Ja, żeby nie tracić czasu, przygotowałem już sobie wcześniej ... (Adam Słodowy)
Obliczenia ab-initio użyte tutaj trwały długie miesiące. Samo ich przygotowanie i przeprowadzenie jest zabiegiem technicznym i nie jest szczególnie kształcące. Potrzebne wyniki tych obliczeń - w postaci kolekcji trajektorii - znajdują się w kartotece data w materiałach do ćwiczeń.
In [6]:
# Wczytanie danych wyliczonych przez program VASP przygotowanych wczesniej
md={}
for k,tr in pickle.load(gzip.open('data/md_PtFePt.p.gz','rb')).items():
T=int(k.split('/')[-1][1:])
md[T]=tr
print(sorted(md.keys()))
W materiałach ćwiczeniowych znajduje się też kilka dalszych przykładów analizy uzyskanych trajektorii. Tutaj zajmiemy się uzyskaniem wykresu kalorymetrycznego ($E_p - E_k$ jako funkcja temperatury $T$), który pozwoli nam rozpoznać temperaturę topnienia nanocząstki o platynowej powierzchni.
In [7]:
def calc_velocities(a1,a2,dt):
'''
Wyliczenie prędkości atomów na podstawie pozycji a1, a2
oraz odstępu czasu dt
'''
dx=(a2.get_positions()-a1.get_positions())
c=Atoms(a2)
tv=sum(c.get_cell(),axis=0)/2
c.set_positions(dx)
c.translate(tv)
c.set_scaled_positions(c.get_scaled_positions())
c.translate(-tv)
return c.get_positions()/dt
In [8]:
def calc_kinen(tr,dt):
'''
Wyliczenie energii kinetycznej atomów w trajektorii tr
przy kroku czasowym dt.
'''
m=tr[0].get_masses()
Ek=[m*sum(calc_velocities(a1,a2,dt)**2,axis=1)/2
for a1,a2 in zip(tr[:-1],tr[1:])]
return Ek
def calc_temp(tr,dt):
'''
Wyliczenie chwilowej temperatury w trajektorii tr
przy kroku czasowym dt
'''
eka=calc_kinen(tr,dt)
ek=sum(eka,axis=-1)
ekd=std(eka,axis=-1)
n=len(tr[0].get_masses())
c=n/(n-1)
return 2*c*ek/(3*n*un.kB), 2*c*ekd/(3*un.kB)
In [9]:
plt.figure(figsize=(8,5))
for T in [300, 1000]:
dt,tr=md[T]; temp,_=calc_temp(tr,dt*un.fs)
avgT=mean(temp) ; devT=std(temp)
plt.plot(dt*np.arange(len(temp))/1000,temp,'-')
plt.axhline(avgT,ls='--',color='r',lw=2)
plt.axhspan(avgT-devT,avgT+devT,color='C7',alpha=0.2)
plt.xlabel('Czas (ps)'); plt.ylabel('Temperatura (K)');
In [10]:
dat=[]
for T in sorted(md):
dt,tr=md[T]
epot=array([s.get_potential_energy() for s in tr])
epot=((epot[:-1]+epot[1:])/2)
ekin=sum(array(calc_kinen(tr,dt*un.fs)),axis=-1)
etot=epot+sum(ekin,axis=-1)
temp_v, temp_d=calc_temp(tr,dt*un.fs)
atemp=mean(temp_v)
dat.append([atemp, mean(epot), mean(ekin),
mean((etot-mean(etot))**2)/(un.kB*(atemp**2))])
dat=array(dat).T
In [11]:
plt.figure(figsize=(8,5))
plt.plot(dat[0],dat[1]-dat[2],'o')
x=(dat[0])[np.arange(dat.shape[1])!=16];
y=(dat[1]-dat[2])[np.arange(dat.shape[1])!=16]
s=inter.UnivariateSpline(x, y, s=1.3)
xx=np.arange(100,2100,5)
plt.plot(xx,s(xx),'k--',lw=2,alpha=0.25);
plt.xlabel('Temperatura (K)'); plt.ylabel('$E_p - E_k$ (eV)');
plt.title('Krzywa kalorymetryczna Pt-12Fe-42Pt')
plt.xlim(100,2100); plt.ylim(-325,-310);
Wszystkie obliczenia prowadzące do pokazanych wyników, oraz wszystkie konieczne dane są dostępne w plikach źródłowych oraz materiałach ćwiczeniowych. Wszystkie materiały można znaleźć pod adresem: https://goo.gl/5JKOJK