Luận văn Nghiên cứu sự hình thành pha tinh thể và thủy tinh của hạt nano Feb bằng phương pháp mô phỏng

÷ 106 nguyên tử đã được xây dựng bằng chương trình thống kê hồi phục (TKHP) song song có HPBXT phù hợp với thực nghiệm. Mô hình 105 và 2×105 nguyên tử của hợp kim CoxB100-x và FeyP100-y VĐH, đã được xây dựng để nghiên cứu đặc trưng lỗ trống và đám lỗ trống [11]. Tuy nhiên, với nghiên cứu cấu trúc vi mô của một số vật liệu kim loại VĐH, đặc trưng nút khuyết tự nhiên, lỗ trống và vai trò của chúng trong cơ chế khuếch tán, sự ảnh hưởng của quá trình hồi phục, nhiệt độ đến các tính chất của hệ, số liệu tính toán được từ các mô hình kích thước khác nhau 104 đến 2×105 nguyên tử là gần như nhau. Trên cơ sở khảo sát kỹ lưỡng các yếu tố nêu trên chúng tôi sử dụng phương pháp ĐLHPT, ở nhiệt độ 0K để nghiên cứu cấu trúc và tính chất vật liệu kim loại VĐH. Để thực hiện mục đích đó chúng tôi xây dựng chương trình TKHP (xem chương 2). Chương trình TKHP là chương trình ĐLHPT ở nhiệt độ T = 0K và do đó động năng của hệ bằng 0. 21 1.3. Lý thuyết cổ điển về mầm và sự phát triển mầm Sự hình thành một thực thể tinh thể mới từ sự hòa tan bắt đầu thông qua quá trình mầm. Mầm tinh thể được định nghĩa như là một hệ các quá trình mà những nguyên tử hoặc phân tử của pha vật chất sắp xếp thành đám của pha khác nó đủ lớn để có thể phát triển đến kích thước lớn hơn. Đám được định nghĩa như là các nhân hoặc nhân có kích thước tới hạn [16]. Mầm có thể là đồng nhất với sự vắng mặt của các hạt ngoài hoặc sự tinh thể trong chất hòa tan. Hoặc có thể không đồng nhất trong sự hiện diện của các hạt ngoài trong chất hòa tan. Cả hai loại mầm được biết nhìn chung như các mầm chính. Mầm loại hai có không gian khi mầm phát triển bởi sự hiện diện của các tinh thể trong chất như vật. Trạng thái siêu bão hòa: Lực điều khiển cần cho các mầm và sự phát triển của tinh thể được gọi như là siêu bão hòa và được định nghĩa như sự khác nhau về thế hóa học giữa phân tử trong chất hòa tan và thế hóa học trong khối của pha tinh thể. 𝜇 = 𝜇𝑠 = 𝜇𝑐 (1.8) trong đó μS là thế hóa học của phân tử trong chất hòa tan và μC là thế hóa học của phân tử trong khối tinh thể. Theo nhiệt động lực học phương trình (1.8) có thể được biểu diễn như sau: 𝜇 = 𝑘𝑇𝑙𝑛𝑆 (1.9) Ở đây k là hệ số Boltsmann, T là nhiệt độ tuyệt đối và S là tỉ lệ siêu bão hòa. Khi μ > 0 chất hòa tan được gọi là siêu bão hòa, nó có nghĩa là mầm hoặc sự tăng trưởng là có thể xảy ra, ngược lại μ < 0 chất hòa tan sẽ không thể siêu bão hòa và sự không hòa tan cần không gian. Sự tạo thành tỉ lệ siêu bão hòa sẽ thay đổi phụ thuộc vào điều kiện của hệ (cụ thể là khí/rắn, hòa tan/rắn, nóng chảy/ rắn). Đối với các mầm và sự tăng trưởng từ các chất hòa tan nó có dạng sau: 22 𝑆 = ∏𝑎𝑖 𝑛𝑖 ∏𝑎𝑖,𝑒 𝑛𝑖 (1.10) Ở đây ni là số ion thứ i trong phân tử của tinh thể, và ai và ai,e là sự hoạt động thực tế và cân bằng của phân tử thứ i trong tinh thể. Dạng năng lượng: Theo lý thuyết mầm, công cần để tạo thành đám gồm n phân tử bằng sự khác nhau giữa năng lượng tự do của hệ trong trạng thái cuối và đầu của hệ, hơn một số hạng liên quan đến sự tạo thành bề mặt giữa các nhân và chất hòa tan. Năng lượng tự do này có thể được biểu diễn bởi (giả thiết nhân có dạng cầu) 𝐺𝑇 = −𝑛∆𝜇 + 4𝜋𝑟 2𝜎 (1.11) Ở đây r là bán kính của nhân và σ là năng lượng tự do bề mặt. Nếu mối phân tử trong tinh thể chiếm một thể tích V, thì mỗi nhân sẽ chứa 4/3π.r3/V phân tử. Thì phương trình (1.11) biến thành dạng sau: 𝐺𝑇 = − 4 3 𝜋 𝑟3 𝑉 ∆𝜇 + 4𝜋𝑟2𝜎 (1.12) Hình 1.1a cho thấy hình vẽ GT như là hàm của bán kính r. Có thể thấy hàm GT đạt cực đại như thế nào, nó thể hiện độ cao rào năng lượng cần để vượt qua đối với mầm (GT*). Giá trị r ở cực đại này r* được định nghĩa là bán kính tới hạn hoặc kích thước các nhân. Giá trị này được xác định bởi công thức: 𝑟∗ = 2𝜎𝑉 𝑘𝑇𝑙𝑛𝑆 (1.13) Có thể chứng minh rằng r* giảm (cũng như G*) khi trạng thái siêu bão hòa tăng (1.13), nghĩa là khả năng có mầm trong hệ xét sẽ càng cao thì trạng thái siêu bão hòa càng cao. Tốc độ mầm: Tốc độ của các mầm (cụ thể là số nhân được tạo thành trên một đơn vị thời gian trên một đơn vị thể tích) có thể được tính bởi loại phương trình Arrhenius: 23 Hình 1.1. Tổng năng lượng tự do theo kích thước của đám (a); Tốc độ mầm như là hàm của trạng thái siêu bão hòa (b); μC cho thấy kích thước tới hạn của trạng thái siêu bão hòa. *G J Aexp kT   (1.14) Ở đây A cũng phụ thuộc vào trạng thái siêu bão hòa. Hình 1.1 cho thấy J như là hàm của trạng thái siêu bão hòa S được thể hiện trên hình 1.1.b. Có thể thấy trong hình vẽ này rằng tốc độ mầm hầu như bằng 0 cho đến tận giá trị tới hạn của trạng thái siêu bão hòa đạt được, sau đó tốc độ tăng theo hàm e mũ. Trạng thái siêu bão hòa tới hạn ΔμC xác định ở vùng phát triển tinh thể có thể thực hiện. 24 Chương 2 PHƯƠNG PHÁP NGHIÊN CỨU Hiện nay mô phỏng đang được thừa nhận như là một phương pháp thứ ba để nghiên cứu các hệ vĩ mô cũng như vi mô. Đây là một phương pháp đóng vai trò liên kết chặt chẽ giữa hai phương pháp lý thuyết và thực nghiệm. Trong chương này chúng tôi trình bày: Phương pháp mô phỏng ĐLHPT, phương pháp TKHP; Kỹ thuật tính toán các đặc trưng cấu trúc như: hàm phân bố xuyên tâm (HPBXT), phân bố số phối trí (SPT); Xây dựng hạt nano vô định hình và xác định HPBXT của mẫu hạt nano ở nhiệt độ khác nhau; Cuối cùng trình bày phương pháp xác định mầm tinh thể và thuật toán. 2.1. Phương pháp động lực học phân tử Phương pháp động lực học phân tử là một công cụ cho phép chúng ta xây dựng mô hình vật liệu nano dựa trên hệ phương trình chuyển động của Newton. Phương trình chuyển động được khảo sát với vận tốc chuyển động của hạt tính bằng thuật toán Verlet theo bước thời gian dt. Xét một hệ gồm N nguyên tử được gieo vào khối hình lập phương cạnh L (hoặc vùng không gian hình cầu, R). Tọa độ ban đầu của nguyên tử có thể lấy ngẫu nhiên nhưng phải thỏa mãn điều kiện không có bất kỳ hai nguyên tử nào quá gần nhau. Dưới tác dụng của lực tương tác, các nguyên tử sẽ dịch chuyển dần đến vị trí cân bằng. Trạng thái cân bằng của mô hình được xác định bởi nhiệt độ và áp suất. Chuyển động của các nguyên tử trong mô hình tuân theo định luật cơ học cổ điển Newton. Đối với hệ gồm N hạt, phương trình chuyển động của định luật hai Newton có thể viết như sau: 25 Phương pháp động lực học phân tử dựa trên phương trình chuyển động Newton: �⃗�𝑖 = 𝑚𝑖�⃗�𝑖 (2.1) 𝑚𝑖�⃗�𝑖 = 𝑚𝑖 𝑑2𝑟𝑖 𝑑𝑡2 = �⃗�𝑖(𝑟1, , 𝑟𝑁) (2.2) trong đó, 𝐹𝑖⃗⃗⃗ là lực tổng hợp tác dụng lên nguyên tử thứ 𝑖 từ các nguyên tử còn lại; 𝑚𝑖 , 𝑎𝑖 và 𝑟𝑖⃗⃗ lần lượt là khối lượng, gia tốc và vị trí của nguyên tử thứ 𝑖. Lực 𝐹𝑖⃗⃗⃗ được xác định theo công thức: �⃗�𝑖 = −∑ 𝜕𝑈𝑖𝑗 𝜕𝑟𝑖𝑗 𝑁 𝑗=1 (2.3) trong đó, 𝑈𝑖𝑗 là thế tương tác giữa nguyên tử thứ 𝑖 và nguyên tử thứ 𝑗 và 𝑟𝑖𝑗 là khoảng cách giữa chúng. Trong mô phỏng động lực học phân tử, ta sử dụng thuật toán Verlet để giải hệ phương trình chuyển động của các nguyên tử theo định luật hai Newton. Với thuật toán này, tọa độ của nguyên tử ở thời điểm (𝑡 + 𝑑𝑡) được xác định thông qua tọa độ của nó ở hai thời điểm 𝑡 và (𝑡 − 𝑑𝑡) bằng biểu thức: 𝑟𝑖⃗⃗ (𝑡 + 𝑑𝑡) = 2𝑟𝑖⃗⃗ (𝑡) − 𝑟𝑖⃗⃗ (𝑡 − 𝑑𝑡) + (𝑑𝑡) 2 �⃗�𝑖(𝑡) 𝑚𝑖 (2.4) Vận tốc ở thời điểm 𝑡 được xác định thông qua tọa độ ở thời điểm (𝑡 − 𝑑𝑡) và (𝑡 + 𝑑𝑡) theo biểu thức: 𝑣𝑖⃗⃗⃗ ⃗(𝑡) = 𝑟𝑖⃗⃗ (𝑡 + 𝑑𝑡) − 𝑟𝑖⃗⃗ (𝑡 − 𝑑𝑡) 2𝑑𝑡 (2.5) Lực 𝐹𝑖⃗⃗⃗(𝑡) được phân tích theo ba thành phần tương ứng với các phương Ox, Oy và Oz của hệ tọa độ Đề các: 26 �⃗�𝑖(𝑡) = 𝐹𝑥𝑖 ⃗⃗ ⃗⃗ ⃗ + 𝐹𝑦𝑖 ⃗⃗ ⃗⃗ ⃗ + 𝐹𝑧𝑖 ⃗⃗⃗⃗⃗ = ∑𝐹𝑥𝑖𝑗 ⃗⃗⃗⃗⃗⃗⃗ 𝑗 +∑𝐹𝑦𝑖𝑗 ⃗⃗⃗⃗⃗⃗⃗ 𝑗 +∑𝐹𝑧𝑖𝑗 ⃗⃗ ⃗⃗ ⃗⃗ 𝑗 (2.6) trong đó ijx j F được xác định như sau: ∑𝐹𝑥𝑖𝑗 ⃗⃗⃗⃗⃗⃗⃗ 𝑗 = 𝑥0⃗⃗⃗⃗⃗∑(− 𝜕𝑈(𝑟𝑖𝑗) 𝜕𝑟𝑖𝑗 ) 𝑥𝑖 − 𝑥𝑗 𝑟𝑖𝑗 (2.7) với 𝑥0⃗⃗⃗⃗⃗ là véctơ đơn vị của trục Ox. Các thành phần ijy j F , ijz j F được xác định tương tự như phương trình (2.7). Khi nghiên cứu các mô hình vật liệu bằng phương pháp ĐLHPT, tùy theo mục đích cần nghiên cứu mà người ta chọn một trong các mô hình sau đây: mô hình NVE, NVT, NPH, NTP, 𝜇TV và 𝜇TP. Trong đó: N, E, V, T, P, H và 𝜇 lần lượt là số nguyên tử, năng lượng toàn phần, thể tích, nhiệt độ, áp suất, entanpy và thế hóa học. Đối với mô hình NVE thì các đại lượng N, V và E không đổi trong suốt thời gian mô phỏng. Còn đối với các mô hình khác sẽ có đại lượng tương ứng không thay đổi. Trong quá trình mô phỏng ĐLHPT, U và K lần lượt là thế năng và động năng của hệ và được tính theo biểu thức sau: 𝑈 =∑𝑈𝑖𝑗(𝑟𝑖𝑗) 𝑖>𝑗 (2.8) 𝐾 =∑ 𝑚𝑖𝑣𝑖 2 2 𝑁 𝑖=1 =∑ 𝑚𝑖 2 ( 𝑟𝑖(𝑡 + 𝑑𝑡) − 𝑟𝑖(𝑡 − 𝑑𝑡) 2𝑑𝑡 ) 2 (2.9) 𝑁 𝑖=1 Năng lượng E của hệ có thể tính theo công thức: 𝐸 = 𝐾 + 𝑈 (2.10) 27 Nhiệt độ của mô hình ĐLHPT có thể được xác định thông qua động năng của hệ theo công thức: 𝑇 = 𝐾 2 3𝑁𝑘𝐵 (2.11) Trong đó 𝑘𝐵 là hằng số Boltzman. Động năng của hệ được xác định thông qua vận tốc của các nguyên tử theo công thức: 𝐾 =∑ 𝑚𝑖𝑣𝑖 2 2 𝑁 𝑖=1 =∑ 𝑚𝑖 2 ( 𝑟𝑖(𝑡 + 𝑑𝑡) − 𝑟𝑖(𝑡 − 𝑑𝑡) 2𝑑𝑡 ) 2 (2.12) 𝑁 𝑖=1 Trong mô hình NVT, để giữ nhiệt độ có giá trị không đổi người ta thường sử dụng kỹ thuật điều chỉnh nhiệt độ (Temperature Scaling). Ý tưởng của thuật toán này là điều chỉnh vận tốc của tất cả các hạt bởi một thừa số được xác định bởi tỷ số giữa nhiệt độ mong muốn và nhiệt độ hiện tại được xác định từ phương trình (2.11). Giả sử nhiệt độ được tính từ phương trình là T, nhiệt độ mong muốn của hệ đạt được là 𝑇0, điều chỉnh vận tốc 𝑣𝑖 của tất cả các nguyên tử theo phương trình sau: 𝑣𝑖 ′ = √ 𝑇0 𝑇 . 𝑣𝑖 (2.13) Chúng ta sẽ thu được: 𝑇′ = 1 3𝑘𝐵𝑁 ∑𝑚𝑖(𝑣𝑖 ′)2 𝑁 𝑖=1 = 1 3𝑘𝐵𝑁 ∑𝑚𝑖 𝑇0 𝑇 𝑣𝑖 2 𝑁 𝑖=1 = 𝑇0 (2.14) Chọn áp suất của mô hình ĐLHPT có thể được điều chỉnh thông qua kích thước của mô hình. Mô hình NPT sẽ điều chỉnh áp suất P thông qua việc nhân tọa độ của tất cả các nguyên tử với thừa số điều chỉnh 𝜆. Khi áp suất của hệ nhỏ hơn giá trị cho phép, ta sẽ chọn 𝜆 > 1, và ngược lại nếu áp suất lớn hơn 28 giá trị cho trước ta chọn 𝜆 < 1. Trong chương trình, áp suất được điều chỉnh như sau: Nhập giá trị áp suất 𝑃𝑚ớ𝑖, nếu 𝑃𝑚ớ𝑖 > 𝑃ℎệ thì 𝜆 = 1 − 𝑑𝑃, ngược lại 𝜆 = 1 + 𝑑𝑃, với giá trị 𝑑𝑃 được chọn là 10−4. Do vậy, tọa độ mới của các nguyên tử được xác định: 𝑥𝑎 ′ [𝑖] = 𝑥𝑎[𝑖]. 𝜆; 𝑦𝑎 ′ [𝑖] = 𝑦𝑎[𝑖]. 𝜆; 𝑧𝑎 ′ [𝑖] = 𝑧𝑎[𝑖]. 𝜆 𝑥𝑏 ′ [𝑖] = 𝑥𝑏[𝑖]. 𝜆; 𝑦𝑏 ′ [𝑖] = 𝑦𝑏[𝑖]. 𝜆; 𝑧𝑏 ′ [𝑖] = 𝑧𝑏[𝑖]. 𝜆 (2.15) Khi đó, kích thước mô hình sẽ có giá trị 𝐿′ = 𝐿𝜆, hoặc R’ = R𝜆 Khi xây dựng mô hình ĐLHPT, các thông số nhiệt độ và áp suất ở thời điểm 𝑡 được xác định như sau: 3 2 𝑁𝑘𝑇(𝑡) = 𝐾(𝑡) (2.16) 𝑇(𝑡) = 2 3 𝐾(𝑡) 𝑘𝐵𝑁 = 1 3𝑁𝑘𝐵 ∑𝑚𝑖(𝑣𝑖(𝑡)) 2 (2.17) 𝑁 𝑖=1 𝑃(𝑡) = 𝑁 𝑉 𝑘𝑇(𝑡) + 1 3𝑉 ∑𝑟𝑖𝑗(𝑡)𝐹𝑖𝑗(𝑡) 𝑖<𝑗 (2.18) Các mô hình mô phỏng NVE, NPT, NVT được sử dụng trong mô phỏng hạt nano. Mô hình NVE cô lập với môi trường bên ngoài do vậy hầu như không chịu tác động của ngoại lực, đây là mô hình có thể sử dụng để khảo sát sự dịch chuyển của các nguyên tử mô hình và từ đó có thể tính được hệ số sự khuếch tán của các nguyên tử. Nhược điểm của mô hình NVE là đẻ khảo sát ở nhiệt độ T và áp suất P cho trước ta phải thực hiện một số rất lớn các bước lặp ĐLHPT, do đó thời gian mô phỏng kéo sẽ dài. Để khắc phục nhược điểm trên, ban đầu chúng ta mô phỏng theo mô hình NPT hoặc NVT để đạt được các thông số T và P đã cho. Bước tiếp theo, thực hiện mô phỏng theo NVE, do đó thời gian mô phỏng sẽ giảm đi rất nhiều. Sơ đồ khối thực hiện phương pháp động lực học phân tử được mô tả trong hình 2.1. 29 Bắt đầu - Đọc các hệ số đặc trưng: T, n, 𝜌, dt, - Chọn các tọa độ và vận tốc ban đầu cho các hạt. - Tính lực tác dụng lên toàn bộ các hạt. - Lấy tích phân các phương trình chuyển động của Newton. - Để các nguyên tử chuyển động tự do dưới tác dụng của lực. 𝑘 < 𝑘𝑚𝑎𝑥? Xác định giá trị trung bình của các đại lượng cần khảo sát. Thực hiện kết quả Kết thúc Đúng Không đúng 𝑘 = 𝑘 + 1 Hình 2.1. Sơ đồ khối phương pháp ĐLHPT 30 2.2. Phương pháp thống kê hồi phục Phương pháp TKHP chính là phương pháp xây dựng mô hình vật liệu nguyên tử ở 0K. Ở nhiệt độ 0K động năng của hệ bằng không và tổng năng lượng của hệ bằng thế năng U của hệ. Điều này có nghĩa là hệ đạt trạng thái cân bằng khi thế năng U của hệ đạt tới giá trị ổn định. Trong phương pháp thống kê phục hồi sau khi tính được 𝐹𝑖 thì từng nguyên tử được dịch đi một khoảng 𝑑𝑟 cho trước theo hướng của lực 𝐹𝑖, trong đó khoảng 𝑑𝑟 là bước dịch chuyển. Giá trị của 𝑑𝑟 ≪ 𝑟0, 𝑟0 là bán kính của nguyên tử. Giả sử nguyên tử thứ 𝑖 trong mô hình được dịch chuyển theo phương của lực tổng hợp từ các nguyển tử còn lại tác dụng lên nó. Lực tổng hợp 𝐹𝑖⃗⃗⃗ tác dụng lên nguyên tử thứ 𝑖 được tính theo công thức: �⃗�𝑖(𝑡) = −∑( 𝜕𝜑(𝑟𝑖𝑗) 𝜕𝑟𝑖𝑗 . 𝑟𝑖𝑗 𝑟𝑖𝑗 ⃗⃗ ⃗⃗ ) (2.19) 𝑖 Trong đó 𝑟𝑖𝑗 là khoảng cách giữa hai nguyên tử 𝑖 và 𝑗; 𝜑(𝑟𝑖𝑗) là thế tương tác giữa các nguyên tử. Lực �⃗�𝑖(𝑡) được phân tích theo 3 thành phần tương ứng với 3 trục 𝑂𝑥, 𝑂𝑦, 𝑂𝑧 của tọa độ Đề các: �⃗�𝑖(𝑡) = 𝐹𝑥𝑖 ⃗⃗ ⃗⃗ ⃗ + 𝐹𝑦𝑖 ⃗⃗ ⃗⃗ ⃗ + 𝐹𝑧𝑖 ⃗⃗⃗⃗⃗ =∑(𝐹𝑥𝑖𝑗 ⃗⃗⃗⃗⃗⃗⃗ 𝑗 + 𝐹𝑦𝑖𝑗 ⃗⃗⃗⃗⃗⃗⃗ + 𝐹𝑧𝑖𝑗 ⃗⃗ ⃗⃗ ⃗⃗ ) (2.20) Trong đó ijx j F được xác định như sau: ∑𝐹𝑥𝑖𝑗 ⃗⃗⃗⃗⃗⃗⃗ 𝑗 = 𝑥0⃗⃗⃗⃗⃗∑(− 𝜕𝑈(𝑟𝑖𝑗) 𝜕𝑟𝑖𝑗 ) 𝑥𝑖 − 𝑥𝑗 𝑟𝑖𝑗 (2.21) Với 𝑥0⃗⃗⃗⃗⃗ là véctơ đơn vị của trục Ox. Các thành phần ijy j F , ijz j F được xác định tương tự như phương trình (2.21). Trong phương pháp TKHP, động năng của mô hình đúng bằng năng lượng của hệ. 31 Hình 2.2. Sơ đồ khối phương pháp TKHP Sau khi lực �⃗�𝑖 được xác định, từng nguyên tử trong mô hình được dịch chuyển đi một khoảng cách 𝑑𝑟 cho trước theo hướng lực tác dụng �⃗�𝑖. Để tăng tốc độ tính toán, trong khoảng 150 bước đầu tiên cho 𝑑𝑟 lớn, khoảng 0.4 Å, sau đó giảm dần và giữ ở mức 0.02 Å cho đến khi hệ đạt trạng thái cân bằng. Sau mỗi bước dịch chuyển, tọa độ của mỗi nguyên tử được xác định lại theo công thức: Bắt đầu Ban đầu gán cho các nguyên tử có tọc độ ngẫu nhiên trong khối lập phương có cạnh dài phù hợp với khối lượng của hệ mà ta cần nghiên cứu . 𝑘 = 1 - Tính lực tác dụng của các nguyên tử xung qunh lên nguyên tử thứ 𝑖. - Dịch chuyển nguyên tử thứ 𝑖 theo chiều lực tác dụng lên một bước dich chuyển 𝑑𝑟 khi đó năng lượng của hệ sẽ giảm. 𝑘 < 𝑘𝑚𝑎𝑥? Xác định giá trị trung bình của các đại lượng cần khảo sát. Thể hiện kết quả Kết thúc Đúng Không đúng 𝑘 = 𝑘 + 1 32 𝑥′(𝑖) = 𝑥(𝑖) + 𝑑𝑟. 𝐹𝑥𝑖 𝐹 (2.22) Trong đó 𝐹𝑥𝑖 và 𝐹 là thành phần lực theo trục Ox và module lực tổng hợp được tính theo các công thức (2.20), (2.21). Các thành phần 𝑦(𝑖) và 𝑧(𝑖) được xác định tương tự như (2.22). Tổng năng lượng của hệ sẽ là: 𝑈 = 1 2 ∑𝜑𝑖𝑗(𝑟) 𝑖<𝑗 (2.23) Sơ đồ khối thực hiện chương trình TKHP được minh họa trên hình 2.2 2.3. Xây dựng hạt nano Như chúng ta đã biết thế tương tác là yếu tố quan trọng ảnh hưởng đến độ tin cậy của các kết quả thu được khi thực hiện các nghiên cứu bằng phương pháp mô phỏng. Vì vậy, việc lựa chọn thế tương tác giữa các nguyên tử thích hợp là cần thiết đối với các mô hình mô phỏng. Mặc dù, đã có một số thế được phát triển để nghiên cứu các hệ kim loại và hợp kim dạng kim loại - á kim VĐH. Tuy nhiên, đối với các mô hình với số lượng nguyên tử lớn và đặc biệt đối với kim loại và hợp kim VĐH với cơ sở Fe, Co và Ni thì thế tương tác cặp vẫn là thế chưa thể thay thế. Trong luận văn này, chúng tôi lựa chọn thế tương tác cặp là thế Pak-Doyama vì tính đơn giản và thực tế cho thấy việc thực hiện thế tương tác này cho phép mô phỏng được vật liệu VĐH có nhiều tính chất phù hợp với số liệu thực nghiệm như: cấu trúc vi mô (HPBXT, TSCT và phân bố SPT), mật độ và một số tính chất vật lý ở các trạng thái hồi phục khác nhau. Thế tương tác này có dạng: 𝜑(𝑟𝑖𝑗) = { 𝑎(𝑟𝑖𝑗 + 𝑏) 4 + 𝑐(𝑟𝑖𝑗 + 𝑑) 2 + 𝑒 𝑘ℎ𝑖 𝑟𝑖𝑗 < 𝑟𝑐𝑢𝑡𝑜𝑓𝑓 0 𝑘ℎ𝑖 𝑟𝑖𝑗 > 𝑟𝑐𝑢𝑡𝑜𝑓𝑓 (2.24) Trong đó 𝜑(𝑟𝑖𝑗) (eV) là thế năng tương tác cặp và rij là khoảng cách giữa các nguyên tử tính bằng Å, rcutoff là bán kính ngắt; các hệ số a, b, c, d và e được xác định theo các số liệu thực nghiệm về module đàn hồi, hằng số mạng như thấy trong bảng 2.1. 33 Bảng 2.1. Hệ số của thế tương tác cặp nguyên tử Pak-Doyama đối với vật liệu kim loại Fe100-xBx VĐH Cặp a (eV/ Å4) b (Å) c (eV/ Å2) d (Å) e (eV) rcắt(Å) Fe-Fe - 0,18892 - 1,82709 1,70192 - 0,50849 - 0,19829 3,44 Fe-B - 0,22407 - 1,47709 2,01855 - 2,15849 - 0,23519 3,09 B-B - 0,08772 - 2,17709 0,79028 - 2,85849 - 0,09208 3,79 Các hạt nano được xây dựng như sau: Trước tiên chúng ta đặt ngẫu nhiên 5000 nguyên tử bên trong một quả cầu với bán kính 28 Å. Cấu hình này được phục hồi để đạt được mức thế năng tối thiểu bằng cách sử dụng phương pháp phục hồi số liệu thống kê hồi phục. Cụ thể là, đối với mỗi nguyên tử, chúng tôi xác định lực tác động lên nó từ các nguyên tử còn lại. Sau đó các nguyên tử chuyển hướng năng lượng đã được xác định bởi tỉ lệ khoảng cách với năng lượng. Quy trình này được thực hiện nhiều lần cho đến khi hệ thống đạt mức năng lượng tối thiểu. Sau đó, chúng tôi hồi phục các mẫu thu được bằng phương pháp mô phỏng MD theo mô hình NVE. Chúng tôi tạo mẫu hạt nano có nhiệt độ từ 700 đến 900K bằng cách nung nóng và phục hồi với 2107 bước. Khi mô phỏng mẫu 900K, số hạt tại khu vực bề mặt đôi khi có động năng lớn do đó chúng có thể tách ra khỏi các hạt nano. Để ngăn chặn điều này, chúng tôi theo dõi khoảng cách giữa mỗi hạt với tâm của hạt nano Ri. Nếu khoảng cách Ri lớn hơn giá trị cố định 28 Å, thì động năng của hạt thứ i được tính bằng không. Quy trình này buộc hạt thứ i trở lại đến khu vực bề mặt nếu nó di chuyển xa hạt nano. 2.4. Xác định hàm phân bố xuyên tâm của hạt nano 2.4.1. Hàm phân bố xuyên tâm Trong mô phỏng vật liệu ở trạng thái VĐH, một đại lượng tuân theo quy tắc thống kê được sử dụng để xác định cấu trúc của vật liệu ở mức nguyên tử, đó là HPBXT. HPBXT có thể được xác định từ thực nghiệm thông qua TSCT. Từ HPBXT, phân bố SPT trung bình và TSCT sẽ được tính toán số trong mô phỏng. 34 Theo cơ học thống kê, hàm tương quan cặp hay HPBXT, 𝑔(𝑟) được xác định như sau: 𝑔(𝑟) = 𝑉 𝑁2 ⟨∑ 𝛿(𝑟𝑖𝑗⃗⃗ ⃗⃗ − 𝑟) 𝑖,𝑗≠𝑖 ⟩ (2.26) trong đó V là thể tích của mẫu vật liệu, N là số nguyên tử chứa trong thể tích V. Phương trình (2.27) có thể viết lại một cách tường minh hơn như sau: 𝑔(𝑟) = 𝑉 𝑁2 ∑ 𝑑𝑟1⃗⃗⃗ ⃗𝑑𝑟2⃗⃗⃗⃗ 𝑑𝑟𝑁⃗⃗⃗⃗⃗𝑃𝑁(𝑟 𝑁)𝛿(𝑟𝑖𝑗⃗⃗ ⃗⃗ − 𝑟) (2.27) 𝑁 𝑖,𝑗≠𝑖 ở đây 𝑟𝑖𝑗⃗⃗ ⃗⃗ = 𝑟𝑖⃗⃗ − 𝑟�⃗⃗�, 𝑟 𝑁 = (𝑟1⃗⃗⃗ ⃗, 𝑟2⃗⃗⃗⃗ , , 𝑟𝑁⃗⃗⃗⃗⃗), 𝑃𝑁(𝑟 𝑁) là hàm phân bố xác suất và 𝑟𝑖⃗⃗ , 𝑟�⃗⃗� là véc tơ toạ độ của các hạt thứ i và j. Véc tơ 𝑟 là một thông số xuất hiện như một biến thực ở vế trái của phương trình (giá trị của 𝑟 do chúng ta chọn). Hàm 𝑔(𝑟) có thể hiểu như là xác suất tìm thấy nguyên tử cách nguyên tử đang xét một véc tơ 𝑟. Đối với hệ đẳng hướng, 𝑔(𝑟) chỉ phụ thuộc vào độ lớn của véc tơ 𝑟. Tích phân (2.27) theo thể tích 𝑉(𝑟, ∆𝑟) giữa 𝑟 và 𝑟 + ∆𝑟 (ở đây 𝑟 = |𝑟|) và giả sử lớp vỏ hình cầu là đủ mỏng, chúng ta sẽ thu được biểu thức: ∫ 𝑑𝑟 𝑉(𝑟,∆𝑟) 𝑔(𝑟) ≈ 4𝜋𝑟2∆𝑟𝑔(𝑟) (2.28) Thay (2.28) vào (2.27) thu được: 𝑔(𝑟) ≈ 𝑉 4𝜋𝑟2∆𝑟𝑁2 ∫ 𝑑𝑟 𝑉(𝑟,∆𝑟) ∑ 𝑑𝑟1⃗⃗⃗ ⃗𝑑𝑟2⃗⃗⃗⃗ 𝑑𝑟𝑁⃗⃗⃗⃗⃗𝑃𝑁(𝑟 𝑁)𝛿(𝑟𝑖𝑗⃗⃗ ⃗⃗ − 𝑟) 𝑁 𝑖,𝑗≠𝑖 = 𝑉 4𝜋𝑟2∆𝑟𝑁2 ∫ 𝑑𝑟 𝑉(𝑟,∆𝑟) ∑𝑑𝑟1⃗⃗⃗ ⃗𝑑𝑟2⃗⃗⃗ ⃗ 𝑑𝑟𝑁⃗⃗ ⃗⃗ 𝑃𝑁(𝑟 𝑁)∑ ∫ 𝑑𝑟𝛿(𝑟𝑖𝑗⃗⃗ ⃗⃗ − 𝑟) 𝑉(𝑟,∆𝑟) 𝑁 𝑖,𝑗≠𝑖 𝑁 𝑖,𝑗≠𝑖 (2.29) Tích phân hàm delta, chúng ta tính được số hạt trong lớp hình cầu là: 𝑛𝑖(𝑟, ∆𝑟) = ∑ ∫ 𝑑𝑟𝛿(𝑟𝑖𝑗⃗⃗ ⃗⃗ − 𝑟) 𝑉(𝑟,∆𝑟)𝑖,𝑗≠𝑖 (2.30) 35 Thay (2.30) vào (2.28) sẽ tìm được: 𝑔(𝑟) ≈ 𝑉 4𝜋𝑟2∆𝑟𝑁2 ∫ 𝑑𝑟1⃗⃗⃗⃗ 𝑑𝑟2⃗⃗⃗⃗ 𝑑𝑟𝑁⃗⃗ ⃗⃗ 𝑃𝑁(�⃗� 𝑁 )∑ 𝑛𝑖(𝑟, ∆𝑟) 𝑖 = 𝑉 4𝜋𝑟2∆𝑟𝑁2 ⟨∑ 𝑛𝑖(𝑟, ∆𝑟) 𝑖 ⟩ (2.31) Phương trình (2.31) có thể viết lại gọn hơn như sau: 𝑔(𝑟 ) = 𝜌(𝑟) 𝜌0 (2.32) trong đó 𝜌0 là mật độ nguyên tử trung bình trong thể tích V của mẫu vật liệu, 𝜌(𝑟)là mật độ nguyên tử ở khoảng cách r tính từ nguyên tử trung tâm. { 𝜌0 = 𝑁 𝑉 𝜌(𝑟) = 1 𝑁 ⟨∑ 𝑛𝑖(𝑟, ∆𝑟)𝑖 ⟩ 4𝜋𝑟2∆𝑟 (2.33) HPBXT cặp trong chương trình mô phỏng được tính toán số bởi: 𝑔𝛼𝛽(𝑟) = 𝑁𝛼𝛽 ⟨∑∑𝛿(𝑟𝑖𝑗⃗⃗ ⃗⃗ − 𝑟) 𝑁𝛽 𝑗 𝑁𝛼 𝑖 ⟩ (2.34) trong đó 𝛼, 𝛽 chỉ các loại nguyên tử và 𝑁𝛼𝛽 = { 𝑁 𝜌0𝑁𝛼(𝑁𝛽 − 1) 𝛼 = 𝛽 𝑁 𝜌0𝑁𝛼𝑁𝛽 𝛼 ≠ 𝛽 (2.35) Ở đây N là tổng số nguyên tử trong mô hình, Nα và Nβ lần lượt là số nguyên tử loại α và β, 𝜌0là mật độ nguyên tử trung bình trong thể tích V. HPBXT tổng cộng được xác định bởi: 𝑔(𝑟) = ∑ 𝑐𝑖𝑏𝑖𝑐𝑗𝑏𝑗𝑔𝛼𝛽(𝑟)𝑖𝑗 [∑ 𝑐𝑖𝑏𝑖𝑖 ]2 (2.36) trong đó ci, cj,bi và bj lần lượt là nồng độ và hệ số tán xạ của nguyên tử loại α và β. 36 Thông tin cũng hết sức quan trọng về cấu trúc địa phương của kim loại và hợp kim VĐH thu được từ phân bố SPT trung bình Zαβ, Zαβ được xác định bằng biểu thức tích phân đỉnh thứ nhất HPBXT tương ứng: 𝑍𝛼𝛽 = 4𝜋𝜌𝛼(𝛽)∫𝑔𝛼𝛽(𝑟)𝑟 2𝑑𝑟 𝑅 0 (2.37) trong đó R là bán kính ngắt, thường được chọn là vị trí cực tiểu ngay sau đỉnh thứ nhất của HPBXT 𝑔𝛼𝛽(𝑟), 𝜌𝛼(𝛽)là mật độ số hạt α hoặc β. Giá trị của Zαβ cho ta biết trong hình cầu có tâm ở vị trí của một nguyên tử loại α và bán kính hình cầu là R có bao nhiêu nguyên tử loại β và ngược lại. Hình 2.3. Sơ đồ minh họa 2 chiều (2D) vị trí các đỉnh của HPBXT đối với cấu trúc VĐH 2.4.2. Hàm phân bố xuyên tâm của hạt nano HPBXT đối với mẫu khối được định nghĩa bởi 𝑔(𝑟) = 4𝜋𝑟2𝑑𝑟𝜌0 (2.38) Các lớp nguyên tử lân cận nhất H àm p h ân b ố x u y ên t âm c ặp 37 Ở đây n(r) là số nguyên tử trong lớp cầu có độ dày dr ở khoảng cách r tới nguyên tử xét; dấu có nghĩa là trung bình trên toàn bộ các nguyên tử khác trong hệ; 𝜌0là mật độ số hạt trong mẫu. Trong trường hợp hạt nano thì lớp cầu có hai phần: Phần ở trong hạt nano chứa một số lượng các nguyên tử và phần ở ngoài hạt nano ở đó không có nguyên tử nào (xem hình 2.4A). Hình 2.4. Minh họa sự xác định HPBXT đối với hạt nano (A); Lõi và bề mặt hạt nano (B); Ba vùng trong hạt nano (C). Nguyên tử Phần bên trong của lớp cầu Phần bên ngoài của lớp cầu Hạt nano A) B) Bề mặt Lõi C) Vùng III Vùng II Vùng I 38 Đối với hạt nano cầu chúng tôi xác định hàm gnano(r) như sau: 𝑔𝑛𝑎𝑛𝑜(𝑟) = 〈𝑛𝑛𝑎𝑛𝑜(𝑟)〉 〈𝑉𝑖𝑛〉𝜌𝑛𝑎𝑛𝑜 (2.39) Ở đây nano là mật độ số hạt nano; nnano(r) là mật độ nguyên tử ở phần bên trong của lớp cầu; Vin và Vout tương ứng là thể tích phần bên trong và phần bên ngoài của lớp cầu; Vin+Vout = 4r2dr. Hạt nano chứa lõi với mật độ bằng mật độ mẫu khối và bề mặt với mật độ thấp hơn như có thể thấy trên hình 2.4B. Nếu tất cả các nguyên tử trên bề mặt được bỏ đi, thì chúng tôi thu được lõi hạt nano. Cấu trúc của các mẫu xây dựng được phân tích thông qua HPBXT của lõi hạt nano. Để nghiên cứu sự phát triển của đám tinh thể chúng tôi chia mẫu thành ba vùng như được minh họa trên hình 2.4C. Ở đây vùng I có bán kính 17 Å. Hai vùng còn lại là các lớp đới cầu có độ dày là 4,5 Å và 6,0 Å. Số nguyên tử trong các vùng này là như nhau. Nhưng mật độ số hạt trong các vùng I, II và III lần lượt là 0,0825, 0,0825 và 0,0304. Nhớ rằng vùng III bao quanh bề mặt hạt nano và có cấu trúc xốp hơn. 2.5. Phương pháp xác định mầm tinh thể Pha vô định hình (VĐH) của mẫu được xây dựng của chúng tôi biến đổi thành pha tinh thể BCC. Do đó để xác định nguyên tử là nguyên tử tinh thể hoặc nguyên tử VĐH, một tiêu chuẩn sau đây được áp dụng. Đó là, nguyên tử được định nghĩa có cấu hình BCC nếu nó thỏa mãn hai điều kiện sau: (1) có 14 nguyên tử lân cận, (2) 6 trong 14 nguyên tử lân cận này có 4 lân cận chung với nguyên tử nhất định, 8 nguyên tử còn lại có 6 lân cận chung với nguyên tử nhất định. Bán kính ngắt được dùng để xác định nguyên tử lân cận là bằng 3.35 Å. Nguyên tử và 14 lân cận thỏa mãn tiêu chuẩn mô tả trên thuộc nguyên tử tinh thể. Tập hợp các nguyên tử như vậy tạo thành một mầm cơ bản. Đám tinh thể chứa một số các mầm cơ bản. Một đám gồm các đám tinh thể nhỏ sẽ có ít nhất một mầm có nguyên tử tinh thể chung với nó. Trong mô phỏng của chúng tôi, chúng tôi tìm thấy một vài mầm riêng biệt có ít nhất 30 nguyên tử hoặc các đám tinh thể lớn trong mẫu vật liệu. 39 Thuật toán xác định mầm: - Lần lượt xét từng nguyên tử thứ i có SPT bằng 14. Khoảng cách từ 14 nguyên tử lân cận của nguyên tử thứ i đến ng