I. INTRODUCTION
Ni-based superalloys cannot satisfy the increasing demand for materials with higher-temperature capabilities. Instead, Co-based superalloys display superior hot corrosion resistance at high temperatures compared to Ni-based alloys. Reference Sims, Stoloff and Hagel1 However, the traditional Co-based superalloys strengthened by carbide precipitates have not been widely used because they exhibit a lower strength at high temperatures with respect to Ni-based alloys strengthened by the formation of γ′ Ni3(Al,Ta) phase. Reference Peng, Miura and Mishima2 So far, the γ′ phase has not been applied in Co-based commercial alloys. Interest in Co-based alloys has been fueled because of the discovery of the L12 ordered Co3(Al,W) precipitate. Reference Sato, Omori, Oikawa, Ohnuma, Kainuma and Ishida3 The Co–Al–W-base alloys with γ + γ′ (fcc + L12) structures exhibit higher yield strength at elevated temperatures compared to Ni-based superalloys. Reference Sato, Omori, Oikawa, Ohnuma, Kainuma and Ishida3,Reference Suzuki and Pollock4 In addition, Co3(Ge,W) and Co3(Ga,W) precipitates with L12 structure have also been found in Co–Ge–W and Co–Ga–W systems. Reference Chinen, Omori, Oikawa, Ohnuma, Kainuma and Ishida5,Reference Chinen, Sato, Omori, Oikawa, Ohnuma, Kainuma and Ishida6 These alloys show higher density than the Ni-based superalloys. However, Makineni et al. Reference Makineni, Nithin and Chattopadhyay7 discovered a new W-free L12 γ′ precipitate of (Co,Ni)3(Al,Mo,Nb) with a composition of Co-37.4 at.% Ni-10.8 at.% Al-8.1 at.% Mo-3.7 at.% Nb, and its density is comparable to Ni-based superalloys. In addition, this tungsten-free alloy has a higher specific 0.2% proof strength both at room temperature and at high temperature compared to the Co-9 at.% Al-9.8 at.% W alloy.
In the past decade, first-principles approaches have become an efficient tool for studying phase stability, mechanical, and thermodynamic properties of crystalline materials. Reference Yu, Cheng and Chen8–Reference Liu12 For example, these properties in L12 γ′-Co3(M,W) (M = Al, Ge, Ga) have been investigated using both first-principles calculations and experiments, Reference Joshi, Vamsi and Karthikeyan13–Reference Yao, Xing and Sun20 indicating that these γ′ precipitates are ductile and have a strengthening effect on the Co matrix. First-principles calculations also show that Ta plays a crucial role in stabilizing the γ′ phase for Co–Al–Mo alloys. Reference Makineni, Samanta, Rojhirunsakool, Alam, Nithin, Singh, Banerjee and Chattopadhyay21 In the present work, a first-principles study is hence performed to investigate those properties of the L12 (Co,Ni)3(Al,Mo,Nb) phase.
II. THEORETICAL METHODS
A. Crystal structure
For an A3B-type compound with the L12 structure, A atoms are located at the face centers, and B atoms are at the cube corners. In the present work, the 32-atom (Co4Ni4)3(Al4Mo3Nb1) supercell was used to mimic the Co-37.5 at.% Ni-12.5 at.% Al-9.375 at.% Mo-3.125 at.% Nb γ′ composition. Reference Makineni, Nithin and Chattopadhyay7 A special quasi-random structure (SQS) attempts to represent a fully random alloy by means of supercells of a few atoms to mimic the most relevant pair and multisite correlation functions of the solution. Reference Jiang18,Reference Zunger, Wei, Ferreira and Bernard22 The work of Rhein et al. Reference Rhein, Dodge, Chen, Titus, Pollock and Anton23 showed that the configurational entropy is ideal across Co3Al–Co3W, and the work from Koßmann et al. Reference Koßmann, Hammerschmidt, Maisel, Müller and Drautz24 showed that SQSs can sometimes accurately represent ground state structures in Co3(Al,X) compounds. The SQS was generated for the L12 (Co4Ni4)3(Al4Mo3Nb1) using the mcsqs code in Alloy Theoretic Automated Toolkit (ATAT) software. Reference Van de Walle, Tiwary, De Jong, Olmsted, Asta, Dick, Shin, Wang, Chen and Liu25,Reference van de Walle26 The generated 32-atom SQS for the L12 (Co,Ni)3(Al,Mo,Nb) compound is shown in Fig. 1. In addition, the 32-atom L12 SQS generated by Jiang was adopted for the Co3(Al0.5W0.5) compound Reference Jiang18 and was used for comparison to the (Co4Ni4)3(Al4Mo3Nb1) compound. The generated 64-atom SQS for the L12 (Co8Ni8)3(Al8Mo6Nb2) compound is shown in Supplementary Material Fig. S1. The atomic positions and correlation functions of both the 32-atom and 64-atom supercells are provided in the Supplementary Material.
32-Atom SQS for L12 (Co,Ni)3(Al,Mo,Nb).
B. Computational details
First-principles calculations of total energies, elastic constants, and electronic structures of the L12 γ′ precipitate were performed using the projector augmented wave (PAW) method as implemented in the Vienna ab initio simulation package (VASP). Reference Kresse and Joubert27,Reference Kresse and Furthmüller28 The generalized gradient approximation proposed by Perdew, Burke, and Ernzerhof (GGA-PBE) was adopted for exchange and correlation effects. Reference Perdew, Burke and Ernzerhof29 The k-point meshes of 5 × 5 × 5 and 10 × 10 × 10 were used for (Co,Ni)3(Al,Mo,Nb) and Ni3Al, respectively. The cutoff energy of the plane wave was set as 500 eV. Convergence was judged when the energy difference between each step of the ionic iterations was less than 4 × 10−6 eV/atom. The electronic density of state (eDOS) was obtained using a modified tetrahedron method. Reference Blöchl, Jepsen and Andersen30 Spin polarized calculations were included to account for the magnetic effects. To study the stability of L12 (Co,Ni)3(Al,Mo,Nb), the phonon density of states were calculated by YPHON code Reference Wang, Wang, Wang, Mei, Shang, Chen and Liu31 using force constants calculated by VASP.
C. Elastic properties
Elastic constants were calculated using an efficient stress–strain method. Reference Shang, Wang and Liu32 This methodology involves applying a set of strains (e = e 1, e 2, e 3, e 4, e 5, and e 6) on a crystal. Here, e 1, e 2, and e 3 refer to normal strains, and e 4, e 5, and e 6 to shear strains. The corresponding stresses (σ = σ1, σ2, σ3, σ4, σ5, and σ6) for the deformed crystal are obtained from first-principles calculations. According to the Hooke’s law, elastic stiffness constants are then calculated from the n set of strains and the corresponding stresses. More details are described in Ref. Reference Shang, Wang and Liu32.
The quasistatic approach was used to compute the isothermal elastic constants as a function of temperature. Reference Wang, Wang, Zhang, Manga, Shang, Chen and Liu33 This approach is based on the fact that the temperature dependence of elastic constants is caused mainly by the volume change of crystal at finite temperatures. This methodology can be summarized in three steps as follows: (a) calculating the elastic constants at 0 K as a function of volume; (b) predicting the volume change as a function of temperature; and (c) computing the temperature dependence of isothermal elastic constants. More details are described in Ref. Reference Wang, Wang, Zhang, Manga, Shang, Chen and Liu33. Generally, the measured elastic constants at high temperatures using the resonance method are isentropic, because the system is adiabatic due to the fact that the elastic waves travel much faster than the rate of heat conduction. For cubic crystals, the thermodynamic relations can be simplified as Reference Shang, Kim, Zacherl, Wang, Du and Liu34
$${{C_{11}^S} \over {C_{11}^T}} = {{C_{12}^S} \over {C_{12}^T}} = {{{C_{\rm{P}}}} \over {{C_{\rm{V}}}}}\;{\rm{and}}\;C_{44}^S = C_{44}^T\quad ,$$
where C
V and C
P are heat capacities at constant volume and constant pressure, respectively.
$C_{ij}^S$
and
$C_{ij}^T$
represent the isentropic and the isothermal elastic constants, respectively. The Voigt–Reuss–Hill approach was used to calculate the bulk (B), shear (G), and Young’s (E) moduli for polycrystalline crystals.
Reference Anderson35
The elastic anisotropy parameter is given by
$$A = {{2{C_{44}}} \over {{C_{11}} - {C_{12}}}}\quad .$$
The Young’s moduli along some orientations in the cubic crystal are given by Reference Gutkin, Ishizaki, Kuramoto and Ovid’ko36
$${E_{100}} = {{({C_{11}} - {C_{12}})({C_{11}} + 2{C_{12}})} \over {{C_{11}} + {C_{12}}}}$$
$${E_{110}} = {\left\{ {{{{C_{11}} + {C_{12}}} \over {({C_{11}} + 2{C_{12}})({C_{11}} - {C_{12}})}} + {1 \over 4}\left( {{1 \over {{C_{44}}}} - {2 \over {{C_{11}} - {C_{12}}}}} \right)} \right\}^{ - 1}}$$
$${E_{111}} = {\left\{ {{{{C_{11}} + {C_{12}}} \over {({C_{11}} + 2{C_{12}})({C_{11}} - {C_{12}})}} + {1 \over 3}\left( {{1 \over {{C_{44}}}} - {2 \over {{C_{11}} - {C_{12}}}}} \right)} \right\}^{ - 1}}\quad .$$
D. Thermodynamic properties
The Helmholtz free energy F(V,T) in a solid at volume V and temperature T can be calculated by the quasiharmonic approach, Reference Shang, Wang, Kim and Liu37
$$F\left( {V,T} \right) = {E_0}\left( V \right) + {F_{{\rm{vib}}}}\left( {V,T} \right) + {F_{{\rm{el}}}}(V,T)$$
where E 0(V) is the static energy at 0 K without the zero point vibrational energy, F vib(V,T) the vibrational free energy, and F el(V,T) the thermal electronic contribution to free energy. F vib(V,T) can be estimated by the empirical Debye–Grüneisen model, Reference Moruzzi, Janak and Schwarz38
$${F_{{\rm{vib}}}}\left( {V,T} \right) = {9 \over 8}{k_{\rm{B}}}{\Theta _{\rm{D}}} + {k_{\rm{B}}}T\left\{ {3\,\ln \left[ {1 - \exp \left( { - {{{\Theta _{\rm{D}}}} \over T}} \right)} \right] - D\left( {{{{\Theta _{\rm{D}}}} \over T}} \right)} \right\}$$
where k B is the Boltzmann’s constant, ΘD the volume-dependent Debye temperature, and D(x) the Debye function given by
$$D\left( x \right) = {3 \over {{x^3}}}\int_0^x {{{{z^3}\,{\rm{d}}z} \over {{e^z} - 1}}} \quad .$$
The volume-dependent Debye temperature ΘD is written as
$${\Theta _{\rm{D}}} = sAV_0^{{1 \over 6}}{\left( {{{{B_0}} \over M}} \right)^{{1 \over 2}}}{\left( {{{{V_0}} \over V}} \right)^\gamma }$$
where s is a scaling factor, A is a constant
${\left( {6{\pi ^2}} \right)^{1/3}}\hbar /{k_{\rm{B}}}$
, V
0 the ground state equilibrium volume, M the atomic mass, and B
0 the bulk modulus. The Grüneisen constant γ is defined by
$$\gamma = {1 \over 2}\left( {1 + {{B'}_0}} \right) - x$$
with x = 2/3 for high temperature case and x = 1 for low temperature case. In the present work, x = 2/3 was used.
Thermal electronic contribution to Helmholtz free energy was determined by Mermin statistics
${F_{{\rm{el}}}}\left( {V,T} \right) = {E_{{\rm{el}}}}\left( {V,T} \right) - T{S_{{\rm{el}}}}(V,T)$
, where E
el is the internal energy and S
el the bare electronic entropy.
Reference Shang, Wang, Kim and Liu37
III. RESULTS AND DISCUSSIONS
A. Structural properties
Calculations of total energies as a function of cell volume were carried out to estimate the lattice constants of (Co,Ni)3(Al,Mo,Nb) and Co3(Al,W) in the L12 structure. The equilibrium lattice constant was deduced from a 4-parameter Birch–Murnaghan equation of state fitting. Reference Shang, Wang, Kim and Liu37 The enthalpy of formation (ΔH) is defined as the enthalpy change of the compound with respect to its constituents in their standard element reference (SER) states. Calculated equilibrium lattice constants and ΔH values at 0 K are summarized in Table I, together with other predictions for the Co3(Al,W) phase. Reference Xu, Han, Wang, Wang, Liu and Liu39 It can be seen from Table I that the calculated lattice constant of Co3(Al,W) is 0.3567 nm, which agrees well with the experimental value of 0.3599 nm. Reference Sato, Omori, Oikawa, Ohnuma, Kainuma and Ishida3 The ΔH value is calculated to be −12.16 kJ/mol-atom for Co3(Al,W), agreeing reasonably well with the other theoretical result (−14.96 kJ/mol-atom). Reference Xu, Han, Wang, Wang, Liu and Liu39 The calculated equilibrium lattice constant of (Co,Ni)3(Al,Mo,Nb) is 0.3594 nm, which is slightly larger than that of Co3(Al,W). It is also seen that this value is in good agreement with the experimental value of 0.3587 nm. Reference Liu, Coakley, Seidman and Dunand40 The equilibrium lattice constant and ΔH value for the 64-atom (Co,Ni)3(Al,Mo,Nb) structure at 0 K were also calculated (see Supplementary Material Table SV).
Calculated equilibrium lattice constant a (nm) and enthalpy of formation ΔH (kJ/mol-atom) for (Co,Ni)3(Al,Mo,Nb) and Co3(Al,W).
a This work.
To understand the stability of L12 (Co,Ni)3(Al,Mo,Nb), the phonon density of states (pDOS) for L12 (Co,Ni)3(Al,Mo,Nb) was calculated in terms of the 32-atom supercell. Figure 2 shows the pDOS of the L12 (Co,Ni)3(Al,Mo,Nb), indicating that it is dynamically stable since imaginary phonon modes are absent.
Calculated total phonon density of states (pDOS) for L12 (Co,Ni)3(Al,Mo,Nb).
B. Electronic structures
The electronic density of states (eDOS) for L12 (Co,Ni)3(Al,Mo,Nb) was further calculated, which provided information on the physical basis of phase stability. Figure 3 shows the total and partial density of states of L12 (Co,Ni)3(Al,Mo,Nb). It is seen that the calculated eDOS presents a well-defined low-energy bonding and high-energy antibonding region. A deep valley near the Fermi level is a pseudo-gap, which separates the bonding states from the antibonding states. This pseudo-gap is mainly due to the d–d hybridization among the transition metal atoms and indicates the existence of covalent-like bonding in (Co,Ni)3(Al,Mo,Nb). Obviously, the phase stability of intermetallic compound depends on the location of the Fermi level. Reference Xu, Oguchi and Freeman41 The system will be stable if the Fermi level falls in the bonding states. Figure 3 shows that the Fermi level is located in the bonding region of the eDOS, which makes the L12 (Co,Ni)3(Al,Mo,Nb) phase stable. In addition, the non-zero densities around the Fermi level occur for (Co,Ni)3(Al,Mo,Nb), implying the thermal electronic contribution to free energy should be considered, especially at high temperatures.
Calculated total and partial electronic density of states (DOS) for L12 (Co,Ni)3(Al,Mo,Nb) with the dashed line representing the Fermi level.
C. Thermodynamic properties
To verify the methodology in the present work, Fig. 4 shows the predicted thermodynamic properties at finite temperatures compared with the experimental results of Ni3Al. Reference Touloukian, Kirby, Taylor and Desai42,Reference Ihsan43 A good agreement is observed between calculations and measurements for Ni3Al. Figure 5 illustrates the calculated thermodynamic properties of (Co,Ni)3(Al,Mo,Nb). Figure 5(a) shows that entropy increases with increasing temperature. In Fig. 5(b), it is shown that C p is proportional to T 3 at temperatures near 0 K. It is indicated in Fig. 5(c) that the thermal expansion coefficient increases rapidly with increasing temperature, then gradually increases at high temperatures for (Co,Ni)3(Al,Mo,Nb). Figure 5(d) shows that the formation enthalpy of (Co,Ni)3(Al,Mo,Nb) gradually decreases with the increasing of temperature.
Thermodynamic properties as a function of temperature for Ni3Al: (a) isobaric heat capacity at constant pressure (C p); (b) linear thermal expansion coefficient; and (c) entropy (S). The recommended experimental values (solid symbols) are also shown for Ni3Al. Reference Touloukian, Kirby, Taylor and Desai42,Reference Ihsan43
Temperature-dependent thermodynamic properties of L12 (Co,Ni)3(Al,Mo,Nb): (a) entropy (S); (b) isobaric heat capacity at constant pressure (C p); (c) linear thermal expansion coefficient; and (d) enthalpy of formation (ΔH).
D. Elastic properties
Calculated elastic properties of (Co,Ni)3(Al,Mo,Nb), Ni3Al, and Co3(Al,W) at 0 K are listed in Table II, together with other predictions for Co3(Al,W). Reference Wang and Wang44 Experimental elastic properties of Ni3Al, Reference Kayser and Stassis45 Co3(Al,W), Reference Tanaka, Ohashi, Kishida and Inui19 and Ni3(Al,Ta) Reference Tanaka, Ohashi, Kishida and Inui19 are also listed for comparison. For cubic crystals, lattice stability requires the elastic stiffness constants that satisfy the following conditions Reference Born and Huang46 :
$${C_{44}} > 0,\;{C_{11}} + 2{C_{12}} > 0,\;{C_{11}} - \left| {{C_{12}}} \right| > 0\quad .$$
Calculated and experimental elastic properties, including C ij , B, G and E, B/G, Cauchy pressure C 12–C 44, Young’s modulus along the [hkl] direction E hkl , E 110/E 100, E 111/E 100 and elastic anisotropy parameter A for (Co,Ni)3(Al,Mo,Nb), Ni3Al, and Co3(Al,W) at 0 K. The unit of elastic properties is GPa, and the properties of Ni3(Al,Ta) are listed for comparison.
Table II indicates that L12 (Co,Ni)3(Al,Mo,Nb) is elastically stable based on Eq. (11). A good agreement is observed between calculations and measurements, in particular for Ni3Al. Reference Kayser and Stassis45 All the shear and Young’s moduli of (Co,Ni)3(Al,Mo,Nb) are larger than those of Ni3Al and Ni3(Al,Ta), but smaller than those of Co3(Al,W). The calculated elastic anisotropy parameter A of (Co,Ni)3(Al,Mo,Nb) is smaller than those of Ni3Al, Co3(Al,W) and Ni3(Al,Ta). The ductile/brittle behaviors of materials have been proposed to be related to the B/G ratio. Reference Pugh47 If B/G < 1.75, the material behaves in a brittle manner, otherwise the material is predicted to be ductile. On the other hand, Pettifor Reference Pettifor48 proposed that the bonding character, reflecting the ductile/brittle behaviors of materials, could be described by the Cauchy pressure C 12–C 44. For metallic bonding, the Cauchy pressure is positive. However, the Cauchy pressure is negative for the directional bonding. It can be seen that the B/G ratio is higher than 1.75, and the Cauchy pressure is positive for (Co,Ni)3(Al,Mo,Nb). Therefore, (Co,Ni)3(Al,Mo,Nb) is predicted to be ductile in nature, and the ductility of (Co,Ni)3(Al,Mo,Nb) is predicted to be higher than that of Co3(Al,W). Calculated elastic properties of (Co,Ni)3(Al,Mo,Nb) based on the 64-atom SQS at 0 K are listed in Supplementary Material Table SVI. Although the 32-atom and 64-atom SQSs have different supercell sizes, the calculated elastic properties are very close to each other. Our test shows that elasticity is not sensitive to the supercell size.
The values of E 110/E 100 and E 111/E 100 are also given in Table II. Here E hkl denotes Young’s modulus along the crystal direction of [hkl]. The values of Young’s modulus are largest in the 〈111〉 direction and smallest in the 〈100〉 direction for (Co,Ni)3(Al,Mo,Nb), Ni3Al, Co3(Al,W), and Ni3(Al,Ta). For example, E 111 is larger than E 100 by 210 GPa for (Co,Ni)3(Al,Mo,Nb), implying a strong orientation dependence for Young’s modulus.
Figure 6 shows the calculated isothermal elastic constant
$C_{ij}^T$
and isentropic elastic constants
$C_{ij}^S$
at different temperatures for the Ni3Al phase. An agreement could be found for
$C_{12}^S$
and
$C_{44}^S$
between experiments
Reference Prikhodko, Yang, Ardell, Carnes and Isaak49
and the present calculations. For
$C_{11}^S$
, the calculated decreasing rate with increasing temperature is slightly smaller than the corresponding experimental measurements. The discrepancy found in
$C_{11}^S$
could be attributed to neglecting of the contributions of F
vib and F
el to these second derivatives.
Reference Wang, Wang, Zhang, Manga, Shang, Chen and Liu33
Calculated isothermal elastic constants (
$C_{11}^T$
: □,
$C_{12}^T$
: △, and
$C_{44}^T$
: ☆) and isentropic elastic constants (
$C_{11}^S$
: ○,
$C_{12}^S$
: ▽, and
$C_{44}^S$
: ☆) for Ni3Al at 300, 600, 900, and 1100 K in comparison with experimental data by Prikhodko et al.
Reference Prikhodko, Yang, Ardell, Carnes and Isaak49
(
$C_{11}^S$
: ●,
$C_{12}^S$
: ▼, and
$C_{44}^S$
: ★).
Figure 7 shows the calculated isentropic elastic constant
$C_{ij}^S$
, elastic properties (B, G, and E), Young’s modulus in the [hkl] direction E
hkl
, E
110/E
100, E
111/E
100, and the elastic anisotropy parameter A for L12 (Co,Ni)3(Al,Mo,Nb) at different temperatures. It can be seen that the isentropic elastic constants decrease monotonically with increasing temperature for the (Co,Ni)3(Al,Mo,Nb) phase. The Young’s modulus in the [111] direction decreases the fastest, while the [100] decreases the slowest with increasing temperature from 300 K to 1300 K. B, G, and E values of (Co,Ni)3(Al,Mo,Nb) decrease linearly with increasing temperature. However, the values of E
111/E
100, E
110/E
100, and A stay nearly constant with increasing temperature for (Co,Ni)3(Al,Mo,Nb). In the (Co,Ni)3(Al,Mo,Nb) phase, [100] is the softest direction, which helps to explain the experimentally observed fcc(γ)/L12(γ′) two-phase microstructure with L12 precipitates along 〈100〉 directions.
Reference Makineni, Nithin and Chattopadhyay7,Reference Jiang18
This two-phase structure is expected to lead to excellent mechanical properties at elevated temperatures.
Calculated isentropic elastic constants
$C_{ij}^S$
, elastic properties (i.e., B, G, and E), Young’s modulus in the [hkl] direction E
hkl
, E
110/E
100, E
111/E
100, and elastic anisotropy parameter A of L12 (Co,Ni)3(Al,Mo,Nb).
IV. CONCLUSIONS
First-principles calculations have been performed to investigate phase stability, electronic structure, elastic, and thermodynamic properties of the L12 (Co,Ni)3(Al,Mo,Nb) phase. A good agreement is found between calculated results and experimental data for Ni3Al, validating the present methodologies such as the Debye model and the quasistatic approach. Calculated results indicate that (Co,Ni)3(Al,Mo,Nb) is predicted to be mechanically and dynamically stable. (Co,Ni)3(Al,Mo,Nb) is predicted to be ductile in nature and possesses covalent-like bonding. Shear and Young’s moduli of (Co,Ni)3(Al,Mo,Nb) are larger than those of Ni3Al and Ni3(Al,Ta) but smaller than those of Co3(Al,W). The values of E 111/E 100, E 110/E 100, and A keep nearly constant with increasing temperature for (Co,Ni)3(Al,Mo,Nb). The temperature has significant effects on heat capacity, entropy, and linear thermal expansion coefficient of (Co,Ni)3(Al,Mo,Nb). It is expected that the present work provides data for the design of novel high-temperature tungsten-free Co-based alloys.
Supplementary Material
To view supplementary material for this article, please visit https://doi.org/10.1557/jmr.2017.8.
The atomic positions of both the 32-atom and 64-atom SQS are listed in Table SI and Table SII, respectively. The correlation functions of both the 32-atom and 64-atom SQS are listed in Table SIII and Table SIV, respectively.
ACKNOWLEDGMENTS
This work is funded by the National Natural Science Foundation of China through grant No. 51201074. Q. Yao acknowledges the support from Jiangsu Government Scholarship for Overseas Studies. The efforts at the Pennsylvania State University were partially supported by the National Science Foundation under grants DMR-1006557 and CMMI-1333999. First-principles calculations were carried out on the LION clusters at the Pennsylvania State University supported by the Materials Simulation Center and the Institute for CyberScience. Calculations were also carried out on the CyberStar cluster funded by the NSF through Grant No. OCI-0821527.
