תורת פונקציונל הצפיפות
תיאור מתמטי

הבעיה שאנו רוצים לפתור היא פתרון משוואת שרדינגר עבור גופים מרובים. כמקובל בבעיות אלו ההנחה הראשונה תהיה להתייחס לגרעיני האטומים כקבועים במקומם (לפי קירוב בורן-אופנהיימר) בזמן שרק האלקטרונים נעים. הנחה זאת מוצדקת בגלל מסתם הגדולה של הגרעינים יחסית לאלקטרונים. את האלקטרונים נתאר בעזרת פונקציית גל של הרבה אלקטרונים  \Psi ({\vec {r}}_{1},\dots ,{\vec {r}}_{N}). משוואת שרדינגר שתיארנו כאן נראית כך:
{\hat {H}}\Psi =\left[{\hat {T}}+{\hat {V}}+{\hat {U}}\right]\Psi =\left[\sum _{i}^{N}(-{\frac {\hbar ^{2}}{2m_{i}}}\nabla _{i}^{2})+\sum _{i}^{N}V({\vec {r}}_{i})+\sum _{i<j}^{N}U({\vec {r}}_{i},{\vec {r}}_{j})\right]\Psi =E\Psi
כאשר עבור מערכת עם  \ N אלקטרונים,  {\hat {H}} הוא ההמילטוניאן,  \ E היא האנרגיה הכללית,  {\hat {T}} אופרטור האנרגיה הקינטית,  {\hat {V}} אופרטור האנרגיה הפוטנציאלית של האלקטרונים כתוצאה מהשדה החשמלי שיוצרים הגרעינים הקבועים במקומם ו- {\hat {U}} אופרטור אנרגיית האינטראקציה בין האלקטרונים לעצמם. הבעיה המורכבת שתיארנו כאן לא ניתנת לפתרון אנליטי וגם לא ניתנת לפירוק לבעיות קטנות יותר של חלקיק בודד בגלל האינטראקציה בין האלקטרונים שמתוארת על ידי  {\hat {U}}.
ישנן שיטות רבות שמנסות לפתור משוואה זו על ידי פיתוח של פונקציית הגל  \Psi  בעזרת דטרמיננטת סלייטר. השיטה הפשוטה ביותר היא שיטת הרטרי-פוק וישנן גם שיטות מתקדמות יותר, אך הבעיה העיקרית שלהן היא המורכבות החישובית הגדולה שהן דורשות.
כאן באה לעזרה שיטת ה-DFT, מכיוון שהיא מאפשרת לייצג את הבעיה של הרבה חלקיקים שמתוארת על ידי  {\hat {U}} בעזרת פונקציה של חלקיק בודד ואז אין לנו צורך ב- {\hat {U}} בתיאור הבעיה. בשיטת ה-DFT משתנה המפתח הוא צפיפות האלקטרונים  n({\vec {r}}), שנתון ע"י:
n({\vec {r}})=N\int {\rm {d}}^{3}r_{2}\int {\rm {d}}^{3}r_{3}\cdots \int {\rm {d}}^{3}r_{N}\Psi ^{*}({\vec {r}},{\vec {r}}_{2},\dots ,{\vec {r}}_{N})\Psi ({\vec {r}},{\vec {r}}_{2},\dots ,{\vec {r}}_{N}).
הקשר הזה הוא הפיך ולכן ניתן לתאר בעקרון את פונקציית הגל במצב היסוד כפונקציונל של צפיפות האלקטרונים במצב היסוד:
\,\!\Psi _{0}=\Psi [n_{0}]
ומכאן ניתן לקבל את ערך התצפית במצב היסוד של כל אופרטור כללי  {\hat {O}} על ידי פונקציונל של  \,\!n_{0}:
O[n_{0}]=\left\langle \Psi [n_{0}]\left|{\hat {O}}\right|\Psi [n_{0}]\right\rangle .
במיוחד מעניינת אותנו אנרגיית מצב היסוד של המערכת שניתנת ע"י:
E_{0}=E[n_{0}]=\left\langle \Psi [n_{0}]\left|{\hat {T}}+{\hat {V}}+{\hat {U}}\right|\Psi [n_{0}]\right\rangle
כאשר התרומה של הפוטנציאל החיצוני  \left\langle \Psi [n_{0}]\left|{\hat {V}}\right|\Psi [n_{0}]\right\rangle  יכולה להירשם במפורש במונחים של צפיפות האלקטרונים במצב היסוד  \,\!n_{0}:
V[n_{0}]=\int V({\vec {r}})n_{0}({\vec {r}}){\rm {d}}^{3}r.
בצורה כללית יותר אנו יכולים לתאר את התרומה של הפוטנציאל החיצוני  \left\langle \Psi \left|{\hat {V}}\right|\Psi \right\rangle  במונחים של צפיפות האלקטרונים  \,\!n:
V[n]=\int V({\vec {r}})n({\vec {r}}){\rm {d}}^{3}r.
הפונקציונלים  \,\!T[n] ו- \,\!U[n] נקראים פונקציונלים אוניברסליים, בעוד ש-  \,\!V[n] נקרא פונקציונל לא אוניברסלי, מכיוון שהוא מתאים לבעיה ספציפית שנובעת מסידור הגרעינים במרחב. מכאן נובע שברגע שנתון לנו הפוטנציאל הספציפי בבעיה  {\hat {V}}, כדי למצוא את אנרגיית מצב היסוד עלינו למזער את הפונקציונל:
E[n]=T[n]+U[n]+\int V({\vec {r}})n({\vec {r}}){\rm {d}}^{3}r
כאשר המשתנה שלנו הוא  n({\vec {r}}) ואנו מניחים שיש לנו פונקציונלים אמינים ל  \,\!T[n] ו- \,\!U[n]. אם נצליח לפתור את בעיית המינימיזציה ולמצוא את אנרגיית מצב היסוד למעשה גם מצאנו את  \,\!n_{0} וכמו שראינו קודם בעזרת הצפיפות במצב היסוד נוכל למצוא את כל ערכי התצפית שמעניינים אותנו.
את בעיית המינימיזציה של  \,\!E[n] ניתן לפתור על ידי שימוש בשיטת כופלי לגראנז'. בשלב ראשון, נתייחס רק לפונקציונל האנרגיה שלא תלוי במפורשות באיבר האינטראקציה בין האלקטרונים:
E_{s}[n]=\left\langle \Psi _{s}[n]\left|{\hat {T}}+{\hat {V}}_{s}\right|\Psi _{s}[n]\right\rangle
כאשר  {\hat {T}} מסמן את אופרטור האנרגיה הקינטית,  {\hat {V}}_{s} זהו הפוטנציאל החיצוני (שנובע מהגרעינים) בו האלקטרונים נעים, כך ש  n_{s}({\vec {r}})\ {\stackrel {\mathrm {def} }{=}}\ n({\vec {r}}).
כך, ניתן לפתור את משוואות Kohn–Sham עבור מערכת העזר שהגדרנו, ללא האינטראקציה:
\left[-{\frac {\hbar ^{2}}{2m}}\nabla ^{2}+V_{s}({\vec {r}})\right]\phi _{i}({\vec {r}})=\epsilon _{i}\phi _{i}({\vec {r}})
שמניבה לנו את האורביטלים  \,\!\phi _{i} שמהם ניתן לקבל את הצפיפות  n({\vec {r}}) של הבעיה המקורית של הגופים המרובים.
n({\vec {r}})\ {\stackrel {\mathrm {def} }{=}}\ n_{s}({\vec {r}})=\sum _{i}^{N}\left|\phi _{i}({\vec {r}})\right|^{2}.
הפוטנציאל האפקטיבי של תיאור החלקיק הבודד ניתן לתיאור בצורה יותר מפורשת:
V_{s}({\vec {r}})=V({\vec {r}})+\int {\frac {e^{2}n_{s}({\vec {r}}\,')}{|{\vec {r}}-{\vec {r}}\,'|}}{\rm {d}}^{3}r'+V_{\rm {XC}}[n_{s}({\vec {r}})]
כאשר האיבר השני מייצג את איבר Hartree שמתאר את הדחייה החשמלית בין אלקטרון לאלקטרון והאיבר האחרון  \,\!V_{\rm {XC}} נקרא פוטנציאל ההחלפה-קורלציה. כאן  \,\!V_{\rm {XC}} כולל בתוכו את כל האינטראקציות בין הגופים.
מכיוון שאיבר Hartree ו- \,\!V_{\rm {XC}} תלויים בצפיפות  n({\vec {r}}) שתלויה באורביטלים  \,\!\phi _{i} שתלויים בפוטנציאל  V_{s}, הבעיה של פתרון משוואות Kohn-Sham צריכה להתבצע בצורה איטרטיבית. בדרך-כלל מתחילים מניחוש התחלתי ל- n({\vec {r}}), אז מחשבים את הפוטנציאל  V_{s} המתאים ופותרים את משוואות Kohn-Sham עבור  \,\!\phi _{i}. בעזרת האורביטלים שנמצאו אפשר לחשב את הצפיפות החדשה ולהתחיל את אותו התהליך מחדש. אפשר לחזור על התהליך הזה עד שמגיעים להתכנסות.