Tag Archive: Bulk modulus

Jan 06

VASP tutor: Structure optimization through Equation-of-State fitting

Materials properties, such as the electronic structure, depend on the atomic structure of a material. For this reason it is important to optimize the atomic structure of the material you are investigating. Generally you want your system to be in the global ground state, which, for some systems, can be very hard to find. This can be due to large barriers between different conformers, making it easy to get stuck in a local minimum. However, a very shallow energy surface will be problematic as well, since optimization algorithms can get stuck wandering these plains forever, hopping between different local minima (Metal-Organic Frameworks (MOFs) and other porous materials like Covalent-Organic Frameworks and Zeolites are nice examples).

VASP, as well as other ab initio software, provides multiple settings and possibilities to perform structure optimization. Let’s give a small overview, which I also present in my general VASP introductory tutorial, in order of increasing workload on the user:

  1. Experimental Structure: This the most lazy option, as it entails just taking an experimentally obtained structure and not optimizing it at all. This should be avoided unless you have a very specific reason why you want to use specifically this geometry. (In this regard, Force-Field optimized structures fall into the same category.)
  2. Simple VASP Optimization: You can let VASP do the heavy lifting. There are several parameters which help with this task.
    1. IBRION = 1 (RMM-DIIS, good close to a minimum), 2 (conjugate gradient, safe for difficult problems, should always work), 3 (damped molecular dynamics, useful if you start from a bad initial guess) The IBRION tag determines how ions are moved during relaxation.
    2. ISIF = 2 (Ions only, fixed shape and volume), 4 (Ions and cell shape, fixed volume), 3 (ions, shape and volume relaxed) The ISIF tag determines how the stress tensor is calculated, and which degrees of freedom can change during a relaxation.
    3. ENCUT = max(ENMAX)x1.3  To reduce Pulay stresses, it is advised to increase the basis set to 1.3x the default value, which is the largest ENMAX value for the atoms used in your system.
  1. Volume Scan (Quick and dirty): For many systems, especially simple systems, the internal coordinates of the ions are often well represented in available structure files. The main parameter which needs optimization is the lattice parameter. This is also often the main change if different functional are used. In a quick and dirty volume scan, one performs a set of static calculations, only the volume of the cell is changed. The shape of the cell and the internal atom coordinates are kept fixed. Fitting a polynomial to the resulting Energy-Volume data can then be used to obtain the optimum volume. This option is mainly useful as an initial guess and should either be followed by option 2, or improved to option 4.
  2. Equation of state fitting to fixed volume optimized structures: This approach is the most accurate (and expensive) method. Because you make use of fixed volume optimizations (ISIF = 4), the errors due to Pulay stresses are removed. They are still present for each separate fixed volume calculation, but the equation of state fit will average out the basis-set incompleteness, as long as you take a large enough volume range: 5-10%. Note that the 5-10% volume range is generally true for small systems. In case of porous materials, like MOFs, ±4% can cover a large volume range of over 100 Å3. Below you can see a pseudo-code algorithm for this setup. Note that the relaxation part is split up in several consecutive relaxations. This is done to further reduce basis-set incompleteness errors. Although the cell volume does not change, the shape does, and the original sphere of G-vectors is transformed into an ellipse. At each restart this is corrected to again give a sphere of G-vectors. For many systems the effect may be very small, but this is not always the case, and it can be recognized as jumps in the energy going from one relaxation calculation to the next. The convergence is set the usual way for a relaxation in VASP (EDIFF and EDIFFG parameters) and a threshold in the number of ionic steps should be set as well (5-10 for normal systems is reasonable, while for porous/flexible materials you may prefer a higher value). There exist several possible equations-of-state which can be used for the fit of the E(V) data. The EOSfit option of HIVE-4 implements 3:
    1. Birch-Murnaghan third order isothermal equation of state
    2. Murnaghan equation of state
    3. Rose-Vinet equation of state (very well suited for (flexible) MOFs)

    Using the obtained equilibrium volume a final round of fixed volume relaxations should be done to get the fully optimized structure.

For (set of Volumes: equilibrium volume ±5%){
	Step 1          : Fixed Volume relaxation
	Step 2→n-1: Second and following fixed Volume relaxation (until a threshold is crossed and the structure is relaxed in fewer than N ionic steps) (IBRION = 2, ISIF=4, ENCUT = 1.3x ENMAX, ICHARG=1, LCHARG=.TRUE., NSW=100) 
	Step n : Static calculation (IBRION = -1, no ISIF parameter, ICHARG=1, ENCUT = 1.3x ENMAX, ICHARG=1, LCHARG=.TRUE., NSW=0) 
Fit Volume-Energy to Equation of State.
Fixed volume relaxation at equilibrium volume. (With continuations if too many ionic steps are required.) 
Static calculation at equilibrium volume
EOS-fitting Diamond and Graphite

Top-left: Volume scan of Diamond. Top-right: comparison of volume scan and equation of state fitting to fixed volume optimizations, showing the role of van der Waals interactions. Bottom: Inter-layer binding in graphite for different functionals.

Some examples

Let us start with a simple and well behaved system: Diamond. This material has a very simple internal structure. As a result, the internal coordinates should not be expected to change with reasonable volume variations. As such, a simple volume scan (option 3), will allow for a good estimate of the equilibrium volume. The obtained bulk modulus is off by about 2% which is very good.

Switching to graphite, makes things a lot more interesting. A simple volume scan gives an equilibrium volume which is a serious overestimation of the experimental volume (which is about 35 Å3), mainly due to the overestimation of the c-axis. The bulk modulus is calculated to be 233 GPa a factor 7 too large. Allowing the structure to relax at fixed volume changes the picture dramatically. The bulk modulus drops by 2 orders of magnitude (now it is about 24x too small) and the equilibrium volume becomes even larger. We are facing a serious problem for this system. The origin lies in the van der Waals interactions. These weak forces are not included in standard DFT, as a result, the distance between the graphene sheets in graphite is gravely overestimated. Luckily several schemes exist to include these van der Waals forces, the Grimme D3 corrections are one of them. Including these the correct behavior of graphite can be predicted using an equation of state fit to fixed volume optimizations.(Note that the energy curve was shifted upward to make the data-point at 41 Å3 coincide with that of the other calculations.) In this case the equilibrium volume is correctly estimated to be about 35 Å3 and the bulk modulus is 28.9 GPa, a mere 15% off from the experimental one, which is near perfect compared to the standard DFT values we had before.

In case of graphite, the simple volume scan approach can also be used for something else. As this approach is well suited to check the behaviour of 1 single internal parameter, we use it to investigate the inter-layer interaction. Keeping the a and b lattice vectors fixed, the c-lattice vector is scanned. Interestingly the LDA functional, which is known to overbind, finds the experimental lattice spacing, while both PBE and HSE06 overestimate it significantly. Introducing D3 corrections for these functionals fixes the problem, and give a stronger binding than LDA.

EOS-fitting for MIL53-MOFs

Comparison of a volume scan and an EOS-fit to fixed volume optimizations for a Metal-Organic Framework with MIL53/47 topology.

We just saw that for simple systems, the simple volume scan can already be too simple. For more complex systems like MOFs, similar problems can be seen. The simple volume scan, as for graphite gives a too sharp potential (with a very large bulk modulus). In addition, internal reordering of the atoms gives rise to very large changes in the energy, and the equilibrium volume can move quite a lot. It even depends on the spin-configuration.

In conclusion: the safest way to get a good equilibrium volume is unfortunately also the most expensive way. By means of an equation of state fit to a set of fixed volume structure optimizations the ground state (experimental) equilibrium volume can be found. As a bonus, the bulk modulus is obtained as well.

Feb 01

Doping of CeO2 as a Tunable Buffer Layer for Coated Superconductors: A DFT Study of Mechanical and Electronic Properties

Authors: Danny E. P. Vanpoucke,
Journal: Developments in Strategic Ceramic Materials:
Ceramic Engineering and Science Proceedings 36(8), 169-177 (2016)
(ICACC 2015 conference proceeding)
Editors: Waltraud M. Kriven, Jingyang Wang, Dongming Zhu,Thomas Fischer, Soshu Kirihara
ISBN: 978-1-119-21173-0
webpage: Wiley-VCH
export: bibtex
pdf: <preprint> 


In layered ceramic superconductor architectures, CeO2 buffer layers are known to form micro cracks during the fabrication process. To prevent this crack formation, doping of the CeO2 layer has been suggested. In this theoretical study, the influence of dopants (both tetravalent and aliovalent) on the mechanical and structural properties of CeO2 is investigated by means of density functional theory. Group IVa and IVb dopants show clearly distinct stability, with the former favouring interface and surface doping, while the latter favour uniform bulk doping. This behaviour is linked to the dopant electronic structure. The presence of charge compensating vacancies is shown to complicate the mechanical and structural picture for aliovalent dopants. We find that the vacancies often counteract the dopant modifications of the host material. In contrast, all dopants show an inverse relation between the bulk modulus and thermal expansion coefficient, independent of their valency and the presence of oxygen vacancies. Based on the study of these idealized systems, new dopants are suggested for applications.

Oct 28

Quasi-1D physics in metal-organic frameworks: MIL-47(V) from first principles

Authors: Danny E. P. Vanpoucke, Jan W. Jaeken, Stijn De Baerdemacker, Kurt Lejaeghere
and Veronique Van Speybroeck
Journal: Beilstein J. Nanotechnol. 5, 1738-1748 (2014)
doi: 10.3762/bjnano.5.184
IF(2014): 2.670
export: bibtex
pdf: <Beilstein> (open access)
Graphical Abstract: (left) Spin density of anti-ferromagnetic MIL-47(V) with ferromagnetic chains. (right) Electronic band structure and density of states.
Graphical Abstract: The MIL-47(V) MOF has one unpaired electron per V site. As a result, different spin configurations are possible, several of which lead to an anti-ferromagnetic state. The spin density of an antiferromagnetic state, containing only ferromagnetic chains is shown on the left. On the right, the electronic band structure of the same system is presented.


The geometric and electronic structure of the MIL-47(V) metal-organic framework (MOF) is investigated by using ab initio density functional theory (DFT) calculations. Special focus is placed on the relation between the spin configuration and the properties of the MOF. The ground state is found to be antiferromagnetic, with an equilibrium volume of 1554.70 Å3. The transition pressure of the pressure-induced large-pore-to-narrow-pore phase transition is calculated to be 82 MPa and 124 MPa for systems with ferromagnetic and antiferromagnetic chains, respectively. For a mixed system, the transition pressure is found to be a weighted average of the ferromagnetic and antiferromagnetic transition pressures. Mapping DFT energies onto a simple-spin Hamiltonian shows both the intra- and inter-chain coupling to be antiferromagnetic, with the latter coupling constant being two orders of magnitude smaller than the former, suggesting the MIL-47(V) to present quasi-1D behavior. The electronic structure of the different spin configurations is investigated and it shows that the band gap position varies strongly with the spin configuration. The valence and conduction bands show a clear V d-character. In addition, these bands are flat in directions orthogonal to VO6 chains, while showing dispersion along the the direction of the VO6 chains, similar as for other quasi-1D materials.

Oct 28

Aliovalent Doping of CeO2: DFT study of oxidation state and vacancy effects

Authors: Danny E. P. Vanpoucke, Patrick Bultinck, Stefaan Cottenier, Veronique Van Speybroeck, and Isabel Van Driessche,
Journal: J. Mater. Chem. A 2(33), 13723-13737 (2014)
doi: 10.1039/C4TA02449D
IF(2014): 7.443
export: bibtex
pdf: <JMaterChemA> <arXiv>


The modification of CeO2 properties by means of aliovalent doping is investigated within the ab initio density functional theory framework. Lattice parameters, dopant atomic radii, bulk moduli and thermal expansion coefficients of fluorite type Ce1-xMxO2-y (with M = Mg, V, Co, Cu, Zn, Nb, Ba, La, Sm, Gd, Yb, and Bi) are presented for 0.00 ≤ x ≤ 0.25. The relative stability of the doped systems is discussed, and the influence of oxygen vacancies is investigated. It is shown that oxygen vacancies tend to increase the lattice parameter, and strongly decrease the bulk modulus. Defect formation energies are correlated with calculated crystal radii and covalent radii of the dopants, and are shown to present no simple trend. The previously observed inverse relationship between the thermal expansion coefficient and the bulk modulus in group IV doped CeO2 [J. Am. Ceram. Soc. 97(1), 258 (2014)] is shown to persist independent of the inclusion of charge compensating vacancies.

Oct 28

Tetravalent Doping of CeO2: The impact of valence electron character on group IV dopant influence

Authors: Danny E. P. Vanpoucke, Stefaan Cottenier, Veronique Van Speybroeck, Isabel Van Driessche, and Patrick Bultinck
Journal: J. Am. Ceram. Soc. 97(1), 258-266 (2014)
doi: 10.1111/jace.12650
IF(2014): 2.610
export: bibtex
pdf: <J.Am.Ceram.Soc.> <arXiv>


Fluorite CeO2 doped with group IV elements is studied within the density functional theory (DFT) and DFT + U framework. Concentration-dependent formation energies are calculated for Ce1−xZxO2 (Z = C, Si, Ge, Sn, Pb, Ti, Zr, Hf) with 0 ≤ x ≤ 0.25 and a roughly decreasing trend with ionic radius is observed. The influence of the valence and near valence electronic configuration is discussed, indicating the importance of filled d and f shells near the Fermi level for all properties investigated. A clearly different behavior of group IVa and IVb dopants is observed: the former are more suitable for surface modifications and the latter are more suitable for bulk modifications. For the entire set of group IV dopants, there exists an inverse relation between the change, due to doping, of the bulk modulus, and the thermal expansion coefficients. Hirshfeld-I atomic charges show that charge-transfer effects due to doping are limited to the nearest-neighbor oxygen atoms.