next_inactive up previous


OpenPhonon: an open source computer code for lattice-dynamical calculations

Alessandro Mirone 1 and Matteo d'Astuto

March 21, 2006


Contents


Introduction

This text is a users manual meant to explain the use of the lattice dynamicals calculation package OpenPhonon. It is by no means an introduction to lattice dynamics, which is well described on references[1,2,3]. Nevertheless, a brief introduction is given in the next section in order to state exactly what kind of calculations the present program can perform. A primary goal of this package is to offer a completely open source code, allowing all the people interested in lattice dynamical calculations to use it and adapt it to the needs of their research. In particular this program is born in order to prepare and interpret inelastic x-ray scattering experiment (IXS), therefore its main new feature is the calculation of IXS intensities from the dynamical structure factor obtained from lattice dynamical calculations. This has been implemented using the DABAX data base for atomic form factors compiled at the ESRF. But the programs could be adapted in order to calculate the response function for other kinds of experiments as inelastic neutron and Raman scattering, nuclear quadrupole resonance etc...

The code is written in Python using its numerialc library NumPy. In order to use the program just a very basic knowledge of the language is required. One-half-a-day tutorial is available at www.python.org.


Basic concept of lattice dynamics

Lattice Hamiltonian is developed to second order around the nominal equilibrium position, given by the experimental crystal structure. This corresponds to a simple implementation of the quasi-harmonic approximation (the positions of the atom at rest are choosen at the minimum of the free energy, in thermodynamic equilibrium) instead of the harmonic approximation (the position of the atom at rest are choose at the minimum of the potential energy, in mechanical equilibrium)

The Hamiltonian is given by the sum of the kinetic part plus the sum of interatomic two-body central potentials. The Hamiltonian degrees of freedom are the atomic cores and atomic shells coordinates. The atomic shells have zero mass and does not contribute to the kinetic part.

The potentials implemented up to now in OpenPhonon are:

1)
The Born-von-Karman potential: $V_{BK}$, which describes an empirical central potential by its second derivative (longitudinal and transverse components) for each interacting couple of atoms. One describes a Born von Karman potential by the force constants $F_L$ (longitudinal) and $F_T$ (transverse) which are given by the formula:


\begin{displaymath} F_L=\frac{d^2V_{BK}(\kappa,\ell,\kappa ' \ell ')}{dr^2} \vert _{r=(r_{\kappa' \ell'}^\circ - r_{\kappa \ell}^\circ)} \end{displaymath} (1)

and


\begin{displaymath} F_T=\frac{1}{r} \frac{dV_{BK}(\kappa,\ell,\kappa ' \ell ')}{dr} \vert _{r=(r_{\kappa' \ell'}^\circ - r_{\kappa \ell}^\circ)} \end{displaymath} (2)

2)
The Coulomb potential:


\begin{displaymath} U_c(\vert\vec{R}_{\kappa,\ell}-\vec{R}_{\kappa ',\ell '}\vert) Z_{\kappa} Z_{\kappa '} \end{displaymath} (3)

where $ U_c=\frac{1}{r} $ and $Z_{\kappa }$ is the charge of the $\kappa$ ion.

3)
The Born-Mayer potential:


\begin{displaymath} V_{BM}= V_{\circ}~\exp\biggl(-\frac{r}{R^\circ_{\kappa}+R^\circ_{\kappa'}}\biggr) \end{displaymath} (4)

4)
The Van der Waals potential:


\begin{displaymath} V_{VW}=-V_{\circ}~\biggl(\frac{R^\circ_{\kappa}+R^\circ_{\kappa'}}{r}\biggr)^6 \end{displaymath} (5)

5)
The Lennard-Jones potential:


\begin{displaymath} V_{VW}=V_{\circ}~\biggl[ \biggl(\frac{R^\circ_{\kappa}+R^\ci... ...(\frac{R^\circ_{\kappa}+R^\circ_{\kappa'}}{r}\biggr)^6 \biggr] \end{displaymath} (6)

6)
Other more complex interaction which have been implemented more recently : Tensorial forces ( dipolar fluctuation mode, Jan-Teller, non central forces ), and angle bonds.


Finally, in condensed matter (in particular for metals) the Coulomb potential is renormalized by the screening effect of the surrounding charges. This effect can be calculated in many ways, we have implemented the simplest approach: the Thomas-Fermi one, i.e. replacing the Coulomb potential by the Yukawa potential:


\begin{displaymath} U_{TF}=\frac{\exp\bigl({\frac{-r}{R_{TF}}}\bigr)}{r} \end{displaymath} (7)

.

The motion of the atom labelled $\kappa$ in the $\ell$th cell is described by the following equations:

$\displaystyle m_{\kappa} \frac{d^2 u_{\kappa\ell}}{dt^2} =$ $\textstyle -X_{\kappa} \sum_{(\kappa,\ell)\neq(\kappa',\ell')} V_{(\kappa \ell ... ... X_{\kappa ' } u_{\kappa ' \ell ' } + Y_{\kappa ' } w_{\kappa ' \ell ' } \bigl)$   (8)
  $\textstyle + X_{\kappa} \sum_{(\kappa,\ell)\neq(\kappa',\ell')} V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{C} Z_{\kappa ' } u_{\kappa \ell }$    
  $\textstyle - \sum_{(\kappa,\ell)\neq(\kappa',\ell')}\Bigl[ V_{(\kappa \ell ~,~\... ...~\kappa ' \ell ' )}^{T} w_{\kappa ' \ell ' }\Bigl] + K_{\kappa} w_{\kappa \ell}$    
  $\textstyle + \Bigl[ \sum_{(\kappa,\ell)\neq(\kappa',\ell')} \bigl( V_{(\kappa \... ...kappa \ell ~,~\kappa ' \ell ' )}^{T} \bigl) - K_{\kappa} \Bigl] u_{\kappa \ell}$    


$\displaystyle \mu_{\kappa} \frac{d^2 w_{\kappa\ell}}{dt^2} =$ $\textstyle -Y_{\kappa} \sum_{(\kappa,\ell)\neq(\kappa',\ell')}{ V_{(\kappa \ell... ... X_{\kappa ' } u_{\kappa ' \ell ' } + Y_{\kappa ' } w_{\kappa ' \ell ' } \bigl)$   (9)
  $\textstyle + Y_{\kappa} \sum_{(\kappa,\ell)\neq(\kappa',\ell')} V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{C} Z_{\kappa ' } w_{\kappa \ell }$    
  $\textstyle - \sum_{(\kappa,\ell)\neq(\kappa',\ell')} \Bigl[ V_{(\kappa \ell ~,~... ...\kappa ' \ell ' )}^{T'}\Bigl] u_{\kappa ' \ell ' } + K_{\kappa} u_{\kappa \ell}$    
  $\textstyle + \Bigl[ \sum_{(\kappa,\ell)\neq(\kappa',\ell')} \bigl( V_{(\kappa \... ...appa \ell ~,~\kappa ' \ell ' )}^{T'} \bigl) - K_{\kappa} \Bigl] w_{\kappa \ell}$    

where equation (8) describes the motion of the centre of mass of each atom while equation (9) refers to the motion of the shell associated to each atom, and its mass $\mu_{\kappa} = 0$. The parameter $V^C$ is given by


\begin{displaymath} V^C=-\frac{\partial^2 U_c(\vert\vec{r}_1-\vec{r}_2\vert)}{\partial\vec{r}_1\partial\vec{r}_2} \end{displaymath} (10)

$K$ are the spring constants which tie the shells to the cores. Symbols $V^{SR}$,$V^{T}$,$V^{T'}$,$V^{S}$ denote non-coulombian interaction. The $V^{SR}$ represent interactions between the cores. Although coulombian interactions act between the cores too, they need a special treatment because of their long range, and have thus their own writing $V^C$. $V^{T}$ are the elastic interactions between a shell and a core on another site, $V^{T'}$ between a core and another site shell and $V^{S}$ are the interactions between two different shells. These terms are the sum of analytical potentials ( Born-Mayer or Van de Waals ) and Born-von-Karman potentials. In the case of analytical potentials the $V$'s are derived from the analogous of equation (10) while in the Born-von-Karman case they are derived from the longitudinal and transverse second derivatives $F_L$ and $F_T$, namely :

\begin{displaymath} V_{(\kappa \ell ~,~\kappa ' \ell ' )}= - \hat{r} \otimes \hat{r} (F_L-F_T) - F_T \end{displaymath} (11)

$\hat{r}$ being the unit vector going from $( ~\kappa ' \ell ' )$ to $(\kappa \ell ~ )$

One can cast together $X_{\kappa \ell}$ (core charge) and $Y_{\kappa \ell}$ (shell charge) in equations (8) and (9) to get $Z_{\kappa \ell}$ (ion charge=core+shell) and $Y_{\kappa \ell}$ by using a different set of displacement variables ( $u_{\kappa \ell}$ and $\bar w_{\kappa \ell}$ ) where $\bar w_{\kappa \ell} = w_{\kappa \ell} - u_{\kappa \ell}$. By performing such substitution, by summing together equations (8) and (9) and considering the limit $\mu_{\kappa} \rightarrow 0$, one obtains :


$\displaystyle m_{\kappa} \frac{d^2 u_{\kappa\ell}}{dt^2} = -Z_{\kappa} \sum_{(\... ...kappa ' } u_{\kappa ' \ell ' } + Y_{\kappa ' } \bar w_{\kappa ' \ell ' } \bigl)$      
$\displaystyle + Z_{\kappa} \sum_{(\kappa,\ell)\neq(\kappa',\ell')} V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{C} Z_{\kappa ' } u_{\kappa \ell }$      
$\displaystyle + Y_{\kappa} \sum_{(\kappa,\ell)\neq(\kappa',\ell')} V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{C} Z_{\kappa ' } \bar w_{\kappa \ell }$     (12)
$\displaystyle - \sum_{(\kappa,\ell)\neq(\kappa',\ell')}\Bigl[ V_{(\kappa \ell ~... ...ll ~,~\kappa ' \ell ' )}^{T'} \Bigl] (u_{\kappa ' \ell ' } - u_{\kappa \ell } )$      
$\displaystyle - \sum_{(\kappa,\ell)\neq(\kappa',\ell')}\Bigl[ V_{(\kappa \ell ~... ...S} + V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{T} \Bigl] \bar w_{\kappa ' \ell ' }$      
$\displaystyle + \sum_{(\kappa,\ell)\neq(\kappa',\ell')}\Bigl[ V_{(\kappa \ell ~... ...}^{S} + V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{T'} \Bigl] \bar w_{\kappa \ell }$      

and for the zero-mass shells

$\displaystyle 0= -Y_\kappa \sum_{(\kappa,\ell)\neq(\kappa',\ell')} V_{(\kappa \... ...~\kappa ' \ell ' )}^{C} Z_{\kappa ' } ( u_{\kappa ' \ell ' } - u_{\kappa \ell})$      
$\displaystyle -Y_\kappa \sum_{(\kappa,\ell)\neq(\kappa',\ell')} V_{(\kappa \ell... ...} V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{C} Z_{\kappa ' } \bar w_{\kappa \ell }$      
$\displaystyle -K_x \bar w_{\kappa \ell } - \sum_{(\kappa,\ell)\neq(\kappa',\ell... ...)}^{S} + V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{T'} \Bigl] u_{\kappa ' \ell ' }$     (13)
$\displaystyle - \sum_{(\kappa,\ell)\neq(\kappa',\ell')} V_{(\kappa \ell ~,~\kappa ' \ell ' )}^{S} \bar w_{\kappa ' \ell ' }$      
$\displaystyle \sum_{(\kappa,\ell)\neq(\kappa',\ell')}\Bigl[ V_{(\kappa \ell ~,~... ...l ~,~\kappa ' \ell ' )}^{T'} \Bigl] ( u_{\kappa \ell } + \bar w_{\kappa \ell })$      

Phonon eigenvectors can be written as Bloch waves:

$\displaystyle u_{\kappa \ell } =u_{\kappa} e^{i \bf {Q} . {\bf {R}}_{\kappa \ell} }$     (14)
$\displaystyle \bar w_{\kappa \ell } =\bar w_{\kappa} e^{i \bf {Q} . {\bf {R}}_{\kappa \ell} }$      

where $\vec R_{\kappa \ell}$ is the equilibrium position of the atom $(\kappa \ell )$. The eigenproblem can then be written as :
$\displaystyle -\omega^2 u$ $\textstyle = - Z ~ C ~ Z ~ u +Z ~ C_0 ~ u - V^{SM} ~ u +V_0^{SM} ~ u$   (15)
  $\textstyle +(Y ~ C_0 - Z ~ C ~ Y - V^S - V^T + V_0^S+V_0^{T'}) ~ \bar w$    
$\displaystyle \bar w$ $\textstyle =( V_0^S+V_0^T-V^S-K+Y ~ C_0 - Y ~ C ~ Y )^{-1} .$    
  $\textstyle (Y ~ C ~ Z -Y ~ C_0 + V^S + V^{T'} - V_0^S - V_0^{T'} ) . u$    

where $u$ and $\bar w$ are now $3N$ dimensional vectors ( N being the number of atoms in a cell), Y,Z,K are $3N$ dimensional diagonal matrices whose diagonals are equal to the shell charges, ionic charge and spring constants respectively. The matrices C (coulombian), and the V's are obtained substituting equation 14 into 15 and summing over all the cells up to convergence. The subscript ``$_o$'' in $C_0$ and in $V_0$'s denotes the so said ``self-forces''. The dynamical matrix is found from equation (15) by elimination of $\bar w$. The summation of $C$ and $C_0$ over the cells needs a special treatment because of the Coulomb long range.

Treatment of the Coulomb long-range interaction

Coulomb long-range makes the real-space summation of forces very slow because of the $\frac{1}{r}$ behaviour. On the other hand in reciprocal spaces Colomb interaction behaves as $\frac{1}{k^2}$ for a point charge. Convergence at high $K$ can be accelerated by gaussian-spreading the point charges, and introducing thus a gaussian factor in the Coulomb interaction. This spreading changes the dynamical matrix, by a little amount if the spreading is small, and the correct result can be recovered by subtracting the difference between the field of a gaussian charge and the field of a point charge. This difference tends to zero very rapidly for distances bigger than the charge spreading and thus the real-space summation converges vary fast. The Coulomb contribution to the dynamical matric is thus made of two contribution : a reciprocal space summation for spreaded charges, and a real space summation for the difference between a point charges and the spreaded charges.

Reciprocal space summation for Coulomb interactions

A point charge $Z$ is spreaded as $\rho(r)$ :

\begin{displaymath} \rho(r)=Z ~(\frac{1}{2 \sigma^2_c \pi })^{3/2}exp( - \frac{r^2}{2 \sigma^2_c}) \end{displaymath}

and in reciprocal space

\begin{displaymath} \rho(k)=Z ~exp( - \frac{k^2 \sigma^2_c }{2 }) \end{displaymath}

The dynamical-matrix block between atom $\eta$ and atom $\xi$ is given by :

\begin{displaymath} Z_{\eta} Z_{\xi} \sum_l exp(i \vec{Q} . (r_{\xi l} -r_{\eta... ...partial_{\eta 0 x} \partial_{\xi l y} V(r_{\xi l}-r_{\eta 0} ) \end{displaymath}

and reciprocal space

\begin{displaymath} \frac{Z_{\eta} Z_{\xi}}{(2 \pi)^3} 4 \pi \int d k^3 \frac{e... ...{k^2} \sum_l exp(i \vec{Q} . (R_l + r_{\xi 0} -r_{\eta 0} )) \end{displaymath}

where $R_l$ is a lattice vector. The above expression can be simplified by the help of the following identity :

\begin{displaymath} \sum_{R_i} exp(i K.R_i) = \frac{(2 \pi)^3}{V} \sum_G \delta(K-G) \end{displaymath}

where V is the volume of the lattice unit cell and $G$ runs over the reciprocal lattice. Finally the matrix block contributing to $Z ~C ~Z$ of equation(15)is given by :


\begin{displaymath} \frac{Z_{\eta} Z_{\xi} 4 \pi }{(2 \pi)^3} \sum_G \frac{ exp... ...G+Q)_x (G+Q)_y exp(iG.( r_{\xi 0} -r_{\eta 0} )) }{ (G+Q)^2 } \end{displaymath} (16)

while the self term ($C_0$ of equation(15)) for atom $\eta$ is obtained by setting $Q=0$ and $Z_{\eta}=1$ in the above expression. In calculating the term given by formula (16) we are actually doing an error : as we sum over the whole lattice we are considering also a spurious affect coming from the interaction of the charge with itself. In a case of rigid ions the self term $ Z C_0$ simply subtracts to $Z ~C ~Z$ and the spurious terms in $C_0$ and in $C$ would delete each other. However in the case of polarizable shells, as can be seen in equation (15), $C_0$ and $C$ have a life on their own, being multiplied sometimes by $Z$ and some other time by $Y$. The spurious term has thus to be removed from $C$ and $C_0 ~Z^{-1}$. That can be obtained removing to the diagonal blocks of $C$ and $C_0 ~Z^{-1}$, for each atom $\eta$ the following quantity :

\begin{displaymath} \frac{4 \pi}{3 \sigma_c^3} \sqrt{\frac{1}{(2 \pi)^3}} \end{displaymath}

That corresponds to the derivative of the electric field at the center of a $3D$ gaussian distribution of charge.

The following subsection explains how the difference between a point charge and a gaussian one is recovered

Real-space summation of the Coulomb gaussian tail

TODO

Screening

Screening can be modeled by a dividing by the static dielectric function

\begin{displaymath} \epsilon(G+Q) \end{displaymath}

the addends of the sum in equation (16)

Files structure and use of OpenPhonon

Structure definition

The first step is to run the program preparation.py, followed by the name of a file containing the structure : 

python preparation.py <structure_file_name>
( see section 6 to see how to invoke OpenPhonon at ESRF or with the pre-compiled version). The structure file specifies the atoms and their position in the cell. This is accomplished by setting the variables AtomNames, CellVectors, PositionsList at the beginning of the input file ( note that the input file is simply executed as a Python script, so one has to use Python sintax).

We show here an example input file for preparation.py referring to a literature case [4]. It can be found in the distribution archive ( binary or source) in the directory test_OP: 

AtomNames=['Cu',   'O',  'Nd' ]
cellvectors=Numeric.array( 
         [ 
           [3.95  , 0.0     ,  0.0],
           [0.0   , 3.95    ,  0.0], 
           [0.0   , 0.0     ,  12.06556] 
         ] )
PositionsList=[
          Numeric.array(  
              [ [0.0,       0.0,       0.0  ],
                [0.5,       0.5,       0.5  ]
              ]    
          ),
          Numeric.array(
              [ [0.5,       0.0,       0.0  ],
                [0.0,       0.5,       0.0  ],
                [0.5,       0.0,       0.25 ],
                [0.0,       0.5,       0.25 ],
                [1.0,       0.5,       0.5  ],
                [0.5,       1.0,       0.5  ],
                [1.0,       0.5,       0.75],
                [0.5,       1.0,       0.75],            
              ]
           ),
           Numeric.array(
              [ [0.0,       0.0,       0.352],
                [0.5,       0.5,       0.148],
                [0.5,       0.5,       0.852],
                [1.0,       1.0,       0.648]        
              ]
           )
]
SeekedTetragon=[('Cu',0),('O',1),('Nd',0),('Nd',2)]

CellsCoo=[(0,0,0),(0,0,0),(0,0,0),(0,0,0)]   # optional
dr_shell=0.05                                # optional
num_nei=3                                    # optional
The SeekedTetragon array must specify a non degenerate 3-D solid that is rotated and translated in all possible ways to look for symmetries and equivalent atoms.

CellsCoo defines the cell in which the vertices are. This variable is optional and is set by default to all vertices in the same cell.

dr_shell defines the width of a shell. Two-body bonds are catalogued according to their length and the type of atoms that they concern. This variable is set by default to 0.1.

Finally num_nei, which is set by default to 2, is the number of neighbours.

The distribution archive contains an example file ( preparationinput ).

This file being set, the script preparation.py can be run (python preparation.py structure-file-name). The script does a rough evaluation of the symmetries and of inequivalent bonds and creates two files. Both files need to be stored somewhere for subsequent use by the program that calculates dispersions. One of these file, cella contains all the needed structural informations, while the other,parinput is a template file that can be edited with the parameters of the interatomic potentials as explained in the next subsection.

Definition of Potentials and eigenvectors calculations

The distribution archive contains parinputedited that corresponds to the test case [4]. Such a file describes the interatomic potentials. Units are CGS. Some knowledge of Python and of the Numerical extension (NumPy) is necessary. Interactions are book-keept by a Python dictionary. For example the dictionary entry 

BK_L['O_G1']['Nd_G0']
describes the longitudinal Born-Karman interactions between oxygen atoms of the equivalency group $1$ and the Lantanium atoms of the equivalence group $0$. The spring values are given in a list where each entry corresponds to a O_G1-Nd_G0 distance. This is a typical example extracted from the parameter file ( file parinputedited in the example directory texttttest_OP ): 
# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
# Interactions between 'O_G1' and 'Nd_G0' 
# Shell N0 going from 2.327062 to 2.327062  like O7 in 0 and Nd2 in 0 0 0 
# Shell N1 going from 4.584508 to 4.584508  like O7 in 0 and Nd2 in -1 0 0 
# Shell N2 going from 5.192371 to 5.192371  like O7 in 0 and Nd1 in 0 0 1 

BMpar.A= 2000  * ev_erg 
BMpar.R= 0.316  * A_cm
distances= [ 2.327062 * A_cm ,4.584508* A_cm ]
BMpar.ZZ= Z_['Nd_G0']*Z_['O_G0' ]
 
print " ['Nd_G0']*Z_['O_G0' ]  "

BK_L['O_G1']['Nd_G0']=map(bm_L,distances)
BK_T['O_G1']['Nd_G0']=map(bm_T,distances)

where $bm_l$ is an user defined function (see input file) of distance.

The parameters file specifies also other thing like, atom polarizability, Charges and analytical potentials. A few words should be spent on analytical potentials that can be of the Born-Mayer, Van der Walls, Lenard Jones type. They are specified by atom-dependent parameters, contrarily to the spring constants which are interaction (bond) dependent, .i.e. are given for each couple of interacting atom-types.

However in the literature one often finds analytical parameter that are bond dependent. In this case on can include them as spring constant. this is done in our example perinputedited file, where the function $bm_L$ is mapped on the shell distances to get spring constants.

Then the calculation is performed running dispersionDeb.py under Python with three arguments: 

python dispersionDeb.py  <parinputedited-filename> <cella-filename> 
  <additional-filename>

The significance of parinputedited and cella has been explained before. The last parameter specifies a file containing some parameters concerning the q-points at which phonons are calculated and other parameters governing the treatement of coulomb interactions. This is an examples : 

astar=2*math.pi/3.95
cstar=2*math.pi/12.06556
eps=0.0011
Q=Numeric.array( [ [astar*0.01*i+eps,-astar*0.01*i+eps,0] \
 for i in range(1,100)])
debyerange=[8,8,2]
debye=astar/4
SMcellrange=[2,2,2]
cellrange= [-1,-1,-1]
Kcellrange=[8,8,2]
sigmacellrange=[8,8,2]
sigmacharge=3
IS_IT_SCREENED=1
selectors=[(1.0,1.0,0.0), (1.0,-1.0,0.0), (1.0,0.0,1.0), (-1.0,1.0,0.0)]

The distribution archive contains an examples file (dispersionDebinput).

The q-points are defined by the variable Q, cstar and astar being just auxiliary variables. When summing up the contributions from spring constants a loop is done over all atoms in several cells. The variable SMcellrange gives the extent of this loop over the unit cells in the three directions a,b,c. Giving $SMcellrange=[0,0,0]$, for example, implies considering just the central unit-cell, but that might loose several interactions, particularly in the considered case where we include also second neighbours. For the example considered, $SMcellrange=[2,2,2]$ is enough. The variable $cellrange$ specifies the extent of the summation of analytical potentials, namely Born-Mayer, Van der Waals, Lennard Jones. The variables $IS\_IT\_SCREENED$ determines if the Coulomb interactions are screened or not. In the case of unscreened interactions the summation shows convergence problems that are cured by decomposing it partly in reciprocal space and partly in real space thanks to an appropriate charge smearing ( see Unisoft manual for details [5]). The Gaussian charge smearing is governed by $sigmacharge$, given in Angstroem. The extensions of reciprocal-space and real-space summation are given by the variables $Kcellrange$ and $sigmacellrange$ respectively. Taking a large sigma for charge smearing will make Kcellrange smaller( faster reciprocal space convergence) but will give a bigger sigmacellrange. Convergence has to be checked against these variables.

In the case of a screened potential, two variables have to be given : $debye$ which multiplies the distance in the exponential damping of the Coulomb potential and the extent of the real space summation, debyerange.

Finally the variable $selectors$ gives a set of directions against which the simmetry of the eigenmodes is checked.

The output is contained in three files : storeall, disp_res, disp_ressq. The ASCII files disp_res, disp_ressq contain the squares and the values of the phonon frequencies respectively. The file storeall stores in Python ``pickled'' form the calculated eigenvectors and frequencies. This file is finally used by the program scattering_shifted.py to calculate scattering intensities.

If the variable $selectors$ is set, an additional file, named polar_res, where for each directions the index of the eigenvalues ( referring to the ordering given in the output file) is given following the order of growing eigenvector projections over the $selectors$ direction. This can help to retrieve eigenvectors of specific symmetries, whose signal is enhanced by the experimental geometry.

Dipolar Fluctuation Model with tensorial force constants

The general expression for the dynamical matrix for the dipolar fluctuation model (DFM) is [8]:

\begin{displaymath} D(\vec{q}) = P(\vec{q}) - T^{\dagger}(\vec{q})[K + S(\vec{q})]^{-1} T(\vec{q}) \end{displaymath}

where $P(\vec{q}),T(\vec{q}),S(\vec{q})$ are Fourier transform of tensorial interaction matrices between couples of sites. The DFM is activated by setting the variable USE_DFM in the input file for dispersionDeb.py : 
USE_DFM=1

Then one has to fill in the entries for Tens_P, Tens_S, Tens_T. Tensor K is scaled to $1$ by rescaling S and T [8]. The following is an example for hcp Cobalt : 

# XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
# Tensorial Interactions between 'Co_G0' and 'Co_G0' 
# Shell N0 going from 0.996332 to 0.996332  like Co1 in 0 and Co0 in 1 0 1 
# Shell N1 going from 1.000000 to 1.000000  like Co1 in 0 and Co1 in -1 1 0 
# Shell N2 going from 1.411622 to 1.411622  like Co1 in 0 and Co0 in -1 1 0 
# Shell N3 going from 1.624000 to 1.624000  like Co1 in 0 and Co1 in 0 0 1 
# Shell N4 going from 1.729936 to 1.729936  like Co1 in 0 and Co0 in 0 2 1 
# Shell N5 going from 1.732051 to 1.732051  like Co1 in 0 and Co1 in -2 1 0 
# Shell N6 going from 1.907191 to 1.907191  like Co1 in 0 and Co1 in 0 1 -1 
# Shell N7 going from 2.000000 to 2.000000  like Co1 in 0 and Co1 in 0 2 0 

hcp_a=2.51
hcp_c=1.624*hcp_a

(a1,b1,g1,d1,
 a2,e2,b2,g2, 
 a3,b3,g3,d3)= [0]*12

(a1,b1,g1, 
 a2,b2,g2,e2,
 a3,b3,g3)=                  (    11925.5, 
                                  -1267.4,
                                  36831.4,
                                  -1819.8,
                                  40089.8,
                                   2171.9,
                                   1502.5,
                                    779.8,
                                    455.2,
                                  -3507.4)

Tens_P['Co_G0']['Co_G0']=[  [Numeric.array([     hcp_a/math.sqrt(3.0) ,
                                                     0, hcp_c/2.0  ]    ),
                           Numeric.array([[a1,0 ,d1],
                                  [0 ,b1,0 ],
                                  [d1,0 ,g1]])
                          ], 
                          [Numeric.array([  0 , hcp_a , 0 ]    ),
                           Numeric.array([[a2  ,e2 ,0],
                                  [-e2 ,b2 ,0 ],
                                  [ 0  ,0 ,g2 ]])  
                          ], 
                          [Numeric.array([  -2*hcp_a/math.sqrt(3.0) , 
                                                    0, hcp_c/2.0 ]    ),
                           Numeric.array([[a3  ,0  ,d3],
                                  [0   ,b3 ,0 ],
                                  [d3  ,0  ,g3 ]])  
                           ]
                        ]

one can see that the general form of Tens_P is given, for each shell, by a typical site-to-site vector and the corresponding tensor. The other tensors for the other interaction of the same shell are found by similitude transforming according to the simmetry group of the crystal. Tens_S and Tens_T are specified in analogous way.

Jahn-Teller tensorial force constants

Jahn-Teller tensorial forces can be treated in an analogous way[6,7]. The parameters Tens_JT has to be specified in the same way as Tens_P, Tens_S, Tens_T but in addition one has to add for each interaction two vectors that are used to calculated the JT factor: this factor is given by the product of two sines having as arguments the scalar product of $Q$ and the vectors. The JT model is activated by setting the variable USE_JT in the input file for dispersionDeb.py : 

USE_JT=1
the way to input the vectors for the JT factor is shown as follows : 
B_JT=100000.0
a=3.95
b=a
Tens_JT={}
Tens_JT['O_G0'] = {}
Tens_JT ['O_G0']['O_G0']= [
                           [    
                             Numeric.array([a/2.0,b/2.0,0]),
                             Numeric.array([[0,B_JT,0.0],
                                            [0,0,0.0],
                                            [0,0,0.0]] ),
                             Numeric.array([a/2.0,0,0]),
                             Numeric.array([0,b/2.0,0])
                           ], 
                           [    
                             Numeric.array([a,0,0]),
                             Numeric.array([[-B_JT,0.0,0.0],
                                            [0,0,0.0],
                                            [0,0,0.0]] ),
                             Numeric.array([a/2.0,0,0]),
                             Numeric.array([a/2.0,0,0])
                           ]
                         ]

where the vectors follows the $3X3$ tensor.

Angle Bond Force Constants

The following input example ( file parinputedited ) 

angleBonds={}
angleBonds['Sb_G0']={}
angleBonds['Sb_G0']['Sb_G0']={}


angleBonds['Sb_G0']['Sb_G0']['Co_G0']=[

[
Numeric.array([0.25   , 0.09212,  -0.08537])*9.035,
0.0002 , 
Numeric.array([0.   , -2*0.15788 ,  0.])*9.035,
0.0002,
anglestrenght,
]

establishes an angular spring for the angle having vertex Sb_G0 ( the first key of angleBonds ) and two sides , one made by the segment connecting Sb_G0 to Sb_G0 (second key) and having components dx,dy,dz = Numeric.array([0.25 , 0.09212, -0.08537])*9.035 . The other side goes from the vertex to Co_G0 and has components Numeric.array([0. , -2*0.15788 , 0.])*9.035.

To activate these interaction you have to write the following keywoards in the input file 

USE_ANGLE_BOND=1
AB_cellrange=[1,1,1]

A scan for every possible atomic triplet is done looping over a cell range whose span is given by the three components AB_cellrange=[1,1,1]. All rotations of the structure symmetry group are then applied to every triplet to check if a rotation can bring the triplet to coincide with the one given by angleBonds, with a tolerance given by 0.0002 ( second and fourth entry of an angleBonds contribution ). The contribution to the dynamical matrix is calculated for an energy given by anglestrenght*(angle-angle_0)**8

Use of an arbitrary screening function for Coulomb interactions

You activate the contribution from for the screened Coulomb interactions by setting the variable USE_SCREENING_FUNCTION=1 in the input file for dispersionDeb.py. The contribution is calculated in reciprocal space. As analytical formula in real space are not easy for an arbitrary screening function the contribution is calculated only in reciprocal space, so that a small enough $\sigma_c$ has to be used and the $K$ range of the summation must be large enough to assure convergence. The following is a tipycal input file for a calculation with a user defined screening function. 
astar=2*math.pi/3.95
cstar=2*math.pi/12.06556
eps=0.0011
Q=Numeric.array( [ [astar*0.5*i+eps,-astar*0.5*i+eps,0] for i in range(1,2)]    )
debyerange=[ -12, -12, -6]
debye=astar/4
SMcellrange=[2,2,2]
cellrange= [-1,-1,-1]
Kcellrange=[-6,-6,-12]
sigmacellrange=[-4,-4,-2]
sigmacharge=0.5
IS_IT_SCREENED=0

USE_SCREENING_FUNCTION =1
ScreenKrange = [6,6,12]

selectors=[(1.0,1.0,0.0), (1.0,-1.0,0.0), (1.0,0.0,1.0), (-1.0,1.0,0.0)]

Note the $K$ range that is specified by the variable ScreenKrange. The contribution from Coulomb interaction without screening function is still active (note IS_IT_SCREENED=0) but gives zero contributions as Kcellrange and sigmacellrange are set to zero. one could still calculate both contribution if the inputted screening function is

\begin{displaymath} \frac{1}{\epsilon}-1 \end{displaymath}

In case on lets IS_IT_SCREENED=1 the contribution from Debye potential must be made silent by setting debyerange to negative values.

The screening function must be set in the parinputedit with the following format : 

def SCREENING_FUNCTION(self, Q):
  moduli2 = Numeric.sum(Numeric.transpose(Q*Q))
  astar=2*math.pi/3.95
  debye=astar/4
  debye=debye*debye
  # res=(1.0/(1.0+debye/moduli2) )-1
  # debye=0.0
  res=(1.0/(1.0+debye/moduli2) )
  return res
Thats just an example. You can write whatever function you whish provided that you return either an array of scalars or an array of $3X3$ tensors, the arrays having the same lenght than $Q$, as it is shown in the above example.

Scattering intensity

The dynamical structure factor could be obtained from the inelastic scattering of the X-ray under the following conditions [9]: i) the scattering process is dominated by the Thompson term and there is no electronic excitation in the energy transfers region, ii) adiabatic approximation is valid, so that the center of mass of the electronic clouds could be identified with the nuclear one. The first assumption is simultaneously satisfied with X-ray with energy in the $\approx 10\div 25~keV$ energy range used for IXS experiments, when we look at energy transfers $\leq 0.4 ~eV$ ($\leq 100 ~THz$). The second one is a very general assumption in condensed matter theory. Under these assumptions we can write the inelastic cross section for one-phonon process as follows. First, we decompose the moment as $\bf {Q}=\bf {G}+\bf {q}$ (in the reciprocal lattice vector $\bf {G}$ plus the phonon pseudo-momentum $\bf {q}$; hence we obtain for a given angle $\Omega$ and frequency $\omega$:


$\displaystyle (\frac{\partial^2\sigma}{\partial\Omega\partial\omega}) \propto$ $\textstyle \frac{k_i}{k_f} (\bf {\hat{e}_i} \cdot \bf {\hat{e}_f})^2 \sum_{\bf ... ...bf {R}_{\kappa}} ~\bf {Q} \cdot \bf {e} (\kappa \vert^{\bf {q}}_\jmath)~\vert^2$    
  $\textstyle \times \langle n_\jmath+ \frac{1}{2}\pm \frac{1}{2}\rangle ~\delta(\omega\mp\omega_\jmath (\bf {q})) ~\Delta(\bf {Q}\mp(\bf {G}+\bf {q}))$   (17)

where $m_{\kappa}$ is the mass for the atom $\kappa$ in the unit cell, $e^{W_{\kappa}}$ is the Debye-Waller factor (which is not calculated in the present version of OpenPhonon), and $n_\jmath$ is the Bose factor for the mode $\jmath$. The phonon polarization term is $\bf {e} (\kappa \vert^{\bf {q}}_\jmath)=u_{\kappa}\sqrt{m_{\kappa}}$ of Eq. 14.

Then the the inelastic cross section can be written as:


$\displaystyle (\frac{\partial^2\sigma}{\partial\Omega\partial\omega}) \propto$ $\textstyle \frac{k_i}{k_f} (\bf {\hat{e}_i} \cdot \bf {\hat{e}_f})^2 \sum_{\bf ... ...\kappa}} ~\bf {Q} \cdot \bf {e} (\kappa \vert^{\bf {G}+\bf {q}}_\jmath)~\vert^2$    
  $\textstyle \times \langle n_\jmath+ \frac{1}{2}\pm \frac{1}{2}\rangle ~\delta(\omega\mp\omega_\jmath (\bf {q})) ~\Delta(\bf {Q}\mp(\bf {G}+\bf {q}))$   (18)

This corresponds to the neutron scattering cross section for 1-phonon scattering provided one replaces the nuclear average scattering length $\bar{b}_{\kappa}$ with the atomic form factor $f_{\kappa}(\bf {Q})$ for the atom $\kappa$ in the unit cell, and the polarisation factor of the Thompson scattering $(\bf {\hat{e}_i} \cdot \bf {\hat{e}_f})^2$. The latest term is fixed $=1$ in this version, because most of the experiments considered where in the scattering confguration with $(\bf {\hat{e}_i} \cdot \bf {\hat{e}_f})^2=1$ [4] or $(\bf {\hat{e}_i} \cdot \bf {\hat{e}_f})^2\approx 1$ (low $\bf {Q}$). Moreover for high energy resolution inelastic X-ray scattering we have $\frac{k_i}{k_f} \cong 1$. For a more detailed description of the IXS cross section, see also the review of E. Burkel [10] and references therein.

In OpenPhonon the calculation is made calling the program scattering_shifted.py in the following way : 

python scattering_shifted.py storeall paramfile
where paramfile is the name of a file like this : 
 
BrillShift_List=[(6,0,0)]
scatt_dictio_f0={ 'Cu_G0':'Cu1+', 'O_G0':'O','O_G1':'0',\
     'Nd_G0':'Nd','Nd_G1':'Nd'   }
scatt_dictio_f12={ 'Cu_G0':'Cu', 'O_G0':'O','O_G1':'0',\
     'Nd_G0':'Nd','Nd_G1':'Nd'   }
Lambda=0.783867
Temp = 15.0
A sample input file is delivered in the distribution archive under the name scatt_shinput. The intensities are calculated over the set of points defined by the variable $Q$ defined in the previous input files, at which points the eigenvectors and frequencies have been calculated and stored in the file store all. The variable $BrillShift\_List$ specifies a set of shift. So one can calculate once the dispersion from $\Gamma$ to $X$ (just to to give an example ) and use the resulting eigenvectors to calculate the scattering intensities at several reciprocal vector zones.

The atomic scattering factors are read using the Dabax database. The $f0$ factors are read from WaasKirf compilation and f1,f2 from Windt datafile.

To specify the scattering factors of the atoms one has to specify the wavelength of X-rays and the two dictionaries, one for $f0$ and the other for (f1,f2). The dictionary keywords are the equivalence class names and their values are Dabax record names like "Cu" for neutral copper or "Cu1+" for copper ion. These names must exist in the database otherwise an error occurs. Finally the temperature ( in Kelvin) has to be entered to properly take into account the Bose-Einstein statistics.

The results are stored in the file scatt_results in the form:
Qx,Qy,Qz,frequency(1),intensity_stoke(1),intensity_antistoke(1),..,
frequency(j),intensity_stoke(j),intensity_antistoke(j),..,
frequency(N),intensity_stoke(N),intensity_antistoke(N)
where j indicates the mode in increasing $\Gamma$ frequency order.

Calculation of the partial phonon density of state (PDOS)

PDOS is constructed as a N points histogram, whose frequency coordinates range from zero up to the maximum eigen-frequency of the system. The histogram bins are filled calculating the eigen-frequencies over the Brillouin zone on a $n \times n \times n$ grid. Symmetries are used to reduce the number of needed calculations. To get the partial DOS for a given site group, the contributions are weighted when added on a particular bin. The weight is calculated as the modulus square of the eigenvector projection over the considered degrees of freedom.

The calculation of the frequencies in all the point of the grid is performed running the program generapuntiDeb.py: 

openphonon generapuntiDeb.py parinput cella generapuntiinp

where the files parinput and cella are the same of section 4.1, and the file generapuntiinp is a file as generapuntiDebinput in the examples directory, with the following structure: 

IS_IT_SCREENED=1
Kcellrange=[8,8,2]
sigmacellrange=[8,8,2]
cellrange= [-1,-1,-1]


debyerange=[8,8,2]
SMcellrange=[2,2,2]
astar=2*math.pi/3.95
kdebye=astar/4

eps=0.001


Nsampling=21

Qall=Numeric.array( [ Numeric.array([0.00001,0,0]) +  qx*cella.Brillvectors[0]*0.5/(Nsampling-1)
                      + qy*cella.Brillvectors[1]*0.5/(Nsampling-1) + qz*cella.Brillvectors[2]*0.5/(Nsampling-1)
                      for  qx in range(0,Nsampling)
                      for  qy in range(0 , Nsampling)
                      for  qz in range(0,Nsampling)
                      ] )


thetasym=45

bdx=math.cos(thetasym*math.pi/180.)
bdy=math.sin(thetasym*math.pi/180.)

basedomain=Numeric.array([ bdx, bdy   ,1])
domains=Numeric.array( [
                           [ bdx, bdy   ,1],
.............................
.............................
.............................

where Nsampling correspond to the grid density $n$.

Subsequently, the number of state per bin is calculated by running the program pa_dos.py: 

openphonon pa_dos.py storeall Number_of_histograms

where storeall is the one created by generapuntiDeb.py, and Number_of_histograms = N as described above.

3D rendering

Structure rendering

The cell structure can be visualized as soon as one creates the file <structure_file> described above. Then one has to create just one additional file, let's call it rendering-file, defining a few variables as in the examples produced below : 
Rendering_Radii=[0.5,0.6,0.7 ]
colors=[ [0.0,0.0,1.0], [1.0,0.0,0.0], [0.0,0.5,0.0] ]
nneig=[1,1,1]

the Rendering_Radii are the visualised sphere radii for the atomic species defined in <structure_file>. The colors are defined as RGB triplets, $1.0$ being the brightest value. The variable $nneig$ specifies how many times the unit cell has to be repeated in $x,y,z$ directions, $[0,0,0]$ meaning just the central one. For example, $[1,1,1]$ specifies a range spanning from $-1$ to $+1$ for each direction. The display window can be popped up using the following command : 

<pymol.com>  structureRendering.py <rendering-file-name> \
   <structure_file-name>

where pymol.com is the $pymol$ interpreter.

The example given in the test_OP directory ( binary distribution ) can be runned using the openPhonon wrapper in the following way : 

openphonon structureRendering.py render preparationinput

EigenModes rendering

The visualisation of eigenmodes is done via structure rendering. 

<pymol.com>   eigenRendering.py  <rendering-file-name> \
      <structure_file-name> <storeall_name>

where storeall is the output file of dispersionDeb.py and rendering-file has the additional information of which mode has to be visualised : 

Rendering_Radii=[0.5,0.6,0.7 ]
colors=[ [0.0,0.0,1.0], [1.0,0.0,0.0], [0.0,0.5,0.0] ]
nneig=[1,1,1]
NPscan=48
NPvector=41
fampli= 2.0

where NPscan specifies $q=Q[NPscan]$ ( NPscan starts from 0)

and NPvector is the eigenvector number ( NPvector goes from 0 up to $3*Natoms-1$ is).


Users Notes for ESRF and precompiled version

At ESRF all the scripts are already installed on the amber machines. You can call them using the wrapper openphonon (as you do with a pre-compiled version). For example, to run 
python scattering_shifted.py storeall paramfile

one types : 

openphonon  scattering_shifted.py storeall paramfile

The pymol interpreter is called also through openphonon as : 

openphonon   eigenRendering.py  <rendering-file-name> <structure_file-name>\
   <storeall_name>

Installation

Installation is done by downloading the distribution archive
(www.esrf.fr/computing/scientific/) and following the instructions therein. One needs Python2.X, a C++ compiler, and the NumPy extension package linked with lapack3.0 .

Alternatively one can download the following file containing a self contained tree of binary codes (pre-compiled version, only for linux):
ftp.esrf.fr/pub/scisoft/ESRF_sw/opbin.tgz

This binary distribution should work with Suse or Redhat $>7.0$. One downloads it on the root directory of a Linux machine. It expands as /opt/ESRF_sw. Binaries and libraries are in /opt/ESRF_sw/bin and /opt/ESRF_sw/lib respectively. You have to set PATH and LD_LIBRARY_PATH accordingly. Test files are in /opt/ESRF_sw/test_OP. Then one should be able to call the openphonon wrapper. Then one replaces the command python by openphonon. Pymol is given along with the binary distribution. The pymol home page is at www.pymol.com.

If you are installing from sources be warned that Numpy might need to be linked against lapack3.0. On some machine the numpy provided lapack_lite give problems with Heigenvectors routine that is a kay ingredient of openphonon. It is easy to correct eventual problem by linking against lapack library, one has to modify the numpy setup.py makefile in the following way ( adapt directories and library name to your lapack installation ): 

# delete all but the first one in this list if using your own LAPACK/BLAS
sourcelist = [os.path.join('Src', 'lapack_litemodule.c'),
#              os.path.join('Src', 'blas_lite.c'), 
#              os.path.join('Src', 'f2c_lite.c'), 
#              os.path.join('Src', 'zlapack_lite.c'),
#              os.path.join('Src', 'dlapack_lite.c')
             ]
# set these to use your own BLAS
library_dirs_list = ["/scisoft/ESRF_sw/linux_i386/lib"]
libraries_list = [ "lapack","blas" ]

Bibliography

1
M. Born, K. Huang ``Dynamical theory of crystal lattices'' (1954) Oxford University Press (Oxford).

2
A. A. Maradudin and S.H. Vosko Rev. Mod. Phys. 40 (1968) 1.

3
P. Brüesch ``Phonons: Theory and Experiments'' (1982) Springer-Verlag (Berlin).

4
M. d'Astuto, P.K. Mang, P. Giura, A. Shukla, P. Ghigna, A. Mirone, M. Braden, M. Greven, M. Krisch, F. Sette, Phys. Rev. Lett. 88 (2002) 167002.

5
Eckold G., Stein-Arsic M., Weber H.J. ``Unisoft - a program package for lattice-dynamical calculations: user manual'' (1986) Kernforschungsanlage Juelich (Juelich).

6
D.V. Fil and others Phys. Rev. B 45 (1992) 5633.

7
C.-Z. Wang and others Phys. Rev. B 59 (1999) 9278.

8
N. Wakabayashi, R. H. Scherm, and H. G. Smith , Phys. Rev. B 25 (1982) 5122-5132.

9
F. Sette, G. Ruocco, M. Krisch, C. Masciovecchio, R. Verbeni, Physica Scripta T66 (1996) 48.

10
E. Burkel, Rep. Prog. Phys. 63 (2000) 171.

About this document ...

OpenPhonon: an open source computer code for lattice-dynamical calculations

This document was generated using the LaTeX2HTML translator Version 2002-2-1 (1.70)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split 0 manual.tex

The translation was initiated by Alessandro Mirone on 2006-03-21

Footnotes

... Mirone1
European Synchrotron Radiation Facility ESRF, BP 220, F-38043 Grenoble Cedex, France
next_inactive up previous
Alessandro Mirone 2006-03-21