אלגוריתמים של תרמוסטט

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

תרמוסטטים ואופני הפעולה שלהם

התרמוסטט של אנדרסון- שיטת התנגשות סטוכסטית

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

$H=\sum _{i=1}^{N}(P_{i}^{2}/2m_{i})+\phi (q)$>

$(1)dq_{i}/dt=\partial H/\partial p_{i},dp_{i}/dt=-\partial H/\partial q_{i}$

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

$P(t,f)=fexp(-ft)$

כאשר P(t,f)dt הוא ההסתברות שההתנגשות הבאה תתרחש במקטע הזמן [t,t+∆t].

הסימולציה

בוחרים סט ראשוני של מיקומים מרחביים ותנעים: $q^{N}(0),P^{N}(0)$ ומבצעים אינטגרל על משוואות התנועה (1) עד שמגיעים לזמן בו מתרחשת ההתנגשות הסטוכסטית הראשונה. נסמן את החלקיק שהשתתף בהתנגשות ב- i. אז התנע של החלקיק לאחר ההתנגשות נבחר אקראית מהתפלגות בולצמן בטמפרטורה T. השינוי בתנע מתרחש מיידית. שאר החלקיקים אינם מושפעים מההתנגשות, ולכן יש לבצע אינטגרל על משוואות התנועה המייצגות את סך כל החלקיקים, עד הזמן בו מתרחשת ההתנגשות הבאה. וחוזר חלילה על האלגוריתם המתואר, זהו אלגוריתם של אנדרסון. התוצאה של תהליך זה היא מסלול מוגבל של N חלקיקים בנפחV: $(q^{N}(t),P^{N}(t))$. מסלול זה יכול לשמש לחישוב ממוצעי זמן של כל גודל שהוא $(F(q)^{N}(t),P^{N}(t);V)$ באופן הבא:

$F^{'}=\lim _{T\to \infty }\int _{0}^{T}F((q)^{N}(t),P^{N}(t);V(t))dt$

השילוב של המכניקה של ניוטון והתנגשות סטוכסטית הופכות את הסימולציה הזו לשרשרת מרקוב. אנדרסון הראה שהחלוקה הקנונית במרחב הפאזה אינה משתנה תחת יישום חוזר ונשנה של האלגוריתם הנ"ל. בשילוב עם העובדה ששרשרת מרקוב אינה מחזורית ואינה ניתנת להפחתה, משתמע שהאלגוריתם של אנדרסון אכן מייצר התפלגות קנונית. למעשה, ממוצע הזמן של כל F כלשהו המחושב באמצעות מסלול אנדרסון שווה לממוצע האנסמבל של F עבור האנסמבל הקנוני בו הטמפרטורה היא T. $F^{'}=1/N!Q(N,V,T)\int _{F}((q)^{N}(t),P^{N}(t);V(t))$

$Q(N,V,T)=1/N!\int exp(-\beta H(q^{N},P^{N};V))dq^{N}dP^{N}$ כאשר Q היא פונקציית החלוקה של האנסמבל הקנוני ו- $\beta =1/kT$

החיסרון

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

התרמוסטט של נוייס הובר - שיטה של המערכת המורחבת

הרחבה של משוואות התנועה

האנרגיה של מערכת של N חלקיקים בטמפרטורה קבועה מבצעת תנודות סביב ערך שיווי המשקל. לכן, על מנת לדמות מערכת שכזו צריך מנגנון כלשהו שיגרום לתנודות אלה. לשם כך, נוייס מציע לגראנגיאן מורחב, לאגרנג'יאן המכיל קואורדינטות ומהירויות נוספות ומאלכותיות. למעשה, שיטת הלגרנג'יאן המורחב הוצגה לראשונה על ידי אנדרסון בסימולציית הלחץ הקבוע. השימוש בשיטה זו נפוץ יותר באופן שהציע נוייס ולכן שיטה זו מכונה גם התרמוסטט של נוייס-הובר. נניח שהמערכת המדומה מכילה N חלקיקים עם קואורדינטות $q'_{i}$, מסות $m_{i}$, אנרגיה פוטנציאלית $\phi (q')$ ותנע $p'_{i}$. דרגת חופש נוספת s משמשת כמערכת חיצונית למערכת המדומה. בנוסף, ישנם משתנים וירטואליים $(q',p',t)$ שקשורים למשתנים האמיתיים (q^',p^',t) באופן הבא:

$q'_{i}=q_{i},p'_{i}=p_{i}/s,t'=\int _{0}^{t}dt/s$

והמהירות האמיתית מבוטאת על ידי:

$dq'_{i}/dt'=sdq'_{i}/dt=sdq_{i}/dt$

כעת ניתן להסתכל על הטרנספורמציות הנ"ל בשינוי זמני: $dt'=dt/s$.

הלגרנג'יאן של המערכת המורחבת של N חלקיקים ומשתנה s במונחי המשתנים הווירטואליים:

$L_{N}=\sum _{i=1}^{N}(m_{i}/2)s^{2}{\dot {q_{i}}}^{2}-\phi (q)+(Q/2){\dot {s}}^{2}-gKTln(s)$

כאשר Q זוהי המסה האפקטיבית שקשורה לs ו-g מייצג במהותו את מספר דרגות החופש של המערכת. עם זאת, ערכו המדויק יבחר כך שתהיה התפלגות קנונית בשיווי משקל. בנוסף, לוגריתם תלות הפוטנציאל במשתנה s היא חיונית לייצור האנסמבל הקנוני. הקשר בין התנע לבין $q_{i}$ ו-s ניתן באופן הבא: $P_{i}=\partial L_{N}/\partial {\dot {q_{i}}}=m_{i}s^{2}{\dot {q_{i}}},p_{s}=\partial L_{N}/\partial {\dot {s}}=Q{\dot {s}}$

באמצעות קשרים אלה ניתן להגיע להמילטוניאן של המערכת המורחבת של N חלקיקים והמשתנה s במונחי המשתנים הווירטואליים: $H_{N}=\sum _{i=1}^{N}P_{i}^{2}/2m_{i}s^{2}+\phi (q)+p_{s}^{2}/2Q+gkTln(s)$

לפי הפורמליזם של ההמילטוניאן משוואות התנועה ניתנות על ידי:

$dq_{i}/dt=\partial H_{N}/\partial P_{i}=P_{i}/m_{i}s^{2}$

$dP_{i}/dt=-\partial H_{N}/\partial q_{i}=-\partial \phi /\partial q_{i}$

$ds/dt=-\partial H_{N}/\partial P_{s}=P_{S}/Q$

$dp_{s}/dt=-\partial H_{N}/\partial s=(\sum _{i=1}^{N}P_{i}^{2}/m_{i}s^{2}-gKT)/s$

כעת ניתן לראות שההמילטוניאן המורחב $H_{N}$ נשמר כאשר המערכת המורחבת נוצרה לפי משוואות התנועה הנ"ל. לכן, שיטה זו מייצרת אנסמבל מיקרוקנוני עבור המערכת המורחבת.

שימור האנסמבל הקנוני

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

$Z=\int \delta [H_{0}(p/s,q)+p_{s}/2Q+gKTln(s)-E]dp_{s}dsdpdq$

כאשר E מייצג את האנרגיה שניתנה מראש ו- $H_{0}$ מייצג את ההמילטוניאן הקלאסי:

$H_{0}(p,q)=\sum _{i=1}^{N}p_{i}/2m_{i}+\phi (q)$

באמצעות החלפת משתנים ניתן להמיר את המשתנים הווירטואליים למשתנים האמתיים, ולקבל: $Z=(1/gkT)exp(((3N+1)/g)(E/kT))\int exp((-(3N+1)/g)(p_{s}/2QkT)dp_{s}\int exp((-(3N+1)/g)(H_{0}(p,q)/kT)dp'q'$

נשים לב, שעבור g= 3N+1 פונקציית החלוקה של המערכת המורחבת זהה לפונקציית החלוקה של המערכת המקורית באנסמבל הקנוני מלבד פקטור קבוע: $Z=C\int exp(-\beta H_{0}(p',q')dp'q'$

(כאשר $\beta =1/kT$)

ואז, פונקציית החלוקה בשיווי משקל:

$\rho (p',q')=exp(-\beta H_{0}(p',q'))$

ולבסוף נקבל:

$F'=<A(p/s,q)>_{e}(extended)=<A(p',q')>_{c}(canonical)$

תרמוסטט גאוסיאן- שינוי קנה המידה של המהירות

זו השיטה הראשונה שהוצעה לקיבוע ערך הטמפרטורה (ללא תנודות). בשיטה זו המהירויות נמדדות לפי קנה המידה:

$p_{i}\to ({\sqrt {T_{0}}}/T)p_{i}$

כאשר $T_{0}$ זו הטמפרטורה הרצויה ו-T זו הטמפרטורה שחושבה באמצעות מהירות החלקיקים.

הרחבה של שיטה זו מרמזת על אילוץ משוואות התנועה כדי לשמור על הטמפרטורה קבועה. עקרון גאוסיאן עבור האילוץ הפחות משמעותי קובע כי כח שנוסף להגבלת תנועה החלקיקים צריך להיות מאונך למשטח. לפי עקרון זה, כח מאלץ $\alpha p_{i}$ מתווסף למשוואת הכוחות (1) וכתוצאה מכך המשוואות לא מוצגות בצורה הקנונית. הפרמטר α נקבע לפי הדרישה שהטמפרטורה, או האנרגיה הקינטית הכוללת, קבועה:

$\sum _{i=1}^{N}p_{i}^{2}/2m_{i}=gkT/2$

$\sum _{i=1}^{N}(p_{i}/m_{i})dP_{i}/dt=0$

$\Rightarrow \alpha =-(\sum _{i=1}^{N}(p_{i}/m_{i})d\phi /dq_{i})/\sum _{i=1}^{N}p_{i}^{2}/m_{i}$

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

חסרון

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

לקריאה נוספת

• (Oxford Graduate Texts) Mark E. Tuckerman - Statistical Mechanics_ Theory and Molecular Simulation-Oxford University Press, USA (2010).pdf
• H.C. Anderson, Molecular dynamics simulations at constant pressure and/or temperature, J. Chem. Phys, 72(4), 15 Feb. 1980.

הערך באדיבות ויקיפדיה העברית, קרדיט,
רשימת התורמים
רישיון cc-by-sa 3.0

29079236אלגוריתמים של תרמוסטט