ตอบคำถาม Band-gap DFT

โดย รศ.ดร.ทีปานิส ชาชิโย ภาควิชาฟิสิกส์ คณะวิทยาศาสตร์ มหาวิทยาลัยนเรศวร

มีคนชวนคุย Advance DFT เจาะลึกตัวทฤษฎีขั้นที่ 8 (ตรงกับบทที่ 8 จาก 10 บท ของตำรา DFT โดย Robert Parr และ Weitao Yang "Density-Functional Theory of Atoms and Molecules" Oxford U. Press) คำถามจาก คุณ Nontapat Wanwieng นศ.ฟิสิกส์ ม.เชียงใหม่

โดย Nontapat Wanwieng นักศึกษาฟิสิกส์ ม.เชียงใหม่

ข้อ 1) LDA เป็นการประมาณเเบบ local ในขณะที่ GGA เป็น semi-local คำว่า loแal เเละ semi-local ในที่นี้หมายความว่าอย่างไรครับ

ลองพิจารณา Exchange-Correlation Energy ของ LDA

\$ E_{XC}^{LDA} = \int {d^3 r\rho (\vec{r})\varepsilon _{xc} [\rho (\vec{r})]} \$

จะเห็นว่า เทอม \$\varepsilon _{xc}^{LDA} [\rho (\vec r)]\$ ขึ้นอยู่กับความหนาแน่นของอิเล็กตรอน ณ จุด \$ \vec{r} \$ ที่กำลังพิจารณานี้เท่านั้น เป็นอิสระจากจุดอื่น   นี้เป็นการประมาณแบบหนึ่ง เรียกว่า Local เพราะขึ้นอยู่กับจุดนี้ จุดเดียว

คุณอาจเห็นว่า ไม่แปลกอะไร? ทำไมต้องถือเป็นการประมาณ? ผมจึงขอยกตัวอย่างการคำนวณพลังงานผลักกันเชิง Coulomb ดังนี้

\$ E_J = \frac{1} {2}\iint {d^3 r_1 d^3 r_2 \frac{{\rho (\vec r_1 )\rho (\vec r_2 )}} {{\left| {\vec r_1 - \vec r_2 } \right|}}} = \frac{1} {2}\int {d^3 r} \rho (\vec r)v_J (\vec r) \$   เมื่อ   \$ v_J (\vec r) = \int {d^3 r'} \frac{{\rho (\vec r')}} {{\left| {\vec r - \vec r'} \right|}} \$

จะเห็นว่า พลังงานศักย์ไฟฟ้า \$ v_J (\vec r) \$ เป็นฟังก์ชันที่ขึ้นอยู่กับตำแหน่งอื่น ๆ ทั้งหมด   กล่าวคือ ต้องอินทิเกรตตลอดทั่วปริภูมิ จึงจะได้มาซึ่งค่าของ \$ v_J (\vec r) \$ ณ จุด \$ \vec{r} \$

ต่างจากกรณี \$ \varepsilon _{xc}^{LDA} [\rho (\vec r)] \$   ที่ขึ้นอยู่กับความหนาแน่น ณ ตำแหน่ง \$ \vec{r} \$ เพียงอย่างเดียว

คราวนี้มาถึง GGA

\$ \varepsilon _{xc}^{GGA} [\rho (\vec r),\nabla \rho ,\nabla ^2 \rho ] \$

ซึ่งเป็นฟังก์ชันของอนุพันธ์ด้วย !!! นี้ไม่ใช่ Local ไปเสียทีเดียว เพราะการหาอนุพันธ์ จะโยงความสัมพันธ์ของจุดข้างเคียงเข้ามาด้วย ยกตัวอย่างเช่นในกรณีของ Finite Difference

\$ \left. {\frac{{df}} {{dx}}} \right|_{@{\text{point}}\,\,i} \approx \frac{{f_{i + 1} - f_i }} {h} \approx \frac{{f_{i + 1} - f_{i - 1} }} {{2h}} \$

จะเห็นว่ามันมีทั้งการใช้จุดข้างๆ มาร่วมพิจารณา นี้ไม่ใช่ Local แต่ก็ไม่ใช่ Global (คือใช้ทุกจุด) เรียกว่า Semi-Local แล้วกัน !!! (Semi แปลว่า ก้ำกึ่ง)

ข้อ 2) ในการคำนวณ energy band gap ของ semiconductor หรือ insulator ด้วย LDA หรือ GGA ทำไมถึงให้ค่าที่ underestimated อย่างมาก (คลาดเคลื่อน 30-100%) เทียบกับ experimental value เสมอ **เคยอ่านผ่านๆว่าเป็นผลจาก self-interaction ,derivative discontinuity ขแง xc potential ซึ่งผมไม่เข้าใจครับ

นี้มีคำศัพท์เทคนิคหลายอัน ที่จะม้วน   เข้ามาตอบคำถามในที่สุด

Hartree-Fock ไม่มีปัญหา Self-Interaction

เป็นที่รู้กันว่า อิเล็กตรอนไม่ผลักตัวเอง (No Self-Interaction) หลักการอันนี้ สะท้อนออกมาให้เห็นตอนที่เราคำนวณพลังงานเฉลี่ย (Expectation Value) ของอิเล็กตรอนในอะตอมไฮโดรเจน คือ

\$ \left\langle \psi \right|\hat H\left| \psi \right\rangle = \underbrace { - \int {d\vec r\psi ^* (\vec r)\frac{1} {2}\nabla ^2 \psi (\vec r)} + }_{{\text{Kinetic}}}\underbrace {\int {d\vec r\psi ^* (\vec r)v(\vec r)\psi (\vec r)} }_{{\text{Potential}}} \$

จะเห็นว่า ไม่มีเทอม \$ + \frac{1} {2}\iint {d\vec r_1 d\vec r_2 }\left| {\psi (\vec r_1 )} \right|^2 \frac{1} {{\left| {\vec r_1 - \vec r_2 } \right|}}\left| {\psi (\vec r_2 )} \right|^2 \$ เพราะถ้ามี !!! ก็จะกลายเป็นว่ากลุ่มหมอกอิเล็กตรอนของไฮโดรเจน ผลักกันเอง ซึ่งผิดหลักฟิสิกส์

คราวนี้มาถึงระบบที่มีหลายอิเล็กตรอน   สมมุติว่าเราใช้ทฤษฎี Hartree-Fock ในการคำนวณ คือสร้างฟังก์ชันคลื่นให้เป็นแบบ Slater determinant แทนด้วยสัญลักษณ์ \$ \left| {\tilde \Psi (\vec x_1 ,\vec x_2 ,\vec x_3 , \cdots ,\vec x_N )} \right\rangle \$ เมื่อ \$ \vec{x} \$ คือสมบัติของอิเล็กตรอนทั้งในแง่ตำแหน่ง \$ \vec{r} \$ และสปิน \$ \vec{s} \$

หากไม่เข้าใจว่า Slater determinant คืออะไร โปรดศึกษาในออนไลน์คอร์ส Electronic Structure Theory [อ้างอิง 1]

จะได้พลังงานเฉลี่ยในกรณีนี้คือ

\$ \begin{gathered} E^{({\text{HF}})} = \sum\limits_i^{\text{occ}} {\int {d\vec r\psi _i^* (\vec r)\left[ { - \frac{1} {2}\nabla ^2 + v(\vec r)} \right]\psi _i (\vec r)} } \\ + \frac{1} {2}\sum\limits_{ij}^{{\text{occ}}} {\iint {d\vec x_1 d\vec x_2 \psi _i^* (\vec x_1 )\psi _i^{} (\vec x_1 )\frac{1} {{\left| {\vec r_1 - \vec r_2 } \right|}}\psi _j^* (\vec x_2 )\psi _j^{} (\vec x_2 )}} \\ - \frac{1} {2}\sum\limits_{ij}^{{\text{occ}}} {\iint {d\vec x_1 d\vec x_2 \psi _i^* (\vec x_1 )\psi _j^{} (\vec x_1 )\frac{1} {{\left| {\vec r_1 - \vec r_2 } \right|}}\psi _j^* (\vec x_2 )\psi _i^{} (\vec x_2 )}} \\ \end{gathered} \$

สมการ (1)

ฟังก์ชัน \$ \psi _i (\vec x) \$ มีชื่อเรียกว่า Orbital โดยอนุโลมคือแทนกลุ่มหมอกของอิเล็กตรอน 1 ตัว (ที่ว่าโดยอนุโลมเพราะจริง ๆ แล้วมันจะแยกออกมาเป็นอิเล็กตรอนตัวนี้ ตัวไหนไม่ได้ เพราะอิเล็กตรอนเหมือนกันหมดทุกตัว จนแยกไม่ออก)

ด้านขวามือของสมการ เทอมที่สอง และสาม ต่างกันตรงที่ดัชนี i และ j โดยทั้งสองเทอมมีชื่อเรียกว่า Coulomb และ Exchange ตามลำดับ

เทอมที่สองเรียกว่า Coulomb เพราะสามารถเขียนให้อยู่ในรูป กลุ่มหมอกอิเล็กตรอน 2 กลุ่มกำลังมีอันตรกิริยาผลักกันแบบคูลอมบ์ คือ กลุ่มที่ i^th: \$ \psi _i^* (\vec x_1 )\psi _i^{} (\vec x_1 ) = \left| {\psi _i } \right|^2 \$ และ กลุ่มที่ j^th: \$ \psi _j^* (\vec x_2 )\psi _j^{} (\vec x_2 ) = \left| {\psi _j } \right|^2 \$ จัดรูปได้ดังนี้

\$ \begin{gathered} + \frac{1} {2}\sum\limits_{ij}^{\text{occ}} {\iint {d\vec x_1 d\vec x_2 \underbrace {\psi _i^* (\vec x_1 )\psi _i^{} (\vec x_1 )}_{\left| {\psi _i } \right|^2 }\frac{1} {{\left| {\vec r_1 - \vec r_2 } \right|}}\underbrace {\psi _j^* (\vec x_2 )\psi _j^{} (\vec x_2 )}_{\left| {\psi _j } \right|^2 }}} \\ = \frac{1} {2}\iint {d\vec r_1 d\vec r_2 \underbrace {\left( {\sum\limits_i^{\text{occ}} {\left| {\psi _i } \right|^2 } } \right)}_{\rho (\vec r_1 )}\frac{1} {{\left| {\vec r_1 - \vec r_2 } \right|}}\underbrace {\left( {\sum\limits_j^{\text{occ}} {\left| {\psi _j } \right|^2 } } \right)}_{\rho (\vec r_2 )}} \\ = \frac{1} {2}\iint {d\vec r_1 d\vec r_2 \frac{{\rho (\vec r_1 )\rho (\vec r_2 )}} {{\left| {\vec r_1 - \vec r_2 } \right|}}} \\ \end{gathered} \$

ในบรรทัดที่สอง มี Self-Interaction ซ่อนอยู่ เพราะว่า ในการนับดัชนี มันมีโอกาสที่ i=j   นี้เป็นว่า กลุ่มหมอกของ Orbital นั้น ผลักกันเอง จะทำให้พลังงานของระบบสูงเกินจริง (แรงผลัก มีพลังงานศักย์เป็นบวก)

แต่ในทฤษฎี Hartree-Fock ยังมีเทอม Exchange เป็นเทอมที่สาม ในสมการ (1) คือ

\$ - \frac{1} {2}\sum\limits_{ij}^{\text{occ}} {\iint {d\vec x_1 d\vec x_2 \psi _i^* (\vec x_1 )\psi _j^{} (\vec x_1 )\frac{1} {{\left| {\vec r_1 - \vec r_2 } \right|}}\psi _j^* (\vec x_2 )\psi _i^{} (\vec x_2 )}} \$

สังเกตว่า ซัมเมชั่นครอบคลุมกรณี i=j ซึ่งถ้าเป็นอย่างนี้ จะทำให้เทอม Exchange มีค่าเท่ากับเทอม Coulomb พอดี เป็นการตัด Self-Interaction ออกโดยสมบูรณ์

คราวนี้มาถึงพลังงานของ DFT ที่เขียนขึ้นครั้งแรกด้วย Walter Kohn และ Liu Sham (KS-DFT) ว่า

\$ E^{{\text{(DFT)}}} = \sum\limits_i {\int {d\vec r\psi _i^* (\vec r)\left[ { - \frac{1} {2}\nabla ^2 + \hat v(\vec r)} \right]\psi _i (\vec r)} } + \frac{1} {2}\iint {d\vec r_1 d\vec r_2 \frac{{\rho (\vec r_1 )\rho (\vec r_2 )}} {{\left| {\vec r_1 - \vec r_2 } \right|}}} + E_{XC} \left[ \rho \right] \$

จะเห็นว่า ทางขวามือของสมการ สองเทอมแรกเหมือนกับกรณี Hartree-Fock ยกเว้นเทอมสุดท้าย   คราวนี้ ในกรณี DFT จะออกแบบพลังงาน Exchange-Correlation \$ E_{XC} \left[ \rho \right]\$ ยังไง? ที่จะลบ Self-Interaction ออกโดยสมบูรณ์ !?

ถ้าเป็น Hartree-Fock นี้ทำได้ เพราะเป็นการมองเชิง Orbital คิดทีละ Orbital จึงสามารถ ลบ Self-Interaction ของแต่ละ Orbital ได้ครบ ได้หมดจดสวยงาม

ถ้าเป็น DFT แล้วล่ะก็ ฟังก์ชัน \$ E_{XC}[\rho]\$ โดยหลักการแล้ว ต้องขึ้นอยู่กับตัวแปร \$ \rho(\vec{r})\$ เท่านั้น แยกออกมาทีละ Orbital ไม่ได้ มันจึงไม่ง่ายเลย ที่จะกำจัด Self-Interaction ออกไปในทฤษฎี DFT

Band-gap ประมาณด้วยผลต่าง HOMO-LUMO

(ก) ความหมายทางฟิสิกส์ของ Band-gap คือพลังงานที่ต้องใช้ เมื่ออิเล็กตรอนในระบบ(ที่มีอิเล็กตรอนอื่นๆอยู่ด้วยเป็นจำนวนมาก) กระโดดขึ้นไปในชั้นพลังงานสถานะกระตุ้น จากการทดลอง เมื่อนำผลึกมาวัด Absorption มันจะมีค่าถี่แสง(คลื่นแม่เหล็กไฟฟ้า)ค่าหนึ่ง ที่การดูดกลืนพุ่งสูงขึ้นรวดเร็วจนผิดสังเกต จุดนั้นคือ Band-gap Energy

(ข) ในการคำนวณ เรามองอิเล็กตรอนที่อยู่ในระบบ มีลักษณะเป็น Orbital ของใครของมัน ซึ่งก็จะมีพลังงาน \$ \epsilon_i \$ เป็นของตัวเอง ดังแสดงในภาพที่ 1 คือค่าพลังงานของแต่ละ Orbital ซึ่งเรียงสูงขึ้นไป

ภาพที่ 1 แสดงระดับพลังงานของ Orbital ซึ่งในระบบผลึกจะเป็นฟังก์ชันของ wave vector หรือ \$ \epsilon_i = \epsilon_i(\vec{k})\$ (เมื่อ \$ \vec{k} \$ แปรผันตรงกับ โมเมนตัมเฉลี่ยของอิเล็กตรอนในสถานะนั้นๆ, บางครั้งเรียกว่า Crystal Momentum) ตีความโดยอนุโลมว่า อิเล็กตรอนที่แม้อยู่ในชั้นพลังงานที่ \$ i^{\text{th}} \$ เดียวกัน แต่ถ้ามีโมเมนตัมพุ่งไปในทิศทางต่างๆกันในผลึก 3 มิติ ก็จะมีพลังงาน \$ \epsilon_i(\vec{k}) \$ แตกต่างกันด้วย

โดยทั่วไปในเซมิคอนดักเตอร์ ที่จำนวนอิเล็กตรอนเป็นเลขคู่ แต่ละ Orbital ก็จะมี 2 อิเล็กตรอนบรรจุอยู่ เติมลงในชั้นพลังงานต่ำสุดก่อน ไล่เรียงสูงขึ้นไป จนถึงชั้นพลังงานสูงที่สุด เรียกว่า HOMO ย่อมาจาก Highest-Occupied Molecular Orbital หรือ สูงที่สุด ที่มีอิเล็กตรอนบรรจุอยู่

ถัดขึ้นไปอีก เป็นชั้นว่างๆ ชั้นแรก หรือ LUMO คือ Lowest Un-occupied Molecular Orbital

อย่างหยาบที่สุดคือมองว่า อิเล็กตรอนจากชั้น HOMO ที่เคยอยู่เดิม กระโดดขึ้นไปอยู่ LUMO ทำให้ผลต่างของพลังงานเป็น \$ \varepsilon _{LUMO} - \varepsilon _{HOMO} \$ ดังแสดงในภาพที่ 1 ด้วยเส้นงูเลื้อยสีแดง

แต่ค่าพลังงาน (ข) จะตีความเป็น ปริมาณฟิสิกส์ (ก) โดยตรงไม่ได้ เป็นเพียงการประมาณ เพราะ ข.ไข่ เราพูดถึงพลังงานของ Orbital ในขณะที่ ก.ไก่ พูดถึงพลังงานของระบบทั้งหมด

โดยทั่วไปเมื่ออิเล็กตรอนกระโดด ไปอยู่ที่อื่น ย่อมเกิดช่องว่างของสถานะขึ้น อิเล็กตรอนตัวอื่นๆจะมีการปรับตัว มีการกระจายตัวต่างออกไป ทำให้พลังงานอื่น ๆ ที่เกี่ยวข้อง ทั้ง Coulomb ทั้งศักย์ไฟฟ้า ทั้ง Exchange มีการปรับตัวหมด ไม่ใช่คิดง่ายๆแค่ \$ \varepsilon _{LUMO} - \varepsilon _{HOMO} \$

อย่างไรก็ตาม Band-gap สามารถคำนวณออกมาด้วยทฤษฎีโดยตรง คือถ้าคิดให้ละเอียดจะยากขึ้น โปรดดูเปเปอร์ของ Perdew และ Levy (1983) https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.51.1884

Band-gap ของ DFT แคบกว่าของ Hartree-Fock

เพื่อจะเข้าใจสมบัติของพลังงาน \$\epsilon_i\$ (ซึ่งจะนำไปสู่ Band-gap) เราจะต้องไล่เรียงดูก่อนว่า \$\epsilon_i\$ มาจากไหน?

มันมาจากการแก้สมการ Kohn-Sham หรือสมการ Hartree-Fock แล้วแต่ทฤษฎีที่คุณเลือกใช้ ซึ่งทั้งสอง มีโครงสร้างของสมการเหมือนกันคือ

\$ \left[ { - \frac{1} {2}\nabla ^2 + \hat v(\vec r) + \hat v_J (\vec r) + \hat v_{xc} } \right]\psi _i = \varepsilon _i \psi _i \$

สมการ (2)

เทอม \$\hat v_J(\vec{r}) \$ คือ Coulomb Interaction ระหว่างอิเล็กตรอน ซึ่งมี Self-Interaction ปนอยู่ภายใน

ในกรณี Hartree-Fock \$\hat v_{xc}\$ ก็จะมีเฉพาะ Exchange แต่เป็น Exchange ที่หัก Self-Interaction จากเทอม Coulomb อย่างสมบูรณ์   ส่วนกรณี DFT ก็จะหักออกบ้าง แต่ไม่ 100%

ภาพที่ 2

ดังแสดงในภาพที่ 2   มาดู Orbital \$ \psi _i \in {\text{occ}} \$ เฉพาะที่เป็น Occupied คือมีอิเล็กตรอนบรรจุอยู่ จากสมการ (2) ด้านซ้ายมือ \$\psi_i\$ จะเข้าไปมีอันตรกิริยากับ \$\hat v_J(\vec{r})\$ ซึ่งภายใน มี Occupied Orbital ดักรออยู่ เกิดเป็น Self-Interaction ขึ้น ทำให้มีพลังงานเพิ่มขึ้นมา กรณี DFT ส่วนที่เพิ่มขึ้นนี้จะถูกหักล้างไปไม่หมด ทำให้พลังงาน \$\epsilon_i\$ ถูกดันให้สูงขึ้น เนื่องจาก Self-Interaction

\$\epsilon_i\$ ของ Occupied Orbital ถูกดันให้สูงขึ้น เนื่องจาก Self-Interaction

อธิบายมา 5 หน้าเพื่อให้เข้าใจแค่นี้

จากภาพที่ 2 เห็นว่า พลังงานของ Orbital กลุ่ม Non-Occupied หรือเรียกว่า Virtual จะไม่ถูกดันขึ้น เพราะอะไร?

เพราะ \$\hat v_J(\vec{r})\$ มีเฉพาะ Occupied ดักรออยู่ (ดูซัมเมชั่นของสมการ 1) แม้นว่า Orbital \$\psi _i \in {\text{virtual}}\$ เข้าไปเจอ ก็จะไม่พบกับ Orbital ตัวเองอยู่ภายใน \$\hat v_J(\vec{r})\$   จึงไม่เกิด Self-Interaction ขึ้น

เมื่อเป็นดังนี้ \$\epsilon_{LUMO}\$ จึงไม่ถูกดันขึ้น ส่งผลให้ภาพรวมแล้ว Band-gap ที่ทำนายด้วย DFT มีค่าหดเล็กลง เนื่องจาก Self-Interaction

ภาพที่ 3 อ้างอิงเปเปอร์ โดย Dovesi et.al (2000)

ในภาพมาจากเปเปอร์ “The Periodic Hartree-Fock Method and Its Implementation in the Crystal Code” โดย Dovesi et.al (2000) จะเห็นว่า Occupied Orbital ของ DFT (LDA) ถูกดันขึ้นอย่างแรง

ไม่ได้แปลว่า Band-gap ของ Hartree-Fock จะแม่นยำ

แม้ Hartree-Fock จะไม่มีปัญหา Self-Interaction แต่มันก็ขาดพลังงาน Correlation ไปอย่างสิ้นเชิง หากคุณลองศึกษา Moller-Plesset Perturbation Theory จะพบว่า ผลของ Correlation ก็คือการ Coupling กันระหว่าง Occupied และ Virtual Orbital ซึ่งการ Coupling นี้เองจะคล้ายทำให้ Virtual Orbital ถูกดึงลงเข้าหา Occupied หรือติดลบลงมา (พลังงาน Correlation มีเครื่องหมายติดลบ)

ดูจากกราฟของโปรแกรม Crystal จะเห็นว่า Virtual Orbital ของ Hartree-Fock (ปราศจาก Correlation) มีพลังงานสูงเกินไป ในขณะที่ LDA (ซึ่งมีพลังงาน Correlation เต็มรูปแบบ) ถูกดึงให้ต่ำลง

ด้วยเหตุนี้ ค่า Band-gap ที่ถูกต้อง จะอยู่กลางๆ ระหว่างของ Hartree-Fock และ DFT (LDA หรือ GGA บางอัน) ไม่จำเป็นนะ !!! ที่ว่าเป็น GGA แล้วจะมีปัญหา Self-Interaction เสมอไป ลองดูเปเปอร์ใน X-Archive ของ T.Chachiyo, H. Chachiyo [อ้างอิง 2, PDF v2 FIG 2] จะพบว่าให้ผลคำนวณใกล้กับ HF ที่สุดในกลุ่มของ LDA, GGA, และ meta-GGA (ความแม่นยำเฉลี่ยเมื่อประยุกต์กับอะตอมในตารางธาตุ)

นิยาม Band-gap ที่รัดกุม ในมุมของ DFT

นิยามของ Energy Gap ที่ใช้แพร่หลายใน DFT [อ้างอิง 3] ก็คือ

\$ G = I - A = [E_{N-1} - E_N)] - [E_N - E_{N+1}] \$

สมการ (3)

ซึ่งผมจะได้ขยายความ ดังต่อไปนี้

ในทางฟิสิกส์ เมื่ออิเล็กตรอนกระโดดออกจากชั้นเดิม ไปเติมลงในชั้นใหม่ (ที่ระดับพลังงานสูงขึ้น) สามารถแตกออกเป็น 2 จังหวะ คือ 1) หลุดออกจากชั้นเดิม นี้มีลักษณะเป็น Ionization Potential แปลว่า พลังงานที่ต้องใช้ ในการดึงมันให้หลุดออกมา   ในทางการคำนวณ Ionization Potentential คือผลต่างของพลังงานระหว่างสถานะทั้งสอง \$ I = (E_{N-1} - E_N)\$ เมื่อ \$ E_{N-1}\$ และ \$ E_{N}\$ คือพลังงานของระบบเมื่อมีจำนวนอิเล็กตรอน \$ N-1 \$ และ \$ N \$ ตัว ตามลำดับ

2) คือไป เติม   ลงชั้นใหม่ นี้มีคำศัพท์เฉพาะเรียกว่า Electron Affinity แปลว่า พลังงานที่คายออกมาเมื่อเพิ่มอิเล็กตรอน(เข้าไปหนึ่งตัว) คุณคงเดาได้ว่า คำนวณอย่างไร มันคือผลต่าง \$ A = (E_N - E_{N+1}) \$

เมื่อนำแนวคิด ของ Ionization Potential และ Electron Affinity มาต่อกัน หรือผนวกให้เป็นกระบวนการเดียว จึงได้เป็นสมการ (3)

สาเหตุที่ต้องเอา Ionization Potential และ Electron Affinity มา ลบ กัน (แทนที่จะบวกกัน) เพราะ ธรรมชาติการนิยามของ \$I\$ และ \$A\$ นั้นสลับกันตั้งแต่ต้น   \$I\$ นิยามเป็นพลังงานที่ต้องป้อนให้ระบบ ในขณะที่ \$A\$ เป็นพลังงานที่ระบบ คาย ออกมา อีกนัยหนึ่ง เครื่องหมายบวกลบ มันสลับกัน ดังนั้น เวลาคำนวณพลังงานที่ต้องป้อนสุทธิ ก็ต้องเป็น \$ I + (-A) \$ หรือ \$ I - A \$ ดังกล่าว

Derivative Discontinuity

ในทางคณิตศาสตร์ แปลว่าฟังก์ชันมีการหัก หรือหลอดกาแฟหัก ทำให้ความชันด้านซ้าย และด้านขวาไม่เท่ากัน เขียนเป็นสมการได้ว่า

\$ \mathop {\lim }\limits_{\eta \to 0^ + } \left\{ {\left. {\frac{{df}} {{dx}}} \right|_{x_0 + \eta } - \left. {\frac{{df}} {{dx}}} \right|_{x_0 - \eta } } \right\} \ne 0 \$

ข้างต้นเป็นการคำนวณความชัน ณ ตำแหน่ง \$ x_0 \$ โดยเทอมแรก เยื้องอยู่ฝั่งขวาเป็นระยะ \$ \eta \$ และเทอมสองอยู่ฝั่งซ้าย เมื่อให้ลิมิตเข้าใกล้ศูนย์ แปลว่าทั้งสองฝั่ง ตีขนาบเข้ามาใกล้ \$ x_0 \$ เรื่อยๆ ถ้าฟังก์ชันไม่หักงอ มันจะบรรจบกัน หรือเท่ากัน ทำให้ผลต่างเป็นศูนย์ แต่ ถ้า มี การ หัก งอ หรือ เกิด Derivative Discontinuity ขึ้น ผลต่างจะไม่เป็นศูนย์

ในบริบทของ Band-gap นั้น โยง อยู่ กับ สมการ (3) แต่ผมจะขอสรุปทฤษฎีพื้นฐานของ DFT เสียก่อน จึงจะโยงเข้าหากันได้

Density Functional Theory มองว่าเมื่อระบบของอิเล็กตรอน จำนวน \$ N \$ ตัว ตกอยู่ภายใต้พลังงานศักย์ภายนอก \$ v(\vec{r})\$ มันจะมีการกระจายตัว ปรับเปลี่ยนตามสภาพภูมิประเทศ (ศักย์) ที่มันอาศัยอยู่ การกระจายตัวนี้แทนด้วยฟังก์ชันความหนาแน่นของอิเล็กตรอน \$ \rho(\vec{r}) \$

พลังงานศักย์ที่ว่านี้ เช่น ศักย์คูลอมบ์จากประจุบวกของนิวเคลียสที่ดึงดูดมันไว้ บางครั้งมาจากสนามไฟฟ้าภายนอกเมื่อคลื่นแม่เหล็กไฟฟ้าวิ่งผ่าน ซึ่งหากภูมิประเทศเปลี่ยนไป หรือ ศักย์ \$ v(\vec{r})\$ เปลี่ยนไป อิเล็กตรอนก็จะกระจายตัวด้วยฟังก์ชัน \$ \rho(\vec{r})\$ เปลี่ยนไปจากเดิม

หลักพื้นฐานของ DFT คือว่า ความหนาแน่นของอิเล็กตรอน \$ \rho(\vec{r})\$ เป็นปริมาณพื้นฐานที่สามารถคำนวณ พลังงานของระบบ เขียนเป็นแผนภาพได้ว่า

ทราบ จำนวนอิเล็กตรอน \$ N \$ และพลังงานศักย์ \$ v(\vec{r})\$ \$ \enspace \rightarrow \enspace \$ จะสามารถคำนวณ \$ \rho(\vec{r}) \$ และ \$ E \$

แผนภาพ (4)

จะเห็นว่า พลังงานของระบบ ขึ้นอยู่กับว่ามีอิเล็กตรอน กี่ตัว หรือ \$ E = E[N] \$ ดังนั้นโดยทั่วไปแล้ว \$ E_{N-1} \ne E_N \ne E_{N+1}\$

แน่นอนว่า \$ N \$ เป็นจำนวนเต็ม แต่สมมุติในทางคณิตศาสตร์ เราจินตนาการเล่นๆว่า มันเป็นเศษทศนิยม เรียกโดยทั่วไปว่า Fractional Occupation แปลว่า จำนวนอิเล็กตรอนในระบบ ไม่ใช่จำนวนเต็ม อาทิเช่น มีอยู่ (4+0.5) อิเล็กตรอน แล้วอย่างนี้ พลังงาน \$ E=? \$

เข้าใจได้ง่ายๆว่า   พลังงานของระบบที่มีอิเล็กตรอน \$4.5\$ ตัว ก็คงอยู่ระหว่าง \$ E_4 \$ และ \$ E_5 \$ ส่วนจะมีค่าเท่าใดนั้น? น่าจะเขียนเป็นฟังก์ชัน ชัดๆ ออกมาได้ยาก เพราะระบบมันซับซ้อน แต่มีผลลัพธ์ทางทฤษฎีที่พิสูจน์ได้ด้วยคณิตศาสตร์ที่น่าทึ่งมาก [อ้างอิง 4] โดย J. Perdew et.al ที่บอกว่า

พลังงาน \$ E_{4.5} = 0.5 E_4 + 0.5 E_5\$ กล่าวคือ แปรผันเชิงเส้นตรง หรือ
\$E_{N+\eta} = (1-\eta) E_N + \eta E_{N+1} \$ เมื่อ \$ 0 \lt \eta \lt 1 \$

โปรดศึกษาเปเปอร์เพื่อดูบทพิสูจน์เมื่อมีเวลาว่าง ผลของทฤษฎีข้างต้น จะได้ว่า เส้นกราฟของพลังงาน \$ E_{N+\eta} \$ จะเป็นเส้นตรงที่เชื่อมระหว่าง \$ E_N \$ และ \$ E_{N+1} \$ ดังแสดงในภาพ

ภาพที่ 4 แสดงพลังงานของ Fractional Occupation ของอะตอม Be (Z=4) เฉพาะกรณี \$ N \$ เป็นจำนวนเต็ม คำนวณด้วยโปรแกรม Siam Quantum LDA (Slater exchange + Chachiyo correlation) Basis Set TZP ส่วนกรณี Fractional Occupation ผมลากเส้นตรงเชื่อมเอาเองเพื่อประกอบการอธิบาย

สิ่งที่เกิดขึ้นจริงในทางฟิสิกส์ ก็คงมีเฉพาะกรณี \$ N = 0, 1, 2, 3, 4, 5 , \cdots \$ แต่ในทางคณิตศาสตร์เราสมมุติได้ !!! และกราฟของพลังงาน เป็นเส้นตรงเชื่อมระหว่าง \$ N \$ ที่เป็นจำนวนเต็ม น่าทึ่งมาก !!! ที่มีทฤษฎีบทพิสูจน์ได้ว่าอย่างนี้

พิจารณาจุด \$ N=4 \$   จุดนี้   มี Derivative Discontinuity เกิดขึ้น เพราะกราฟมันหักงอ ความชันด้านขวา ด้านซ้าย ไม่ เท่า กัน   คราวนี้เราโยงเข้าหา Band-gap ยังจำกันได้ จากสมการ (3) ว่า \$ G = [E_{N-1} - E_N)] - [E_N - E_{N+1}] \$ หรือ จัดรูปได้ว่า

\$ G = [E_{N+1} - E_N] - [E_N - E_{N-1} ] \$

คราวนี้ลองพิจารณา อนุพันธ์ \$ \frac{dE}{dN}\$ ณ ซีกขวาของ \$ N = 4\$ อนุพันธ์ มันก็คือ ความชัน แต่เนื่องจากกราฟ เป็นเส้นตรง เราแค่เอาส่วนสูง หาร ฐาน ของสามเหลี่ยมได้เลย ดังนั้น ณ ซีกขวามือ \$ \frac{dE}{dN} = \frac{E_{N+1} - E_N}{\delta N}\$ แต่เรากำลังพูดถึง จำนวนนับของอิเล็กตรอนที่เพิ่มทีละหนึ่ง \$ \delta N = 1 \$ ดังนั้น \$ \frac{dE}{dN} = E_{N+1} - E_N\$

ใช้ตรรกะลักษณะนี้ กับกรณีซีกซ้ายของ \$ N = 4\$ เราเขียน Band-gap ได้ว่า

\$ G = \mathop {\lim }\limits_{\eta \to 0^ + } \left\{ {\left. {\frac{{\partial E}} {{\partial N}}} \right|_{N_0 + \eta } - \left. {\frac{{\partial E}} {{\partial N}}} \right|_{N_0 - \eta } } \right\} \$

หรือ Band-gap มีค่าเท่ากับ Derivative Discontinuity ของพลังงาน (เทียบกับการเปลี่ยนแปลงจำนวนอิเล็กตรอน) นั่นเอง

ทีนี้ เกี่ยวอะไรกับ "Derivative Discontinuity ของ XC Potential" ในคำถามข้อ 2?   ก็เพราะว่า พลังงานดังกล่าว ประกอบด้วยหลายเทอม ทั้งพลังงานจลน์ พลังงานาศักย์ พลังงานผลักกันเองระหว่างอิเล็กตรอน และที่สำคัญ ในทฤษฎี DFT คือ Exchange-Correlation นั่นเอง ดังนั้น Derivative Discontinuity ของ XC จึงมีผลต่อการคำนวณ Band-gap ด้วยประการฉะนี้

ข้อ 3) มี Exchange-correlation มากมายให้เลือกใช้ ทำไมต้องมีหลายแบบ มีแบบที่ "ดีที่สุด" หรือไม่ ขึ้นอยู่กับอะไร เเละจะเลือกใช้เเบบใดมีข้อพิจารณาอย่างไรครับ

ตอบเล่น ๆ คือ เลือกใช้ Chachiyo GGA ตอบจริง ๆ คือ

ต้นเหตุของการมี X-C เยอะ ๆ เพราะ ยังไม่มีใครหาสูตร "กรณีใด ๆ " เจอ   ดังนั้นก็ต้องอาศัยการประมาณ หรือตีกรอบเป็นกรณีเฉพาะ ไอ้การประมาณนี่แหละ เป็นบ่อเกิดแห่งความหลากหลาย

คือ ถ้าจะประมาณ มันก็ต้องมี "โมเดล" ในการคิดก่อน   ทำนองว่า นักฟิสิกส์จะวิเคราะห์การเคลื่อนที่ของแม่โคในอากาศ ก็อาจโมเดลว่า แม่โคเป็น ทรงกลมตันรัศมีหนื่ง /// หรือ Einstein โมเดลผลึกว่าเป็นของแข็งสั่นแบบสปริงเพื่อคำนวณความจุความร้อนเป็นต้น /// วกมาสูตร X-C ก็ต้องมีโมเดล   บางคนใช้การวิเคราะห์อิเล็กตรอนด้วย Perturbation บางคนใช้ Planewave+Quasi Particle หลัง ๆ มานี้ มีบางกลุ่มใช้ Neural Network หรือ Machine Learning (ตีพิมพ์ลงเนเจอร์ด้วยนะอันนี้)

โมเดลใคร โมเดลมัน มันก็เลยมีหลายสูตร ทีนี้   แต่ละสูตร มักมีตัวแปรที่ไม่ทราบค่า ลอยๆไว้   บางทีเขาเรียก parameters ที่จะต้องนำไป fit กับระบบที่ใช้เป็นกรณีศึกษา (สูตรชาชิโยไม่มี fitting parameters แต่ใช้เงื่อนไขที่ลู่เข้ากรณีที่ได้จากทฤษฎีสนามควอนตัม [อ้างอิง 6])

fitting parameters ก็ต้องเอามา fit กับระบบที่ใช้เป็นกรณีศึกษา ทีนี้แหละ เกิด ความ เฉพาะ เจาะจงมากขึ้น /// ถ้าคุณเอาไป fit กับโลหะ เอ้า สูตร XC คุณก็จะใช้ได้ เฉพาะกับโลหะ เอาไปใช้กับเซมิคอนดักเตอร์ก็จะไม่ดี /// หรือถ้าเอาไป fit กับ โมเลกุล เอ้า ไปใช้กับผลึกที่เป็น periodic ก็จะเสีย นี่แหละ สรุปแล้ว

สาเหตุหลักที่มีหลายสูตร เพราะ ต้องใช้โมเดลในการคิด ถ้าเริ่มจากโมเดลต่างกันจะได้สูตรที่มีโครงร่างต่างๆกัน   และเมื่อนำไป fit curve กับระบบ ตัวอย่าง   ต่างกัน ก็จะได้สูตรสุดท้ายที่ต่างออกไปอีก

วิธีการเลือกใช้ ให้ค้นดูด้วย keyword: band-gap DFT benchmark ในกูเกิล เพื่อค้นหาเปเปอร์ไว้ประกอบการเลือก XC ที่เหมาะกับระบบของคุณ   ยกตัวอย่างเช่น [อ้างอิง 5]

ภาพที่ 5 อ้างอิง Jochen Heyd et.al "Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional"

ด้วยรักและผูกพัน, ทีปานิส ชาชิโย

Keyword: ควอนตัม, Band-gap, Periodic Potential, DFT

อ้างอิง

• [1] ทีปานิส ชาชิโย, ออนไลน์คอร์ส Electronic Structure Theory https://sites.google.com/view/siamphysics/electronic-structure
• [2] T. Chachiyo, and H. Chachiyo. "Simple and accurate exchange energy for density functional theory." arXiv preprint arXiv:1706.01343 (2017). or peer-reviewed version https://doi.org/10.3390/molecules25153485
• [3] J. Perdew et.al "Understanding Band Gaps of Solids in Generalized Kohn-Sham Theory"
• [4] J. Perdew et.al "Density-Functional Theory for Fractional Particle Number: Derivative Discontinuities of the Energy"
• [5] The Journal of Chemical Physics 123(17):174101
• [6] อุดมศิลป์ ปิ่นสุข. "สมการมหัศจรรย์:สมการชาชิโย." Thai Journal of Physics 36.4 (2019): 100-106.