VanVleckCalculator is a Python package for the calculation of the van Vleck distortion modes of octahedra, along with various other parameters including octahedral volume, angular shear, and various distortion parameters such as bond length distortion index, quadratic elongation, and the Van Vleck octahedral distortion modes.
It is developed by Liam Nagle-Cocco, a PhD student in physics at the University of Cambridge. For questions or suggestions, email [email protected].
- Python 3.8 or above
- numpy
- pymatgen
Download the van_vleck_calculator.py file, or better yet, clone the repository using Git or GitHub Desktop so you can pull any updates.
Then, in a Jupyter notebook (or whatever you use to write Python), write the following:
import sys
sys.path.insert(1,r"C:\Users\User\Documents\GitHub\VanVleckCalculator\code")
from van_vleck_calculator import *
where the path should be the path containing VanVleckCalculator's code. This will enable VanVleckCalculator to be accessed.
For more in-depth information, see the documentation or the code itself. See also worked examples.
To initialise an Octahedron object, it is necessary to first create a Pymatgen Structure object for your unit cell. The easiest way to do this is by importing a CIF. It is also possible to manually create a Pymatgen.core.structure object using known lattice parameters and atomic coordinates.
The following code can be used to generate a Pymatgen.core.structure object using a CIF file.
from pymatgen.core.structure import Structure
struc = Structure.from_file("cif_name.cif")
If you are working with a very large CIF file (such as from a molecular dynamics simulation or big box analysis of pair distribution function data), it is recommended that you set the optional frac_tolerance argument to be zero, to avoid shifting some atomic sites to high-symmetry positions. This requires Pymatgen version 2023.01.20 or later. For example:
from pymatgen.core.structure import Structure
struc = Structure.from_file("cif_name.cif",frac_tolerance=0)
This approach can also be used to import a VASP CONTCAR file, and several other crystal structure file types as documented on the Pymatgen website.
It is recommended you use a CIF, but if this is not feasible/desirable for your particular case, the following approach may be used. This example will generate a pymatgen.core.structure object for α-MnO2 in the I4/m space group.
from pymatgen.core.structure import Structure
from pymatgen.core.lattice import Lattice
latt = Lattice(matrix=([9.85,0,0],[0,9.85,0],[0,0,2.86]))
local_coords = (
[0.35049,0.16700,0],
[0.15137,0.19876,0],
[0.54139,0.16782,0]
)
struc = Structure(lattice=latt,species=["Mn","O","O"],coords=local_coords)
struc = struc.from_spacegroup(lattice=latt,species=["Mn","O","O"],coords=local_coords,sg="I4/m")
For more information on this, see the documentation for Pymatgen.
To generate an Octahedron object, it is necessary to tell the code where to consider the "centre" of your octahedron and which pymatgen.core.structure object to search in. There are other optional arguments, giving the ability to forbid certain species from being ligands or mandate ligands can only be a certain species, giving the range over which to search for ligands, and other options.
The only two mandatory arguments are a pymatgen.core.sites.PeriodicSite object (see here) and the pymatgen.core.structure object. Sites within a structure can be accessed by indexing, like an array. In the α-MnO2 example given previously, an Octahedron object can easily be created for the MnO6 octahedra as follows:
oct1 = Octahedron(struc[0],struc)
where struc[0] will refer to the 0th site in struc, which in this case will be a Mn site. Here, default values will be used for the range over which to search for ligands, and any site within this range will be considered a ligand. Often, the default values will work well enough. If we want to manually set a maximum distance over which to search for ligands, and exclude Mn ions from having other Mn ions as ligands, we could write the following:
oct1 = Octahedron(struc[0],struc,forbidden_ligands=["Mn"],ligands_max_distance=2.5)
Now, only atoms that are not Mn can be assigned as ligands to the oct1 object, and only atoms less than 2.5 Ångstroms will be considered. This is fine for α-MnO2, but what if we have a larger number of elements? We can force O to be the ligand as follows:
oct1 = Octahedron(struc[0],struc,possible_ligands=["O"])
For the sake of argument, if we had a system containing many different species, and we wanted O and Cl as possible ligands, we could use the following:
oct1 = Octahedron(struc[0],struc,possible_ligands=["O","Cl"])
Van Vleck defined a set of six modes based on the symmetry of the octahedra. These are labelled Q1 to Q6, where Q2 and 3 are the modes sensitive to Jahn-Teller distortion, and Q4 to Q6 are angular shear modes. To calculate these modes in the way van Vleck proposed in his original paper requires a Cartesian basis with the origin as the centre of the octahedron and the core-ligand bonds along the axes of the basis. Determination of the basis is complicated when there is angular distortion, as there is no way to orient the octahedron in the basis such that all core-ligand bonds are along axes. VanVleckCalculator is a tool to overcome this challenge and calculate the van Vleck distortion modes of an octahedron by rotating the octahedron to minimise the deviation of each ligand from the axes.
In addition to the van Vleck modes, this same rotation algorithm is used to calculate the angular shear, anti-shear, and angular shear fraction defined in the method paper accompanying this code package (under review).
For details on the application and interpretation of the modes, we refer the reader to the method paper accompanying this code package (under review) and the literature. This section will explain the details of how the algorithm works.
The first step, before performing the octahedral rotation, is to fix the origin. By default, this position is the atom in the centre of the octahedron. However, for some purposes this is not appropriate. For instance, if one is analysing the output of a big box Pair Distribution Function refinement or a Monte Carlo simulation, the central ion may exhibit some thermal motion from its ideal position. This will have a significant effect on the calculated Van Vleck modes, and so in this situation the optional argument should be used to fix the centre of the octahedron to be the crystallographic site. Likewise, the pseudo or second-order Jahn-Teller distortion may factor into the decision about where to situate the origin.
Below are examples of the use of the method calculate_van_vleck_distortion_modes for various possible origins. Here the origin is fixed as the position of the central cation (the default):
VanVleckModes = oct1.calculate_van_vleck_distortion_modes()
where VanVleckModes will by a Python list containing [Q1,Q2,Q3,Q4,Q5,Q6].
If you wish to fix the coordinates which are to be the origin, this is done as follows:
VanVleckModes = octahedron.calculate_van_vleck_distortion_modes(octahedral_centre=[x,y,x])
where x, y, and z are the absolute coordinates of the centre of the octahedron.
Additionally, it is also possible to fix the origin as the average position of the 6 ligands, by fixing the octahedral_centre argument as "average_ligand_position". An example of this:
VanVleckModes = octahedron.calculate_van_vleck_distortion_modes(octahedral_centre="average_ligand_position")
The mathematics of the Van Vleck modes assumes a perfect octahedron with no angular distortion, i.e. where all ligand-core-ligand angles are an integer number of 90°. In practice, coordination octahedra in crystal systems typically exhibit some angular distortion. For the sake of calculating the van Vleck modes, this gives rise to two possible approaches:
- Ignore the angular distortion and perform the calculation as if the vector from the centre of the octahedron to each ligand is the axis.
- Attempt to rotate the octahedron such that each core-ligand bond has an assigned axis, and the average angle between each core-ligand bond and the assigned axis is as low as possible. Then perform the calculation using the coordinate of each ligand in the axis to which the ligand is assigned.
VanVleckCalculator defaults to taking the second approach, using an algorithm which automatically rotates the octahedron to attempt to minimise the angle between each core-ligand bond and the crystallographic axes. However, there is the functionality to take both approaches.
This approach ignores the angular distortion entirely. It can be performed by supplying the argument ignore_angular_distortion to the method, for instance:
VanVleckModes = oct1.calculate_van_vleck_distortion_modes(ignore_angular_distortion=True)
This ignore_angular_distortion argument is compatible with the octahedral_centre argument for fixing the origin. It is also compatible with the specified_axes argument described in the next section. It is not compatible with the output_pairs argument.
Note that this approach will always give values of zero for the Q4, Q5, and Q6 modes.
In order to perform the Van Vleck calculation along orthogonal axes, the octahedron must be rotated such that the three axes of the octahedron correspond to the axes in space. In practice, this is not possible to achieve exactly if an octahedron has angular distortion, and so a match must be made as closely as possible. In VanVleckCalculator, there are three approaches to this problem:
- By default, the following approach is taken. In each of the xy, xz, and yz planes, the angle to rotate 4 (the 4 which should end up being within that plane) of the 6 atoms such that within that plane, all atoms are on an axis, is calculated, and the octahedron is then rotated by the average of the four angles. The result is a set of axes optimised for the shape of the octahedron. This approach may not be appropriate for analysis of a supercell from big-box PDF analysis or Monte-Carlo methods, because the random motion of the ligands may lead to some random rotation of the calculated axes of the octahedron relative to the crystallographic axes. If the algorithm fails, a warning will be displayed to the user.
- The second approach is achieved by manually inputting a set of vectors to the
calculate_van_vleck_distortion_modes()method using thespecified_axesargument, and fixing the argumentautomatic_rotation=False. By this approach, three orthogonal vectors are given, and CrystalPolyhedra will then rotate the octahedron such that these three vectors correspond to the x-, y-, and z-axes respectively. - The third approach is a combination of the first two. A set of axis are input using the
specified_axesargument, as in option 2, but the argumentautomatic_rotation=True. The code will rotate the octahedron such that these correspond to the axes in space. However, then the first option is used, but with a starting point closer to the actual axes.
Here the first method is used, i.e. the axes are automatically determined:
VanVleckModes = octahedron.calculate_van_vleck_distortion_modes()
Here, the second approach is taken:
VanVleckModes = octahedron.calculate_van_vleck_distortion_modes(
specified_axes=[ [1,1,0], [1,-1,0], [0,0,1] ],
automatic_rotation=False
)
where [ [1,1,0], [1,-1,0], [0,0,1] ] are presumed to be the axes of the octahedron.
Finally, option 3 is as follows:
best_guess_axes = [[1,-np.sqrt(2),1],[1,np.sqrt(2),1],[-1,0,1]]
VanVleckModes = octahedron.calculate_van_vleck_distortion_modes(specified_axes=best_guess_axes)
To check that this operation was successful, three methods can be used to check the position of ligands around the centre following the axis rotation. These three methods take all the same arguments as the methods which calculate the van Vleck modes.
visualise_sites_for_van_vleck()method can be used to produce a 3D plot showing the positions of the pairs (each shown in a different colour) around the centre of mass; to check the rotation worked, make sure that the ligand positions are as close as possible to the axes.output_sites_for_van_vleck()method can be used to return a Python list with shape (3,2,3) where the first axis contains three pairs of opposite ligands, the second axis refers to each site in a pair, and the third axis is the x-, y-, and z- positions of the site.output_and_visualise_sites_for_van_vleck()method performs both the above operations.
An alternative approach to checking the pair positions would be to use the output_pairs argument, which gives a second output identical to the output_sites_for_van_vleck() method. For example:
VanVleckModes, pair_positions = oct1.calculate_van_vleck_distortion_modes(output_pairs=True)
Finally, a rotational tolerance parameter can be set, which determines whether the rotation is judged to have been successful. The default value is 0.25. This can be set with the argument rotation_tolerance given to any of the van Vleck methods described in this section. The boolean argument omit_failed_rotations can be set to True if any calculations where the rotation is deemed to have failed (as defined by the rotation_tolerance) will return None for all parameters. This may be useful when analysing many octahedra within a supercell.
The calculate_van_vleck_distortion_modes() method calculates the 6 modes, and can be implemented as follows:
modes = oct1.calculate_van_vleck_distortion_modes()
The calculate_van_vleck_jahn_teller_params() method returns two parameters based on the magnitude √(Q2 + Q3) and the angle arctan(Q2/Q3) of the Jahn-Teller modes Q2 and Q3. This angle will be in the range 0 to 2π. This will return the angle in radians, but using an optional argument degrees=True will give an angle in radians. I.e.:
params = oct1.calculate_van_vleck_jahn_teller_params()
will give [mag,angle] where angle is in radians. Degrees can be obtained as follows:
params = oct1.calculate_van_vleck_jahn_teller_params(degrees=True)
Another Van Vleck calculation method is calculate_degenerate_Q3_van_vleck_modes() which returns a list containing [Q3,-0.5Q3+0.5√3 Q2,-0.5Q3-0.5√3 Q2], which are degenerate and equivalent modes for the three possible axes of elongation/compression.
The van Vleck modes were introduced in Van Vleck, J. H. "The Jahn‐Teller Effect and Crystalline Stark Splitting for Clusters of the Form XY6." The Journal of Chemical Physics 7.1 (1939): 72-84. The first known use of the approximation where angular distortion is ignored was in Kanamori, Junjiro. "Crystal distortion in magnetic compounds." Journal of Applied Physics 31.5 (1960): S14-S23.
Once an Octahedron object is initialised, various parameters can be calculated.
The volume of a tetrahedron is calculated from the distances between its vertices using the Cayley-Menger Determinant. For an Octahedron object, an octahedron is split into 8 tetrahedra and the volume of each of these is calculated. Any other polyhedron has its volume determined by splitting the surface into triangles, and calculating the volume of the tetrahedron bound by the vertices of the triangle and the centre of the polyhedron.
To obtain the volume, the calculate_volume() method should be used. This can be done as follows:
oct1.calculate_volume()
Note volume is given in Ångstroms cubed.
This approach to the calculation was inspired by the work in Swanson, Donald K., and R. C. Peterson. "Polyhedral volume calculations." The Canadian Mineralogist 18.2 (1980): 153-156.
Octahedral surface area is defined by splitting up the surface area into a set of triangles and summing the area of these triangles. It can be done using the calculate_surface_area() method. This can be done as follows:
oct1.calculate_surface_area()
where the returned value is in units of Ångstroms squared.
This approach to the calculation was inspired by the work in Swanson, Donald K., and R. C. Peterson. "Polyhedral volume calculations." The Canadian Mineralogist 18.2 (1980): 153-156.
The bond length distortion index parametrises the distortion in bond length from the centre of the octahedron. It takes the average of the differences between each bond length (from centre to vertex) and the average bond length. It is defined as follows:
where
It can be calculated using the calculate_bond_length_distortion_index() parameter as follows:
oct1.calculate_bond_length_distortion_index()
Another variation involves calculating this same parameter but using the average position of the 6 ligands as an alternative to using the central atom. This may be useful for instance for the second-order Jahn-Teller distortion, when the central atom is offset. This can be achieved using the octahedral_centre string argument as follows:
polyhedron.calculate_bond_length_distortion_index(octahedral_centre = "average_ligand_position")
It was first defined in Baur, W. H. "The geometry of polyhedral distortions. Predictive relationships for the phosphate group." Acta Crystallographica Section B: Structural Crystallography and Crystal Chemistry 30.5 (1974): 1195-1215.
Bond angle variance parameterises the degree of angular distortion from a perfect octahedron. The vertex-core-vertex angles in a perfect octahedron 90° respectively. In the presence of external constraints, these angles may be distorted. It is defined as follows:
where
It can be calculated using the calculate_bond_angle_variance() parameter as follows:
oct1.calculate_bond_angle_variance()
Output can be calculated in radians by using the optional argument degrees set to False, i.e. calculate_bond_angle_variance(degrees=False).
It was first defined in Baur, W. H. "The geometry of polyhedral distortions. Predictive relationships for the phosphate group." Acta Crystallographica Section B: Structural Crystallography and Crystal Chemistry 30.5 (1974): 1195-1215.
Quadratic elongation parameterises the elongation of a polyhedron. It is defined as follows:
where
It can be calculated using calculate_quadratic_elongation() as follows:
oct1.calculate_quadratic_elongation()
It was first defined in Robinson, Keith, G. V. Gibbs, and P. H. Ribbe. "Quadratic elongation: a quantitative measure of distortion in coordination polyhedra." Science 172.3983 (1971): 567-570
Effective coordination number parameterises the bond length distortion of a polyhedron in an equivalent manner to the bond length distortion index. It is defined as follows:
where,
It can be calculated using calculate_effective_coordination_number() as follows:
oct1.calculate_effective_coordination_number()
It was first defined by Rudolf Hoppe in a 1979 work, see later paper here: Hoppe, Rudolf, et al. "A new route to charge distributions in ionic solids." Journal of the Less Common Metals 156.1-2 (1989): 105-122.
For materials exhibiting a second-order Jahn-Teller distortion, it is useful to parameterise the degree of this distortion using the off-centering metric. Two examples of this are given in the literature.
An off-centering distance has been defined in the literature. This is given by the following equation:
It can be calculated in VanVleckCalculator as follows:
off_centre_dist = oct1.calculate_off_centering_distance()
An alternative metric was defined by PS Halasyamani in 2004. First, three angles
Given these values, the off-centering metric can then be calculated as follows:
It can be calculated in VanVleckCalculator as follows:
off_centre_param = oct1.calculate_off_centering_metric()
If you use this code in your work, please cite both this GitHub repository and the paper L. A. V. Nagle-Cocco and S. E. Dutton, Journal of Applied Crystallography (2024), 57, 1.