ניתוח המסות ונוזלים מסימולציות דינמיקה מולקולרית Ab Initio עם חבילת UMD

Summary

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

Abstract

פיתחנו חבילת קוד פתוח מבוססת פייתון כדי לנתח את התוצאות הנובעות מסימולציות של דינמיקה מולקולרית של AB initio של נוזלים. החבילה מתאימה ביותר ליישומים במערכות טבעיות, כמו סיליקט ותחמוצות נמסות, נוזלים על בסיס מים ונוזלים סופר-ביקורתיים שונים. החבילה היא אוסף של סקריפטים של פייתון הכוללים שתי ספריות מרכזיות העוסקות בפורמטי קבצים ועם קריסטלוגרפיה. כל קבצי ה- Script מופעלים בשורת הפקודה. אנו מציעים תבנית פשוטה לאחסון המסלולים האטומיים והמידע התרמודינמי הרלוונטי של הסימולציות, שנשמר בקבצי UMD, עומד עבור דינמיקה מולקולרית אוניברסלית. חבילת UMD מאפשרת חישוב של סדרה של תכונות מבניות, הובלה ותרמודינמית. החל בפונקציית חלוקת הזוגות הוא מגדיר אורכי קשר, בונה מטריצת קישוריות בין-אטומית, ובסופו של דבר קובע את ההפרטציה הכימית. קביעת חיי המין הכימי מאפשרת לבצע ניתוח סטטיסטי מלא. לאחר מכן סקריפטים ייעודיים מחשבים את ההעתקים הריבועיים הממוצעים עבור האטומים, כמו גם עבור המינים הכימיים. ניתוח המתאם העצמי המיושם של המהירויות האטומיות מניב את מקדמי הדיפוזיה ואת ספקטרום הרטט. אותו ניתוח המיושם על הלחצים מניב את הצמיגות. החבילה זמינה דרך אתר GitHub ובאמצעות דף ייעודי משלה של פרויקט ERC IMPACT כחבילת גישה פתוחה.

Introduction

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

עם הופעת כוח המחשוב והזמינות של מחשבי-על, שתי טכניקות אטומיסטיות מספריות עיקריות משמשות כיום לחקר המצב הדינמי של מערכת אטומית לא גבישית, מונטה ^קרלו1 ודינמיקה מולקולרית (MD)^1,2. בסימולציות מונטה קרלו נדגם באופן אקראי המרחב התצורה; שיטות מונטה קרלו מראות קנה מידה ליניארי במקביל אם כל תצפיות הדגימה אינן תלויות זו בזו. איכות התוצאות תלויה באיכות מחולל המספרים האקראיים ובנציגות הדגימה. שיטות מונטה קרלו מראות קנה מידה ליניארי במקביל אם הדגימה אינה תלויה זו בזו. בדינמיקה מולקולרית (MD) המרחב התצורלי נדגם על ידי מסלולים אטומיים תלויי זמן. החל מתצורה נתונה, המסלולים האטומיים מחושבים על ידי שילוב משוואות התנועה הניוטון. ניתן לחשב את הכוחות הבין-אטומיים באמצעות מודלים של פוטנציאלים בין-אטומיים (ב-MD קלאסי) או בשיטות של עקרונות ראשונים (ב-AB initio, או בעקרונות הראשונים, MD). איכות התוצאות תלויה באורך המסלול וביכולתו לא להימשך למינימה מקומית.

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

כתוצאה מכך, הרבה מאמץ הוקדש להדמיית התוצאות, וחבילות שונות זמינות היום, בפלטפורמות שונות, בקוד פתוח או לא [^Ovito3, VMD4, ^Vesta5, ^Travis6 וכו ']. כל כלי ההדמיה הללו מתמודדים ביעילות עם מרחקים בין-אטומיים, וככאלה, הם מאפשרים חישוב יעיל של פונקציות התפלגות זוגות ומקדמי דיפוזיה. לקבוצות שונות המבצעות סימולציות דינמיקה מולקולרית בקנה מידה גדול יש תוכנה קניינית לנתח מאפיינים שונים אחרים הנובעים מהסימולציות, לפעמים ב- shareware או בצורות אחרות של גישה מוגבלת לקהילה, ולעתים מוגבלות בהיקפן ובשימושן בחבילות מסוימות ספציפיות. אלגוריתמים מתוחכמים לחילוץ מידע על מליטה בין-אטומית, דפוסים גיאומטריים ותרמודינמיקה מפותחים ומיושמים בחלק מהחבילות הללו3,4,5,6,7 וכו '.

כאן אנו מציעים את חבילת UMD - חבילת קוד פתוח שנכתבה בפייתון כדי לנתח את הפלט של סימולציות דינמיקה מולקולרית. חבילת UMD מאפשרת חישוב של מגוון רחב של מאפיינים מבניים, דינמיים ותרמודינמיים (איור 1). החבילה זמינה דרך אתר GitHub (https://github.com/rcaracas/UMD_package) ובאמצעות דף ייעודי (http://moonimpact.eu/umd-package/) של פרויקט ERC IMPACT כחבילת גישה פתוחה.

כדי להפוך את הטיפול לאוניברסלי וקל יותר, הגישה שלנו היא קודם כל לחלץ את כל המידע הקשור למצב התרמודינמי ולמסלולים האטומיים מקובץ הפלט של הריצה המולקולרית-דינמית בפועל. מידע זה מאוחסן בקובץ ייעודי, שתבניתו אינה תלויה בחבילת MD המקורית שבה בוצעה הסימולציה. אנו קוראים לקבצים אלה קבצי "umd", אשר מייצג דינמיקה מולקולרית אוניברסלית. בדרך זו, חבילת UMD שלנו יכולה לשמש בקלות על ידי כל קבוצת ab initio עם כל תוכנה, כל עם מאמץ מינימלי של הסתגלות. הדרישה היחידה להשתמש בחבילה הנוכחית היא לכתוב את הניתח המתאים מהפלט של תוכנת MD מסוימת לתבנית הקובץ umd, אם זה עדיין לא קיים. לעת עתה, אנו מספקים מנתחים כאלה עבור ^VASP8 ואת חבילות ^QBox9 .

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

קבצי umd הם קבצי ASCII; הרחבה אופיינית היא "umd.dat" אך לא חובה. כל רכיבי הניתוח יכולים לקרוא קבצי ASCII בתבנית umd, ללא קשר לסיומת השם בפועל. עם זאת, חלק מהסריפטים האוטומטיים שנועדו לבצע סטטיסטיקה מהירה בקנה מידה גדול על פני מספר סימולציות מחפשים במיוחד קבצים עם הסיומת .dat umd. כל מאפיין פיזי מבוטא בשורה אחת. כל שורה מתחילה במילת מפתח. בדרך זו התבנית ניתנת להתאמה רבה ומאפשרת להוסיף מאפיינים חדשים לקובץ umd, תוך שמירה על הקריאות שלו לאורך גירסאות. 30 השורות הראשונות של קובץ umd של סימולציה של פירוליט ב 4.6 GPa ו 3000 K, בשימוש להלן בדיון, מוצגים באיור 2.

איור 2: תחילתו של קובץ umd המתאר סימולציה של פירוליט נוזלי ב 4.6 GPa ו 3000 K.
לאחר הכותרת מופיע התיאור של כל תמונה. כל מאפיין כתוב בשורה אחת, המכילה את שם המאפיין הפיזי, הערך והיחידות, כולם מופרדים באמצעות רווחים. אנא לחץ כאן כדי להציג גירסה גדולה יותר של איור זה.

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

אז כל שלב של הסימולציה מפורט. ראשית, הפרמטרים התרמודינמיים המיידיים ניתנים, כל אחד על קו אחר, המציין (i) את שם הפרמטר, כמו אנרגיה, לחצים, לחץ הידרוסטטי שווה ערך, צפיפות, נפח, פרמטרי סריג וכו ', (ii) הערך שלו, ו -(iii) יחידותיו. שולחן המתאר את האטומים מגיע הבא. קו כותרת מעניק את האמצעים השונים, כגון מיקומים קרטזיים, מהירויות, מטענים וכו ', והיחידות שלהם. ואז כל אטום מפורט בשורה אחת. לפי קבוצות של שלושה, המתאימים לשלושת צירי x, Y, z , הערכים הם: המיקום המופחת, מיקומי הקרטזי מקופלים לתא הסימולציה, עמדות הקרטזיות (שלוקחות בחשבון כראוי את העובדה שאטומים יכולים לחצות מספר תאי יחידה במהלך סימולציה), המהירויות האטומיות והכוחות האטומיים. שני הערכים האחרונים הם סקלרים: מטען ורגע מגנטי.

שתי ספריות עיקריות מבטיחות את תפקודה התקין של החבילה כולה. ספריית umd_process.py עוסקת בקבצי umd, כגון קריאה והדפסה. ספריית crystallography.py עוסקת בכל המידע הקשור למבנה האטומי בפועל. הפילוסופיה הבסיסית של הספרייה crystallography.py היא להתייחס לסריג כאל מרחב וקטורי. הפרמטרים של תא היחידה יחד עם הכיוון שלהם מייצגים את וקטור הבסיס. "החלל" יש סדרה של תכונות סקלריות (נפח ספציפי, צפיפות, טמפרטורה, ומספר ספציפי של אטומים), תכונות תרמודינמיות (אנרגיה פנימית, לחץ, קיבולת חום, וכו '), וסדרת תכונות טנזוריות (מתח וגמישות). אטומים מאכלסים את החלל הזה. כיתת "הסריג" מגדירה הרכב זה, לצד חישובים קצרים שונים, כגון נפח ספציפי, צפיפות, קבלת הסריג ההדדי מהישיר וכו '. כיתת "אטומים" מגדירה את האטומים. הם מאופיינים בסדרה של תכונות סקלריות (שם, סמל, מסה, מספר אלקטרונים וכו ') וסדרת תכונות וקטוריות (המיקום בחלל, ביחס לבסיס הווקטורי המתואר במחלקת הסריג, או ביחס לקואורדינטות קרטזיות אוניברסליות, מהירויות, כוחות וכו '). מלבד שתי מחלקות אלה, ספריית crystallography.py מכילה סדרה של פונקציות לביצוע מגוון בדיקות וחישובים, כגון מרחקים אטומיים או כפל תאים. הטבלה המחזורית של הרכיבים כלולה גם כמילון.

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

Protocol

1. ניתוח הדינמיקה המולקולרית פועל

הערה: החבילה זמינה דרך אתר GitHub (https://github.com/rcaracas/UMD_package) ובאמצעות דף ייעודי (http://moonimpact.eu/umd-package/) של פרויקט ERC IMPACT כחבילת גישה פתוחה.

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

דגל	משמעות	סקריפט המשתמש בו	ערך ברירת מחדל
-h	עזרה קצרה	כל
-f	שם קובץ UMD	כל
-אני.	שלבי תרמיזציה שיש לבטל	כל	0
-אני.	קובץ קלט המכיל את איגרות החוב הבין-אטומיות	דגימה	איגרות חוב.קלט
-s	דגימה של התדירות	msd, speciation	1 (כל צעד נחשב)
-a	רשימת אטומים או אניונים	דגימה
-c	רשימת ציטוטים	דגימה
-אני.	אורך אג"ח	דגימה	2
-t	טמפרטורה	תנודות, ריולוגיה
-v	דיסקרטיזציה של רוחב חלון הדגימה של המסלול לניתוח תזוזה מרובעת ממוצעת	msd	20
-z	דיסקרטיזציה של תחילת חלון הדגימה של המסלול לניתוח תזוזה מרובעת ממוצעת	msd	20

טבלה 1: הדגלים הנפוצים ביותר המשמשים בחבילת UMD והמשמעות הנפוצה ביותר שלהם.

2. התחל בהפיכת הפלט של סימולציית MD המבוצעת בקוד עקרונות ראשונים, כגון ^VASP8 או ^QBox9, לקובץ UMD.
   1. אם סימולציות MD נעשו ב- VASP, אז בסוג שורת הפקודה:
      VaspParser.py -f -i <ביטוליטיאלStep>
      כאשר -f flag מגדיר את שם קובץ VASP OUTCAR ואת אורך התרמיזציה של - i.
      הערה: השלב הראשוני, המוגדר על-ידי -i מאפשר להשליך את השלבים הראשונים של הסימולציות, המייצגות את התרמיזציה. בריצת דינמיקה מולקולרית טיפוסית, החלק הראשון לחישוב מייצג את התרמיזציה, כלומר את הזמן שלוקח למערכת את כל האטומים לתאר התפלגות דמוית גאוס של הטמפרטורה, ולמערכת כולה להפגין תנודות של הטמפרטורה, הלחץ, האנרגיה וכו 'סביב ערכי שיווי המשקל. חלק תרמיזציה זה של הסימולציה לא צריך להילקח בחשבון בעת ניתוח המאפיינים הסטטיסטיים של הנוזל.
3. שנה את . קבצים ב - . קבצי xyz כדי להקל על הדמיה על חבילות שונות אחרות, כמו ^VMD4 או ^Vesta5. בסוג שורת הפקודה:
   umd2xyz.py -f -i -s
   כאשר -f מגדיר את שם . קובץ umd , -i מגדיר את תקופת התרמיזציה שיש למחוק, ו - s תדירות הדגימה של המסלול המאוחסן ב. קובץ umd . ערכי ברירת המחדל הם -i 0 –s 1, כלומר בהתחשב בכל שלבי הסימולציה, מבלי להיפטר מהם.
4. הפוך את קובץ umd לקבצי POSCAR מסוג VASP באמצעות קובץ ה- Script umd2poscar.py; ניתן לבחור תמונות של הסימולציות בתדירות מוגדרת מראש. בסוג שורת הפקודה:
   umd2poscar.py -f -i -l -s
   כאשר –l מייצג את השלב האחרון שיש להפוך לקובץ POSCAR. ערכי ברירת המחדל הם -i 0 -l 10000000 -s 1. ערך זה של – גדול מספיק כדי לכסות מסלול שלם טיפוסי.

2. בצע את הניתוח המבני

1. הפעל את קובץ ה- script gofrs_umd.py כדי לחשב את פונקציית התפלגות הזוג (PDF) _gᴀʙ(r) עבור כל זוגות הסוגים האטומיים A ו- B (איור 3). הפלט כתוב בקובץ ASCII אחד, מופרד באמצעות טאבים, כאשר הסיומת gofrs.dat. בסוג שורת הפקודה:
   gofrs_umd.py -f -s < Sampling_Frequency > -d -i
   הערה: ברירות המחדל הן Sampling_Frequency (התדירות לדגימה של המסלול) = שלב אחד; דיסקרטיזציהבחטרוול (להתוויית g(r)) = 0.01 Å; InitialStep (מספר השלבים בתחילת המסלול שנמחקים) = 0. ה- PDF הרדיאלי, _gᴀʙ(r) הוא המספר הממוצע של אטומים מסוג B במרחק d_ᴀʙ בתוך קליפה כדורית של רדיוס r ועובי dr המתמקדת באטומים מסוג A (איור 3):
   
   
   עם ρ הצפיפות האטומית, _NA ו- _NB מספר האטומים מסוג A ו- B, ו- δ(r−_rᴀʙ) פונקציית הדלתא השווה ל- 1 אם האטומים A ו- B נמצאים במרחק בין r ל- r + dr. האבסיסה של המקסימום הראשון של _gᴀʙ(r) מעניקה את אורך הקשר בעל ההסתברות הגבוהה ביותר בין האטומים מסוג A ו- B, שהוא הקרוב ביותר למרחק קשר ממוצע שאנו יכולים לקבוע. המינימום הראשון מפריד את היקף ספירת התיאום הראשונה. מכאן שהאינטגרל מעל מסמך PDF עד למינימום הראשון מעניק את מספר התיאום הממוצע. הסכום של הפיורה הופך את _הג'י(r) לכל זוגות הסוגים האטומיים A ו- B מניב את תבנית עקיפה של הנוזל, כפי שהושג באופן ניסיוני עם דיפרקטור. עם זאת, במציאות, כפי שלעתים קרובות ספירות התיאום מסדר גבוה חסרות מן _gᴀʙ(r), דפוס עקיפה לא ניתן להשיג בשלמותו.

איור 3: קביעת פונקציית התפלגות הזוג.
(א) עבור כל אטום של מין אחד (לדוגמה אדום), כל האטומים של המינים המתאמים (לדוגמה אפור ו/או אדום) נספרים כפונקציה של מרחק. (ב) גרף התפלגות המרחק המתקבל עבור כל תמונה, שבשלב זה הוא רק אוסף של פונקציות דלתא, הוא בממוצע על פני כל האטומים וכל התמונות ומשוקלל על ידי התפלגות הגז האידיאלית כדי ליצור (ג) את פונקציית התפלגות הזוגות הרציפה. המינימום הראשון של ה- g(r) הוא הרדיוס של ספירת התיאום הראשונה, המשמשת בהמשך בניתוח הדגימה. אנא לחץ כאן כדי להציג גירסה גדולה יותר של איור זה.

2. לחלץ את מרחקי הקשר הבין-אטומיים הממוצעים כמו radii של ספירות התיאום הראשונות. עבור זה, לזהות את המיקום של המקסימום הראשון של פונקציות _gᴀʙ(r): התווה את הקובץ .dat gofrs ביישום גיליון אלקטרוני ולחפש את מקסימה ומינימה עבור כל זוג אטומים.
3. זהה את הרדיוס של ספירת התיאום הראשונה, כמינימום הראשון של מסמך PDF, _gᴀʙ(r), באמצעות תוכנת גיליון אלקטרוני. זהו הבסיס לכל הניתוח המבני של הנוזל; מסמך PDF מניב את מצב הקשר הממוצע של האטומים בנוזל.
4. לחלץ את המרחקים של המינימה הראשונה, כלומר, את המורסה, ולכתוב אותם בקובץ נפרד, שנקרא, למשל, bonds.input. לחלופין, הפעל את אחד מקבצי ה- Script analyze_gofr של חבילת UMD כדי לזהות את מקסימה ואת המינימה של פונקציות _gᴀʙ(r). בסוג שורת הפקודה:
   analyze_gofr_semi_automatic.py
5. לחץ על המיקום של המקסימום ואת המינימום של הפונקציה _gᴀʙ(r) המוצגת בגרף שנפתח על ידי התוכנית. קובץ ה- Script סורק באופן אוטומטי את התיקיה הנוכחית, מזהה את כל קבצי gofrs.dat ומבצע את הניתוח עבור כל אחד מהם. לחץ שוב על המקסימום והמינימום בחלון בכל פעם שהתסריט זקוק לניחוש ראשוני משכיל.
6. פתח והסתכל על הקובץ שנוצר באופן אוטומטי הנקרא bonds.input המכיל את מרחקי הקשר הבין-אטומיים.

3. בצע את ניתוח הדגימות

1. לחשב את הטופולוגיה של מליטה בין האטומים, באמצעות המושג קישוריות בתוך תורת הגרפים: האטומים הם הצמתים והקשרים הבין-אטומיים הם הנתיבים. הסקריפט speciation_umd.py זקוק למרחקי הקשר הבין-אטומיים המוגדרים בקובץ bonds.input .
   הערה: מטריצת הקישוריות בנויה בכל שלב: שני אטומים הנמצאים במרחק קטן יותר מהרדיוס של ספירת התיאום הראשונה המתאימה שלהם נחשבים מקושרים, כלומר מחוברים. רשתות אטומיות שונות בנויות על ידי התייחסות לאטומים כאל צמתים בגרף שחיבוריהם מוגדרים על ידי קריטריון גיאומטרי זה. רשתות אלה הן המין האטומי, וההרכב שלהן מגדיר את ההשתתפות האטומית בנוזל המסוים הזה (איור 4).

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

2. הפעל את סקריפט ההמחקה כדי להשיג את מטריצת הקישוריות ולקבל את הפוליהדרה הקואורדינציה או את הפולימר. בסוג שורת הפקודה:
   speciation_umd.py -f -s -i -l -c -a -m -r
   כאשר הדגל -i נותן את הקובץ עם מרחקי הקשר הבין-אטומיים, שהופקו למשל בשלב הקודם. לחלופין, הפעל את התסריט באורך יחיד אחד עבור כל איגרות החוב המוגדרות על-ידי הדגל -l.
   הערה: הדגל -c מציין את האטומים המרכזיים ואת -דגל ה- ligands. הן אטומים מרכזיים והן ליגנדים יכולים להיות מסוגים שונים; במקרה זה, הם חייבים להיות מופרדים על ידי פסיקים. דגל מ' מעניק את הזמן המינימלי שמין חייב לחיות כדי להיחשב בניתוח. כברירת מחדל זמן מינימלי זה הוא אפס, כל המופעים נספרים בניתוח הסופי.
   1. הפעל את קובץ ה- script speciation_umd.py עם הדגל – r 0, אשר מדגים את גרף הקישוריות ברמה הראשונה כדי לזהות את הפוליהדרה של התיאום. לדוגמה, אטום מרכזי, המוזכר כקטטה , עשוי להיות מוקף באנונים אחד או יותר (איור 4). סקריפט ההפרטה מזהה כל אחד ואחד מהפוליהדרה של התיאום. הממוצע המשוקלל של כל התיאום polyhedra נותן את מספר התיאום, זהה לזה המתקבל משילוב של PDF. בסוג שורת הפקודה:
      speciation_umd.py -f -i -c -a -r 0
      הערה: מספרי התיאום הממוצעים בנוזלים הם מספרי שברים. שבר זה נובע מהמאפיין הממוצע של הקואורדינציה. ההגדרה המבוססת על דגימה מניבה ייצוג אינטואיטיבי ואינפורמטיבי יותר של מבנה הנוזל, שבו הפרופורציות היחסיות של המינים השונים, כלומר, קואורדינציות, מכומתות.
   2. הפעל את קובץ ה- script speciation_umd.py עם הדגל – r 1, אשר מדגים את גרף הקישוריות בכל רמות העומק כדי להשיג את הפילמור. לרשת דרך הגרף האטומי יש עומק מסוים, שכן אטומים נקשרים רחוק יותר לאג"ח אחרות (למשל, ברצפים של סימונים ואניונים מתחלפים) (איור 4).
3. פתח את שני הקבצים . פופול.dat ו. Stat.dat ברציפות; אלה מהווים את הפלט של סקריפט ההפרטה. כל אשכול כתוב בשורה אחת, המציין את הנוסחה הכימית שלו, את הזמן שבו הוא נוצר, את הזמן שבו הוא מת, את חייו, מטריצה עם רשימת האטומים היוצרים את האשכול הזה. שרטט את חייו של כל אשכול אטומי של כל המינים הכימיים שנמצאו בסימולציה כפי שנמצא בקובץ .dat .popul (איור 5).
4. שרטט את ניתוח האוכלוסייה עם השפע של כל מין, כפי שנמצא ב . קובץ stat.dat . ניתוח זה, הן מוחלט והן יחסי, מתאים לסטטיסטיקה בפועל של הפוליהדרה התיאום עבור המקרה -r 0; במקרה של פילמור, עם -r 1 זה צריך להיות מטופל בזהירות כמו כמה נורמליזציה על המספר היחסי של אטומים ייתכן שיהיה צורך ליישם. השפע מתאים לאינטגרל לאורך החיים. ה. הקובץ .dat מפרט גם את הגודל של כל אשכול, כלומר, כמה אטומים יוצרים אותו.

4. מקדם דיפוזיה ממוחשב

1. לחלץ את תזוזות ריבוע ממוצע (MSD) של האטומים כפונקציה של זמן כדי לקבל את הדיפוזיה העצמית. הנוסחה הסטנדרטית של MSD היא:
   
   שבו הטרומיות הן נחירות מחדש. בעזרת הכלי MSD, ישנן דרכים שונות לנתח את ההיבטים הדינמיים של הנוזלים.
   הערה: T הוא הזמן הכולל של הסימולציה ו- N _α הוא מספר האטומים מסוג α. הזמן הראשוני _t0 הוא שרירותי ומשתרע על פני המחצית הראשונה של הסימולציה. _Ninit הוא מספר הפעמים ההתחלתיות. τ הוא רוחב מרווח הזמן שעליו מחושב MSD; הערך המרבי שלו הוא חצי מאורך הזמן של הסימולציה. ביישום MSD טיפוסי, כל חלון מתחיל בסוף החלון הקודם. אבל דגימת ספסר יכולה להאיץ את החישוב של MSD, מבלי לשנות את השיפוע של MSD. לשם כך, החלון i-th מתחיל בזמן _t0(i), אבל החלון (i+1)-th מתחיל בזמן _t0(i) + τ + v, שבו הערך של v מוגדר על-ידי המשתמש. באופן דומה, רוחב החלון גדל בשלבים נפרדים שהוגדרו על-ידי המשתמש, ככזה: τ(i) = τ(i-1) + z. הערכים של z ("צעד אופקי") ו- v ("צעד אנכי") הם חיוביים או אפסיים; ברירת המחדל עבור שניהם היא 20.
2. חשב את MSD באמצעות סדרת קבצי script msd_umd . הפלט שלהם מודפס ב. msd.dat הקובץ, שבו MSD של כל סוג אטומי, אטום או אשכול מודפס בעמודה אחת כפונקציה של זמן.
   1. חשב את ה- MSD הממוצע של כל סוג אטומי. ה- MSD מחושב עבור כל אטום ולאחר מכן ממוצע עבור כל סוג אטומי. קובץ הפלט מכיל עמודה אחת עבור כל סוג אטומי. בסוג שורת הפקודה:
      msd_umd.py -f -z -v <ומפ ההגיוני> -b <בוליסטי>
   2. חשב את ה- MSD של כל אטום. ה- MSD מחושב עבור כל אטום ולאחר מכן ממוצע עבור כל סוג אטומי. קובץ הפלט מכיל עמודה אחת עבור כל אטום בסימולציה ולאחר מכן עמודה אחת עבור כל סוג אטומי. תכונה זו מאפשרת לזהות אטומים המתפזרים בשתי סביבות שונות, כגון נוזל וגז, או שני נוזלים. בסוג שורת הפקודה:
      msd_all_umd.py -f -z -v -b <בוליסטי>
   3. לחשב את MSD של המינים הכימיים. השתמש באוכלוסיית האשכולות המזוהים עם קובץ ה- Script של הדגימה ומודפסים ב . קובץ .dat פופול . ה- MSD מחושב עבור כל אשכול בודד. קובץ הפלט מכיל עמודה אחת עבור כל אשכול. כדי להימנע מלשקול פולימרים בקנה מידה גדול, להציב מגבלה על גודל האשכול; ברירת המחדל שלו היא 20 אטומים. בסוג שורת הפקודה:
      msd_cluster_umd.py -f -p -s -b <בליסטי> -c
      הערה: ערכי ברירת המחדל הם: -b 100 –s 1 – c 20.
3. התווה את ה- MSD באמצעות תוכנה מבוססת גיליון אלקטרוני (איור 6). בייצוג יומן רישום של MSD לעומת זמן, זהה את שינוי השיפוע. הפרד את החלק הראשון, בדרך כלל קצר, המייצג את המשטר הבליסטי , כלומר שימור מהירות האטומים לאחר התנגשויות. החלק השני הארוך יותר מייצג את המשטר הדיפוזי , כלומר פיזור מהירות האטומים לאחר התנגשויות.
4. חשב את מקדמי הדיפוזיה מהשיפוע של MSD כ:
   
   כאשר Z הוא מספר דרגות החופש (Z = 2 עבור דיפוזיה במישור, Z = 3 עבור דיפוזיה בחלל), ו- t הוא שלב הזמן.

5. פונקציות מתאם זמן

1. חשב את פונקציות מתאם הזמן כמדד לאינרציה של המערכת באמצעות הנוסחה הכללית:
   
   A יכול להיות מגוון משתנים תלויי זמן, כגון עמדות אטומיות, מהירויות אטומיות, לחצים, קיטוב וכו ', כל אחד מניב – באמצעות קשרי גרין-קובו12,13 – תכונות פיזיות שונות, לפעמים לאחר טרנספורמציה נוספת.
2. לנתח את המהירויות האטומיות כדי להשיג את ספקטרום הרטט של הנוזל ואת הביטוי החלופי של מקדמי דיפוזיה עצמית אטומית.
   1. הפעל את קובץ ה- script vibr_spectrum_umd.py כדי לחשב את פונקציית המתאם האוטומטי של המהירות האטומית (VAC) עבור כל סוג אטומי ולבצע את התמרת פורייה המהירה שלו. בסוג שורת הפקודה:
      vibr_spectrum_umd.py <>
      כאשר –t היא הטמפרטורה שהמשתמש חייב להגדיר. התסריט מדפיס שני קבצים: . vels.scf.dat הקובץ עם הפונקציה VAC עבור כל סוג אטומי ו . קובץ .dat vibr עם ספקטרום הרטט נרקב על כל מין אטומי ואת הערך הכולל.
   2. פתח וקרא את .dat vels.scf. התווה את הפונקציה VAC מהקובץ .dat vels.scf באמצעות תוכנה דמוית גיליון אלקטרוני.
   3. שמור את החלק האמיתי של פורייה VAC. זה מה שמניב את ספקטרום הרטט, כפונקציה של תדירות:
      
      איפה מ' המסות האטומיות.
   4. התווה את ספקטרום הרטט מהקובץ .dat הגיבר באמצעות תוכנה דמוית גיליון אלקטרוני (איור 7). זהה את הערך הסופי ב - ω = 0 המתאים לאופי הדיפוזי של הנוזל ולפסגות השונות של הספקטרום בתדר סופי. זהה את ההשתתפות של כל סוג אטומי בספקטרום הרטט.
      הערה: הפירוק של סוגי האטום מראה כי לאטומים שונים יש תרומות שונות של ω= 0, המתאימות למקדמי הדיפוזיה שלהם. הצורה הכללית של הספקטרום היא הרבה יותר חלקה עם פחות תכונות מאשר עבור מוצק מתאים.
   5. בקליפה, קרא את האינטגרל מעל ספקטרום הרטט, אשר מניב את מקדמי הדיפוזיה עבור כל מין אטומי.
      הערה: תכונות תרמודינמיות ניתן להשיג על ידי שילוב מן הספקטרום הרטט, אבל התוצאות יש להשתמש בזהירות בגלל שני קירובים: השילוב תקף בתוך קירוב מעין הרמוני, אשר לא בהכרח להחזיק בטמפרטורות גבוהות; ויש להשליך את החלק דמוי הגז בספקטרום המתאים לפיזור. לאחר מכן האינטגרציה צריכה להיעשות רק על פני החלק דמוי הסריג של הספקטרום. אבל הפרדה זו דורשת בדרך כלל כמה שלבים וחישובים נוספים לאחר ^{העיבוד14}, שאינם מכוסים על ידי חבילת UMD הנוכחית.
3. הפעל את סקריפט viscosity_umd.py כדי לנתח את המתאם העצמי של הרכיבים מתח טנזור כדי להעריך את הצמיגות של להמיס. בסוג שורת הפקודה:
   viscosity_umd.py -f -i -s -o -l
   הערה: תכונה זו היא גישוש ויש לנקוט בכל התוצאות בזהירות. ראשית, בדוק ביסודיות את ההתכנסות של הצמיגות ביחס לאורך הסימולציה.
   1. לגזור את הצמיגות של הנוזל מן המתאם העצמי של טנסור ^הלחץ15 כמו:
      
      כאשר V ו- T הם הנפח והטמפרטורה בהתאמה, _κB הוא קבוע בולטמן σ _ij המרכיב האלכסוני של הלחץ-טנסור, לידי ביטוי בקואורדינטות קרטזיות.
   2. השתמש בהתאמה נאותה יותר כדי לקבל הערכה חזקה יותר של ^{הצמיגות15,16} ולהימנע מהרעש של פונקציית המתאם האוטומטי של טנסור הלחץ שעשוי לנבוע מהגודל הסופי ומשך הזמן הסופי של הסימולציות. עבור פונקציית המתאם האוטומטי של טנסור הלחץ, השתמש בטופס הפונקציונלי ^הבא15,16 המניב תוצאות טובות:
      
      כאשר A, B, τ1, τ2 ו- ω הם פרמטרים מתאימים. לאחר האינטגרציה, הביטוי עבור הצמיגות הופך:
      

6. פרמטרים תרמודינמיים הנובעים מההדמיות.

1. הפעל averages.py כדי לחלץ את הערכים הממוצעים ואת ההתפשטות (כסטיית תקן) עבור לחץ, טמפרטורה, צפיפות ואנרגיה פנימית מקבצי umd . בסוג שורת הפקודה:
   averages.py
   עם -s 0 כברירת מחדל.
2. חשב את השגיאה הסטטיסטית של הממוצע, באמצעות שיטות חסימה.
   הערה: ישנם טעמים שונים של שיטה זו. בעקבות עבודתם של אלן ו Tildesley2, זה נפוץ בממוצע על רצפים של בלוקי זמן, של אורך ארוך יותר ויותר, ולהעריך את סטיית התקן ביחס לממוצע האריתמטי17. ניתן להגיע להתכנסות בגבול של גדלי בלוק רבים וארוכים מספיק, כאשר הדגימה אינה מאומתת. למרות שערך הסף בפועל עבור ההתכנסות בדרך כלל צריך להיבחר באופן ידני.
   1. השתמש בשיטת ^{החצייה18}: החל מדגם הנתונים הראשוני, בכל שלב κ, חצי את מספר הדגימות על ידי ממוצע על כל שתי דגימות תואמות רצופות מהשלב הקודם κ−1:
      
   2. הפעל את קובץ ה- Script fullaverages.py כדי לבצע את הניתוח הסטטיסטי המלא, כולל השגיאה של הממוצע. בסוג שורת הפקודה:
      fullaverages.py -u <>
      הערה: קובץ ה- Script הופך לאוטומטי עד כדי חיפוש כל קבצי .umd.dat בספריה הנוכחית וביצוע הניתוח עבור כולם. ברירות המחדל הן -s 0 –u 0. עבור -u 0 הפלט הוא מינימלי, עבור -u 1 פלט הוא במלואו, עם מספר יחידות חלופיות מודפסות החוצה. קובץ Script זה דורש תמיכה גרפית, מכיוון שהוא יוצר תמונה גרפית לבדיקת ההתכנסות להערכת השגיאה בממוצע.

תוצאות

פירוליט הוא דגם רב-רכיבים של סיליקט נמס (0.5Na2O 2CaO 1.5Al2O3 4FeO 30MgO 24SiO2) שמשווה בצורה הטובה ביותר את הרכב כדור הארץ של סיליקט בתפזורת – הממוצע הגיאוכימי או כל כוכב הלכת שלנו, למעט הליבה מבוססת הברזל ^שלו19. כדור הארץ המוקדם נשלט על ידי סדרה של אירועי התכה בקנה מידה ^גדול20, האחרון אולי בלע את כדור הארץ כולו, לאחר עיבויו לדיסק הפרוטולונארי21. פירוליט מייצג את הקירוב הטוב ביותר להרכב של אוקיינוסי מאגמה בקנה מידה פלנטרי כזה. כתוצאה מכך, חקרנו בהרחבה את התכונות הפיזיות של פירוליט להמיס בטווח טמפרטורות 3,000\u20125,000 K ו 0\u2012150 טווח לחץ GPa מן סימולציות ab initio מולקולרית-דינמיקה ביישום VASP. תנאים תרמודינמיים אלה מאפיינים לחלוטין את תנאי האוקיינוס המגמה הקיצוניים ביותר בכדור הארץ. המחקר שלנו הוא דוגמה מצוינת לשימוש מוצלח בחבילת UMD לכל ניתוח מעמיק של ^המסות22. חישבנו את ההתפלגות ואת אורכי הקשר הממוצעים, עקבנו אחר השינויים בתיאום המזון הקטיונלי והשווינו את התוצאות שלנו למחקרים ניסיוניים וחישוביים קודמים על סיליקטים אמורפיים של קומפוזיציות שונות. הניתוח המעמיק שלנו עזר לפרק מספרי תיאום סטנדרטיים לבוחרים הבסיסיים שלהם, לתאר את נוכחותו של פוליהדרה קואורדינציה אקזוטית בהמסה, ולחלץ תקופות חיים לכל הפוליהדרה הקואורדינציה. הוא גם תיאר את החשיבות של דגימה בסימולציות הן מבחינת אורך המסלול והן מבחינת מספר האטומים הקיימים במערכת המעוצבת. באשר לאחר העיבוד, ניתוח UMD אינו תלוי בגורמים אלה, עם זאת, יש לקחת בחשבון בעת פרשנות התוצאות המסופקות על-ידי חבילת UMD. כאן, אנו מראים כמה דוגמאות כיצד חבילת UMD יכול לשמש כדי לחלץ כמה תכונות אופייניות של נמס, עם יישום pyrolite מותך.

פונקציית התפלגות הזוג Si-O המתקבלת מהתסריט gofrs_umd.py מראה כי הרדיוס של ספירת התיאום הראשונה, שהוא המינימום הראשון של פונקציית g(r), טמון סביב 2.5 אנגסטרומים ב- T = 3000 K ו- P = 4.6 GPa. המקסימום של g(r) טמון ב 1.635 Å - זהו הקירוב הטוב ביותר לאורך העיקול. הזנב הארוך נובע מהטמפרטורה. באמצעות מגבלה זו כמרחק הקשר של Si-O, ניתוח ההפרטציה מראה שיחידות _SiO4, שיכולות להימשך עד כמה פיקוסקונדים, שולטות בהמסה (איור 5). יש חלק חשוב של המסה המציג פילמור חלקי, כפי משתקף על ידי נוכחות של dimers כמו _Si2O7, ושלישים כמו יחידות _Si3Ox. החיים המקבילים שלהם הם בסדר של picosecond. לפולימרים מסדר גבוה יותר יש אורך חיים קצר בהרבה.

איור 5: חיי המינים הכימיים של סי-או.
הדגימה מזוהה בהמסה מרובת רכיבים ב-4.6 GPa ו-3,000 אלף. התוויות מסמנות את המונומרים _SiO3, _SiO4 ו- _SiO5 ואת הפולימרים השונים של _SixOy. אנא לחץ כאן כדי להציג גירסה גדולה יותר של איור זה.

הערכים השונים של השלבים האנכיים והאופקיים, המוגדרים על-ידי הדגלים -z ו- -v לעיל, מניבים דגימות שונות של ה- MSD (איור 6). אפילו ערכים גדולים של z ו- v מספיקים כדי להגדיר את המדרונות ובכך את מקדמי הדיפוזיה של האטומים השונים. הרווח בזמן עבור לאחר העיבוד הוא מדהים כאשר הולכים לערכים גדולים של z ו- v. ה- MSD מציע קריטריון אימות חזק מאוד לאיכות הסימולציות. אם חלק הדיפוזיה של MSD אינו ארוך מספיק, זהו סימן לכך שהסימולציה קצרה מדי, ואינה מצליחה להגיע למצב הנוזלי במובן הסטטיסטי. הדרישה המינימלית לחלק המתפזר של MSD תלויה מאוד במערכת. אפשר לדרוש שכל האטומים ישנו את האתר שלהם לפחות פעם אחת במבנה ההיתוך כדי שהוא ייחשב ^כנוזל10. דוגמה מצוינת עם יישומים במדעים פלנטריים היא סיליקט מורכב נמס בלחצים גבוהים קרוב או אפילו מתחת לקו הנוזל ^שלהם11. האטומים של Si, הציטוטים העיקריים היוצרים רשת, מחליפים אתרים אחרי יותר משני תריסר פיקוסקודים. סימולציות קצרות מסף זה יהיו במידה ניכרת מתחת לדגימה של המרחב התצורה האפשרי. עם זאת, כמו אניונים תיאום, כלומר אטומי O, לנוע מהר יותר מאשר אטומי Si המרכזיים, הם יכולים לפצות על חלק הניידות האיטית של Si. ככזה, המערכת כולה אכן יכולה לכסות דגימה טובה יותר של המרחב התצורלי מאשר להניח רק מן עקירות Si.

איור 6: תזוזות מרובעות ממוצעות (MSD).
ה- MSD מאויר עבור מספר סוגים אטומיים של המסת סיליקט מרובה רכיבים. הדגימה עם שלבים אופקיים ואנכיים שונים, z ו- v, מניבה תוצאות עקביות. עיגולים מוצקים: -z 50 –v 50. עיגולים פתוחים: -z 250 –v 500. אנא לחץ כאן כדי להציג גירסה גדולה יותר של איור זה.

לבסוף, פונקציות VAC האטומיות מניבות את ספקטרום הרטט של ההמסה. איור 7 מראה את הספקטרום באותו תנאי לחץ וטמפרטורה כמו לעיל. אנו מייצגים את התרומות של אטומי Mg, Si ו- O, כמו גם את הערך הכולל. בתדר אפס יש ערך סופי של הספקטרום, אשר מתאים לאופי הדיפוזיה של ההיתוך. הפקת המאפיינים התרמודינמיים בספקטרום הרטט צריכה להסיר את התו הדיפוזי דמוי הגז הזה מאפס, אך גם לקחת בחשבון כראוי את הריקבון שלה בתדרים גבוהים יותר.

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

Discussion

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

הקבצים .dat umd מתאימים להרכבים המשמרים את מספר החלקיקים ללא שינוי לאורך הסימולציה. חבילת UMD יכולה לקרוא קבצים הנובעים מחישובים שבהם הצורה והנפח של תיבת הסימולציה משתנים. אלה מכסים את החישובים הנפוצים ביותר, כמו NVT ו- NPT, שם מספר החלקיקים, N, טמפרטורה T, נפח, V, ו / או לחץ, P, נשמרים קבועים.

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

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

Materials

Name	Company	Catalog Number	Comments
getopt library			open-source
glob library			open-source
matplotlib library			open-source
numpy library			open-source
os library			open-source
Python software	The Python Software Foundation	Version 2 and 3	open-source
random library			open-source
re library			open-source
scipy library			open-source
subprocess library			open-source
sys library			open-source