Luận văn Nghiên cứu cấu trúc và cơ chế khuếch tán trong sio2 lỏng bằng phương pháp mô phỏng

chuyển của nguyên tử được mô tả bởi công thức: = a(t) + 6D.t (1.11) 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). Đối với chất lỏng đơn giản, hệ số khuếch tán vào cỡ 10-5 cm2s-1, chỉ cần vài trăm bước là đủ để nhận được kết quả tin cậy. Với mẫu chất lỏng nhớt 18 (SiO2, Al2O3), hệ số tự khuếch tán gần điểm nóng chảy (Tm) khoảng 10-6 cm2s-1, số bước phải tăng lên đến cỡ hàng trăm ngàn mới đạt được kết quả tốt. Do đó thời gian đạt tới đường tiệm cận sẽ tăng lên. Nếu hệ số khuếch tán vào cỡ 10-7 cm2s-1 hoặc nhỏ hơn thì về mặt thực hành không thể xác định hệ số này bằng phương pháp động lực học phân tử trong giới hạn tốc độ máy tính hiện nay. Tuy nhiên, phương pháp động lực học phân tử lại cho khả năng nghiên cứu ảnh hưởng của áp suất, nhiệt độ và khối lượng nguyên tử lên hệ số khuếch tán. Trong trường hợp đơn giản nhất, mô hình chất lỏng có dạng quả cầu cứng, độ linh động của nguyên tử được xác định bởi hệ số kích cỡ D(m/kT)1/2/ (trong đó m là khối lượng hạt,  là đường kính hạt, T là nhiệt độ tuyệt đối, k là hằng số Boltzmann). Ở nhiệt độ không đổi, hệ số này được tính toán bởi phương pháp động lực học phân tử, nó giảm dần với thể tích. Vì thế hệ số tự khuếch tán gần như không thể xác định được cho trường hợp V/Vo < 1.45, với Vo là thể tích tinh thể xếp chặt tương ứng ở nhiệt độ 0K. Nếu để mô tả sự thay đổi của hệ số khuếch tán theo nhiệt độ theo phương trình Arrhenius thì kết quả thu được bị lệch so với thực tế. Lí do là trên thực tế cơ chế khuếch tán trong chất lỏng không phụ thuộc kích hoạt tự nhiên, mà phụ thuộc vào nhiệt độ theo công thức 2D aT . Khi mô phỏng hệ phi tinh thể với thế tương tác gần, sự gia tăng áp suất làm giảm hệ số tự khuếch tán. Trong trường hợp sử dụng thế cặp Gauss: 2 0u(r) U .exp( br )  (1.12) Stilling và Weber đã phát hiện sự gia tăng bất bình thường của hệ số tự khuếch tán phụ thuộc áp suất của hệ. Hiệu ứng này xuất hiện vì khi khoảng cách giảm thì hàm lực –du(r)/dr vượt qua giá trị cực đại và sau đó giảm theo. Đối với chất lỏng làm lạnh nhanh thì giá trị thu được thường lớn hơn thực tế. Ví dụ trong mẫu của hệ thống sử dụng thế Lennad - Jones, giá trị của hệ số 19 khuếch tán tại nhiệt độ 0,1Tm là 10-7 cm2s-1, nghĩa là quá lớn so với thực tế. Với mẫu KCl và SiO2 cũng có kết quả tương tự. Sự sai lệch này có thể do kích thước quá nhỏ của mẫu hoặc do hiệu ứng lực tương tác xa, thực tế tương tác này là rất nhỏ rất khó xác định nên thường bỏ qua. Trong tương lai, với việc ứng dụng những thuật toán mới và máy tính hiện đại chúng ta hy vọng có thể nghiên cứu hệ với hàng triệu hạt. Các mẫu vật liệu khác nhau và các phương pháp phân tích khác nhau đã đạt được thành công về việc mô tả cơ chế khuếch tán tập thể trong chất lỏng siêu lạnh. Người ta dự đoán rằng các nguyên tử có thể chia làm hai loại: một loại như chất rắn không thể dịch chuyển, loại còn lại như chất lỏng có thể tham gia vào quá trình di chuyển khuếch tán. Thể tích hiệu dụng của ô lưới mà trong đó hạt đang lưu trú và tỷ lệ giữa thể tích này và thể tích tới hạn cho phép xác định nguyên tử thuộc loại nào. Ví dụ trong hệ với tương tác cặp hoàn toàn là đẩy thuộc dạng: 120u(r) (r / r)  , trong đó r là khoảng cách nguyên tử,  và ro là các hằng số, vai trò quyết định thuộc về hệ số kích cỡ: 1/4 3 * 0 B Nr . V k T              (1.13) Với N là số hạt trong thể tích V. Khi * của vật liệu phi tinh thể tăng từ 1,1 đến 1,5 thì thành phần hạt như trong chất lỏng giảm từ 1 đến 0, nghĩa là ảnh hưởng lên khuếch tán. Tiêu chuẩn này không đơn thuần là thể hiện tính cấu trúc mà còn thể hiện ý tưởng về mối tương quan giữa độ linh động của hạt và thể tích tự do. Trong vùng ổn định của trạng thái lỏng, phương pháp động lực học phân tử khiến ta dễ dàng và chắc chắn trong việc tính hệ số khuếch tán nếu thế tương tác giữa các hạt được biết tới. Mặc dù vậy, trong nhiều năm gần đây, dường như không có nhiều sự phát triển thực sự trong việc phân tích cơ chế khuếch tán tập thể cho chất lỏng ở vùng lạnh sâu. Quá trình khuếch tán của Si và O thông qua cơ chế chuyển động tập thể của cả khối SiOx (x=4, 5, 6). 20 Chương 2 PHƯƠNG PHÁP TÍNH TOÁN Trong chương này, kỹ thuật mô phỏng hệ ôxít SiO2 ở trạng lỏng sẽ được trình bày chi tiết. Trước tiên, chúng tôi trình bày phương pháp mô phỏng động lực học phân tử, thế Born-Mayer-Huggins, gần đúng Ewald-Hansen. Tiếp theo chúng tôi trình bày các kỹ thuật phân tích cấu trúc như : Hàm phân bố xuyên tâm, phân bố số phối trí, phân bố góc liên kết. Cuối cùng, chúng tôi trình bày cơ chế khuếch tán thông qua độ dịch chuyển bình phương trung bình và các loại phản ứng trong đám thông qua phản ứng SiOxSiOx' với x, x'=4, 5 và OSiy OSiy' với y, y'=2. 2.1. Phương pháp động lực học phân tử Phương pháp ĐLHPT là một công cụ cho phép chúng ta xây dựng mô hình vật liệu 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. Tọa độ ban đầu của các nguyên tử có thể lấy ngẫu nhiên nhưng phải thoả 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: Phương pháp ĐLHPT dựa trên phương trình chuyển động Newton: F m a i i i  (2.1) 2 i i i i i 1 N2 d r m a =m =F(r ...,r ) dt (2.2) 21 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 theo công thức: iji 1 ij U F r N j      (2.3) 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. 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 hai 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ằng biểu thức: 2 i i i i i F (t) r (t + dt) = 2r (t) - r (t - dt) + (dt) m (2.4) 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 ii r (t+dt)-r (t-dt) v (t) = 2dt (2.5) 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.6) trong đó ijx j F được xác định như sau: 0 ( ) . . ij ij i j x j ij ij U r x x F x r r           (2.7) 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.7). 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, 22 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: ij ij i>j U= U (r ) (2.8) 2N N 2i i i i i i=1 i=1 m v m r (t+dt)-r (t-dt) K= ( ) 2 2 2dt   (2.9) Năng lượng E của hệ có thể tính theo công thức: E=K+U (2.10) 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.11) 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: 2N 2i i i i i=1 1 m v r (t+dt)-r (t-dt) K= ( ) 2 2 2dt N i i m    (2.12) 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 (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: 23 ' 0 i i T v = v T (2.13) 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.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ẽ λ > 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 vậy, toạ độ mới của các nguyên tử được xác định: 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.15) 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: 3 . ( ) ( ) 2 NkT t K t (2.16) N 2 i i=1B B 2 K(t) 1 T(t)= = m (v (t)) 3 k N 3Nk i (2.17) 1 ( ) ( ) ( ) ( ) 3 ij ij i j N P t kT t r t F t V V     (2.18) Các mô hình mô phỏng NVE, NPT, NVT được sử dụng trong mô phỏng . 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 24 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 đến, 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. Thế tương tác Cho đến nay, đã có một một số hàm thế tương tác được phát triển để mô phỏng ô xít SiO2 ở cả trạng thái lỏng, tinh thể và VĐH [2, tr.5]. Tuy nhiên, yêu cầu cần thiết đối với các hàm thế đó là phải bảo đảm sự ổn định của các tứ diện SiO4 và phải thể hiện được bản chất của liên kết Si-O. Một phần của liên kết này có bản chất cộng hoá trị và một phần có bản chất ion. Thêm nữa, hàm thế cũng phải có giá trị thích hợp khi mô phỏng các dạng thù hình khác nhau của ô xít SiO2. Một số hàm thế thoả mãn điều kiện trên là thế tương tác cặp, một số hàm thế khác kết hợp cả thế tương tác cặp và thế ba thành phần. Nhiều hàm thế được rút ra từ dạng thế Born-Mayer-Huggins như sau: 2 6 8( ) exp( ) i j ij ij ij ij ij ij ij ij ij q q e U r A B r C r D r r       (2.19) Trong phương trình (2.19), ( )ijU r là thế tương tác giữa nguyên tử thứ i và nguyên tử thứ j; ijr là khoảng cách giữa chúng; qi, qj là điện tích của các nguyên tử i, j; Aij, Bij, Cij và Dij là các thông số tương ứng của hàm thế. Hai thành phần sau cùng ở vế phải của phương trình (2.19) liên quan đến tính đa cực của các nguyên tử và có thể bỏ qua nếu như các nguyên tử không có sự phân cực. Để nghiên cứu cấu trúc và cơ chế khuếch tán của ô xít SiO2, chúng tôi đã sử dụng thế tương tác cặp BKS. Tuy loại thế này khá đơn giản nhưng thực tế cho thấy việc sử dụng nó đã cho phép mô phỏng được nhiều tính chất của SiO2 ở trạng lỏng (như cấu trúc, mật độ và sự giãn nở nhiệt, v.v... ) Các kết quả mô phỏng phù hợp tốt với thực nghiệm. Thế tương tác này có dạng: 25 2 6( ) exp( ) i j ij ij ij ij ij ij ij q q e U r A B r C r r     (2.20) trong đó, rij là khoảng cách giữa hai nguyên tử thứ i và j. Các hệ số của hàm thế được cho trong bảng sau: Bảng 2.1. Các hệ số của thế BKS đối với hệ SiO2 Cặp nguyên tử Aij (eV) Bij (Å-1 ) Cij (eV Å6) Điện tích (e) O-O 1388,773 2,760 175,000 qO = − 1,2 Si-O 18003,757 4,873 33,538 qSi = + 2,4 Si-Si 0,0 0,0 0,0 2.3. Gần đúng Ewald-Hansen Các hàm thế tương tác sử dụng trong mô phỏng thường có hai phần: tương tác gần và tương tác xa. Thế tương tác U(r) ~ 1 sr được xem là tương tác gần khi s >3. Khi s  3 thế V(r) được xem như là thế tương tác xa. Xét một hệ gồm các nguyên tử được chứa trong một hình cầu bán kính R, năng lượng tương tác giữa các nguyên tử tỷ lệ với 1 sr ( r là khoảng cách giữa hai nguyên tử). Thế năng tác dụng nên nguyên tử thứ i là: ij( ) N s i j E i r   (2.21) trong đó N là số nguyên tử trong hình cầu. Giả sử hệ là đẳng hướng với mật độ , phương trình (2.21) có thể viết lại như sau: 2 0 ( ) 4 R s E i r dr r    (2.22) Trong luận văn này, chúng tôi đã sử dụng kỹ thuật gần đúng Ewald đối với tương tác xa (tương tác coulomb) và kỹ thuật ngắt tương tác đối với tương tác gần. 26 Xét một hệ gồm N nguyên tử đặt trong không gian tính toán là một hình lập phương có kích thước L. Thế tương tác của các nguyên tử trong hệ và các ảnh của nó với điều kiện biên tuần hoàn được xác định như sau: 0 1 1 ij, 1 U 2 N N i j n i j n q q r       (2.23) Trong đó qi , qj lần lượt là điện tích của nguyên tử thứ i và j; n là vectơ toạ độ tâm của các hình hộp ảnh và n = (n1, n2, n3) = n1.Lx + n2.Ly +n3.Lz với x, y, z là vectơ đơn vị trong hệ toạ độ Đề các. Tâm hình hộp chứa hệ N nguyên tử mô phỏng có tọa độ n = (0, 0, 0) còn các hình hộp ảnh có tâm ở toạ độ L.n theo ba chiều với n tiến đến vô cùng. Hình 2.1 mô tả mô hình tính toán gần đúng Ewald –Hansen trong không gian 2 chiều, mạng tuần hoàn 3x3 được dựng lên từ ô cơ sở có tâm n(0,0) n(-1,1) n(0,1) n(1,1) n(-1,0) n(0,0) n(1,0) n(-1,-1) n(0,-1) n(1,-1) Hình 2.1. Mô hình tính toán gần đúng Ewald –Hansen trong không gian 2 chiều, mạng tuần hoàn 3x3 được dựng lên từ ô cơ sở có tâm n(0,0). Trong biểu thức, tổng đầu tiên chỉ xét điều kiện i  j với n = 0; rij,n là khoảng cách giữa một nguyên tử với các nguyên tử khác trong hệ hoặc trong các hình hộp ảnh. ij,n jn i i jr r r r r nL     (2.24) Trong mô phỏng ĐLHPT, việc xác định thế tương tác Coulomb chiếm phần lớn thời gian mô phỏng. Tuy nhiên, sử dụng gần đúng Ewald-Hansen cho phép giảm đáng kể thời gian tính toán. 2.4. Xác định các đặc trưng vi cấu trúc và tính chất của mô hình Nội dung của phần này tập trung vào việc trình bày lần lượt cách xác định các đại lượng đặc trưng cho vi cấu trúc của vật liệu như hàm phân bố xuyên tâm, phân bố số phối trí, phân bố góc. Ngoài ra, cách xác định hệ số khuếch tán, một đại lượng quan trọng liên quan mật thiết đến vi cấu trúc SiO2 lỏng cũng sẽ được trình bày chi tiết. 27 2.4.1. Hàm phân bố xuyên tâm Trong mô phỏng vật liệu ở trạng thái lỏng 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á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 sẽ đượ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.25) 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.25) có thể viết lại một cách tường minh hơn như sau: 1 2 ij2 , ( ) ... ( ) ( ) N N N N i j j V g r drdr dr P r r r N      (2.26) ở đây rij = ri - rj và ri, rj là véc tơ toạ độ của các hạt thứ i và thứ j. Véc tơ 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 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éc tơ r. Đối với hệ đẳng hướng, g(r) chỉ phụ thuộc vào độ dài của véc tơ r. Lấy tích phân qua 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 ( , ) ( ) 4. ( ) V r r dr g r r rg r    (2.27) Thay phương trình (2.26) vào (2.27) chúng ta có: 1 2 ij2 2 ,( , ) 2 ij2 2 , ( , ) ( ) ... ( ) ( ) 4 ... ( ) ( ) 4 N N N N i j iV r r N N N N i j i V r r V g r dr dr dr dr p r r r r rN V dr dr p r dr r r r rN                  (2.28) 28 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à ij ( , ) ( , ) ( )i i j V r r n r r dr r r       (2.29) Thay (2.29) vào (2.28) chúng ta tìm được: 1 22 2 2 2 ( ) ... ( ) ( , ) 4 ( , ) 4 N N N i i i i V g r dr dr dr P r n r r r rN V n r r r rN           (2.30) Phương trình (2.30) có thể viết lại một cách đơn giản như sau: 0 ( ) ( ) r g r    (2.31) với 0 chính 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 2 1 ( , ) ( ) 4 i i N V n r r N r r r             (2.32) HPBXT cũng có thể được xác định từ thực nghiệm. Đại lượng có thể đo được trực tiếp từ thực nghiệm nhiễu xạ là cường độ nhiễu xạ I(). Trong đó,  là góc giữa tia tới và tia tán xạ. Gọi kin và kout tương ứng là véc tơ sóng tới và véc tơ sóng tán xạ. Bởi vì tán xạ là đàn hồi, |kin|=|kout|, với k=kin- kout chúng ta có: 4 K= k = sin( / 2) in    (2.33) Cường độ tán xạ có thể được tách thành hai phần: thừa số dạng nguyên tử f(K) và thừa số cấu trúc S(K) như sau: I(θ)=f(K)NS(K) (2.34) 29 Thừa số hình dạng đặc trưng cho loại nguyên tử và phụ thuộc vào việc hiệu chỉnh thiết bị đo. Thừa số cấu trúc được xác định bởi phương trình (2.35) và chứa tất cả các thông tin về vị trí của các nguyên tử. N l,m 1 S(K)= exp[iK(r -r )] N l m (2.35) Liên hệ giữa thừa số cấu trúc với hàm HPBXT, chúng ta dùng định nghĩa chuẩn về hàm phân bố xuyên tâm (2.25) và biểu diễn chuyển đổi Fourier của hàm Dirac delta như trong phương trình sau: iKx 1 ( ) 2 x dK e      (2.36) Tách tổng trong phương trình (2.35) thành 2 phần ứng với l=m và lm N N l l m N l m 1 1 S(K)= exp[iK.0] exp[iK.(r -r )] N 1 1 exp[iK.(r )] N l m lm N         (2.37) Áp dụng chuyển đổi fourier trong không gian ba chiều đối với (2.35) chúng ta được: iK.r 3 1 1 dK e ( ) ( ) ( ) (2π) N lm l m S K r r r N       (2.38) Thay phương trình (2.25) vào (2.38) chúng ta được: iK.r 3 1 N dK e S(K)=δ(r)+ g(r) (2π) V (2.39) Giản ước hàm delta trong phương trình (2.39) chúng ta được: iK.r 3 1 N dK e [S(K)-1]= g(r) (2π) V (2.40) Từ phương trình (2.40) chúng ta thấy HPBXT, g(r) có thể được xác định từ thực nghiệm thông qua thừa số cấu trúc. Trong luận văn này các HPBXT thành phần của hệ SiO2 lỏng được tính trong chương trình mô phỏng như sau: 30   βα NN , αβ ij i j 0 αβ 0 α β ( ) N δ(r -r) α, β Si (Ge), O α=β ( 1) N N α β ρ N N g r N N N                (2.41) Ở đâ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à loại ; 0 là mật độ nguyên tử trung bình trong thể tích V. Hàm phân bố xuyên tâm tổng cộng được xác định như sau: , ij 2 ( ) ( ) i i j j i i i c b c b g r g r c b            (2.42) Trong đó ci, cj là nồng độ của nguyên tử loại α và β, bi và bj là hệ số tán xạ cho nguyên tử loại α và β. 2.4.2. Xác định số phối trí và độ dài liên kết Số phối trí trung bình Z được xác định bằng biểu thức tích phân píc thứ nhất HPBXT tương ứng: 2 0 4 ( ) cr jZ g r r dr   (2.43) trong đó, rc là bán kính ngắt, thường được chọn là vị trí cực tiểu ngay sau píc thứ nhất của HPBXT g(r). 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 là rc, có bao nhiêu nguyên tử loại . Từ vị trí của đỉnh thứ nhất trong các hàm phân bố xuyên tâm thành phần cho phép ta xác định được độ dài liên kết giữa các cặp nguyên tử. Cụ thể, từ vị trí đỉnh thứ nhất của HPBXT thành phần gSi-Si(r) ta suy ra được khoảng cách lân cận gần nhất giữa hai nguyên tử Si-Si. Tương tự, ta có thể tính được độ dài liên kết giữa các cặp nguyên tử O-Si và O-O từ vị trí đỉnh thứ nhất của HPBXT thành phần gO-Si(r) và gO-O(r). 31 2.4.3. Xác định phân bố góc Như chúng ta đã biết, cấu trúc của ôxít SiO2 ở trạng thái lỏng được tạo thành từ các đơn vị cấu trúc cơ bản SiO4, SiO5, SiO6. Phân bố góc Si-O-Si sẽ cho biết liên kết của các nguyên tử trong một đơn vị cấu trúc cơ bản SiOx (thông tin về trật tự gần). Còn phân bố góc O-Si-O cho biết liên kết giữa các đơn vị cấu trúc cơ bản. Như vậy, từ kết quả tính phân bố góc có thể xác định được sự thay đổi trật tự gần trong các đơn vị cấu trúc cũng như sự thay đổi trật tự trong khoảng trung, liên quan đến sự thay đổi liên kết giữa các đơn vị cấu trúc. Trong các công trình [22,9] các kết quả tính toán phân bố góc chưa cung cấp được thông tin chi tiết về trật tự gần của các đơn vị cấu trúc cơ bản. Phân bố góc O-Si-O trong các công trình trên nhận được bằng cách tính trung bình cho tất cả các đơn vị cấu trúc trong hệ. Do vậy, trong luận văn này, việc tính toán phân bố góc được thiết lập cho từng đơn vị cấu trúc SiOx một cách riêng biệt. Vì thế so với kết quả trong các công trình trên, kết quả tính toán trong luận văn này cung cấp nhiều thông tin chi tiết hơn về vi cấu trúc. Thuật toán sử dụng xác định phân bố góc trong luận văn này được thiết lập như sau: Trước tiên, chúng tôi tiến hành xác định tập hợp các đơn vị cấu trúc cơ sở SiOx. Sau đó, ứng với mỗi loại đơn vị cấu trúc, chúng tôi xác định phân bố góc O-Si-O. Tuy nhiên, đối với phân bố góc Si-O-Si, chúng tôi cũng tiến hành xác định tất cả các đơn vị cấu trúc OSiy (y=2, 3, 4) trước, sau đó xác định phân bố góc Si-O-Si. Để xác định góc O-Si-O hoặc Si-O-Si, khi biết toạ độ của các nguyên tử tương ứng, chúng tôi làm như sau: giả sử chúng ta xét một tập gồm ba nguyên tử với toạ độ tương ứng: O1(x1,y1,z1); Si(x2,y2,z2); O2(x3,y3,z3). Góc O-Si-O được xác định bằng biểu thức: 32            2 2 2 2 2 2 2 1 2 1 2 1 212121 . ... arccos nmlnml nnmmll  (2.44) trong đó: 1 2 1 1 2 1 1 2 1 2 2 3 2 2 3 2 2 3 ; ; ; ; ; ; l x x m y y n z z l x x m y y n z z             (2.45) Ở đây, trật tự của các toạ độ x1, x2, x3 (tương tự với toạ độ y và z) đóng vai trò quan trọng trong việc xác định sự phân bố góc. Góc Si-O-Si được xác định hoàn toàn tương tự như đối với góc O-Si-O. 2.4.4. Xác định cơ chế khuếch tán Các mẫu vật liệu khác nhau và các phương pháp phân tích khác nhau đã đạt được thành công về việc mô tả cơ chế khuếch tán tập thể trong chất lỏng siêu lạnh. Quá trình khuếch tán trong vật liệu ô xít SiO2 ở trạng thái lỏng là chủ đề của những nghiên cứu chuyên sâu trong các thập kỷ qua. Kể từ năm 1965, Adam và Gibbs [3] đề xuất một lý thuyết để mô tả sự chuyển động của chất lỏng siêu lạnh. Lý thuyết này cho rằng các chất lỏng thay đổi cấu trúc của nó theo các vùng độc lập. Các vùng có chứa các nguyên tử lại có xu hướng hình thành các đám, từ đó phát triển các khái niệm về động lực học không gian không đồng nhất. Đối với các chất lỏng động không đồng nhất cũng được xây dựng bởi các thí nghiệm và mô phỏng. Kết quả thấy rằng hầu hết các nguyên tử linh động và không linh động đã hình thành các đi chuyển chậm và nhanh trong không gian. Chúng tôi đã nghiên cứu sự khuếch tán trong SiO2 lỏng sau một cách tiếp cận rằng sự di chuyển các hạt được đánh giá thông qua phản ứng SiOxSiOx' với x, x'=4, 5 và OSiy OSiy' với y, y'=2. 33 Hình 2.2. Mô hình hóa các loại phản ứng trong SiO2 lỏng A) là phản ứng SiO4SiO5, SiO5SiO4 và liên kết Si 12- O2 được thay thế bởi liên kết Si12- O6; B) là phản ứng SiO5SiO4, SiO4SiO5 và các liên kết Si 10- O5 bị hỏng, sau đó phục hồi.