Setting up a Point Charge Model which yields given Crystal Field Parameters

Sometimes some or all crystal field parameters are available, either from literature or determined from the high temperature expansion of inverse susceptibility etc. In these cases the point charge model should reproduce these crystal field parameters and we describe here a procedure which has been tested with some success on a number of cases.

As an example we consider the ternary compound TmNiC$_2$, denoting point charges on Tm, Ni and C by $C_{Tm}$, $C_{Ni}$ and $C_{C}$, respectively. The full calculation with all input files and logbook calc.html is available in examples/tmcni2 .

From the high temperature expansion of the inverse susceptibility the second order crystal field parameters of this orthorhombic phase have been determined by [38] to be $B_2^0=0.22$ meV and $B_2^2=-0.33$ meV. We make now use of the fact, that every crystal field parameter is linear dependent on the point charge, compare equations (103) and (104). Thus we get the following set of linear equations:


$\displaystyle B_2^0$ $\textstyle =$ $\displaystyle {\rm Ni}_{20} C_{Ni}+{\rm C}_{20} C_{C}+{\rm Tm}_{20} C_{Tm}$ 
$\displaystyle B_2^2$ $\textstyle =$ $\displaystyle {\rm Ni}_{22} C_{Ni}+{\rm C}_{22} C_{C}+{\rm Tm}_{22} C_{Tm}$ 
$\displaystyle 0$ $\textstyle =$ $\displaystyle C_{Ni}+2 {\rm C}_{C}+ C_{Tm}$(109)

The last equation ensures the charge neutrality and should also be taken into account. The coefficients ${\rm Ni}_{lm}$, ${\rm C}_{lm}$ and ${\rm Tm}_{lm}$ can be obtained by a point charge calculation with setting all charges zero except one (which is set to $1\vert e\vert$) and evaluating the computed crystal field parameter $B_l^m$.

The following code snippet shows how to calculate the coefficients C ${\rm C}_{20}$ and ${\rm C}_{22}$ using cif2mcphas in linux:

cif2mcphas -pc 20 -sp -nm Ni -ch Tm=0,Ni=0,C=1 TmNiC2.cif
#  extract coefficient 
. getvariable B20 Tm1_1.sipf
# stream editor sed is used to ensure if it is zero ("not found") a 0.0 is stored in C20 instead of "not"
export C20=$(echo $MCPHASE_GETVARIABLE_VALUE | sed s/not/0.0/ )
. getvariable B22 Tm1_1.sipf
export C22=$(echo $MCPHASE_GETVARIABLE_VALUE | sed s/not/0.0/ )

...and in windows:

call cif2mcphas -pc 20 -sp -nm Ni -ch Tm=0,Ni=0,C=1  TmNiC2.cif 
echo B20=0.0 > dd
echo B22=0.0 >> dd
REM then push the sipf file onto dd - if it contains B20 or B22 these will be read
type Tm1_1.sipf >> dd
call getvariable B20 dd
set C20=%MCPHASE_GETVARIABLE_VALUE%
call getvariable B22 dd
set C22=%MCPHASE_GETVARIABLE_VALUE%

Figure 19: Convergence of the point charge calculation of the crystal field parameter $B_2^0$ in TmNiC$_2$.
Image B20convergence

Convergence is reached in the current system when point charges up to a distance of $r_{max}=20$ Å are taken into account, see fig. 19. For the radial integrals $\langle r^l \rangle$ of Tm$^{3+}$ we used the valus given by [39] including the shielding factors $\sigma$. Thus we arrived at a system of 3 equations (117) with 3 unknown charges. This system may be readily solved. However, it turns out, that the computed charges are unreasonable large and have unreasonable signs.

Figure 20: Point charge model for TmNiC$_2$, large spheres correspond to atomic positions, small spheres to electronic bond charges $E^1, \dots , E^6$.
Image model

Therefore, we introduced some additional point charges and placed these on the lines connecting a pair of atoms. This choice was motivated by the intuition, that the covalent bonding between the elements will lead to some anisotropic electronic charge density. We had to introduce six such point charges, these are indicated by $E^1, \dots , E^6$ and refer to electronic bond charges for C-C, C-Ni, C-Ni, Tm-C, Tm-Ni and Tm-Tm pairs, respectively (see fig 20) Here we show the corresponding part of the input cif file:

loop_
  _atom_site_label
  _atom_site_type_symbol
  _atom_site_fract_x
  _atom_site_fract_y
  _atom_site_fract_z
  _atom_site_U_iso_or_equiv
  _atom_site_adp_type
  _atom_site_occupancy
  _atom_site_site_symmetry_order
  _atom_site_calc_flag
  _atom_site_refinement_flags_posn
  _atom_site_refinement_flags_adp
  _atom_site_refinement_flags_occupancy
  _atom_site_disorder_assembly
  _atom_site_disorder_group
 Tm1 Tm 0.5000 0.0000 0.3850(2) 0.0078(2) Uani 1 4 d S T P . .
 Ni1 Ni 0.0000 0.0000 -0.0034(3) 0.0089(3) Uani 1 4 d S T P . .
 C1 C 0.0000 0.349(2) 0.1851(18) 0.0105(16) Uani 1 2 d S T P . .
; E1 C-C bond charge as 0 0.5 z as C1
 E1 E 0.0000 0.500(2) C1        0.0105(16) Uani 1 2 d S T P . .
; atomic radii  Tm 227 pm C 170pm  Ni 163pm
; E2 E3 C-Ni bonds at 
E2 E Ni1+(C1-Ni1)*163/(163+170) Ni1+(C1-Ni1)*163/(163+170) Ni1+(C1-Ni1)*163/(163+170)   0.0105(16) Uani 1 2 d S T P . .
;  ... and
E3 E Ni1+(C1-Ni1)*163/(163+170) Ni1+0.5+(C1-Ni1-0.5)*163/(163+170) Ni1+0.5+(C1-Ni1-0.5)*163/(163+170)	  0.0105(16) Uani 1 2 d S T P . .
; E4 Tm - C bond charge at
 E4 E C1+(Tm1-C1)*170/(227+170)  C1+(Tm1-C1)*170/(227+170)   C1+(Tm1-C1)*170/(227+170)    0.0105(16) Uani 1 2 d S T P . .
; E5 Tm - Ni bond charge at
 E5 E N1+(Tm1-Ni1)*163/(227+163) N1+(Tm1-Ni1)*163/(227+163) N1+(Tm1-Ni1)*163/(227+163)   0.0105(16) Uani 1 2 d S T P . .
; Tm - Tm bond charge at
E6 E  Tm1 Tm1+0.25  Tm1+0.25 0.0105(16) Uani 1 2 d S T P . .

The set of equations (117) has to be extended for these bond electron point charges and we get


$\displaystyle B_2^0$ $\textstyle =$ $\displaystyle {\rm Ni}_{20} C_{Ni}+{\rm C}_{20} C_{C}+{\rm Tm}_{20} C_{Tm}+\sum_{i=1}^6 E^i_{20} C_{E^i}$ 
$\displaystyle B_2^2$ $\textstyle =$ $\displaystyle {\rm Ni}_{22} C_{Ni}+{\rm C}_{22} C_{C}+{\rm Tm}_{22} C_{Tm} +\sum_{i=1}^6 E^i_{22} C_{E^i}$ 
$\displaystyle 0$ $\textstyle =$ $\displaystyle C_{Ni}+2 {\rm C}_{C}+ C_{Tm} +C_{E^1}+2C_{E^2}+2C_{E^3}+4C_{E^4}+2C_{E^5}+2C_{E^6}$(110)

By variation of the charges $C_{E^i}$ a set of point charges and crystal field parameters could be found, which reproduces closely the crystal field effects reported in [38], see table 3.


Table 3: Point Charge Model for TmNiC$_2$
atom bond charge($\vert e\vert$) da db dc
Tm   2.159 0.5 0 0.385
Ni   -0.339 0 0 -0.0034
C   3.230 0 0.349 0.1851
$E^1$ C-C -2.6945 0 0.5 0.1851
$E^2$ C-Ni -1.5506 0 0.17083 0.08887
$E^3$ C-Ni -0.3955 0 0.42609 0.34412
$E^4$ Tm-C -0.2008 0.21411 0.19955 0.2707
$E^5$ Tm-Ni -0.2084 0.20897 0 0.16233
$E^6$ Tm-Tm -0.2369 0.5 0.25 0.635
$B_2^0$ 0.22 meV        
$B_2^2$ -0.328 meV        
$B_4^0$ -0.00105 meV        
$B_4^2$ 0.00296 meV        
$B_4^4$ -0.00994 meV        
$B_6^0$ -3.604e-06 meV        
$B_6^2$ 7.969e-06 meV        
$B_6^4$ 1.574e-05 meV        
$B_6^6$ 2.11e-05 meV        

This set of parameters was used to calculate the thermodynamic and spectroscopic physical properties using mcphasit and mcdispit, respectively. Note that in such calculations the bonding charges $E^1, \dots , E^6$ have to be removed, because electrons are correlated and cannot move freely in a solid (compare Pauli principle, failure of the Drude model). Technically we have to remove the E charges from the list of atoms in mcphas.j to prevent mcphas from calculating too large energy u and mcdisp from calculating optical modes due to vibrations of electronic bonding charges E. This can be done with reduce_unitcell with option -delatom. In examples/tmnic2/calc.bat the procedure is described in detail.