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

בתחום הכימיה החישובית, **אופטימיזציה גאומטרית** (באנגלית: **Molecular geometry optimization**) הוא תהליך לקביעת סידור אטומים במרחב כך שהכוח הנקי הפועל על כל אטום יתקרב לאפס, בהתאם למודל חישובי של קשרים כימיים, והנקודה במשטח האנרגיה הפוטנציאלית (PES) תהיה נקודת קיצון. קבוצת האטומים יכולה להיות מולקולה בודדת, יון, מצב מעובה, מצב מעבר או אפילו שילוב של אלה. לדוגמה, במודל חישובי מסוים ייתכן שהבסיס יהיה מכניקת הקוונטים.

לדוגמה, כאשר מנסים לבצע אופטימיזציה לגאומטריה של מולקולת מים, המטרה היא לקבוע את אורכי הקשרים בין המימן לחמצן ואת זווית חמצן-מימן-חמצן באופן שממזער את הכוחות שמושכים את האטומים אחד אל השני או דוחפים אותם זה מזה.

הסיבה העיקרית לביצוע אופטימיזציה גאומטרית נובעת מהמשמעות הפיזיקלית של המבנה המתקבל: המבנים שעברו אופטימיזציה לעיתים קרובות תואמים את המבנים שבהם החומר נמצא בטבע, והגאומטריה שלהם משמשת במחקרים ניסויים ותאורטיים במגוון תחומים כמו מבנה כימי, תרמודינמיקה, קינטיקה כימית, ספקטרוסקופיה ועוד.

ברוב המקרים, אך לא בהכרח תמיד, התהליך שואף למצוא סידור של האטומים שמייצג מינימום אנרגיה מקומי או גלובלי. לעיתים, במקום לחפש מינימום גלובלי, יעד האופטימיזציה יכול להיות מצב מעבר – נקודת אוכף על משטח האנרגיה הפוטנציאלית. בנוסף, במהלך התהליך ייתכן שקואורדינטות מסוימות (כגון אורך קשר כימי) יקבעו כקבועות במהלך האופטימיזציה.

## גאומטריה מולקולרית ופרשנות מתמטית
### תיאור המיקום של האטומים
ניתן לתאר את גאומטריית קבוצה של אטומים באמצעות וקטור המכיל את מיקומם. וקטור זה יכול לכלול את הקואורדינטות הקארטזיות של האטומים, או, במקרה של מולקולות, קואורדינטות פנימיות המורכבות מאורכי קשר, זוויות קשר וזוויות דיהדרליות.

אנרגיה כפונקציה של המיקומים:

בהינתן קבוצה של אטומים ווקטור r שמתאר את מיקומיהם, מגדירים את האנרגיה כפונקציה של המיקומים, E(r). כך, בעיית אופטימיזציה גאומטרית היא בעיית אופטימיזציה מתמטית שמטרתה למצוא את הווקטור r שבו E(r) נמצא במינימום מקומי. כלומר, אנו מחפשים נקודה בה:
* הנגזרת הראשונה של E(r) ביחס למיקומי האטומים,  \frac{\vartheta E (r)}{\vartheta r}, היא וקטור האפס.
* מטריצת הנגזרות השניות,  (\frac{\vartheta ^2E (r)}{\vartheta r_i \vartheta r_j})_{ij} הידועה גם כמטריצת הסיאן, המתארת את העקמומיות של משטח האנרגיה הפוטנציאלית ב-r, מכילה ערכים עצמיים חיוביים בלבד.

חיפוש מצבי מעבר:

מקרה מיוחד של אופטימיזציה גאומטרית הוא חיפוש אחר גאומטריית מצב מעבר – מצב שבו המערכת נמצאת בנקודת אוכף על משטח האנרגיה.

## מודלים חישוביים
המודל שמספק הערכה מקורבת ל-E(r) יכול להתבסס על: מכניקת הקוונטים (באמצעות תאוריית פונקציונליות הצפיפות או שיטות חצי-אמפיריות), שדות כוח, או שילוב של שיטות אלה (כמו במודלים QM/MM).

## תהליך האופטימיזציה
עם מודל חישובי מתאים וניחוש התחלתי לגאומטריה הנכונה, מתבצע תהליך אופטימיזציה איטרטיבי, לדוגמה:

# **חישוב הכוחות**: מחשבים את הכוח הפועל על כל אטום, כלומר,  \frac{\vartheta E (r)}{\vartheta r}.
# **בדיקת תנאי עצירה**: כאשר הכוח הפועל על האטומים קטן מערך סף שנקבע מראש, תהליך האופטימיזציה מסתיים.
# **עדכון המיקומים**: אם לא, מעדכנים את מיקום האטומים על ידי צעד חישובי  \vartriangle r שמיועד להקטין את הכוחות.
# **חזרה על התהליך**: מחזירים את התהליך להתחלה עד שהכוחות מתחת לסף המוגדר.

באופן זה, באמצעות סדרת חישובים ועדכונים, מתקרבים בהדרגה לגאומטריה האופטימלית שבה האנרגיה היא מינימלית והכוחות כמעט מאוזנים.

## היבטיים מעשיים של האופטימיזציה
כפי שתואר לעיל, ניתן להשתמש בשיטות כגון מכניקת הקוונטים כדי לחשב את האנרגיה, E(r), יחד עם הגרדיאנט של משטח האנרגיה הפוטנציאלית – כלומר, הנגזרת של האנרגיה ביחס למיקום האטומים,  \frac{\vartheta E}{\vartheta r} וכן את מטריצת הנגזרות השניות של המערכת  (\frac{\vartheta ^2E (r)}{\vartheta r_i \vartheta r_j})_{ij} הידועה גם כמטריצת הסיאן, המתארת את עקמומיות המשטח בנקודה r.

אלגוריתם לאופטימיזציה יכול לנצל חלק או את כל הערכים של  \frac{\vartheta E}{\vartheta r} E(r), ו-  (\frac{\vartheta ^2E (r)}{\vartheta r_i \vartheta r_j})_{ij} במטרה למזער את הכוחות הפועלים על האטומים. תאורטית, ניתן להשתמש בכל שיטת אופטימיזציה – כגון ירידת גרדיאנט, גרדיאנט מצומד או שיטת ניוטון – אך למעשה נוטים להעדיף אלגוריתמים המשתמשים במידע על עקמומיות משטח האנרגיה, כלומר במטריצת הסיאן. עם זאת, עבור רוב המערכות המעשיות חישוב המטריצה הזו יכול להיות יקר מאוד, ולכן נהוג לאמוד אותה בעזרת הערכות שמבוססות על ערכים עוקבים של הגרדיאנט, כפי שקורה בתהליך אופטימיזציית קוויזי-ניוטון.

## בחירת מערכת הקואורדינאטות
הבחירה במערכת הקואורדינטות המתאימה היא קריטית להצלחת האופטימיזציה. לדוגמה, קואורדינטות קארטזיות עשויות להיות מיותרות כיוון שמולקולה לא ליניארית עם N אטומים מכילה  3N - 6 דרגות חופש של ויברציה, בעוד שקואורדינטות קארטזיות מציעות 3N ממדים. יתרה מזו, קואורדינטות קארטזיות מציגות תלות הדדית גבוה, מה שמתבטא בכך שבמטריצת הסיאן קיימים הרבה איברים מחוץ לאלכסון שאינם שואפים לאפס. מצב זה עלול לגרום לבעיות נומריות בתהליך האופטימיזציה, משום שקשה לקבל הערכה מדויקת למטריצה, וחישובה מצריך משאבים רבים. לעומת זאת, כאשר האנרגיה מתוארת באמצעות שדות כוח סטנדרטיים, פותחו שיטות יעילות לגזירת מטריצת הסיאן בקואורדינטות קארטזיות תוך שמירה על סיבוכיות חישובית בסדר גודל דומה לזו של חישוב הגרדיאנט. מנגד, קואורדינטות פנימיות נוטות להראות תלות הדדית נמוכה יותר, אך בנייתן וקביעתן עשויה להיות מורכבת, ובחלק מהמערכות למשל, במערכות עם סימטריה גבוהה, תיאורן יכול להיות מאתגר. בנוסף, תוכנות כימיה חישובית רבות ומודרניות כוללות כלים אוטומטיים ליצירת מערכות קואורדינטות הולמות לצורך האופטימיזציה.

### הגבלת דרגות חופש
ניתן לצמצם את מספר דרגות החופש בתהליך האופטימיזציה. לדוגמה, אפשר להקפיא מיקומים של אטומים מסוימים או לקבוע מראש ערכים עבור אורכי קשר וזוויות. לעיתים קרובות, דרגות חופש אלו מכונות "קפואות".

### אופטימיזציית מצבי מעבר
מבני מצב מעבר (Transition state) יכולים להיקבע על ידי חיפוש נקודות אוכף על פני משטח האנרגיה הפוטנציאלית (PES) של החומר הכימי הרלוונטי. נקודת אוכף מסדר ראשון היא מיקום על פני השטח המתאים למינימום בכל הכיוונים פרט לאחד; נקודת אוכף מסדר שני היא מינימום בכל הכיוונים פרט לשניים, וכן הלאה. בהגדרה מתמטית, נקודת אוכף מסדר *n* מאופיינת כך:  {\partial E \over\partial r} ומטריצת הסיאן (Hessian),  {\partial \partial E\over\partial r_i\partial r_j}, מכילה בדיוק *n* ערכים עצמיים שליליים.

האלגוריתמים לאיתור גאומטריות של מצבי מעבר מתחלקים לשתי קטגוריות עיקריות: שיטות מקומיות ושיטות חצי-גלובליות. שיטות מקומיות מתאימות כאשר נקודת ההתחלה לאופטימיזציה קרובה מאוד למצב המעבר האמיתי (הגדרה של "קרובה מאוד" תינתן בהמשך), ושיטות חצי-גלובליות משמשות כאשר מחפשים לאתר את מצב המעבר עם ידע מוקדם מועט מאוד לגבי הגאומטריה שלו. ישנן שיטות, כגון שיטת הדימר (Dimer method), שנכללת בשתי הקטגוריות.

### חיפושים מקומיים
אופטימיזציה מקומית דורשת ניחוש ראשוני של מצב המעבר שקרוב מאוד למצב המעבר האמיתי. "קרוב מאוד" משמעותו בדרך כלל שהניחוש הראשוני חייב לכלול מטריצת הסיאן (Hessian) עם ערך עצמי שלילי אחד, או שהערך העצמי השלילי המתאים לקואורדינטת התגובה גדול בערכו (ערך מוחלט הכי גדול) מהערכים העצמיים השליליים האחרים. בנוסף, הווקטור העצמי עם הערך העצמי השלילי ביותר חייב להתאים לקואורדינטת התגובה, כלומר, עליו לייצג את השינוי הגאומטרי הקשור בתהליך שמצב המעבר שלו מחפשים.

בהינתן התנאים ההתחלתים הללו, אלגוריתם אופטימיזציה מקומית יכול לנוע "במעלה" לאורך הווקטור העצמי עם הערך העצמי השלילי ביותר ו"במורד" לאורך שאר דרגות החופש, באמצעות שיטה דומה לשיטת קוואזי-ניוטון (quasi-Newton).

### שיטת דימר
שיטת דימר (Dimer method) יכולה לשמש למציאת מצבי מעבר אפשריים ללא ידע מקדים על המבנה הסופי או לשיפור ניחוש ראשוני של מבנה מעבר. הרעיון המרכזי הוא ליצור "דימר" – זוג של נקודות קרובות מאוד זו לזו על ה-PES. האלגוריתם מזיז את הדימר כלפי מעלה (לכיוון אנרגיה גבוהה יותר) תוך כדי חיפוש אחר הכיוון שבו העקמומיות היא הנמוכה ביותר (כלומר, הערך העצמי השלילי ביותר של מטריצת הסיאן). הדימר מסתובב בהתמדה כדי לוודא שהתנועה אכן מתבצעת לאורך כיוון זה. כתוצאה מכך, השיטה מאפשרת להתקדם אל עבר נקודת אוכף (מצב מעבר) ביעילות, גם ללא מידע מוקדם על המבנה המדויק של מצב המעבר.

### טכניקת הרפיה באמצעות אקטיבציה (ART)
טכניקת הרפיה באמצעות אקטיבציה (Activation Relaxation Technique, ART) היא גם שיטה פתוחה למציאת מצבי מעבר חדשים או לשיפור נקודות אוכף ידועות על ה-PES. השיטה פועלת על ידי מציאת הכיוון שבו העקמומיות היא השלילית ביותר (המחושבת באמצעות אלגוריתם Lanczos) והתנועה בכיוון זה, בדומה לשיטת הדימר. עם זאת, ההתקדמות אינה ישירה: לאחר כל תנועה בכיוון זה (הנקראת "אקטיבציה"), המערכת מתאזנת מחדש על גבי המישור המאונך לכיוון התנועה ("הרפיה"). מחזורי האקטיבציה וההרפיה מאפשרים למערכת לנוע ביעילות לכיוון נקודת האוכף תוך מניעת סטיות מיותרות. בשל כך, ART מתאימה במיוחד למצבים שבהם רוצים למצוא מסלולים חדשים או לשפר חישובים קודמים של מצבי מעבר.

### שיטות שרשרת של מצבים (Chain-of-State Methods)
שיטות "שרשרת של מצבים" נועדו למצוא קירוב גאומטרי למצב המעבר בהתבסס על הגאומטריות של המגיבים והתוצרים. לאחר יצירת הקירוב, ניתן להשתמש בו כנקודת התחלה טובה יותר לשיפור האופטימיזציה באמצעות חיפוש מקומי, כפי שתואר קודם.

בשיטות אלו משתמשים בסדרה של וקטורים – כלומר, נקודות על פני משטח האנרגיה הפוטנציאלית (PES) – המחברות בין המגיב לתוצר ומחלקות את מסלול התגובה לחלקים קטנים (דיסקרטיים). נקודות אלו מכונות לעיתים "חרוזים" בשל האנלוגיה לשרשרת חרוזים המחוברים ביניהם בחוטים או קפיצים. לעיתים קרובות, החרוזים נוצרים תחילה באמצעות אינטרפולציה פשוטה בין המגיב לתוצר.  r_i = \frac{i}{N} r_{product} + 1 - \frac{i}{N} r_{reactant} 

כאשר  i \in 0, 1, ..., N. לכל חרוז  r_i יש אנרגיה E(r_i) וכוחות  -\partial E\over\partial r_i הפועלים עליו, כאשר מתבצעת אופטימיזציה תחת אילוצים על מנת למצוא ייצוג מדויק של מסלול התגובה. כדי למצוא את הייצוג המדויק, התהליך מצריך אילוצי מיקום כדי למנוע מהחרוזים להתכנס ישירות חזרה אל גאומטריות המגיב או התוצר.

לעיתים קרובות, אילוץ זה מושג על ידי projecting- הקרנת הרכיבים של הכוח הפועל על כל חרוז, או לחלופין על ידי לאלץ את תנועת כל חרוז במהלך האופטימיזציה, אשר מקבילים למסלול התגובה. לדוגמה, אם לצורך נוחות מגדירים את גרדיאנט האנרגיה בתור  g_i = {\partial E \over\partial r_i} אז גרדיאנט האנרגיה בכל חרוז, לאחר החסרת הרכיב גרדיאנט האנרגיה אשר מקביל למסלול התגובה, נתון על ידי: g_i^\perp = g_i-\tau_i(\tau_i \cdot g_i) = (I-\tau_i \tau_i^T)g_i

כאשר *I* מטריצת היחידה, ו-  \tau_i וקטור יחידה המייצג את המשיק למסלול התגובה בנקודה  r_i.

על ידי הקרנת הרכיבים של גרדיאנט האנרגיה או של צעד האופטימיזציה אשר מקבילים למסלול התגובה, האלגוריתם מצמצם באופן משמעותי את הנטייה של כל אחד מהחרוזים להתכנס ישירות אל מינימום מקומי, מה שמאפשר תיאור מדויק יותר של מסלול התגובה.

### שיטת מעבר סינכרוני ליניארי
השיטה הפשוטה ביותר בקטגוריה זו היא מעבר סינכרוני ליניארי (Linear Synchronous Transit, LST). בשיטה זו נלקחות נקודות בין גאומטריות המגיב והתוצר, והנקודה בעלת האנרגיה הגבוהה ביותר נבחרת להמשך אופטימיזציה בחיפוש מקומי. שיטה מתקדמת יותר היא מעבר סינכרוני ריבועי (Quadratic Synchronous Transit, QST), אשר מרחיבה את LST בכך שהיא מאפשרת מסלול תגובה פרבולי במקום ליניארי. במקום לחבר את מבני המגיב והתוצר בקו ישר ולבחור את הנקודה בעלת האנרגיה הגבוהה ביותר להמשך אופטימיזציה (כפי שנעשה ב-LST), QST משתמשת במסלול מעוקל יותר כדי לתאר טוב יותר את מסלול המעבר.

לאחר יצירת הפרבולה, הנקודה בעלת האנרגיה הגבוהה ביותר לא רק נבחרת אלא גם **מותאמת באמצעות אופטימיזציה בכיוון אנכי לפרבולה**. השיטה מחפשת את מעבר האנרגיה האמיתי על פני משטח האנרגיה הפוטנציאלי (PES) ולא רק את הנקודה הגבוהה ביותר על המסלול המחושב מראש.

### שיטת הרצועה האלסטית הדחופה (Nudged elastic band)
בשיטת הרצועה האלסטית הדחופה (NEB), החרוזים המייצגים את מסלול התגובה מקבלים כוחות קפיציים בנוסף לכוחות הכימיים  {-\partial E \over\partial r_i} , מטרת הכוחות הקפיציים היא לשמר מרווחים אחידים בין החרוזים. הכוחות הקפיציים עבור נקודה i מוגדרת:  f_i = f_i^ \parallel - g_i^ \perp

כאשר:  f_i^ \parallel = k[(( r_{i+1} - r_1) - (r_i - r_{i-1})) \cdot \tau_i]\tau_j

הוא כוח הקפיץ המקביל למסלול בנקודה  r_i, k קבוע הקפיץ ו \tau_i כמו שמוגדר למעלה.

ביישום מסורתי של השיטה, הנקודה בעלת האנרגיה הגבוהה ביותר משמשת להמשך חיפוש מקומי של מצב המעבר. קיימות וריאציות שונות של NEB, כגוןClimbing Image NEB, שבה הנקודה בעלת האנרגיה הגבוהה ביותר נדחפת עוד יותר כלפי מעלה כדי להשיג קירוב טוב יותר למצב המעבר. שיפורים נוספים כוללים שילוב רגרסיה בתהליך גאוסוני (Gaussian Process Regression) כדי להפחית את מספר החישובים הנדרשים. למערכות בעלות גאומטריה לא-אוקלידית (שאינן במרחב  R^2), כמו מערכות מגנטיות, קיימת הרחבה של השיטה בשם Geodesic Nudged Elastic Band, המותאמת למבנים גאומטריים מסובכים יותר.

### שיטת המיתר (String Method)
שיטת המיתר (String Method) משתמשת בקווים רציפים המחברים בין הנקודות, כדי למדוד ולהטיל את אילוצי המרחק ביניהן ולחשב את המשיק למסלול בכל נקודה. בכל שלב בתהליך האופטימיזציה, ניתן להזיז את הנקודות בהתאם לכוח הפועל עליהן בניצב לנתיב. אם אילוץ המרחק השווה בין הנקודות אינו מתקיים לאחר ההזזה, ניתן לבצע חלוקה מחדש של הנקודות, תוך שימוש בייצוג הקו הרציף של הנתיב ליצירת וקטורים חדשים עם מרווחים מתאימים.

גרסה מתקדמת של שיטה זו היא שיטת המיתר המתפתח (Growing String Method), שבה מסלול התגובה נבנה באופן הדרגתי מהקצוות (המגיבים והתוצרים) תוך כדי התקדמות האופטימיזציה.

## השוואה עם טכניקות אחרות
אופטימיזציה גאומטרית שונה באופן יסודי מסימולציית דינמיקה מולקולרית. בסימולציית דינמיקה מולקולרית, מתבצע סימולציה של תנועת המולקולות בזמן, תוך התחשבות בטמפרטורה, בכוחות כימיים, במהירויות התחלתיות, בתנועת בראונית של ממס ועוד, באמצעות יישום חוקי תנועה של ניוטון. המשמעות היא שנתיבי התנועה של האטומים שמתקבלים בחישוב נושאים משמעות פיזיקלית. לעומת זאת, באופטימיזציה גאומטרית לא מתקבל "נתיב" בעל משמעות פיזיקלית – כאן המטרה היא למזער את הכוחות הפועלים על כל אטום בקבוצת האטומים, והדרך בה מושג מינימום זה אינה משמעותית בפני עצמה. אלגוריתמים שונים לאופטימיזציה עשויים להגיע לאותה תוצאה, המבנה בעל האנרגיה המינימלית, אך עשויים לעשות זאת בדרכים שונות.

## ראו גם
*כימואינפורמטיקה
*פיזיקה חישובית
*כימיה חישובית
*דינמיקה מולקולרית
*מודל כימי קינטי

## הערות שוליים

קטגוריה:כימיה פיזיקלית
חישובית
קטגוריה:כימיה תאורטית
ואן דר ואלס