Herein, spin-orbit effect for describing heavy-element systems was included in composite methods using Kramers symmetry, and the isomerization pathway of TlH dimers was successfully predicted.
Recently, composite methods have been developed to obtain accurate thermochemical properties while retaining efficient computational costs, but most of them have focused on organic molecules. However, unlike light elements, the relativistic effect becomes significantly large when we describe heavy-element systems, and the electron density property on the atomic and molecular levels greatly changes. Therefore, the Dirac equation is used to describe heavy-element systems, but it is complex to use. Effective core potential (ECP) is one of the most successful tools for heavy-element systems because it contains a scalar relativistic effect while it has the non-relativistic form of the Hamiltonian. Computational efficiency is also expected because the number of electrons that must be described is reduced by the substitution of core electrons to potential.
A composite method was summarized that separates calculations to low level and obtains higher-level total energy by combining each contribution of the basis set effect.
The two-component theory-based composite method uses DFT/SO-DFT for geometry optimization, and 2nd-order Muller-Plesset/cc-pVDZ is used as a base energy. The total energy is obtained by addition of the polarization effect, diffuse function effect, an additional correlation correction term, and a spin-orbit correction term to base energy. The composite method is divided to a SP-set and a KRSP-set. One-component calculation is used for the SP-set, and an additional spin-orbit correction term is added, but two-component relativistic calculation is used for the KRSP-set, and no more spin-orbit correction term is required.
The 5th-, 6th-, and 7th-row elements of group 13, np1configuration were targeted to check the increased spin-orbit effect. In the total energy of the SP-set and KRSP-set, the spin-orbit correction term was dominant for the 6th and 7th row (In: 24%, Tl: 62%, E113: 46%), and polarization effect and diffuse function effect follow. The contribution of the additional electron correlation correction term was reduced going down the periodic table.
AREP level geometry optimization is recommended for 5th-row elements, but SOREP is strongly recommended for 7th-row elements. For 6th-row elements, the optimization level can be selected depending on the relative importance of accuracy and efficiency. Even for a closed-shell system, the spin-orbit correction term mainly contributes to the total energy, and describing the exact spin-orbit correction is important.
TCB2 is suggested for extending TCB1 to open-shell systems. Starting from the one-electron relativistic theory based SP-set of TCB1, Kramers unrestricted Hartree-Fock (KUHF) and Kramers unrestricted Moller-Plesset second-order perturbation (KUMP2) was tested for spin-orbit correction terms. First of all, the best combinations of methods and basis set for ionization potential (IP) for 6th-row elements was studied. cc-pVDZ, cc-pVTZ, aug-cc-pVDZ, aug-cc-pVTZ, uncontracted cc-pVDZ, and uncontracted cc-pVTZ were used. KUMP2/uncontracted cc-pVTZ matched with the experimental data with 2.16 kcal/mol errors. This result was applied to IP of Tl, and SO-CASSCF of the ccCA method showed 7.74kcal/mol of errors, while TCB2 showed only 2.03 kcal/mol.
The basis set of the energy-adjusted effective core potential should be uncontracted for exactly describing spin-orbit splitting; however, it increases computational costs. Therefore, a new approach is required for describing spin-orbit splitting with the required accuracy while keeping as many of the basis sets contracted as possible. Spin-orbit splittings of Tl atoms are analyzed by partially uncontracts p- and d- orbital parts. When two more p-parts were uncontracted, while others remained contracted, it showed a good match with reference data. However, d- orbital part uncontraction led to worse results.
The complete basis set limit (CBS limit) was tested for uncontracted basis sets. For the TlH molecule with 60 electron core ECP, CBS limits are fitted by Y=Y_0+ A/N for Hartree-Fock method and by Y=Y_0+ Be^cN for electron correlation methods. Herein, N has two cases, X, the cardinal number of the basis set and Lmax, and L, the maximum principal quantum number. Here, X3, X4, L3, and L4 formation was applied for above the equation, and X3 showed the best reliability to fitting equations. Based on this, the CBS limit is still available when basis sets for energy-adjusted ECP are uncontracted.
In this study, the isomerization pathway of TlH dimers was confirmed for the first time. HHTlTl, Tl(μ-H)TlH, HTlTlH, and Tl(μ-H)2Tl were obtained at the SO-B3LYP/aug-cc-pVTZ level. Each transition state structure was also verified by vibrational frequencies and IRC tests. Isomerization is continuously process from Tl(μ-H)2Tl (standard, 0kcal/mol) toHTlTlH (18.667 kcal/mol),Tl(μ-H)TlH (13.012 kcal/mol), and HHTlTl ( 20.307 kcal/mol). H atoms in bent HTlTlH structure are flipped through the linear transition structure and same isomerization is process from this flipped structure. Contrast to previous researches, they predicted that H atoms in HHTlTl structure move to inside and change to Tl(μ-H)2Tl, but Tl(μ-H)2Tlis only obtained from HTlTlH based on this research. Also second order saddle point(SOSP) works as a transition state of the path between Tl(μ-H)2Tl and HTlTlH.
Compared with Tl(μ-H)2Tl and Tl(μ-H)Tl, the dissociation energy of HTlTlH is lowered when the spin-orbit effect is included. As the Tl-Tl bond is only dominant for the TlH-TlH bond in HTlTlH, a pure σ-bond has a disadvantage because of the mixing between spinors that have different symmetry. For this case, a pure σ-bond consists of 1/3 p1/2 and 2/3 p3/2.
Fuzzy bond order analysis shows reliable results that HHTlTl has a tighter bond than a single bond; otherwise, HTlTlH, whose Tl is connected to H, has a weaker bond than a single bond. All bonds in Tl(μ-H)2Tl are in the same order, 0.538, which makes it hard to dissociate it to monomers.
본 연구에서는 6주기 이상의 중원소를 포함한 화합물을 묘사하는 경우 필요되는 스핀-오비트 효과를 크라머스 대칭을 사용하여 기존의 복합방법론에 포함하였으며 6주기 13족 원소인 수소화탈륨 이합체에서의 구조이성질체 변화 경로를 스핀-오비트 효과를 포함하여 예측하였다.
현대 화학에서 화합물의 열역학적 성질을 묘사하는 경우, 유효한 정확도를 유지하면서 효율적인 계산을 하고자 하는 목적으로 복합 방법론이 많이 개발되었으나 낮은 주기의 원소들을 포함하는 유기 화합물을 대상으로 하고 있다. 5,6 주기 이상의 무거운 원소들을 포함하는 경우, 핵의 전하량이 증가함에 따라 상대론 효과가 크게 나타나게 되고 이에 의해 원자 및 분자 수준에서의 전자 밀도가 크게 달라져 스핀-오비트 효과를 정확히 묘사하는 것이 필요하다. 이를 위해 디락 방정식의 해를 구하는 방법이 사용되고 있으나 수식의 복잡성으로 인해 쉽게 활용하기에는 어려움이 있다. 동일한 목적으로 개발된 유효중심 퍼텐셜의 경우, 비 상대론적 해밀토니안의 형태를 가지면서도 스칼라 상대론적 효과들이 포함되어 있을 뿐만 아니라 전자의 일부를 퍼텐셜로 치환함으로써 계산 효율성까지 도모할 수 있다.
복합 방법론은 상대적으로 낮은 수준의 계산을 통해 계산 방법과 기저함수의 기여도를 그 속성별로 나누어 더함으로 높은 수준의 계산에 해당하는 결과를 얻는 방법이라 요약할 수 있다. 상대론적 복합 방법론에서는 DFT/SO-DFT를 통하여 최적화된 구조를 얻고 2nd order Muller-Plesset 과 cc-pVDZ 기저함수를 가장 낮은 수준의 계산 값으로 삼아 이를 바탕으로 기저함수의 Polarization effect, Diffuse function effect, 그리고 계산방법론에서 추가적인 전자 상관효과를 추가적으로 더 포함함으로써 분자의 최종적인 에너지를 얻게된다. 복합 방법론은 대상 분자에 따라 SP-set과 KRSP-set 으로 나뉘게 되며 SP-set의 경우 AREP계산을 기본으로 하고 추가적인 스핀-오비트 보정항이 추가되며, KRSP-set의 경우 이-성분 상대론적 계산을 바탕으로 하므로 추가적인 스핀-오비트 보정항을 필요로 하지 않는다.
대상 분자로는 분자수준에서까지 상대론 효과를 크게 보이는 np1 원소인 13족을 대상으로 하였으며 효과의 차이를 보기 위하여 5,6,7주기 원소인 In, Tl, E113원소를 대상으로 하였다.
SP-set 과 KRSP-set 에서 각 효과들의 기여도를 살펴보면, 6,7주기의 경우 스핀-오비트 보정항이 가장 크게 기여를 하고 ( In: 24%, Tl: 62%, E113: 46%) 그 뒤로 기저함수 분극성 효과, 분산성 효과가 전체 에너지에 기여를 하며 추가적인 전자 상관 효과에 의한 보정은 주기가 아래로 내려갈 수록 기여도가 크게 감소하였다. TCB1에서 시작이 되는 구조는 5주기에서는 AREP수준에서의 계산이 권장되는 반면, 7주기에서는 SOREP 수준의 계산이 반드시 요구된다. 6주기의 경우 계산 효율성과 정확성 측면에서 선택적으로 결정할 수 있다. 닫힌 껍질 계에서도 스핀오비트 항의 기여도가 가장 높기 때문에 무거운 원소를 묘사하기 위해서는 정확한 수준의 스핀-오비트 효과가 반영되는 것이 중요하다.
닫힌 껍질 계를 묘사하는 TCB1 방법론의 열린 껍질 계로의 확장 버전으로 TCB2가 제안되었다. 일-성분 상대론적 계산을 기반으로 하는TCB1의 SP-set을 바탕으로하여 기존의 방법론들이 스핀-오비트 완전 활성 공간 자체 일관 성 장 이론 (SO-CASSCF)을 사용한 것과 달리 크라머스 제한되지 않은 하트리-폭 방법(KUHF), 크라머스 제한되지 않은 완전 활성 묄러-플레셋 2차 섭동론(KUMP2)을 바탕으로 하여 스핀-오비트 효과의 기여도를 살펴보았다. TCB1의 확장에 앞서, 6주기 원소의 이온화 에너지가 가장 잘 묘사되는 방법론과 기저함수의 조합을 살펴보았다. cc-pVDZ, cc-pVTZ, aug-cc-pVDZ, aug-cc-pVTZ, uncontracted cc-pVDZ, uncontracted cc-pVTZ 기저함수가 사용되었고, KUMP2/uncontracted-cc-pVTZ 조합이 실험치와의 편차를 2.16 kcal/mol정도로 보이며 가장 근사한 결과를 보여주었다. Tl과 Pb의 이온화 에너지 계산에 이를 적용해 본 결과, 기존의 SO-CASSCF방법으로 포함한 경우 Tl에서 실험치와 7.74kcal/mol의 차이를 보인 반면, TCB2는 2.03 kcal/mol로 훨씬 나은 효율을 보여주었다.
에너지 조정 유효중심퍼텐셜을 사용한 경우 스핀-오비트 갈라짐을 정확히 묘사하기 위해서는 기저함수의 수축을 풀어주는 것이 권장되나, 수축성이 제거된 기저함수의 경우 기저함수 행렬의 크기가 커지게 되어 적분 시간이 크게 소모된다. 따라서 최대한 기저함수를 수축된 상태로 유지하면서 스핀-오비트 갈라짐을 높은 정확도로 묘사하는 방법이 있을 경우 계산 시간 측면에서 큰 장점을 얻게 된다. 이를 위해6p1 전자배치를 갖는Tl원자의 스핀-오비트 갈라짐을 수축된 기저함수의 수축성을 p-,d-오비탈에 대해 부분적으로 제거해주며 살펴보았고 그 결과 p-오비탈의 수축성을 두 개 이상 추가로 풀어준 경우 기준 결과와 좋은 상관성을 보여주나 스핀-오비트 갈라짐을 분석하고자하는 p-오비탈이 아닌 외부의d-오비탈의 경우 수축성을 풀어주게 되면 오히려 p-오비탈 모사에 나쁜 영향을 끼쳐 기준 결과에서 벗어나는 경향성을 보여주었다. s-,f-오비탈 까지 전부 수축성을 제거한 경우, p-오비탈만 수축성을 제거한 것보다 기준치에서 벗어나는 결과를 보이는 것으로 보아 스핀-오비트 갈라짐을 묘사하고자 하는 오비탈에 대해서만 수축성을 제거해주는 것이 가장 효율적이고 정확한 결과를 보임을 알 수 있다.
마지막으로 완전한 기저 함수 한계값 (complete basis set limit; CBS limit)에 대한 수식을 수축되지 않은 기저 함수(uncontracted basis set)에 대해서도 살펴보았다. 60개의 핵전자를 치환한 스핀-오비트 유효 퍼텐셜을 가진 TlH분자에서, 하트리-폭 수준에서는 Y=Y_0+ A/N 으로, 전자 상관 방법론에서는 전체 에너지와 전자 상관 에너지 각각을 나눠 Y=Y_0+ Be^cN의 식에 대해 어느 정도 근사하게 묘사되는지 살펴보았다. 이 때, N은 주어진 기저함수에서의 최대 자기 양자수 값인 동시에 기저 함수의 기수 값인 X와 최대 주 양자수의 값을 의미하는 L (통상적으로 L값은 X+1과 일치한다) 두 변수에 대해 각각 두 가지로 나누어 X3, X4, L3, 그리고 L4 의 형태로 수식에 적용하였다. 각각을 위의 수식에 적용하여 CCSD(T)/aug-cc-pV5Z 계산 결과와 비교해 보았을 때, 수식과의 일치도 및 에너지 값의 타당성 면에서 X3 형태로 수식에 적용 되었을 때 가장 신뢰할만 한 결과를 보여주었다. 이를 바탕으로, 유효 중심 퍼텐셜을 사용한 경우 기저함수의 수축성을 제거한 경우에도 완전한 기저 함수 한계값을 통해 효율적으로 높은 수준의 계산 결과를 기대할 수 있음을 볼 수 있었다.
다음으로, TlH 이합체의 구조전환 메커니즘을 최초로 살펴보았다. HHTlTl, Tl(μ-H)TlH, HTlTlH, Tl(μ-H)2Tl 의 네 가지 안정한 상태를SO-B3LYP/aug-cc-pVTZ수준에서 얻었으며, 각 구조간의 전이 상태 구조를 확인, 최종적인 구조 전환 메커니즘을 획득하였다. 구조는 가장 안정한Tl(μ-H)2Tl (기준 상태, 0kcal/mol)를 기준으로 하여HTlTlH (18.667 kcal/mol), Tl(μ-H)TlH (13.012 kcal/mol), HHTlTl ( 20.307 kcal/mol) 순으로 전이가 일어나며 뒤틀린HTlTlH 구조의 경우 선형구조를 전이구조로 삼아 수소의 위치가 젖혀지며 거울 대칭된 구조를 기준으로 동일한 구조 전환이 진행된다. 기존의 결과와 대비하여 특이한 점은, 기존의 낮은 수준 계산에서는 HHTlTl구조에서 H원자의 위치가 안쪽으로 이동하며Tl(μ-H)2Tl구조로 전환될 것이라 예측한 반면, Tl(μ-H)2Tl구조는 HTlTlH구조에서 얻어짐을 볼 수 있었다. 또한, Tl(μ-H)2Tl에서 HTlTlH로 전환되는 경우, 2차 안장점 (second order saddle point)이 전이구조의 역할을 하게 된다.
Tl(μ-H)2Tl, Tl(μ-H)Tl, HTlTlH 세 가지 구조의 결합 에너지를 비교해 보았을 때, 앞의 두 구조와는 달리 HTlTlH의 경우, 스핀-오비트 효과가 포함되었을 때 결합 에너지가 더 낮아지는 결과가 나타났다. 이는 세 구조의 TlH-TlH간 결합을 살펴보았을 때, 앞의 두 가지 구조가 H원자를 매개체로 하여 결합이 일부 형성되는 반면, HTlTlH는 Tl간의 결합으로만 결합이 이루어지며 6p 오비탈에 전자가 하나 있는 Tl 의 경우 스핀-오비트 효과에 의해 스피너간의 섞임이 일어나 순수한σ 결합을 형성하기 위해서는 낮은 에너지를 가지는p1/2 만이 결합에 사용되는 것이 아니라1/3 p1/2 과2/3 p3/2 의 혼합으로 나타나게 되어σ 결합을 형성하는 경우 불리하게 작용하기 때문이다.
전자 위상 분석 기법을 사용하여, 퍼지 결합 차수 분석을 해 본 결과, 신뢰성 있는 결과를 얻을 수 있었다. HHTlTl에서의 Tl-Tl결합은 1차보다 약간 강한 상태로 결합 된 반면, 양쪽 Tl에 각각 수소가 붙어 있는 경우 1차 보다 약한 결합 상태를 보이고 있다. 또한 Tl(μ-H)2Tl 의 경우 네 개의 결합 세기가 동일하게 0.538로 나타나고 있어 TlH 단위체로 분해하기 어려울 것으로 예상된다.
그로스 전자 분포 분석을 통해 각 원소별 전하 차이를 AREP와 SOREP 수준에서 살펴보았고, HHTlTl구조의 경우 스핀-오비트 효과가 포함되면서 전자밀도의 쏠림 현상이 더 크게 나타났다. 이 경우, 6s와6pz가 H과의 결합에 기여하고6py는 Tl간의 결합에 관계되는 것으로 보인다. 반면, HTl(μ-H)Tl의 경우 스핀-오비트 효과가 포함되면서 전자밀도의 쏠림현상이 완화되는 것이 관찰되었다. HTlTlH 과 Tl(μ-H)2Tl 구조의 경우, 본 연구에서 진행한 기저함수 조합보다 작은 기저함수를 쓴 기존의 QCI수준의 계산결과와 비교해 보았을 때, 상대적으로 전자 밀도 쏠림이 더 적게 나타났다.
