< MuSR.Wien2k | Index | a collection of python hints >

Dipsum trial with ASE

Define the structure and view it
Calculate dipolar sums
Check Fe3O4
View muon sites
Save results with pickle
Point-charge muon potential

The example in this page is Fe₃O₄ above the Verwey transition.

Define the structure

See the example fe3o4.py. To run
import fe3o4.py
reload(fe3o4) # to upgrade, after first import and subsequent modifications
To show the parameters use:

fe3o4.a
fe3o4.fe3o4.get_initial_magnetic_moments()
fe3o4.fe3o4._get_atomic_numbers()
fe3o4.fe3o4._get_cell()  # shows numeric zeros of 5.14e-16
fe3o4.fe3o4._get_positions() # what units? __doc__ is not informative
fe3o4.fe3o4.positions # the same, so this is a numerical matrix
fe3o4.fe3o4.pop() # WARNING! it kills the last atom

Calculate the dipolar sums

This is the content of a file (bdipp.py) that does it:

import ase
import numpy as np
from ase.lattice.spacegroup import crystal
def dipsum(lattice,musite,a,n):
   """This routine is a new one that calculates dipolar sums
      in the finite lattice made of (2*n+1)^3 cells, 
      within a Lorentz sphere of radius R(Angstroems) 
      centered at musite. R = a, 2a, 3a,..., na and for each
      a row of Bdip components is added. Bdip is calculated 
      in each shell and finally cumsum-med
   """
   number_sites=lattice.get_number_of_atoms() 
   latticemu=lattice.extend(musite) # adds the muon site at the end
   # whose index is number_sites +1 for mu, -1 since indices start from 0
   Bdip=np.zeros((n,3))
   print latticemu.__getitem__(number_sites)
   for i in xrange(number_sites-1):
      site=lattice.__getitem__(i)
      if site.get_symbol() == 'Fe': 
                                 # this is for Fe3O4
         r = latticemu.get_distance(i,number_sites) # number_sites as an index
                                               # is the last, i.e. the muon
         ur = (site.get_position() - musite.get_position()) / r 
          # unit vector in the direction joining i-th site and muon site
         j=int(np.floor_divide(r,a)) # j is the index of the Lorentz shell 
         # of width a (0 is the index of the first Lorentz sphere) 
         if j < n: # within the largest Lorentz sphere inside lattice
              Bdip[j,:] = \
              Bdip[j,:] + \
              ( 3.0 * np.dot(site.get_initial_magnetic_moment(),ur) * ur - \
               site.get_initial_magnetic_moment() ) / r**3
   Bdip=Bdip*0.9274009  # conversion in Tesla
# m_mu dot Bdip_i = (gamma_mu hbar I) dot (mu_0/4pi) gamma_e hbar 
#                    (3 (S dot ur)ur - S)/r^3
# (mu_0/4pi) mu_B 1e-30= 0.9274009
   return Bdip
n_quadrant=7
n_cells=2*n_quadrant+1
a = 8.3940 
magmoms=np.array([(0., 0., -4.44),(0., 0., -4.44),(0., 0., -4.44),(0., 0., -4.44),
(0., 0., -4.44),(0., 0., -4.44),(0., 0., -4.44),(0., 0., -4.44),
(0., 0., 4.17),(0., 0., 4.17),(0., 0., 4.17),(0., 0., 4.17),
(0., 0., 4.17),(0., 0., 4.17),(0., 0., 4.17),(0., 0., 4.17),
(0., 0., 4.17),(0., 0., 4.17),(0., 0., 4.17),(0., 0., 4.17),
(0., 0., 4.17),(0., 0., 4.17),(0., 0., 4.17),(0., 0., 4.17),
(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),
(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),
(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),
(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),
(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),
(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),
(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),
(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0),(0.0, 0.0, 0.0)]) 
fe3o4=crystal(['Fe','Fe','O'], 
magmoms=magmoms, basis=[(0.12500,  0.12500,  0.12500), 
                (0.5, 0.5, 0.5), (0.25480,  0.25480,  0.25480)], 
                setting=2,  
                spacegroup=227, cellpar=[a, a, a, 90, 90, 90],
                size=(n_cells,n_cells,n_cells),pbc=False)
#mu=ase.Atom('H',((0.1309+n_quadrant)*a,(0.2856+n_quadrant)*a,
#                 (0.2838+n_quadrant)*a),mass=1.008/9)#0.1309 0.2856 0.2838
#mu=ase.Atom('H',((0.2856+n_quadrant)*a,(0.1309+n_quadrant)*a,
#                 (0.2838+n_quadrant)*a),mass=1.008/9)#0.2856 0.1309 0.2838
mu=ase.Atom('H',((0.2850+n_quadrant)*a,(0.2850+n_quadrant)*a,
                 (0.1302+n_quadrant)*a),mass=1.008/9)#0.2850 0.2850 0.1302
# run the three sites independently the muon above 250 is tunneling among these sites
# other runs for all three must be performed with moments along [1 -1 1], [-1 1 1] and [1 1 -1]
# M=2*0.9274009*sum(magmoms,0)/a**3 # bulk magnetization 
                                # = moment of a cell/volume of a cell
B=dipsum(fe3o4,mu,a,n_quadrant)
Bc=np.cumsum(B,0) # each row contains the sum for a Lorentz radius. 
Row 0 in B corresponds to column 0 in R.

Check Fe₃O₄

T>250K: run the script with three sites and take the average: it is nearly zero. Same turning moments

along [-1 1 1], [1 -1 1] and [1 1 -1]. In each case the Lorentz field is roughly 0.21T and a contact term of 0.21T justifies the observed field of 0.42T.

T_V<T<250K: same spin orientation, local average

View crystallographycally equivalent muon sites

This does it

import ase
import numpy as np
from ase.lattice.spacegroup import crystal
a = 8.3940
fe3o4=crystal(['Fe','Fe','O','H'], 
         basis=[(0.12500,  0.12500,  0.12500), 
                (0.5, 0.5, 0.5), 
                (0.25480,  0.25480,  0.25480),
                (0.2850, 0.2850, 0.1302)], 
                setting=2,  
                spacegroup=227, cellpar=[a, a, a, 90, 90, 90],
                size=(1,1,1),pbc=False)
ase.visualize.view(fe3o4)

To improve this simple view set Show bonds from the dropdown 'View' (or ctrl-b)

Save results with pickle

f=open('fe3o4.pickle','w')
data1={'Bc1111':Bc1111,'Bc1112':Bc1112,'Bc1113':Bc1113,'M111':M111} # this is a dictionary
pickle.dump(data1,f)
f.close()

Now try more fe3o4.pickle'. To load

f=open('fe3o4.pickle','r')
data=pickle.load(f) # now data contains the dictionary

Calculate muon potential

Ewald trick for convergence. Let's consider just point-charge ions (DFT will correct this). The muon sees the electrostatic potential {$V$} of these charges and in first approximation it localizes on a constrained minimum, such as, e.g. at a distance 1.1 A from Oxygen. Therefore a very simple strategy is to look for the minimum of {$V$} on this sphere. However the sum of the Coulomb potentials of lattice ions does not converge well, if at all. The Ewald trick is to sum over a cell in real space, to grant charge neutrality, and then to sum in Fourier space exploiting lattice periodicity. If {$\mathbf r$} is the muon site

{$ V({\mathbf r}) = \int d{\mathbf r}^\prime \rho_T({\mathbf r}^\prime) \phi(|{\mathbf r}-{\mathbf r}^\prime|) $}

Here

{$\phi(r)=1/4\pi\varepsilon_0 r$}

is the Coulomb function and, if {$\rho_c({\mathbf r})=\sum_l q_k\delta({\mathbf r}-{\mathbf r}_k)$} is the charge density inside a cell, the total density

{$\rho_T({\mathbf r})=\sum_{n,m,l}\sum_k q_k\delta({\mathbf r}-{\mathbf r}_k-n{\mathbf a}_1-m{\mathbf a}_2-l{\mathbf a}_3)$}

turns out to be the convolution of {$\rho_c$} with the lattice function

{$L(\mathbf{r}) =\sum_{m,n,l} \delta({\mathbf r} - m {\mathbf a}_1 - n {\mathbf a}_2 - l {\mathbf a}_3)$},

where {${\mathbf a}_i$} are the lattice base vectors. Then the first equation states that {$V$} itself is the convolution of {$\rho_T$} and {$\phi$}. Exploiting the convolution theorem the Fourier transform {$V(q)$} is therefore just the product of the Fourier transform of these functions

{$V({\mathbf q})=L({\mathbf q})\rho_c({\mathbf q})\phi({\mathbf q})$}

with

{$L(\mathbf{q}) =\frac{2\pi^3}{{\mathbf a}_1\times{\mathbf a}_2\cdot{\mathbf a}_3}\sum_{m,n,l} \delta({\mathbf q} - m {\mathbf b}_1 - n {\mathbf b}_2 - l {\mathbf b}_3)$}

where {${\mathbf b}_i$} are the reciprocal lattice base vectors. This allows us to compute

{$V({\mathbf q})=\frac{2\pi^3}{{\mathbf a}_1\times{\mathbf a}_2\cdot{\mathbf a}_3}\rho_c({\mathbf Q})\phi({\mathbf Q})\delta({\mathbf q}-{\mathbf Q})$}

restricted to reciprocal lattice vectors {${\mathbf Q}=m {\mathbf b}_1 + n {\mathbf b}_2 + l {\mathbf b}_3$} and to Fourier transform it back to obtain {$V({\mathbf r})$}. The last ingredient is the Fourier transform of the Coulomb function {$\phi$}. In order to compute it start with the Yukawa function {$\phi_\lambda(r)=e^{-\lambda r}/r$} and choose spherical coordinates with the {$\hat z$} axis along {$\mathbf q$}, so that the plain wave becomes {$e^{iqr\cos\theta}$}. The integral is straightforward. Then take the limit for {$\lambda\rightarrow 0$}. The Fourier transform of the Yukawa function is

{$\phi_\lambda(q)=\frac{4\pi}{\lambda^2+q^2}$}

Here is the python implementation (with a bug!)

Top

< MuSR.Wien2k | Index | a collection of python hints >

ASEStart

Dipsum trial with ASE

Table of Contents

Define the structure

Top

Calculate the dipolar sums

Top

Check Fe₃O₄

Top

View crystallographycally equivalent muon sites

Top

Save results with pickle

Top

Calculate muon potential

Top

ASEStart

Dipsum trial with ASE

Table of Contents

Define the structure

Top

Calculate the dipolar sums

Top

Check Fe3O4

Top

View crystallographycally equivalent muon sites

Top

Save results with pickle

Top

Calculate muon potential

Top

Check Fe₃O₄