In [1]:
%%HTML
<link rel="stylesheet" type="text/css" href="custom.css">
<style>
 .rise-enabled {background-color: #ffc;}
 .rise-enabled .text_cell {font-size: 120%;}
</style>
<script>
 MathJax.Hub.Config({displayAlign: 'left'});
</script>



In [2]:
# inicjalizacja
import ase
import ase.io
import numpy
from ase.visualize import view
from ase.build import bulk

Warsztaty modelowania w nanofizyce


Drgania struktury krystalicznej - fonony

Jan Łażewski

Zakład Komputerowych Badań Materiałów, Instytut Fizyki Jądrowej PAN, Kraków



In [3]:
# przygotowanie rysunku
t=ase.io.read('FIGs/MnAs_phM06.traj',index=':')
v=view(t,viewer='nglview')

v.custom_colors({'Mn':'green','As':'blue'})
v.view._remote_call("setSize", target="Widget", args=["500px", "500px"])
#v.view.center_view()
v.view.background='#ffc'
v.view.parameters=dict(clipDist=-200)

# wyświetlenie rysunku
v


The installed widget Javascript is the wrong version.

Dygresja na marginesie


Wykład został przygotowany w postaci tzw. żywej prezentacji w środowisku IPython Notebook



szczególne podziękowania dla prof. Pawła Jochyma za zachętę i okazaną pomoc



In [4]:
# przygotowanie rysunku
t=ase.io.read('FIGs/MnAs_softM48_20.traj',index=':')
v=view(t,viewer='nglview')

v.custom_colors({'Mn':'green','As':'blue'})
v.view._remote_call("setSize", target="Widget", args=["500px", "500px"])
#v.view.center_view()
v.view.background='#ffc'
v.view.parameters=dict(clipDist=-200)

# wyświetlenie rysunku
v


The installed widget Javascript is the wrong version.

Opis drgań kryształu w przybliżeniu harmonicznym - metoda bezpośrednia


K. Parlinski, Z.Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997).

K. Parlinski, PHONON software (Cracow, Poland, 1997-2017),
http://wolf.ifj.edu.pl/phonon/, http://www.computingformaterials.com/

Rozwinięcie energii potencjalnej w przybliżeniu harmonicznym

$$ \newcommand{\mb}{\boldsymbol} \newcommand{\h}{\hspace{-0.05cm}} $$

W przybliżeniu harmonicznym energię potencjalną sieci krystalicznej $E\bigl(\bigl\{ {\mb R}({\mb n},\mu)\bigr\}\bigr)$ można rozwinąć w szereg względem wychyleń ${\mb U}({\mb n},\mu)$ atomów
z ich położeń równowagi:

$$ \begin{eqnarray} E\bigl(\{ {\mb R}\}\bigr)\; =\; E_0\;\; + & & \h \sum_{{\mb n},\mu}\frac{\partial E\bigl(\bigl\{{\mb R}({\mb n},\mu)\bigr\}\bigr) } {\partial {\mb U}({\mb n},\mu)}\, {\mb U}({\mb n},\mu) \nonumber \\ + & & \h \frac{1}{2}\sum_{\substack{{\mb n},\mu \\ {\mb m},\nu}} \frac{\partial^2E\bigl(\bigl\{ {\mb R}({\mb n},\mu)\bigr\}\bigr) } {\partial {\mb U}({\mb n},\mu)\,\partial {\mb U}({\mb m},\nu)} \,{\mb U}({\mb n},\mu)\,{\mb U}({\mb m},\nu)\, \text{,} %\label{eq:harmE} \end{eqnarray} $$

gdzie ${\mb R}({\mb n},\mu)$ oznacza pozycję atomu $\mu$ w komórce prymitywnej ${\mb n}$, a sumowanie przebiega po całym krysztale.

Macierz stałych siłowych


W rozwinięciu $\,E\bigl(\bigl\{ {\mb R}({\mb n},\mu)\bigr\}\bigr)$ pierwsze pochodne znikają z warunków minimum energii w równowadze, zaś drugie pochodne definiują macierz stałych siłowych:

$$ \begin{eqnarray} \Phi_{ij}({\mb n},\mu;{\mb m},\nu ) = \frac{\partial^2 E\bigl(\bigl\{ {\mb R}({\mb n},\mu)\bigr\}\bigr) } {\partial U_i({\mb n},\mu)\,\partial U_j({\mb m},\nu)}\, \text{,} \end{eqnarray} $$

reprezentującą oddziaływanie między atomami w krysztale.


Siły Hellmanna-Feynmana (HF)


Każde wychylenie ${\mb U}({\mb n},\mu)$, nawet pojedynczego atomu, powoduje powstanie sił HF na wszystkich atomach w układzie:

\begin{eqnarray} {\mb F}({\mb n},\mu) = - \frac{\partial E}{\partial {\mb R}({\mb n},\mu)}\, . \end{eqnarray}

In [5]:
import ase.io
from ase.data import covalent_radii 
from ase.data.colors import cpk_colors
from ipywidgets import IntSlider, FloatSlider, interactive, fixed, HBox
from glob import glob


#scene.caption= """Right button drag or Ctrl-drag to rotate "camera" to view scene.
#To zoom, drag with middle button or Alt/Option depressed, or use scroll wheel.
#  On a two-button mouse, middle is left + right.
#Touch screen: pinch/extend to zoom, swipe or two-finger rotate."""

#clr={'Mn':color.green, 'As':color.blue}
clr={25:(0,0.51,0), 33:(0,0,0.95)}


def broom():
    for o in scene.objects:
        o.visible=False
        
class crystal:
    def __init__( self, trj, fscale=1, rscale=1 ):
        self.trj=trj
        self.rscale=rscale
        self.fscale=10**fscale
        self.n=0
        scene.fov=0.05
        scene.background=vec(1,1,0.8)
        scene.center=vec(*tuple(trj[0].get_center_of_mass()))
        self.atoms=[sphere(pos=vec(*tuple(r)), 
                           radius=self.rscale*covalent_radii[Z], 
                           color=vec(*clr[Z])) 
                     for Z, r in zip(trj[0].get_atomic_numbers(), trj[0].get_positions())]
        
        self.forces=[arrow(pos=vec(*tuple(r)), 
                           axis=vec(*tuple(f)), 
                           color=color.red) 
                     for r, f in zip(self.trj[0].get_positions(), self.fscale*trj[0].get_forces())]
   
    def update_frame( self, n=None, fscale=None, rscale=None ):
        if n is not None:
            self.n=n
        if fscale is not None:
            self.fscale=10**fscale
        if rscale is not None:
            self.rscale=rscale    
        for a, fv, r, f, Z in zip(self.atoms, self.forces, 
                               self.trj[self.n].get_positions(), 
                               self.fscale*self.trj[self.n].get_forces(),
                               self.trj[self.n].get_atomic_numbers()):
            a.pos=vec(*tuple(r))
            a.radius=self.rscale*covalent_radii[Z]
            fv.pos=vec(*tuple(r))
            fv.axis=vec(*tuple(f))
       
    def show(self, maxr=5, maxf=2.5):
        w1=interactive(self.update_frame, 
                       n=IntSlider(min=0,max=len(self.trj)-1,step=1,value=0,
                                   description='n:'), 
                       fscale=fixed(None), 
                       rscale=fixed(None));
        w2=interactive(self.update_frame, 
                       n=fixed(None), 
                       fscale=FloatSlider(min=1,max=maxf,value=1.0,
                                          description='siła:'), 
                       rscale=fixed(None));
        w3=interactive(self.update_frame, n=fixed(None), 
                       fscale=fixed(None), 
                       rscale=FloatSlider(min=0.1,max=maxr,value=1.0,
                                          description='promień:'));
        return HBox([w1, w2, w3])

In [6]:
from vpython import scene, vec, sphere, arrow, color



In [7]:
broom()
            
#trj=ase.io.Trajectory('data/md_T_1000.traj')
#trj=ase.io.read('FIGs/OUTCAR1', index=':')
trj=[ase.io.read(fn) for fn in sorted(glob('FIGs/OUTCAR_??'))]

c=crystal(trj)
c.show(1,2.5)

Zależność sił HF od stałych siłowych


W przybliżeniu harmonicznym z rozwinięcia energii otrzymujemy zależność:

$$ \begin{eqnarray} F_i({\mb n},\mu) = -\sum_{{\mb m},\nu,j} \Phi_{ij}({\mb n},\mu;{\mb m},\nu ) \,U_j({\mb m},\nu) . \end{eqnarray} $$

Macierz dynamiczna


Macierz dynamiczną definiujemy jako transformatę Fouriera macierzy stałych siłowych:

$$ \begin{eqnarray} {\mb D}({\mb k};\mu,\nu) = \frac{1}{\sqrt{M_\mu M_\nu}} \sum_{\mb m} {\mb \Phi}(0,\mu;{\mb m},\nu)\; \text{e}^{-2\pi i{\mb k}\cdot[ {\mb R}(0,\mu)-{\mb R}({\mb m},\nu)]} \text{,} \end{eqnarray} $$

gdzie ${\mb k}$ jest wektorem falowym, $M_\mu$ i $M_\nu$ masami atomów $\mu$ i $\nu$, zaś sumowanie przebiega po wszystkich atomach kryształu.

Równanie własne


Rozwiązanie równania własnego:

$$ \begin{eqnarray} {\mb D}({\mb k})\,{\mb e}({\mb k},j) = \omega^2({\mb k},j)\,{\mb e}({\mb k},j) \end{eqnarray} $$

poprzez diagonalizację macierzy dynamicznej daje częstości drgań $\,\omega({\mb k},j)$
i wektory polaryzacji $\,{\mb e}({\mb k},j)$.

Modelowane układy


Modelowane układy


Modelowane układy


Modelowane układy


Modelowane układy


Modelowane układy


Modelowane układy


Modelowane układy


Modelowane układy


Wyliczane własności dynamiczne


Metodą bezpośrednią można wyliczyć między innymi następujące własności materiałowe:

  • fononowe relacje dyspersji
  • intensywności pików fononowych
  • parametry Grüneissena
  • wektory polaryzacji + animacja
  • całkowite i cząstkowe widma DOS
  • funkcje termodynamiczne: E, S, F, G, c$_\text{v}$
  • czynniki Debye'a-Waller'a
  • rozszerzalność cieplna
  • rozpraszanie neutronów i X
  • diagramy fazowe, granice stabilności struktur

(file:///home/lazewski/pubs_txt/others/orals/warsztaty2017/warsztatyJL17.pdf)


In [8]:
!jupyter nbconvert FononyDFT.ipynb --to slides


[NbConvertApp] Converting notebook FononyDFT.ipynb to slides
[NbConvertApp] Writing 276174 bytes to FononyDFT.slides.html