using ic1ion as a stand-alone program to do crystal field calculations

ic1ion can be started by issuing the command ic1ion (reading input from single ion parameter file mcphas.ic). Alternatively, the single ion parameter file can be specified on the commandline as

ic1ion filename

ic1ion also accepts these options:

which outputs in the file results/ic1ion.out only the header containing the free ion (Coulomb and spin-orbit) and single ion (crystal field and magnetic fields) parameters.
The single ion parameter file has to be of the form (note the crystal field parameters in Wybourne notation can be calculated by pointc):
# single ion property file example for module ic1ion


# Configuration *********************************************
# (i) either specify single ion type by IONTYPE 


# (ii) or by giving parameters explicitly:
# from Griffith 1971 p 487  B=1115  C=4366 ZETA=533
# p379 CoIII=Co2+ is 3d7 configuration
conf = d7
# from equations p 83 A=F0-49*F4  B=F2-5*F4 C=35*F4
# we get F4=C/35  F2=B+5*F4
# however F0 is not determined from this table ...

# for charges and chrgplt to work you need to specify the number of electrons in the shell



# some keywords, each in a separate line (each keyword sets a flag to true and holds
# for the rest of the file to be read)

# --- add a magnetic field (Tesla) ---

# -- choose the basis to be used for output of eigenstates in ic1ion.out ----
# Either:
# basis=mSmL
# For |L,S,mS,mL> basis. Or

# For |L,S,J,mJ> basis (default if none specified).

# to allow output of more than 4 components of the eigenvector in ic1ion.out, 
# you need to specify the number X of eigenvectors here:

# eigenvectors=X

# --- some faster algorithms for matrix diagonalisation ...

# uncomment the next line to save matrices which are often used the first time they are
# calcultated in the directory results/mms, and then reread for subsequent instances.

# save_matrices

# if the next line is uncommented then it will use the Arnoldi method to calculate only 
# a few of the lowest eigenvalues in Icalc() [those for which the thermal occupation
# factor at 300K is nonzero]. Note: this method may sometimes give erroneous eigenvectors
# due to numerical instabilities. In addition, it may be disabled in Mac OS X versions.

# arnoldi

# Another more reliable method to find only a few of the lowest eigenvalue uses the 
# relatively robust representation algorithm in Lapack and can be accessibly by 
# uncommenting the keyword: "partial" below. This method takes about 2/3 the
# time of the full computation of Icalc(), whereas the arnoldi routine takes about
# 1/5 the time.

# partial

# Finally the full Hamiltonian matrix may be truncated at some levels, leaving only
# the lowest lying energy levels to be used in the calculation of the moments.
# This option is activated by uncommenting the line below, specifying a ratio X of
# levels to keep. For example, if X=0.2, for f2, 20% of the full 91 levels == 18
# levels will be used for the Icalc() calculation.
# This routine first calculates and diagonalises the singleion Hamiltonian without 
# the Zeeman or mean-field terms, and then uses these eigenvectors to make a unitary 
# transformation matrix to rotate the multipolar and L and S operators into this
# new basis. Then it truncates the matrix, leaving only the X fraction lowest 
# levels.
# This method is much faster (for low values of X) even than the arnoldi method, but
# is much more approximate. Note The most computationally extensive operations are 
# carried out at the start - the diagonalisation of the single ion Hamiltonian, and
# matrix-matrix multiplication to rotate the operator matrices. Once done the rotated 
# matrices are stored in memory and then reused when the <J> expectation values have 
# to be evaluated. the routine uses the Icalc_parstorage matrix for storing the 
# rotated matrices.
# One should experiment with the value of X and compare it with using the full matrix 
# diagonalisation before relying on the results. 
# This approximation method is probably not necessary for d-electron matrices (as 
# those are usually quite small), but for f-electrons it is quite useful!
# Uncomment the next line to use the truncate method with 1/5th of all levels

# truncate_matrix=0.2

# a list of crystal field parameters (Llm denote Wybourne parameters)

# Blm denote Stevens normalised parameters. They are related to Wyborne
# parameters by: Blm = l(lm) t(lm) Llm  where l(lm) is a l- and m-dependent
# factor and t(lm) is the Stevens operator equivalent factor, usually 
# denoted alpha, beta, gamma for l=2,4,6 respectively. Please note that
# by default the factor <J||alpha||J> etc. for a J-multiplet is used for
# f-electron ions whilst <L||alpha||L> etc. ignoring spin is used for 
# d-electron ions as per convention in the literature. You can force the
# use of one or the other type of factors by uncommenting the following:



# Note that the parameters output in the header of the ic1ion.out file
# is given in the type of the last input parameter read (e.g. Stevens
# normalisation in this case).

# magnetisation calculation can be triggered by the following parameters
# default magnetisation output units is Bohr magneton per ion

# if the following is uncommented, ic1ion will output magnetisation in emu/mol

# if the following is uncommented, ic1ion will output magnetisation in Am^2/mol

 xT   = 1
 xHa  = 0
 xHb  = 0
 xHc  = 0
 xmin = 1
 xstep= 1
 xmax = 300
 yT   = 0
 yHa  = 1
 yHb  = 0
 yHc  = 0
 ymin = 1
 ystep= 1
 ymax = 1

# the following parameters are read, when ic1ion is used as a module for other
# programs in the mcphase program suite, such as mcphas, mcdiff, mcdisp ...

# Neutron Scattering Length (10^-12 cm) (can be complex)
#  ... note: - if an occupancy other than 1.0 is needed, just reduce 
#              the scattering length linear accordingly
# Neutron Magnetic Form Factor coefficients - thanks to J Brown
#   d = 2*pi/Q      
#   s = 1/2/d = Q/4/pi   
#   sin(theta) = lambda * s
#   r= s*s = Q*Q/16/pi/pi
#   <j0(Qr)>=   FFj0A*EXP(-FFj0a*r) + FFj0B*EXP(-FFj0b*r) + FFj0C*EXP(-FFj0c*r) + FFj0D
#   <j2(Qr)>=r*(FFj2A*EXP(-FFj2a*r) + FFj2B*EXP(-FFj2b*r) + FFj2C*EXP(-FFj2c*r) + FFj2D
#   <j4(Qr)>=r*(FFj4A*EXP(-FFj4a*r) + FFj4B*EXP(-FFj4b*r) + FFj4C*EXP(-FFj4c*r) + FFj4D
#   <j6(Qr)>=r*(FFj6A*EXP(-FFj6a*r) + FFj6B*EXP(-FFj6b*r) + FFj6C*EXP(-FFj6c*r) + FFj6D
#   Dipole Approximation for Neutron Magnetic Formfactor:
#        -Spin Form Factor       FS(Q)=<j0(Q)>
#        -Angular Form Factor    FL(Q)=<j0(Q)>+<j2(Q)>
#        -Rare Earth Form Factor F(Q) =<j0(Q)>+<j2(Q)>*(2/gJ-1)
FFj0A=0.4085 FFj0a=23.8526 FFj0B=0.6091 FFj0b=8.2456 FFj0C=-0.1676 FFj0c=0.0415 FFj0D=0.1496     
FFj2A=3.4386 FFj2a=16.5303 FFj2B=1.9638 FFj2b=6.1415 FFj2C=0.2997 FFj2c=2.2669 FFj2D=0.0009     
FFj4A=-0.4759 FFj4a=14.0462 FFj4B=0.2747 FFj4b=3.7306 FFj4C=0.2458 FFj4c=1.2504 FFj4D=0.0057 

# radial wave function parameters (for use in programs pointc, chrgplt, charges) 
# R_Np,XIp(r)= r^(Np-1) . exp(-xi r) . (2 XIp)^(Np+0.5) / sqrt(2Np!)  
# values tabulated in clementi & roetti Atomic data and nuclear data tables 14 (1974) 177-478
# Co2+ is isoelectronic to Fe+, looking at page  422 of Clemente & Roetti 
# the 3D radial wave function is expanded as R(r)=sum_p C_p R_Np,XIp(r)
N1=3 XI1=4.95296 C1=0.36301 
N2=3 XI2=12.2963 C2=0.02707 
N3=3 XI3=7.03565 C3=0.14777
N4=3 XI4=2.74850 C4=0.49771 
N5=3 XI5=1.69027 C5=0.11388 
# for program pointcharge alternatively, the radial matrix elements may be given:

# Radial Matrix Elements (e.g. Abragam Bleaney 1971 p 399)
#<r^2> in units of a0^2 a0=0.5292 Angstroem
#<r^4> in units of a0^4 a0=0.5292 Angstroem
#<r^6> in units of a0^6 a0=0.5292 Angstroem

#pointcharges charge[|e|]  x[A] y[A] z[A]
pointcharge=  0.3            3    4    2

# Debye-Waller Factor: sqr(Intensity)~|sf|~EXP(-2 * DWF *s*s)=EXP (-W)
#                      with s=sin(theta)/lambda=Q/4pi
# relation to other notations: 2*DWF=Biso=8 pi^2 <u^2>
# unit of DWF is [A^2]
The output files are written to directory results and contain energy eigentstates (ic1ion.out) , polycrystal neutron transitions and intensities (ic1ion.trs) and, optional, the magnetisation (ic1ion.mag). Note in these output files the parameters alpha/beta/gamma in the header are not the Steven's operator equivalent factors, but rather the parameters of the configuration interactions. The operator equivalent factors are listed as <L||alpha||L> etc. or <J||alpha||J> etc. depending on whether $J$- or $L$-operator equivalents are used. $L$- ($J$-)operators are default for $d$- ($f$-)electron ions. In practice, as the $J$- or $L$- operators are valid only within an $L$- or $J$-manifold, whilst ic1ion calculates the splitting of the full $l^\nu$ configuration with many manifolds, the operator equivalent factors only serve to convert the input Stevens parameters into Wybourne normalisation, which is used internally by the program to calculate the Hamiltonian, for details on crystal field parameter conventions see appendix E.

Note that if you wish to compare the output of ic1ion with published energies/wavefunctions using Stevens parameters for rare earths, you should set the Slater F2, F4, F6 and spin-orbit zeta integrals to very high values - such as 10$^7$ meV, in order to obtain the $LS$-coupling limit. For comparison with Stevens parameters within an $L$-manifold for $d$-electron ions, the Slater integrals should be set to high values but the spin-orbit integral to zero. In this case the ic1ion output will have twice the degeneracy of a pure $L$-manifold calculation because the spin is taken into account in ic1ion and generates an additional two-fold degeneracy. If a non-zero (e.g. literature specified) spin-orbit integral zeta is specified, then this degeneracy will be lifted.

Martin Rotter 2017-01-10