Virtual crystal approximation study of the complex refractory carbides based on Ti-Nb-Mo-V-C system with CASTEP computer code Pavlo Prysyazhnyuka, Roman Bishchaka, Sergiy Korniyb, Myroslav Panchuka and Volodymyr Kaspruk c a Ivano-Frankivsk National Technical University of Oil and Gas, Karpatska str., 15, Ivano-Franokivsk, 76019, Ukraine b Karpenko Physico-Mechanical Institute National Academy of Sciences of Ukraine, Address, Lviv, Index, Ukraine c Ternopil Ivan Puluj National Technical University, Ruska str., 56, Ternopil, 46001, Ukraine Abstract In the present study virtual crystal approximation with first principles CASTEP computer code was performed for investigation of multicomponent carbides in Ti-Nb-Mo-V-C system. The analysis of a computed properties data for fifteen carbides with equimolar metals components ratios showed, that elastic constants have significant positive deviation from rule of mixtures. Systematization of computed quantitative (elastic and lattice constants, hardness and fracture toughness) together with qualitative (charge density) data allowed us to distinguish ternary (Nb,Ti,Mo)C carbide as the most suitable compound for development wear-resistant materials and coatings. Keywords 1 Computing from first-principles, CASTEP code, multicomponent carbides, elastic constants, wear-resistance 1. Introduction Progress in engineering material science, structure-based materials and drug design highly depends on efficiency of implementation well revised theories and models into high-throughput computer codes. In the recent decades first-principles calculations based on density functional theory [1] has been the dominant technique for the quantum mechanical modelling of 3D periodic solids due to its high robustness, predictive power and reproducibility of results, when validating against experimental data. Computing properties of disordered solid solutions using so-called virtual crystal approximation (VCA) technique [2] within DFT framework is an effective approach to reach sufficient accuracy at relatively low computational cost, when compared to supercell [3], special quasirandom structure (SQS) [4] and effective cluster interaction (ECI) [5] methods. Similar approaches based on rectangular and hexagonal cells [6, 7] for modelling complex cyber-physical systems are used in analysis of medical and biological processes [8]. The general limitation of the VCA is the low accuracy in calculation of excess properties such as mixing enthalpy or volume for the interstitial solid solutions [9], which are characterized by significant deviation from Vegard’s law. However, mechanical properties, such as elastic constants, in many cases can be determined with high precision ITTAP’2021: 1nd International Workshop on Information Technologies: Theoretical and Applied Problems, November 16–18, 2021, Ternopil, Ukraine EMAIL: pavlo1752010@gmail.com (P. Prysyazhnyuk); bishchakr@gmail.com (R. Bishchak); kornii@ipm.lviv.ua (S. Korniy); myroslav_panchuk@ukr.net (M. Panchuk); kasprukv@tntu.edu.ua (V. Kaspruk). ORCID: 0000-0002-8325-3745 (A. 1); 0000-0001-6876-7142 (A. 2); 0000-0003-3998-2972 (A. 3); 0000-0002-4898-2707 (A. 4); ©️ 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) even for multicomponent systems like high entropy alloys (HEAs) and ceramics (HECs). S. Wang at el. [10] calculated the lattice constants (LC) using VCA implemented in Cambridge Serial Total Energy Package (CASTEP) [11] for the and SQS method for the FeCoNiCr – based HEAs and showed that the accuracy of both methods is very close to each other and to available experimental data, when calculated systems of elements with similar atomic radius and electronic configuration. Results of successful VCA computations using CASTEP code for HECs of Ti-Zr-Hf-V-Nb-Ta-C system, based on data analysis of computed properties, including heats of formations, elastic constants and melting points reported in [12]. It was found that TiHfVTaC quaternary carbide has the best combination of properties among fifteen investigated compounds. Experimental investigations of complex refractory compounds require special high temperature equipment, pure fine or ultrafine powders of initial components, which makes such type of experiments very expensive. Using first-principal computing codes for prediction properties of complex compounds allows to decrease number of experiments and determine most significant trends in composition – properties relationship, using data analysis information technologies. 2. Computational details and models The general concept underlying the VCA for a (A1-xBx)C-type compounds is the compositionally averaging of AC and BC parent phases pseudopotentials and obtaining resulting external potential (VVCA(r)) according to the formula : 𝑉𝑉𝐶𝐴 (𝑟) = (1 − 𝑥)𝑉𝐴𝐶 (𝑟) + 𝑥𝑉𝐵𝐶 (𝑟), (1) where x – composition. Within the DFT approach for a system with Nv valence electrons total energy (Etot) of a can be written as 1 1 𝑛 (𝑟)𝑛(𝑟 ′ ) (2) 𝐸𝑡𝑜𝑡 [{𝜙𝑖 }, {𝑅𝐼 }] = 𝑈({𝑅𝐼 }) + ∑ 𝜙𝑖 |− ∇2 + 𝑉𝑒𝑥𝑡 | + ∬ 𝑑𝑟𝑑𝑟 ′ + 𝐸𝑋𝐶 [𝑛], 2 2 |𝑟 − 𝑟 | ′ 𝑖 𝐼 where RI shows location of the lattice site I and 𝑉𝑝𝑠 is the corresponding pseudopotentials, 𝑛(𝑟) is the electron density, 𝐸𝑋𝐶 is the exchange-correlation energy, 𝑈({𝑅𝐼 }) is the energy of ions interaction and 𝑉𝑒𝑥𝑡 is the external pseudopotential, which in terms of VCA for a MC-type multicomponent carbides can be represented as 𝑉 (𝑟, 𝑟 ′ ) = ∑ ∑ 𝑤 𝐼 𝑉 𝑀 (𝑟 − 𝑅 , 𝑟 ′ − 𝑅 ), (3) 𝑒𝑥𝑡 𝑀 𝑝𝑠 𝐼𝑀 𝐼𝑀 𝐼 𝑀 𝑀 where 𝑉𝑝𝑠 is a pseudopotential for a given metal component M, which located on a lattice site I (fig. 𝐼 1) and 𝑤𝑀 represents its relative amount, so the total value on each site is equal to 1. Lattice models were constructed for MC-type carbides with FCC structures (spacegroup - Fm3m), where M was pure Ti, Nb, V, Mo and their combinations in equiatomic ratios, while other sites were filled with carbon. Figure 1: Crystal structure of the MC-type carbide with FCC lattice In the present study calculations were performed by VCA method, implemented CASTEP 19.11 code based on DFT. The equilibrium lattice geometry was computed by BFGS [13] algorithm within general gradient approximation (GGA) for the Perdew-Burke-Ernzerhof exchange-correlation functional [14]. The electron-ion interactions were treated by ultrasoft Vanderbilt pseudopotentials (C19 version) generated on the fly. The integration over Brillouin zone was performed using 8×8×8 Monkhorst-Pack [15] k-points mesh at the 520 eV plane wave energy cutoff. Convergence tolerance of self-consistent field (SCF) was set to 10-5 eV/atom, 3⨯10-2 eV/Å and 5⨯10-2 GPa for total energy, max ionic force and max stress, respectively. After the reaching of equilibrium geometry, “stress- strain” method at strain amplitude of 0.03 Å was used to obtain components of elastic tensor for further prediction of mechanical properties. 3. Results and discussion To cover all MC-type carbides of the Ti-Nb-Mo-V-C system with equimolar ratio of metal components 15 structures (four unary, six binary, four ternary and one quaternary carbide) were optimized (relaxed). Their lattice parameters are listed in table 1. The comparison calculated results against available experimental data reported in [16, 17, 18] for some unary stochiometric carbides shows that calculation error in all cases does not exceed 10 % (typical error of DFT calculations). This indicates good accuracy of structure relaxation with chosen parameters and allows to extend calculations for the mechanical properties determination using “stress-strain” method. Table 1 Lattice constants of the calculated structures Structure/a, Å Structure/a, Å Structure/a, Å NbC/4.4826 (Nb,V)C/4.2411 (Nb,Ti,V)C/4.2598 TiC/4.3291 (Nb,Mo)C/4.4111 (Nb,Ti,Mo)C/4.1647 VC/ 4.1363 (Ti,V)C/4.2038 (Nb,V,Mo)C/4.2727 MoC/4.3641 (Ti,Mo)C/4.1907 (Ti,V,Mo)C/4.2336 (Nb,Ti)C/4.2463 (V,Mo)C/4.2179 (Nb,Ti,V,Mo)C/4.2491 At the next stage, components Cij of a 6×6 stress tensor were computed for all structures in order to calculate elasticity characteristics: bulk modulus (B), shear modulus (G) and Young's modulus (E) by Voight-Reuss-Hill (VRH) averaging scheme using following equations: 1 2 (4) 𝐵 = (𝐶11 + 𝐶22 + 𝐶33 ) + (𝐶12 + 𝐶13 + 𝐶23 ), 9 9 1 1 (5) 𝐺= (𝐶11 + 𝐶22 + 𝐶33 − 𝐶12 − 𝐶13 − 𝐶23 ) + (𝐶11 + 𝐶22 + 𝐶33 ), 15 5 9𝐺𝐵 (6) 𝐸= . 𝐺 + 3𝐵 The calculated elastic properties are shown in fig. 2. As can be seen, elastic properties demonstrate significant deviations from rule of mixtures. For example, complex carbides (Nb,Ti)C and (Nb,Ti,Mo)C show extremely high values of all elastic moduli, which is much more higher than for pure unary carbide components as well as for calculated according to the rule of mixtures or statistical prediction [19], while their lattice parameters well agree with the Vegard’s law. In the general case, it can be seen that increasing of the alloying system complexity leads to the increasing of all elastic moduli. According to the Oganov and Mazhnick [20] there is a strong correlation between Young’s modulus, Poisson’s ratio (ν) and hardness, measured by Vickers method (HV), which is described by following formula: 0,585 (7) 9𝐸 (1 − 2𝜈)2 𝐻𝑉 = 2 ( ) − 3 , 8(1 + 𝜈)3 where Poisson’s ratio can be calculated for known B and G using (6) and following equations: 𝐸 𝐸 (8) 𝐵= , 𝐺= . 3(1 − 2𝜈) 2(1 + 𝜈) Another important materials characteristic, which can be derived using computed components of a stress tensor, is a fracture toughness represented through the critical tension stress intensity coefficient (KIc) as 1 1 − 3 (9) 𝐾𝐼𝑐 = 𝑓𝐸𝑁 𝛼0 2 𝑉06 [𝜁(𝜈)𝐸]2 where 𝑓𝐸𝑁 is the electronegativity factor, α0 is the constant depending on the chemical bonding characteristics in material (α0 is equal to 8840 GPa for covalent and ionic crystals), V0 is the volume per 1 atom and 𝜁 (𝜈) is the empirical fitting function of Poisson’s ratio. Figure 2: Computed elastic properties for MC-type carbides with different metal components ratio The Vickers hardness as well as fracture toughness, calculated by Eqs. (7) and (9), respectively, using utilities, which include machine learning approach, developed by USPEX laboratory [21] are shown in Fig. 3. As can be seen from figure, increasing overall trends for HV and KIc, when alloying system becomes more complex are very similar to the elastic moduli changes in the same systems. (Ti,Nb)C and (Nb,Ti,Mo)C also show the most favorable combination of properties in terms of hardness-toughness relationship, which is most desirable for many practical applications [22]. Moreover, despite high hardness of the (Nb,Ti,Mo)C its fracture toughness is highest among all investigated structures, which is unusual for carbide-like materials with high amount of covalent bonds. Such combination of properties makes this compound very good candidate for using as a component of impact wear resistant materials and coatings [23] for exploitation in extreme operating conditions, caused by high specific loads. Figure 3: Computed Vickers hardness and fracture toughness for MC-type carbides with different metal components combinations in equiatomic ratios Analysis of electron density data maps (EDM) (fig. 4) shows that areas with low electron density (highlighted in orange) become smaller, when number of mixed metal core atoms increasing. For the (Nb,Ti,Mo)C carbide isocharge lines indicating strength of interatomic bonding show that electron density between in metal-carbon (M-C) as well as in metal – metal (M-M) pairs is relatively low. This can be described as a formation covalent bonds provided by M-C couples and metallic bonding provided through M-M interactions. Other calculated structures are also characterized by mixture covalent-metal bonding, but its manifestation is relatively weaker. Weak interatomic bonding was detected for (Nb,Ti)C carbide, despite its mechanical properties is rather high. This discrepancy needs additional investigations, so it does not allow to classify this structure as an optimal. However, computed differences in electron density distributions gives only qualitative proof of the nonlinearity composition – properties relationships in given systems. Figure 4: Computed charge density maps of unary and complex carbides with high elastic modulus 4. Conclusions The effect of alloying of MC-type (FCC structure) refractory carbides with Nb, Ti, Mo and V in equimolar rations and all possible combinations on the elastic moduli, Vickers hardness and fracture toughness has been studied via first principles computing. CASTEP computer code with implemented virtual crystal approximation approach. The calculated values of lattice constants and elastic moduli, checked against literature experimental data for a known carbide structures of a given system are in good agreement, because the maximum error does not exceed 10 %. It was found, that ternary carbide of (Nb,Ti,Mo)C has the best combination of predicted mechanical characteristics, such as Young modulus (~ 650 GPa), Vickers hardness (~ 29 GPa) and fracture toughness (~ 6.5 MPa·m1/2), besides electron density map for this compound indicates presence strong covalent-metallic interatomic bonding, which corresponds to the high mechanical properties. The general limitation of the proposed approach is the rather low accuracy of the thermodynamic stability regions prediction. This requires additional calculation or experimental investigation of formation energies for selected compounds with promising mechanical properties. So, it can be expected that using complex (Nb,Ti,Mo)C carbides in pure state or as a reinforcing phases in composites and composite coatings can be useful for increasing their abrasion wear resistance. 5. References [1] Kohn, Walter, and Lu Jeu Sham. "Self-consistent equations including exchange and correlation effects." Physical review 140.4A (1965): A1133. [2] Bellaiche, L., and David Vanderbilt. "Virtual crystal approximation revisited: Application to dielectric and piezoelectric properties of perovskites." Physical Review B 61.12 (2000): 7877. [3] Prysyazhnyuk, Pavlo, et al. "Analysis of the Effects of Alloying With Si and Cr on the Properties of Manganese Austenite Based on Ab Initio Modelling." Eastern-European Journal of Enterprise Technologies 6.12 (2020): 108. doi: 10.15587/1729-4061.2020.217281 [4] Zunger, Alex, et al. "Special quasirandom structures." Physical review letters 65.3 (1990): 353. [5] Gonis, A., et al. "Configurational energies and effective cluster interactions in substitutionally disordered binary alloys." Physical Review B 36.9 (1987): 4630. [6] Martsenyuk, Vasiliy P. et al. “On Application of Latticed Differential Equations with a Delay for Immunosensor Modeling.” Journal of Automation and Information Sciences 50.6 (2018): 55–65. doi: 10.1615/jautomatinfscien.v50.i6.50 [7] Martsenyuk, V., A. Sverstiuk, and I. S. Gvozdetska. “Using Differential Equations with Time Delay on a Hexagonal Lattice for Modeling Immunosensors.” Cybernetics and Systems Analysis 55.4 (2019): 625–637. doi: 10.1007/s10559-019-00171-2 [8] Martsenyuk V., Klos-Witkowska A., Sverstiuk A., Bahrii-Zaiats O., Bernas M., Witos K. Intelligent Big Data system based on scientific machine learning of cyber-physical systems of medical and biological processes. The Fourth International Workshop on Computer Modeling and Intelligent Systems (CMIS-2021). 27 of April 2021, Zaporizhzhia, Ukraine. p. 34 - 48. [9] Levchuk, K. H., T. M. Radchenko, and V. A. Tatarenko. "High-Temperature Entropy Effects in Tetragonality of the Ordering Interstitial–Substitutional Solution Based on Body-Centred Tetragonal Metal." Metallofiz. Noveishie Tekhnologii 43 (2021): 1-26. doi: 10.15407/mfint.43.01.0001. [10] Wang, Shen, et al. "Comparison of two calculation models for high entropy alloys: Virtual crystal approximation and special quasi-random structure." Materials Letters 282 (2021): 128754. doi: 10.1016/j.matlet.2020.128754. [11] Clark, Stewart J., et al. "First principles methods using CASTEP." Zeitschrift für Kristallographie-Crystalline Materials 220.5-6 (2005): 567-570. [12] Liu, Shi-Yu et al. “Phase Stability, Mechanical Properties and Melting Points of High-Entropy Quaternary Metal Carbides from First-Principles.” Journal of the European Ceramic Society 41.13 (2021): 6267–6274. doi: 10.1016/j.jeurceramsoc.2021.05.022. [13] Perdew, John P., Kieron Burke, and Matthias Ernzerhof. "Generalized gradient approximation made simple." Physical review letters 77.18 (1996): 3865. [14] Pfrommer, Bernd G., et al. "Relaxation of crystals with the quasi-Newton method." Journal of Computational Physics 131.1 (1997): 233-240. [15] Monkhorst, Hendrik J., and James D. Pack. "Special points for Brillouin-zone integrations." Physical review B 13.12 (1976): 5188. [16] Lutsak, D. L., et al. "Formation of structure and properties of composite coatings TiB2-TiC-Steel obtained by overlapping of electric-arc surfacing and self-propagating higherature synthesis." Metallofizika i Noveishie Tekhnologii 38.9 (2016): 1265-1278. [17] Prysyazhnyuk, Pavlo et al. “Development of the Composite Material and Coatings Based on Niobium Carbide.” Eastern-European Journal of Enterprise Technologies 6.12 (96) (2018): 43– 49. doi:10.15587/1729-4061.2018.150807. [18] Kurlov, A. S., and A. I. Gusev. “Effect of Nonstoichiometry on the Lattice Constant of Cubic Vanadium Carbide VCy.” Physics of the Solid State 59.8 (2017): 1520–1525 doi:10.1134/s1063783417080133. [19] Trembach, Bohdan et al. “Application of Taguchi Method and ANOVA Analysis for Optimization of Process Parameters and Exothermic Addition (CuO-Al) Introduction in the Core Filler During Self-Shielded Flux-Cored Arc Welding.” The International Journal of Advanced Manufacturing Technology 114.3-4 (2021): 1099–1118. doi:10.1007/s00170-021-06869-y [20] Mazhnik, Efim, and Artem R. Oganov. "A model of hardness and fracture toughness of solids." Journal of Applied Physics 126.12 (2019): 125109. doi: 10.1063/1.5113622. [21] USPEX Tools and Utilities, 2021. URL: https://uspex-team.org/en/uspex/tools. [22] Ropyak, L. Ya, T. O. Pryhorovska, and K. H. Levchuk. "Analysis of materials and modern technologies for PDC drill bit manufacturing." Progress in Physics of Metals 21.2 (2020): 274- 301. doi: 10.15407/ufm.21.02.274 [23] Ivanov, Oleksandr, et al. "Improvement of abrasion resistance of production equipment wear parts by hardfacing with flux-cored wires containing boron carbide/metal powder reaction mixtures." Management Systems in Production Engineering 28 (2020): 178–183. doi: 10.2478/mspe-2020-0026