Luận văn Nghiên cứu cấu trúc và động học của Al2O3 - 2SiO2 lỏng bằng phương pháp mô phỏng

nhất trong hóa học vô cơ là họ cấu trúc Zeolite. Cấu trúc này, thể hiện nhiều hiệu ứng đặc biệt khi nén. Giống như silica, sự giảm thể tích của zeolite khi nén là do sự giảm thể tích tự do của các ống kênh thông qua sự quay của các khối tứ diện và sự uốn của các góc liên kết Si-O-Si. Khi nghiên cứu sự thay đổi cấu trúc của zeolite natrolite tự nhiên khi nén, kết quả cho thấy thể tích của đơn vị cấu trúc cơ bản sẽ giảm liên tục khi nén và không có sự chuyển pha cấu trúc ở áp suất nén nhỏ hơn 7 GPa. Ở áp suất nén lớn hơn 7 GPa, mẫu vật liệu sẽ chuyển sang pha có cấu trúc VĐH, tương tự như quá trình chuyển pha cấu trúc quan sát được trên silica. Quá trình chuyển pha từ pha cấu trúc tinh thể sang pha cấu trúc VĐH cho thấy, có sự phá vỡ cấu trúc mạng tuần hoàn (gồm các tứ diện sắp xếp tuần hoàn) để chuyển sang mạng gồm các tứ diện sắp xếp ngẫu nhiên. Nếu các ống kênh (channels) và các lỗ trống (holes) trong framework đủ lớn chúng ta có thể nhận được hệ số 12 nén có giá trị âm. Hiện tượng hệ số nén có giá trị âm xảy ra khi kích thước phân tử của môi trường nén đủ nhỏ để có thể đi vào các ống kênh và lỗ hổng. Bởi vì áp lực bên trong của cấu trúc zeolite dường như lớn hơn áp lực môi trường nén, vì thế tổng thể tích của hệ (tinh thể có cấu trúc zeolite và môi trường nén) giảm. Tuy nhiên, thể tích của tinh thể có cấu trúc zeolite tăng khi có sự nén. Bằng cách tương tự, Hazen đã quan sát được thấy các hiện tượng dị thường khác thường quan sát thấy trong quá trình nén vật liệu có cấu trúc framework là sự tăng kích thước theo một chiều hoặc hai chiều khi nén thủy tinh (hệ số nén có giá trị âm trên một hoặc hai chiều). Hiện tượng này cũng được đề cập như hiện tượng mà tỷ số Poisson có giá trị âm. Một ví dụ khác về họ, các vật liệu có cấu trúc framework mà sự chuyển pha của nó được đặc trưng bởi sự uốn của các liên kết trong mạng là vật liệu perovskites. Vật liệu này có công thức chung là ABX3 có thể được xem như bắt nguồn từ cấu trúc loại WO3. WO3 (BX3) tạo thành mạng ba chiều từ các khối bát diện WO6, W(B) ở các vị trí đỉnh của các đơn vị cấu trúc hình lập phương, cation A ở vị trí tâm của hình lập phương. Trong các cấu trúc lý tưởng góc B-X-B ở đỉnh của các bát diện là 180°. Cấu trúc framework của các khối bát diện này dễ bị méo bởi sự quay của các bát diện cứng. Vì vậy, sẽ dẫn đến sự uốn của cấu trúc framework. Hiện tượng méo này có thể xảy ra ở điều kiện áp suất thường tùy thuộc vào bản chất của các cation A và B. Ở điều kiện nhiệt độ và áp suất thường, cấu trúc hình học của vật liệu phụ thuộc vào kích thước tương đối của các cation và liên kết giữa chúng. Tương tự như vậy, ảnh hưởng của quá trình nén tới cấu trúc của vật liệu phụ thuộc vào kích thước tương đối ở điều kiện thường và hệ số nén tương đối của các cation. Trong hầu hết các trường hợp thì cation A dễ bị nén hơn. Do đó, dẫn đến sự nghiêng của các khối đa diện. Một ví dụ về uốn của cấu trúc framework là sự méo của cấu trúc orthorhombic trong MgSiO3 khi áp suất tăng. Nhiều tinh thể có cấu trúc được hình thành bởi sự xếp chặt của các anion, còn các cation thì ở vị trí xen kẽ của các anion. Các anion có kích thước lớn và 13 mềm tiếp xúc với nhau, các cation có kích thước nhỏ và cứng nằm ở vị trí tách biệt nhau và chỉ tiếp xúc với các anion xung quanh. Số anion bao quanh một cation được xác định bởi tỷ số của kích thước của cation và anion. Trong cấu trúc xếp chặt của các anion, sự giảm thể tích khi tăng áp suất không phải do sự quay của các đa diện cứng gây ra mà sự giảm thể tích là do sự xếp chặt hơn của các anion, do đó nó làm giảm khoảng cách anion-anion trong cấu trúc. Sự giảm “kích thước” của các anion không ảnh hưởng đến kích thước của cation tương ứng , điều này dẫn đến tỷ số kích thước của cation/anion tăng lên, kết quả là số phối trí tăng. Sự tăng số phối trí thường kèm theo sự tăng khoảng cách liên kết cation-anion. Đây là một hiệu ứng gây ngạc nhiên đối với một chuyển pha áp suất cao. Tuy nhiên, sự mở rộng của hình cầu chứa các lân cận gần nhất được bù trừ bởi sự co lại của hình cầu chứa các lân cận thứ hai. Do đó, điều quan trọng là không chỉ nhìn vào các lân cận gần nhất của các cation. Khi nhìn vào cả hình cầu lân cận thứ nhất và thứ hai của cation, có thể thấy sự tăng số phối trí của các cation làm cho môi trường liên kết của các cation ít bị méo hơn. Cần nhớ rằng, sự tăng số phối trí có thể được giải thích hợp lý về mặt định lượng với sự méo môi trường liên kết của các cation. Sự nén một cấu trúc xếp chặt của các anion trước hết dẫn đến sự giảm khoảng cách anion-anion nhưng không sắp xếp lại. Điều này làm giảm tất cả các khoảng cách liên kết cation-anion làm cho các liên kết trở nên bền vững hơn (gọi là sự tăng liên kết). Trong khuôn khổ phương pháp hóa trị liên kết của Brown, sự tăng liên kết cation-anion xung quanh một nguyên tử nói chung gây ra sự “quá liên kết” của các nguyên tử, điều này dẫn đến sự không ổn định của cấu trúc, do mối quan hệ hàm mũ giữa cường độ liên kết và độ dài liên kết. Một tập hợp các độ dài liên kết cân bằng (tập hợp các liên kết riêng có cùng độ dài) sẽ có tổng hóa trị liên kết thấp hơn do với một tập các liên kết khác (có cùng số liên kết, cùng giá trị trung bình của độ dài liên kết) nhưng có các độ dài liên kết riêng khác nhau (sự sắp xếp bị 14 méo). Sự tăng số phối trí sẽ làm giảm sự quá liên kết (overbonding) của các nguyên tử bằng cách làm cho môi trường liên kết xung quanh nguyên tử ít bị méo hơn. Một ví dụ về hiệu ứng này là sự chuyển pha của silica (SiO2) từ coesite (nguyên tử Si có số phối trí là 4, nguyên tử O có số phối trí là 2) thành stishovite (nguyên tử Si có số phối trí là 6, nguyên tử O có số phối trí là 3). Trong khi cấu trúc của coesite là sự xếp chặt của các tứ diện SiO4 (các tứ diện SiO4 liên kết với nhau thông qua các đỉnh của tứ diện), cấu trúc của stishovite giống cấu trúc rutile được đặc trưng bởi các khối tám mặt SiO6 liên kết với nhau thông qua cạnh chung (edge sharing). Một số ví dụ khác về sự tăng số phối trí khi tăng áp suất là sự chuyển pha từ cấu trúc loại NaCl (số phối trí 6) sang cấu trúc loại CsCl (số phối trí 8) trong alkaly chroride hoặc chuyển pha áp suất cao từ cấu trúc loại zincblende (số phối trí 4) sang cấu trúc loại NaCl hoặc β-Sn (số phối trí 6) trong một số vật liệu bán dẫn. Trong các công trình nghiên cứu gần đây, có hai cơ chế được đưa ra để giải thích cho sự tăng số phối trí của các cation và đảm bảo sự cân bằng điện của hệ. Cơ chế thứ nhất liên quan đến sự hình thành cấu trúc “triclusters” (cấu trúc mà một nguyên tử O liên kết với 3 cation), cơ chế thứ hai là một nguyên tử trong liên kết cầu đi vào không gian hình cầu phối trí của khối AlO4 để hình thành Al có số phối trí cao. Cả hai cơ chế đều tạo ra [3]O, nhưng trong cơ chế thứ nhất thì [3]O liên kết với ba tứ diên (AlO4 hoặc SiO4), còn cơ chế thứ hai thì [3]O liên kết với một hoặc nhiều nguyên tử Al có số phối trí cao (AlO5 hoặc AlO6). Kết quả nghiên cứu cho thấy cấu trúc của nhôm silicate ở áp suất cao phụ thuộc vào nồng độ của Al và phụ thuộc vào tỷ số NBO/T (số nguyên tử O không tham gia liên kết cầu trên một tứ diện). Nó dường như chỉ ra rằng, nồng độ NBO trong nhôm silicat ở áp suất thường là một thông số quan trọng trong việc hình thành cấu trúc của nó ở áp suất cao. Tuy nhiên, cho đến nay vẫn còn rất ít các công trình nghiên cứu về vấn đề này. 15 1.4. Tính chất động học của các hệ ôxít 1.4.1. Cơ chế khuếch tán Khuếch tán nói chung, đặc biệt là khuếch tán trong vật liệu ôxít dạng lỏng luôn là đề tài được nghiên cứu bằng cả lý thuyết và thực nghiệm. Tuy nhiên, để hiểu rõ được cơ chế khuếch tán không phải là một vấn đề đơn giản. Nhiều nghiên cứu đã đưa ra nhiều cách giải thích cho cơ chế khuếch tán trong vật liệu ôxít như: Cơ chế vacancy (cơ chế nút khuyết), cơ chế khuếch tán tập thể, cơ chế khuếch tán theo các khe trống trong mạng tinh thể (khuếch tán xen kẽ) - Cơ chế khuếch tán xen kẽ Các nguyên tử khuếch tán (có kích thước rất nhỏ hơn các nguyên tử nền) chuyển động qua khe của các nguyên tử nền và định sứ ở lỗ hổng tạo ra bởi các nguyên tử nền (hình 1.1). Hình 1.1. Cơ chế khuếch tán xen kẽ - Cơ chế khuếch tán Vacancy (Nút khuyết) Các nguyên tử nhẩy từ vị trí cân bằng của nó sang vị trí nút khuyết lân cận (hình 1.2). Theo cơ chế khuếch tán vacancy, người ta cho rằng quá trình khuếch tán của nguyên tử là sự trao đổi vị trí giữa các nguyên tử nằm tại nút mạng với các vacancy bên cạnh và hệ số khuếch tán sẽ tỷ lệ với nồng độ cân bằng vacancy. Nguyên tử nền Nguyên tử khuếch tán 16 Nguyên tử dung môi Nguyên tử chất tan Hình 1.2. Cơ chế khuếch tán qua nút khuyết Năng lượng kích hoạt bằng tổng năng lượng tạo vacancy và năng lượng dịch chuyển vacancy. - Cơ chế khuếch tán tập thể Trong dung dịch thay thế, kích thước nguyên tử của chất tan tương tự kích thước nguyên tử dung môi. Cơ chế khuếch tán của các nguyên tử chất tan và các nguyên tử dung môi là sự đổi chỗ trực tiếp của các nguyên tử lân cận hoặc có thể đổi chỗ theo đường vòng (hình 1.3). Hình 1.3. Cơ chế khuếch tán tập thể Đối với hệ phi tinh thể, người ta sử dụng phương pháp ĐLHPT để tính hệ số tự khuếch tán. Trong quá trình hồi phục của phương pháp ĐLHPT, thời gian mỗi bước nhảy thực hiện độ dịch chuyển nhỏ có liên quan đến sự dịch chuyển của các nguyên tử lân cận. Điều này dẫn đến sự xuất hiện cơ chế Nút khuyết (vacancy) 17 khuếch tán tập thể. Bình phương trung bình của độ dịch chuyển của nguyên tử được mô tả bởi công thức: = a(t) + 6D.t (1.1) Sự có mặt của số hạng thứ nhất là do chuyển động dao động của nguyên tử, khi thời gian lớn thì thành phần này không thay đổi và sự phụ thuộc vào thời gian của sẽ tiệm cận với đường thẳng. Dựa vào độ dốc của đường thẳng này ta tính được hệ số khuếch tán (D). 1.4.1. Động học không đồng nhất Hiểu được cấu trúc và động học của silica lỏng là rất quan trọng từ những quan điểm cơ bản và thực tiến. Theo các thực nghiệm nghiên cứu, silica lỏng có tính chất riêng biệt. Đối với trường hợp giàu silica lỏng không phụ thuộc vào áp suất của độ nhớt. Tuy nhiên, đối với trường hợp nghèo silica lỏng lại phụ thuộc vào áp suất. Sự gia tăng của số phối trí trong nguyên tử hoặc giảm trong góc liên kết Si–O–Si được cho là không phụ thuộc vào áp suất, hơn nữa silica lỏng còn thể hiện động học không đồng nhất. Điều này có nghĩa là chất lỏng bao gồm các vùng linh động và không linh động riêng biệt, ở đó nguyên tử chuyển động chậm hoặc nhanh. Sự không đồng nhất của cấu trúc liên quan đến động học không đồng nhất. Để làm rõ cơ chế khuếch tán, mô phỏng ĐLHPT là một công cụ thích hợp. Có rất nhiều mô phỏng động học không đồng nhất (dynamic heterogeneity) của silica lỏng không phụ thuộc vào áp suất. Các kỹ thuật số như: Hàm tương quan, trực quan hóa, phân tích đám được sử dụng rộng rãi để tiến hành động học không đồng nhất (ĐHKĐN). Gần đây, mô phỏng silicate chủ yếu tập trung vào việc sắp xếp các đơn vị TOx, ở đây T là mạng lưới ban đầu hoặc mạng mới. Các đơn vị TOx được kết nối với nhau thông qua cầu O. Mức độ polymer hóa được tính bằng tỷ lệ các nguyên tử NBO (cầu không có nguyên tử O tham gia) trên mỗi tứ diện được sử dụng để thiết lập mối quan hệ giữa động học và cấu trúc. 18 SiO2 và hợp chất của chúng với các ôxít khác nhau (Al2O3, MgO, PbO, CaO, Na2O, ) là thành phần chính của lớp vỏ trái đất. Vì vậy, AS2 được sự quan tâm cả khoa học vật liệu và khoa học trái đất. AS2 là một silic hai nguyên đơn giản, trong đó Al hoạt động như mạng mới và mạng lưới ban đầu. Cụ thể có một mạng SiO4 chia góc và các cation linh động. Theo các nghiên cứu về phổ Raman và nhiễu xạ tia X, nồng độ của silicate nóng chảy được đi kèm theo áp suất gây polymer hóa. NBO và tri-cluster đóng một vai trò quan trọng đối với tính chất vận chuyển khối lượng. Mô phỏng ĐLHPT cho thấy tỷ lệ tương đối của TOx phụ thuộc mạnh vào chất hóa học của silicate nóng chảy. Đặc biệt số phối trí cao của Al có khuynh hướng tăng với khối lượng Al2O3. Các kênh giàu Al được tìm thấy cung cấp bằng chứng cho các trật tự hóa học trong chất tan. Mô phỏng trên Na2O-SiO2-Al2O3 và CaO-SiO2-Al2O3 được hỗ trợ trong việc phân tích thực nghiệm của khuếch tán dị thường. Hơn nữa, nó được thể hiện một sự tương tác thú vị giữa động học và môi trường địa phương. Cụ thể là sự tương quan rõ ràng giữa ĐLHPT và mật độ xác suất tìm thấy. Các vùng có khuynh hướng chuyển động nhanh thì nồng độ của Ca và Al cao, trong khi các vùng có khuynh hướng chuyển động chậm thì chứa silicon. Trong nghiên cứu này, chúng tôi sử dụng phân tích về sự khuếch tán và cấu trúc AS2 lỏng trong một phạm vi áp suất rộng. Chúng tôi đề xuất một mô hình đám SC mà không tính xác định kết nối và cho biết thêm chi tiết về môi trường địa phương của ion trong chất lỏng hình thành mạng, phân tích tập trung vào cơ chế khuếch tán và động học không đồng nhất cho AS2 lỏng. Để thu được thông tin đầy đủ hơn về cấu trúc cũng như động học trong hệ không đồng nhất, người ta nghiên cứu các hàm tương quan theo cả biến không gian và thời gian. Những hàm tương quan này có thể thu được từ thực nghiệm tán xạ nơtron. 19  Hàm tương quan hai điểm Xét hàm tương quan mật độ phụ thuộc cả theo không gian và thời gian:      1 1G r, r, r r r , t r ,0 .     (1.2) Biểu diễn theo các mật độ vi mô ta thu được:         N N 1 1 j i i 1 j 1 G r, r , t r r r t r r 0          . (1.3) Đối với các hệ không đồng nhất, hàm  1G r, r , t chỉ phụ thuộc vào khoảng cách. Lấy tích phân theo thể tích ta thu được:        N N j i i 1 j 1 1 G r, t r r t r 0 N       (1.4) Ở thời điểm ban đầu t = 0 hàm van – Hove  G r, t được giản ước và ta thu được:            N N j i i 1 j 1 1 G r,0 r r 0 r 0 N r g r           (1.5) Hàm van – Hove có thể tách làm hai thành phần: phần đồng nhất  SG r,0 và thành phần không đồng nhất  DG r,0 :      S DG r,0 G r,0 G r,0  . (1.6) Trong đó,        N S i i i 1 1 G r,0 r r 0 r t N      . (1.7)        N D j j i j 1 G r,0 r r 0 r t N      . (1.8) 20 Hàm tương quan van – Hove có ý nghĩa vật lý là mật độ xác suất tìm hạt thứ i trong lân cận của r ở thời điểm t biết hạt thứ j ở trong lân cận của vị trí ban đầu ở thời điểm t = 0.  SG r, t là mật độ xác suất tìm hạt thứ j khác hạt thứ i ở thời điểm t khi biết vị trí ban đầu của hạt thứ i ở thời điểm t = 0. Từ điều kiện chuẩn hóa  SG r, t ta thu được:  SdrG r, t 1 . (1.9) Tương tự, ta có:  DdrG r, t N 1  . (1.10) Khi xét giới hạn thời gian t  , hệ không nhớ được cấu hình ban đầu nữa và các hàm tương quan khong phụ thuộc vào khoảng cách r nữa:    S S r t 1 limG r, t limG r, t 0 V     . (1.11)    D D r t N 1 limG r, t limG r, t V       . (1.12) 21 Chương 2 PHƯƠNG PHÁP NGHIÊN CỨU Trong chương này, kỹ thuật mô phỏng các hệ ôxít sẽ được trình bày chi tiết. Trước hết chúng tôi trình bày về mô hình động lực học phân tử. Tiếp theo, chúng tôi tạo ra mẫu AS2 bằng cách gieo ngẫu nhiên 5500 nguyên tử trong đó có 1000 nguyên tử Al, 3500 nguyên tử O và 1000 nguyên tử Si vào khối hình lập phương, kích thước khối hình lập phương được tính toán từ mật độ của các vật liệu thực. Sau đó, chúng tôi tạo ra các mẫu vật liệu ở các nhiệt độ khác nhau và áp suất phòng. Cuối cùng chúng tôi xác định các đặc trưng và tính chất của mẫu. 2.1. Xây dựng mô hình động lực học phân tử Trước hết, chúng tôi xét một hệ gồm N nguyên tử gieo vào khối hình lập phương cạnh L. Tọa độ ban đầu của các 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 của Newton. Đối với hệ gồm N hạt, hệ chuyển động tuân theo định luật II Newton có dạng như sau: 2 i i i i i 1 N2 d r m a =m =F (r ...,r ) dt (2.1) Trong đó, Fi là lực tổng hợp tác dụng lên nguyên tử thứ i từ các nguyên tử còn lại; mi và ai lần lượt là khối lượng và gia tốc của nguyên tử thứ i. Lực Fi được xác định bởi công thức: N ij i j 1 ij U F r      (2.2) Trong đó, ijU là thế tương tác giữa nguyên tử thứ i và nguyên tử thứ j và rij là khoảng cách giữa chúng. 22 Trong mô phỏng ĐLHPT, 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 II Newton. Trong thuật toán này, toạ độ của nguyên tử i ở thời điểm (t + dt) được xác định thông qua tọa độ của nó ở hai thời điểm t và (t - dt) bởi biểu thức:          2 i i i i i F t r t + dt = 2r t - r t dt + dt m  (2.3) Vận tốc ở thời điểm t được xác định thông qua tọa độ ở thời điểm (t - dt) và (t + dt) theo biểu thức:      i i i r t+dt -r t dt v t = 2dt  (2.4) Lực Fi(t) đượ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:   i i i ij ij ij y zi x x y z j j j F t F F F F F F        (2.5) Trong đó ijx j F được xác định như sau: ij ij i j x 0 j ij ij U(r ) x x F X . . r r           (2.6) Với 0X là véctơ đơn vị của trục OX. Các thành phần ij ijy zF , F được xác định tương tự như phương trình (2.6). Khi nghiên cứu các mô hình vật liệu bằng phương pháp ĐLHPT, tuỳ theo mục đích cần nghiên cứu mà người ta thường 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ế hoá 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ó cá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: 23  ij ij i>j U= U r (2.7)     2 2N N i ii i i i=1 i=1 r t dt -r t dtm v m K= 2 2 2dt           (2.8) Năng lượng E của hệ có thể tính theo công thức: E=K+U (2.9) 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: B 2 T=K 3Nk (2.10) Trong đó kB 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 2N N i ii i i i=1 i 1 r t dt -r t dtm v m K= 2 2 2dt           (2.11) 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.10). Giả sử nhiệt độ được tính từ phương trình (2.10) là T, nhiệt độ mong muốn của hệ đạt được là T0, điều chỉnh vận tốc vi của tất cả các nguyên tử theo phương trình sau: ' 0 i i T v = v T (2.12) Chúng ta sẽ thu được:   N N 2 ' 20 i i i i 0 i=1 i=1B B T1 1 T'= m v = m v =T 3k N 3k N T   (2.13) 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 24 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 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 Pmới, nếu Pmới >Phệ thì λ = 1- dP, ngược lại λ = 1+ dP với giá trị dP được chọn là 10-4. Do đó, toạ độ mới của các nguyên tử được xác định như sau: a a a a a a b b b b b b x' [i] = x [i]. ; y' [i] = y [i]. ; z' [i] = z [i]. x' [i] = x [i]. ; y' [i] = y [i]. ; z' [i] = z [i].       (2.14) Khi đó, kích thước mô hình sẽ có giá trị L’ = Lλ Khi xây dựng mô hình ĐLHPT, các thông số nhiệt độ và áp suất ở thời điểm t được xác định như sau:        2N i i i=1B B K t2 1 T t = = m v t 3 k N 3Nk  (2.15)        ij ij i j N 1 P t kT t r t F t V 3V     (2.16) 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ố tự 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 sẽ kéo 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 mô hình NVE, do đó thời gian mô phỏng sẽ được giảm đi rất nhiều. 2.2. Xây dựng mẫu Al2O3-2SiO2 lỏng Hệ Al2O3-2SiO2 (AS2) bao gồm 5500 nguyên tử (1000 nguyên tử Al, 1000 nguyên tử Si và 3500 nguyên tử O) ở các nhiệt độ và áp suất khác nhau được 25 xây dựng bằng phương pháp động lực học phân tử (ĐLHPT) sử dụng thế tương tác cặp Born – Mayer - Huggins với hàm tương tác được xác định theo biểu thức:     2 i j 6 8 ij ij ij ij ij ij ij ij ij q q e U r A exp B r C r D r r       (2.17) Ban đầu gieo ngẫu nhiên các nguyên tử vào hình lập phương, để loại bỏ cấu hình nhớ ban đầu, mẫu được nung nóng đến nhiệt độ 6000 K. Quá trình cân bằng tại 6000 K thu được bằng cách chạy hồi phục với 50 ps. Sau đó, mẫu được làm lạnh xuống tới 2000 K. Mẫu được chạy thống kê hồi phục NPT (áp suất và nhiệt độ không đổi) ở áp suất phòng. Sau đó, mẫu được chạy NVE (thể tích và năng lượng không đổi) để mẫu đạt trạng thái cân bằng, nghĩa là năng lượng không đổi (trên 50 ps). Tiếp theo, từ mẫu 2000 K chúng tôi tạo ra 7 mẫu khác nhau ở các nhiệt độ 2000, 2500, 3000, 3500, 4000, 4500, 5000 K và chạy NVE như mẫu 2000 K. Tương tự như vậy các mẫu áp suất khác nhau được xây dựng bằng cách lấy mẫu AS2 ở nhiệt độ 3500 K nén đến các áp suất khác nhau từ 0, 5, 10, 15, 20 GPa sau đó được chạy hồi phục NVE để đạt trạng thái cân bằng. Thông qua phân tích các mẫu trên, chúng tôi nhận được các thông số vi cấu trúc của AS2. 2.3. Xác định các đặc trưng cấu trúc mô hình 2.3.1. Xác định hàm phân bố xuyên tâm Trong mô hình vật liệu ở trạng thái lỏng và VĐH, một đại lượng tuân theo quy tắc thống kê được sử dụng để xác định vi cấu trúc của vật liệu ở mức nguyên tử, đó là hàm phân bố xuyên tâm (HPBXT). HPBXT cũng có thể được xác định từ thực nghiệm thông qua thừa số cấu trúc. Thông qua HPBXT có phân bố về số phối trí trung bình, khoảng cách liên kết trung bình và góc liên kết trung bình cũng được xác định. Theo cơ học thống kê, hàm tương quan cặp hay HPBXT g(r) được xác định như sau:    ij2 i, j i V g r r r N     (2.18) 26 Trong đó V là thể tích của mẫu vật liệu và N chính là số nguyên tử chứa trong thể tích V đó. Phương trình (2.18) có thể viết lại một cách như sau:       N N 1 2 N N ij2 i, j j V g r dr dr ...dr P r r r N     (2.19) Ở đây, ri,j =ri – rj và ri , rj là véctơ tọa độ của các hạt thứ i và thứ j. Véctơ r 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 r là do chúng ta chọn). Hàm g(r) có thể hiểu là tỷ lệ thuận với xác suất tìm thấy nguyên tử cách nguyên tử trung tâm một véctơ r. Đối với hệ đẳng hướng, g(r) chỉ thuộc vào độ dài của véctơ r. Lấy tích phân của thể tích V(r, r) giữa r và r+dr và giả sử rằng lớp vỏ hình cầu là đủ mỏng chúng ta sẽ thu được:      2 V r, r dr g r 4 r g. r r     (2.20) Thay phương trình (2.20) vào (2.19) chúng ta có:             N N 1 2 N N ij2 2 i, j iV r, r N N 2 N N ij2 2 i, j i V(r, r) V g r dr dr dr ...dr p r r r 4 r rN V dr ...dr p r dr r r 4 r rN                  (2.21) Tích phân hàm delta, chúng ta sẽ tính được số hạt trong lớp hình cầu là ni(r,r) và       N i ij i, j i V r, r n r, r dr r r        (2.22) Thay (2.22) vào (2.21) chúng ta tìm được:         N 1 2 N N i2 2 i i2 2 i V g r dr dr ...dr P r n r, r 4 r rN V n r, r 4 r rN           (2.23) Phương trình (2.23) có thể viết lại một cách đơn giản như sau:     0 r g r    (2.24) 27 Với 0 là mật độ nguyên tử trung bình trong thể tích V của mẫu vật liệu và (r) là mật độ nguyên tử ở khoảng cách r tính từ nguyên tử trung tâm.     0 i i 2 N