Xác định các tính chất nhiệt động cân bằng pha lỏng-Lỏng, lỏng-hơi của các hệ hóa học bằng phương pháp phân tích nhiệt động thống kê, hóa lượng tử và mô phỏng toàn cục Monte-Carlo

1 BỘ GIÁO DỤC VÀ ĐÀO TẠO TRƯỜNG ĐẠI HỌC ĐÀ LẠT ----oOo---- BÁO CÁO ĐỀ TÀI NGHIÊN CỨU KHOA HỌC CẤP BỘ GIAI ĐOẠN 2007-2009 XAÙC ÑÒNH CAÙC TÍNH CHAÁT NHIEÄT ÑOÄNG CAÂN BAÈNG LOÛNG-HÔI, LOÛNG-LOÛNG CUÛA CAÙC HEÄ HOÙA HOÏC BAÈNG PHÖÔNG PHAÙP PHAÂN TÍCH NHIEÄT ÑOÄNG THOÁNG KEÂ, HOÙA LÖÔÏNG TÖÛ VAØ MOÂ PHOÛNG TOAØN CUÏC MONTE CARLO MÃ SỐ: B2007-14-10 CHỦ NHIỆM ĐỀ TÀI: TS. PHẠM VĂN TẤT ĐÀ LẠT-2009 2 DANH SAÙCH NHÖÕNG NGÖÔØI THAM GIA ThS. Trần Kim Cương Khoa Hóa học, Đại học Đà Lạt ThS. Trần Kim Oanh Khoa Hóa học, Đại học Đà Lạt ThS. Nguyễn Hải Hà Khoa Hóa học, Đại học Đà Lạt ThS. Phạm Hầu Thanh Việt Khoa Hóa học, Đại học Đà Lạt CN. Lê Thị Thanh Trân Khoa Hóa học, Đại học Đà Lạt 3 MUÏC LUÏC Tóm tắt kết quả 5 Đặt vấn đề 7 Tổng quan 9 I. Lý thuyết hóa lượng tử 9 I.1. Cơ sở lượng tử 9 I.2. Phương pháp tính toán lượng tử 9 II. Phương pháp nhiệt động thống kê 12 II.1. Các phương trình trạng thái 12 II.2. Các mô hình hệ số hoạt độ 15 II.3. Lý thuyết mô hình COSMO 16 II.4. Lý thuyết mô phỏng Monte Carlo 23 III. Các nhóm hợp chất 32 Phương pháp 34 I. Các điều kiện thực hiện 34 I.1. Dữ liệu cấu trúc 34 I.2. Phần mềm 34 I.3. Sơ đồ thực hiện chung 34 II. Phương trình trạng thái 35 III. Mô hình hệ số hoạt độ 36 IV. Mô hình COSMO 36 IV.1. Tối ưu hóa cấu hình 36 IV.2. Tính toán điện tích và năng lượng pha lỏng 36 IV.3. Tính toán điện tích bề mặt sigma 37 IV.4. Tính toán hệ số hoạt độ 37 IV.5. Tính toán sai số 37 V. Mô phỏng Monte Carlo 38 V.1. Xác định hàm thế tương tác 38 V.2. Xác định điều kiện thực hiện mô phỏng Monte Carlo GEMC 39 V.3. Tính toán tính chất nhiệt động cân bằng pha 40 4 Kết quả và thảo luận 41 A. Cân bằng lỏng - hơi 41 I. Mô hình hệ số hoạt độ 41 1. Hệ bậc 2 của các hệ ete vòng với 1-clorohexan 41 2. Hệ bậc 2 của etyl 1,1-dimetyletylete(1)+propan-1-ol(2) 45 3. Hệ bậc 2 axít acetic(1) + dietyl ete (2) 47 4. Hệ bậc 2 etanol (1)+nước(2); methanol(1)+benzene(2); toluene(1)+clorobenzen(2) 51 5. Hệ bậc 2 2,3-butandiol với n-butanol, n-butylacetat và etylacetat 53 II. Mô hình COSMO-SAC 55 II.1. Hệ bậc 2 gồm các hỗn hợp bậc 2 tạo thành từ các hợp chất trong nhóm ankan+ankan; ankan+ancol; ankan+ceton, ancol+nước trong khoảng nhiệt độ (183,15 – 623,15 K) và áp suất (0,1-19,0MPa) 55 II.2. Hệ bậc 2 gồm phenol(1)+styren(2); etylmercaptan(1)+n-butan(2); ter-butylmercaptan(1)+propan(2) và dimetylete(1)+propan(2) 58 III. Mô phỏng Monte Carlo 61 III.1. Sử dụng hàm thế kiểu Lennard-Jones 62 III.2. Sử dụng hàm thế tương tác ab initio 71 B. Cân bằng lỏng - lỏng 77 I. Độ tan của aspirin trong dung môi 77 II. Cân bằng lỏng – lỏng của hệ bậc ba 79 Kết luận 82 Tài liệu tham khảo 83 Các kết quả đạt được 86 Phụ lục 87 5 TOÙM TAÉT KEÁT QUAÛ NGHIEÂN CÖÙU ÑEÀ TAØI KHOA HOÏC VAØ COÂNG NGHEÄ CAÁP BOÄ Tên đề tài: Xác định các tính chất nhiệt động cân bằng pha lỏng-lỏng, lỏng-hơi của các hệ hóa học bằng phương pháp phân tích nhiệt động thống kê, hóa lượng tử và mô phỏng toàn cục Monte-Carlo. Mã số: B2007-14-10. Chủ nhiệm đề tài: Phạm Văn Tất Tel: 0974655560 E-mail: phvtat@yahoo.co.uk Cơ quan chủ trì đề tài: TRƯỜNG ĐẠI HỌC ĐÀ LẠT Cơ quan và cá nhân phối hợp thực hiện: Đại học Cologne, Cộng hòa Liên bang Đức Thời gian thực hiện: 24 tháng Mục tiêu - Xây dựng cơ sở dữ liệu của một nhóm khoảng 20-30 hợp chất hóa học bao gồm các đơn chất vô cơ, chất hữu cơ và các chất lỏng thông thường. - Xây dựng phương trình thế cho các hệ hóa học từ tính toán nhiệt động thống kê và hóa học lượng tử. Tính toán cân bằng pha từ phương trình thế mới và từ mô hình hệ số hoạt độ, mô hình COSMO và kỹ thuật mô phỏng Monte Carlo toàn cục. Nội dung chính Các mô hình hệ số hoạt độ NRTL, Wilson, UNIQUAC và UNIFAC, mô hình lượng tử COSMO-SAC và phương pháp mô phỏng Monte Carlo toàn cục (GEMC) được sử dụng để tính toán cân bằng lỏng-lỏng, lỏng-hơi của đơn chất, hỗn hợp bậc hai của các chất và độ tan của hợp chất hữu cơ trong các dung môi. Phương pháp COSMO được thực hiện từ tính toán lượng tử phiếm hàm ab initio. Mật độ điện tích bao quanh phân tử được tính từ phương pháp phiếm hàm DFT-BP với bộ hàm cơ sở DNP. Cân bằng pha của hệ đơn và hệ bậc hai cũng như độ tan của các chất được dự đoán từ hệ số hoạt độ của chất nhận được từ phương pháp tính toán lượng tử và bằng kỹ thuật mô phỏng Monte Carlo GEMC khi sử dụng các phương trình thế kiểu Lennard-Jones, hay phương trình thế mới nhận được từ các tính toán hóa học lượng tử ab initio ở mức độ lý thuyết chính xác cao. Các kết quả cân bằng pha và độ tan của hợp chất trong các dung môi nhận được từ các kỹ thuật tính toán lý thuyết phù hợp tốt với số liệu tính toán và thực nghiệm từ các tài liệu khác. 6 SUMMARY Project Title: Evaluation of thermodynamic properties of liquid-liquid, vapor- liquid equilibria for chemical systems by using the approaches of statistical thermodynamic analysis, quantum calculation and Monte Carlo Simulation. Code number: B2007-14-10. Coordinator: Pham Van Tat Tel: 0974655560; E-mail: phvtat@yahoo.co.uk Implementing Institution: University of Dalat Cooperating Institution: Physical Institute, University of Cologne, Germany Duration: 24 months (from 6/2007 to 6/2009) Objectives: - Constructing the database of a set of 20-30 compounds. It includes of pure inorganic and organic substances and normal liquids. - Establishing the pair potentials of chemical systems from statistical thermodynamic analysis and quantum calculation. The phase equilibria of them are calculated by using new pair potentials, activity coefficient models, COSMO-SAC model and Gibbs Ensemble Monte Carlo Simulation. Main contents: In this work, we used the activity coefficient models NRTL, Wilson, UNIQUAC, UNIFAC, and the conductor-like screening model COSMO-SAC (segment activity coefficient) resulting from ab initio calculations, and the Gibbs Ensemble Monte Carlo simulation (GEMC) to predict the vapor-liquid equilibria of pure components, the binary mixtures and the solubility of organic substances in the different solvents. The models NRTL, Wilson, UNIQUAC and UNIFAC were constructed with the theory of statistical thermodynamic analysis. The method COSMO was carried out by using the ab initio calculations of density functional theory (DFT). The segment profiles of organic molecules were calculated by the density functional theory VWN-BP with the basis set DNP. The activity coefficients of substances were obtained from these sigma profiles. The vapor-liquid equilibrium data of pure components and the binary systems were predicted using the activity coefficients derived from the quantum computations. The solubility data of compounds in the different solvents were also calculated using these similar ways. The Gibbs Ensemble Monte Carlo simulation technique was also used to generate the vapor-liquid equilibrium data for the pure components and the binary mixtures using the Lennard-Jones-style intermolecular potential and the site-site intermolecular potentials ab initio. These were resulted from the high level of theory CCSD(T) with the basis sets aug-cc-pVmZ (m = 2, 3). The data of vapor-liquid equilibria and solubility for several systems resulting from the theoretical calculations were compared with the experimental data and those from literatures. The discrepancies between them were insignificant. 7 ÑAËT VAÁN ÑEÀ Hướng nghiên cứu xác định các tính chất nhiệt động cân bằng lỏng-hơi và lỏng- lỏng của các hệ hóa học khác nhau bằng phương pháp nhiệt động thống kê, hóa lượng tử và mô phỏng toàn cục Monte Carlo hiện nay còn rất mới ở Việt Nam và có rất ít tài liệu nghiên cứu đi sâu vào lĩnh vực nghiên cứu này. Hiện nay chỉ có một vài kết quả nghiên cứu liên quan đến các tính chất cơ chế phản ứng của các hợp chất dựa trên các nguyên lý cơ bản của hóa học lượng tử. Ở một số nước trên thế giới việc nghiên cứu và dự đoán các tính chất nhiệt động cân bằng pha dựa trên lý thuyết nhiệt động thống kê, tính toán lượng tử hiện đại và các kỹ thuật mô phỏng Monte Carlo vẫn đang là vấn đề chính được nhiều nhà khoa học thuộc nhiều cơ sở của các Đại học trên thế giới như châu Âu, Nhật Bản, Hàn Quốc và Mỹ khai thác nghiên cứu. Xuất phát từ nhu cầu thực tế không thể thiếu đối với các lĩnh vực khoa học ứng dụng và khoa học cơ bản cần có những số liệu về tính chất nhiệt động của các hệ cân bằng lỏng-hơi và lỏng- lỏng của các hợp chất và hỗn hợp bậc 2 và bậc 3 của các chất. Các lĩnh vực có thể sử dụng các dữ liệu của đề tài này bao gồm lĩnh vực nghiên cứu độc tố môi trường, kỹ thuật xử lý môi trường, hóa nông nghiệp, hóa dầu, nhiên liệu lỏng, vật liệu tổng hợp, hóa dược và hóa thực phẩm. Đặc biệt lĩnh vực môi trường rất cần đến các dữ liệu thông tin về tính chất hóa lý của các hợp chất hữu cơ, như tính chất phân tán, cân bằng lỏng hơi, cân bằng lỏng- lỏng liên quan đến độ tan của các chất trong môi trường nước. Dựa vào tính chất hóa lý của các hợp chất hữu cơ, các nhà kỹ thuật xử lý môi trường có thể đánh giá, quản lý và đưa ra biện pháp xử lý các chất gây ô nhiễm theo một kỹ thuật xử lý như kỹ thuật chiết, kỹ thuật chưng cất. Trong các ngành dược phẩm, hóa thực phẩm và đặc biệt là kỹ thuật lọc dầu việc tách chất và trộn chất luôn là vấn đề cần thiết liên quan đến các kỹ thuật thiết kế thiết bị sản xuất thích hợp. Vì có những ứng dụng tế như vậy cho nên đề tài này nhằm nghiên cứu đưa ra các phương pháp tính toán mới trên cơ sở của lý thuyết nhiệt động thống kê, hóa 8 học lượng tử và kỹ thuật mô phỏng Monte Carlo cho phép giải quyết các vấn đề về thông tin hóa lý của các hệ chất lỏng, giải quyết được nhiều vấn đề khó khăn đo đạc trong thực nghiệm khi mà phương tiện, trang thiết bị phòng thí nghiệm phục vụ đo đạc hiện nay chưa thể đáp ứng được đầy đủ ở mọi cơ sở nghiên cứu khoa học và cơ sở sản xuất. Chính từ những nhu cầu thực tế như vậy mà đề tài này cố gắng nghiên cứu đưa ra các kỹ thuật tính toán tính chất nhiệt động của một số hợp chất dựa trên lý thuyết nhiệt động thống kê, hóa lượng tử và các phương tiện hiện có của trường Đại học Đà Lạt để giải quyết các vấn đề liên quan đến các hợp chất hữu cơ. Mục tiêu chính của đề tài cần thực hiện bao gồm: - Xây dựng cơ sở dữ liệu của một nhóm khoảng 20 hợp chất hóa học bao gồm các đơn chất vô cơ, chất hữu cơ, các chất lỏng ion thông thường. - Xây dựng phương trình thế cho các hệ hóa học từ tính toán nhiệt động thống kê và hóa học lượng tử. - Tính toán các tính chất nhiệt động cân bằng pha lỏng-hơi, lỏng-lỏng từ phương trình thế được xây dựng từ lý thuyết nhiệt động thống kê, mô hình hệ số hoạt độ, mô hình COSMO và kỹ thuật mô phỏng Monte Carlo toàn cục. - So sánh các kết quả nghiên cứu với các phương pháp tính toán khác và các nguồn dữ liệu thực nghiệm khác nhau. 9 TOÅNG QUAN I. LÝ THUYẾT HÓA LƯỢNG TỬ I.1. Cơ sở lượng tử [3,4] Nền tảng của hóa lượng tử là mô hình sóng về nguyên tử, coi nguyên tử được tạo thành từ một hạt nhân mang điện tích dương và các điện tử quay xung quanh. Các điện tử trong mô hình sóng là các đám mây điện tử chuyển động trên các quỹ đạo và vị trí của chúng được đặc trưng bởi một phân bố xác suất chứ không phải là một điểm rời rạc. Để biết được phân bố xác suất, và tính chất của hạt đơn trong không gian 3 chiều người ta giải phương trình Schrödinger: (1) Ở đây ψ là hàm sóng ở vị trí r (x,y,z) và thời điểm t; V là thế năng; m là khối lượng của hạt. Phương trình Schrödinger không phụ thuộc thời gian được sử dụng chủ yếu: (2) Ở đây E là năng lượng và V là thế năng. Các thông tin về hệ lượng tử nhận được từ việc giải phương trình Schrödinger. Đối với hệ nhiều điện tử phương trình Schrödinger được giải gần đúng. I.2. Phương pháp tính toán lượng tử [3] Để giải gần đúng phương trình Schrödinger, các phương pháp khác nhau được đưa ra dựa trên phương trình Roothaan, tất cả để giải quyết vấn đề thế năng tương tác giữa các điện tử với nhau trên cơ sở các phương trình chứa tích phân Coulomb và các tích phân xen phủ giữa các electron. Hiện nay có ba nhóm phương pháp tính gần đúng khác nhau: - Phương pháp ab initio (phương pháp không kinh nghiệm): CI, CISD, QCISD(T), MPn (n = 2-4) - Phương pháp lý thuyết hàm mật độ (DFT) 10 - Phương pháp bán kinh nghiệm sử dụng các tham số thực nghiệm: CNDO, NDDO, MINDO, AM1, PM3, ZINDO … 1. Phương pháp ab initio [3] Hiện nay phương pháp này được coi là phương pháp chính xác nhất, phương pháp này sử dụng các tham số lý thuyết ban đầu của hạt và thực hiện tính tất cả các giai đoạn ước tính của tất cả các điện tử đối với phương trình Schrödinger. Các bộ hàm cơ sở được sử dụng cùng các phương pháp trong các tính toán ab initio có dạng Gauss: 1 N i icμ μ μ φ χ = =∑ (3) Ở đây cμi là hệ số giãn nở orbital phân tử. Các hàm cơ sở χ1,…, χN cũng được chọn là chuẩn hóa. Các phương pháp chính được dùng chủ yếu hiện nay: - Phương pháp tương quan điện tử (Electron Correlation Method): o Tương tác cấu hình (Configuration Interaction (CI)) o Tương tác cấu hình đầy đủ (Full CI) o Tương tác cấu hình giới hạn (Limited Configuration Interaction) - Phương pháp nhiễu loạn Möller-Plesset (Möller– Plesset Perturbation (MP)): o Phương pháp nhiễu loạn bậc 2: MP2 o Phương pháp nhiễu loạn bậc 3: MP3 o Phương pháp nhiễu loạn bậc 4: MP4 Trong phương pháp nhiễu loạn, hàm sóng bị nhiễu loạn biểu diễn dưới dạng một chuỗi lũy thừa: 0 2 1 2 2 3 3 ...ψ ψ λ ψ λ ψ λ ψ= + + + + (4) 2. Phương pháp hàm mật độ (Density Functional Theory (DFT)) [12] Phương pháp này dựa trên phương pháp cơ học lượng tử mô hình của Thomas- Fermi-Dirac và dựa trên sự tương quan điện tử thông qua hàm mật độ điện tử. Các phương pháp tương quan được sử dụng: - VWN, Perdew-Zunger 81, LYP, Perdew-Wang 91, PBE96, HCTH98 - Dùng với tổ hợp: B3-LYP, B3-PW91, EDF1, Becke97 11 3. Phương pháp bán thực nghiệm [3] Để giảm khó khăn trong tính toán, các phương pháp bán thực nghiệm khác nhau được đưa ra: - Thay thế tích phân đa tâm bằng tích phân ít tâm. - Thay các tích phân và các đại lượng vật lý trong phương trình lượng tử bằng các tham số thực nghiệm. - Các phương pháp bán thực nghiệm thường dùng : o CNDO (Complete Neglect of Differential Overlap) Phương pháp này dựa trên phương pháp ZDO bỏ qua hoàn toàn sự xen phủ vi phân. Đây là phương pháp kinh nghiệm đơn giản nhất, tính theo giải thuật trường tự hợp, được dùng để tính toán các tính chất điện tử ở trạng thái cơ bản. o INDO (Intermediate Neglect of Differential Overlap) Phương pháp bỏ qua sự xen phủ trung gian vi phân, được sử dụng để tính toán cấu tạo lớp electron và mật độ spin của các phân tử thuận từ có electron độc thân. o MINDO (Modified Intermediate Neglect of Differential Overlap) Quá trình tính toán của phương pháp này thực hiện theo thuật toán trường tự hợp, phương pháp này tính tích phân đẩy 2 tâm theo công thức của Klopmaan. Phương pháp này áp dụng để tính toán độ dài liên kết, nhiệt tạo thành, hằng số lực và thế ion hóa đối với các hợp chất hữu cơ. o AM1 (Austin Model 1) Phương pháp này được cải tiến từ phương pháp MNDO áp dụng cho các nguyên tử thuộc chu kỳ thứ nhất, thứ 2 và thứ 3 của bảng hệ thống tuần hoàn trừ kim loại chuyển tiếp. Phương pháp sử dụng thuật toán trường tự hợp. Phương pháp có nhiều ưu điểm tính toán các tính chất phân tử, tối ưu hóa hình học, năng lượng toàn phần và nhiệt tạo thành. o PM3 (Parametric Model 3) Phương pháp PM3 sử dụng thuật toán trường tự hợp, tương tự phương pháp AM1 với tất cả tham số tối ưu đầy đủ, về ý nghĩa nó có một tập tham số tốt nhất áp dụng cho nhiều nhóm nguyên tố thuộc nhóm chính, trừ kim loại chuyển tiếp. o ZINDO (Zerner’s INDO) Phương pháp ZINDO cải tiến từ phương pháp CNDO và INDO, phương pháp cũng sử dụng thuật toán trường tự hợp, có thể áp dụng cho hầu hết các hợp chất sinh học, các chất cơ kim và các chất hữu cơ. Thuật toán của ZINDO được tham số 12 hóa cho nhiều nguyên tố thuộc chu kỳ một, hai và ba và có thể mở rộng cho nhiều hợp chất. 4. Các bộ hàm cơ sở dùng trong các tính toán ab initio và DFT [3,51,53] - Bộ hàm cơ sở kiểu Slater: STO-nG (n = 2-6), những bộ hàm này có kích thước nhỏ, là hàm mũ của của PGSTO được xác định bằng cách khớp cho STO. Bộ hàm cơ sở STO-3G được sử dụng rộng rãi. - Bộ hàm cơ sở kích thước trung bình: 3-21G, 6-31G và 6-311G là các bộ hàm cơ sở thuộc nhóm phân chia hóa trị. - Bộ hàm cơ sở thuộc nhóm tương quan thích hợp (Correlation consistent Basis Set): aug-ccpVDZ, aug-cc-pVTZ, aug-cc-pVQZ, cc-pVDZ, cc-pVTZ, cc-pVQZ, cc-pV5Z và cc-pV6Z). Thường dùng để thực hiện tính toán cấu trúc tương quan điện tử. Các bộ hàm cơ sở aug-cc-pVDZ, cc-pVDZ và cc-pVTZ hay được sử dụng. II. PHƯƠNG PHÁP NHIỆT ĐỘNG THỐNG KÊ II.1. Các phương trình trạng thái Các phương trình trạng thái có thể được thành lập từ các dữ liệu thực nghiệm hay có thể được thành lập từ số liệu tính toán lượng tử. Các tham số của các phương trình trạng thái a và b được xác định từ các giá trị Tc và Pc cho phép tính toán được đầy đủ các tính chất của hệ. 1. Phương trình trạng thái van-der-Waals [15] Phương trình trạng thái này đơn giản nhất, thường được sử dụng để mô tả khả năng và không có khả năng thực hiện của phương trình trạng thái đơn giản. Phương trình này cần giá trị nhập Tc và Pc tương ứng với các giá trị a và b cho mỗi thành phần tương ứng: (5) Các tham số sử dụng cho hợp chất đơn: và (6) Các tham số đối với hệ chất lỏng hỗn hợp được tính từ quy tắc trộn. 2. Phương trình trạng thái Peng-Robinson [19] Phương trình này được sử dụng rộng rãi trong lĩnh vực nhiệt động kỹ thuật hóa học. Phương trình cho các kết quả dự đoán tỷ trọng chất lỏng tốt hơn phương trình 13 Soave-Redlich-Kwong. Phương trình cần 3 giá trị nhập cho mỗi hợp chất là Tc, Pc và yếu tố ω: (7) Các tham số sử dụng cho hợp chất đơn: với (8) với 3. Phương trình trạng thái Redlich-Kwong [20] Phương trình được cải tiến và phát triển từ phương trình trạng thái van-der- Waals; phương trình này cũng chỉ được sử dụng để tính toán mô tả khả năng thực hiện của một phương trình trạng thái. Phương trình này cần các giá trị nhập Tc và Pc cho mỗi hợp chất. (9) Các tham số cho đơn chất và (10) 4. Phương trình trạng thái Soave-Redlich-Kwong [20] Phương trình này phát triển từ phương trình Redlich-Kwong, ở đây tham số a được tính toán phụ thuộc nhiệt độ. Phương trình cũng cần có 3 tham số nhập Tc, Pc và yếu tố ω: (11) Các tham số của đơn chất với (12) 14 với 5. Phương trình trạng thái Peng-Robinson cải tiến [30] Phương trình này được thay đổi từ phương trình Peng-Robinson với sự dịch chuyển thể tích không ảnh hưởng đến việc tính toán áp suất hơi. Hàm α đưa ra với sự phụ thuộc của nhiệt độ của tham số a(T). Phương trình này khác biệt chủ yếu ở hàm alpha và tham số hiệu chỉnh: (13) với và với (14) 6. Phương trình trạng thái Deiters [47] Phương trình này được xây dựng trên cơ sở lý thuyết chuỗi cầu cứng nhiễu loạn với mục đích nhận được các giá trị nhiệt độ, áp suất và tỷ trọng tới hạn của các hệ đơn chất. Năm 1982 Deiters tính toán các giá trị cân bằng pha của hỗn hợp bậc hai từ phương trình chuỗi đa thức. Phương trình cần 3 tham số nhập a, b và c: (15) Với và pijk hằng số nhận được từ Deiters 15 Phương trình trạng thái của Deiters đã tính thành công cân bằng lỏng-hơi của nhiều hệ hóa học đơn chất và hỗn hợp bậc 2 và bậc 3 khác nhau. II.2. Các mô hình hệ số hoạt độ [8,9,10,11,17,31,56] Các mô hình NRTL, UNIQUAC và Wilson sử dụng để tính toán hệ số hoạt độ và cân bằng lỏng - hơi của các hệ chất lỏng ở các nồng độ và nhiệt độ khác nhau. Mô hình Wilson đưa ra để tính năng lượng tự do Gibbs ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ Λ−= ∑∑ N j ijj N i i E xx RT g ln (16) Hệ số hoạt độ được tính theo phương trình ∑∑∑ Λ Λ−+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ Λ−= N k N j kjj kik N j ijji x xx 1lnlnγ (17) Ở đây ]/)(exp[)/( RTVV iiijijij λλ −−=Λ ; Vi là thể tích mol lỏng của cấu tử i; λij - λii là năng lượng xác định bằng thực nghiệm. xi là phân số mol của cấu tử i. T là nhiệt độ, K. Mô hình NRTL được phát triển từ phương trình Wilson sử dụng cho hệ bậc hai ∑ ∑∑= N i ki N i ijijiN i i E G xG x RT g τ (18) Hệ số hoạt độ được tính theo phương trình ⎟⎟ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎜⎜ ⎝ ⎛ −+= ∑ ∑∑∑∑ ∑ N k kkj N k kjkjk ij N j N k kkj ijj k N k ki N j jjiji i xG Gx xG Gx xG xG τ τ τ γln (19) Các tham số tương tác bậc hai của mô hình ;)ln(/ TDTCTBA jijijijiji +++=τ )exp( jijijiG τα−= và αij = αji. Mô hình UNIQUAC được mở rộng từ mô hình Wilson sử dụng cho nhiều cấu tử 16 ∑ ∑∑ ∑ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−+Φ= i j jijii i i i i ii i i ii E xqxqz x xq RT g τθφ θ lnln 2 ln (20) Hệ số hoạt độ tính theo phương trình ∑∑∑∑ −+⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−Φ−+Φ+ Φ= N j N k kjk ijj ii N j jiji N j jj i i i i i i i i i qqqlxx lqz x τθ τθτθθγ lnln 2 lnln (21) Ở đây ])ln(/)(exp[ TDTCRTuuA ijijjjijijij ++−−=τ và 1))(2/( +−−= iiii rqrzl ; z = 10. Giá trị iΦ và iθ được xác định ∑ =Φ N k kk ii i xr xr ; ∑ = N k kk ii i xq qxθ (22) Ở đây qi = Awi / (2,5.109) với Awi là diện tích van der Waals và ri = Vwi/ 15,17 với Vwi là thể tích van der Waals; xi là phần mol. II.3. Lý thuyết mô hình COSMO Trong phần này gồm hai phần chính liên quan đến lý thuyết tính toán cân bằng pha của các hệ dung dịch chất lỏng trên cơ sở của lý thuyết lượng tử: mô hình COSMO (conductor-like screening model), COSMO-RS (real solvents), COSMO- SAC (segment activity coefficient) và các giá trị bề mặt sigma solvat hóa (gọi là giá trị sigma profile). Mô hình COSMO được Klamt [14] đưa ra đầu tiên là sự mở rộng của mô hình sovat hóa điện ly; tiếp theo là mô hình COSMO-SAC được phát triển bởi Lin và Sandler [15], các mô hình này với mục đích là dự đoán tính chất nhiệt động của pha lỏng. Phương pháp nhiệt động solvat hóa đòi hỏi bề mặt sigma được tính toán chỉ với cấu trúc phân tử và cơ học lượng tử. Giá trị bề mặt sigma là quan trọng trong các tính toán tính chất cân bằng pha đối với các hệ bậc hai. Cơ sở của mô hình COSMO là bề mặt sử dụng của dung môi đối với phân tử chất tan [14]. Về cơ bản COSMO đặt phân tử vào trong một hốc được tạo ra trong khuôn khổ môi trường giống nhau là dung môi. Hình 1 mô tả quá trình solvat hóa lý tưởng. 17 Hình 1. Mô tả quá trình solvat hóa lý tưởng Hình 2. So sánh quá trình solvat hóa lý tưởng Klamt và bình thường Trong Hình 1 và Hình 2, năng lượng tự do solvat hóa *soliGΔ đặc trưng cho sự thay đổi năng lượng tự do Gibbs lên quan đến sự dịch chuyển phân tử i ở vị trí cố định ở trạng thái khí lý tưởng tới vị trí cố định ở trong dung dịch S. Năng lượng tự do tạo thành hốc dung môi *caviGΔ đặc trưng bằng sự thay đổi năng lượng tự do Gibbs cần để tạo thành một hốc dung môi trong phạm vi dung dịch S của kích thước phân tử i chính xác. Năng lượng tự do điện tích *chgiGΔ đặc trưng bằng năng lượng tự do Gibbs cần để loại bỏ điện tích bề mặt khỏi bề mặt của hốc phân tử. Solvat hóa lý tưởng (Klamt) Solvat hóa bình thường 18 Phương trình biểu diễn năng lượng solvat hóa như là tổng năng lượng tạo thành hốc trong dung dịch và năng lượng tự do điện tích: * * * * * * * / / sol cav chg i i i sol cav res iS i i i S i S G G G G G G G Δ = Δ + Δ Δ = Δ + Δ + Δ (23) Ở đây *iSiGΔ và */resi SGΔ là năng lượng tự do solvat hóa lý tưởng và phục hồi 1. Giá trị bề mặt sigma (sigma profiles) Giá trị bề mặt điện tích sigma (sigma profiles) là phân bố xác suất mật độ điện tích bề mặt của một phân tử hay một hỗn hợp. Các mô hình COSMO xây dựng kích thước hình thể hốc phân tử trong khuôn khổ độ dẫn theo một nhóm các quy tắc và các hướng nguyên tử. Lưỡng cực phân tử và môment cao hơn nhận được điện tích từ môi trường bao quanh đến bề mặt của hốc. Bề mặt điện tích sigma pi của phân tử i là xác suất xác định một phần điện tích bề mặt [14,15,30]: ( ) ( ) / ( ) /i i i i ip n n A Aσ σ σ= = ( ) /i i i effn n A a σ σ= =∑ (24) ( )i iA A σ σ=∑ Ở đây ni(σ) là số phần bề mặt mật độ điện tích bề mặt σ; ni là số phần bề mặt tổng cộng bao quanh một hốc phân tử; Ai là diện tích bề mặt của hốc phân tử; Ai(σ) là diện tích bề mặt tổng cộng của tất cả các phần bề mặt với một mật độ điện tích bề mặt σ; ni(σ) và Ai(σ) có mối quan hệ Ai(σ) = aeff ni(ni(σ)) xác định bởi Lin [15]. Giá trị bề mặt điện tích sigma của hỗn hợp cũng là giá trị trung bình trọng số của bề mặt điện tích thành phần đơn và không giới hạn cho số cấu tử: ( ) ( ) ( ) ( ) i i i i i i i i i i i S i i i i i i i i x n p x A p x p p x n x x A σ σ σ σ = = = ∑ ∑ ∑ ∑ ∑ ∑ (25) Lin và Sandler đã xác định mật độ điện tích bề mặt trung bình của phần bề mặt σ* để tính mật độ điện tích bề mặt mới σ. 19 2 2 2 * 2 2 2 2 2 2 2 2 2 2 2 exp exp n eff mn n n n eff n eff m n eff mn n n eff n eff r r d r r r r r r d r r r r σ σ ⎛ ⎞−⎜ ⎟⎜ ⎟+ +⎝ ⎠= ⎛ ⎞−⎜ ⎟⎜ ⎟+ +⎝ ⎠ ∑ ∑ (26) Ở đây bán kính hiệu dụng /eff effr a π= là một tham số hiệu chỉnh của mô hình; dmn là khoảng cách giữa hai phần bề mặt điện tích m và n. Các phần bề mặt điện tích cặp đôi có mật độ điện tích σm và σn. Bán kính trung bình rav = 0,81764 sử dụng để tính bề mặt điện tích. Mật độ điện tích trung bình σm được tính theo phương trình: 2 2 2 * 2 2 2 2 2 2 2 2 2 2 2 exp exp n av mn n n n av n av m n av mn n n av n av r r d r r r r r r d r r r r σ σ ⎛ ⎞−⎜ ⎟+ +⎝ ⎠= ⎛ ⎞−⎜ ⎟+ +⎝ ⎠ ∑ ∑ (27) Một mẫu biểu diễn diện tích bề mặt trung bình nhận được từ tính toán lượng tử: Hình 3. Mật độ điện tích bề mặt của: nước, butylamin, 1-butanol và butyronitril 2. Mô hình COSMO-SAC Mô hình COSMO-SAC được Lin và Sandler phát triển [15,20]. Đầu tiên hệ số hoạt độ được xác định từ sự khác biệt năng lượng tự do Gibbs của cấu tử đơn trong Nước Butylamin 1-butanol Butyronitril 20 dung dịch; Sau đó Lin va Sandler đã sử dụng thêm số hạng tổ hợp Staverman- Guggenheim để cải thiện việc tính toán năng lượng tự do tạo thành hốc dung môi. * * / / / /ln ln chg chg SGi S i i i S i S G G RT γ γΔ −Δ= + (28) Số hạng tổ hợp Staverman-Guggenheim /ln SGi Sγ tính cho ảnh hưởng của kích thước và hình dạng của chất tan và dung môi. /ln ln ln2 SG i i i i S i i j j ji i i z q l x l x x φ θ φγ φ= + + − ∑ (29) Ở đây / ; / ,i i i j j i i i j j j j x q x q x r x rφ θ= =∑ ∑ và ( / 2)( ) ( 1)i i i il z r q r= − − − (30) Với xi là phân số mol của cấu tử i; z là số tọa độ được lấy là 10; thể tích chuẩn hóa ri = Vi/ r và diện tích bề mặt qi = Ai/ q; Ai và Vi là diện tích bề mặt hốc và thể tích hốc dung môi; r và q là tham số mô hình COSMO-SAC. Trong mô hình của Klamt [14,30], môi trường được coi là đồng nhất và độ dẫn hoàn hảo, trong khi Lin và Sandler lại xác định năng lượng tự do điện tích như là tổng của số hạng năng lượng tự do solvat hóa lý tưởng và năng lượng tự do phục hồi cần thiết để loại bỏ điện tích thuộc bề mặt điện tích khỏi bề mặt của hốc dung môi sau khi đặt chất tan trong hốc dung môi. Vì năng lượng solvat hóa lý tưởng đối với chất tan trong dung môi S, hay trong chất lỏng tinh khiết i là tương đương, phương trình hệ số hoạt độ đơn giản: * * / / / /ln ln res res SGi S i i i S i S G G RT γ γΔ −Δ= + (31) Lin và Sandler xác định năng lượng tự do phục hồi như là tổng của các kết quả bề mặt sigma và giá trị hệ số hoạt độ trên tất cả phần bề mặt. ** // ( ) ( ) ln ( )m m m resres Si S i m i i m s m GG n n p RT RT σ σ σ σ σ σ⎡ ⎤ΔΔ = = Γ⎢ ⎥⎢ ⎥⎣ ⎦∑ ∑ (32) Ở đây ΓS(σm) là hệ số hoạt độ đối với một phần bề mặt m của mật độ điện tích σ. Hệ số hoạt độ phần bề mặt trong dung dịch ΓS(σm) và trong chất lỏng tinh khiết Γi(σm) nhận được từ cơ học thống kê [14,15]: 21 ( , )ln ( ) ln ( ) ( ) exp ( , )ln ( ) ln ( ) ( ) exp n n m n s m s n S n m n i m i n i n Wp RT Wp RT σ σ σ σσ σ σ σ σσ σ σ ⎧ ⎫−Δ⎪ ⎪⎡ ⎤Γ = − Γ⎨ ⎬⎢ ⎥⎣ ⎦⎪ ⎪⎩ ⎭ ⎧ ⎫−Δ⎪ ⎪⎡ ⎤Γ = − Γ⎨ ⎬⎢ ⎥⎣ ⎦⎪ ⎪⎩ ⎭ ∑ ∑ (33) Ở đây ( , )m nW σ σΔ là năng lượng trao đổi cần để nhận một cặp của các phần bề mặt từ một cặp phần bề mặt trung hòa điện. Nó bao gồm những đóng góp từ tương tác tĩnh điện hay năng lượng không tương thích Emf, tương tác liên kết hydro Ehb và tương tác phi tĩnh điện Ene của các cặp bề mặt. Phương trình tính ( , )m nW σ σΔ như sau: 2'( , ) ( ) max[0, ]min[0, ] 2m n m n hb acc hb don hb W cασ σ σ σ σ σ σ σ⎛ ⎞Δ = + + − +⎜ ⎟⎝ ⎠ (34) Ở đây 'α là hằng số năng lượng không tương thích và là kết quả tính toán của giá trị yếu tố khả phân cực fpol = 0,64; hằng số α được xác định theo các quan hệ: 3/2 0 0,3 ' eff pol a f α α α = ∈ = (35) Phương trình cuối cùng để tính hệ số hoạt độ: / /ln ( )[ln ( ) ln ( )] ln m SG i S i i m S m i m i Sn p σ γ σ σ σ γ= Γ − Γ +∑ (36) Mô hình COSMO-SAC có thể viết như sau: ( ) ( , )ln ( ) ln ( )expi n m ni m i n i i p W A RT σ σ σσ σ⎧ ⎫−Δ⎡ ⎤Γ = − Γ⎨ ⎬⎢ ⎥⎣ ⎦⎩ ⎭∑ (37) Hay / / 1ln ' ( )[ln ( ) ln ( )] ln m SG i S i m S m i m i S eff p a σ γ σ σ σ γ= Γ − Γ +∑ (38) 3. Tính toán giản đồ pha Tính toán giản đồ pha đối với các hệ bậc 2 và bậc 3 ở các nhiệt độ khác nhau, sau khi thực hiện các tính toán mật độ điện tích bề mặt của phân tử và giá trị hoạt độ tương ứng. Áp suất tổng cộng được sử dụng trong các tính toán giản đồ pha: 22 0 tot i i i i p p x γ=∑ (39) Ở đây 0ip là các áp suất hơi của cấu tử tinh khiết đơn; ix và iγ là phân số mol của các cấu tử i trong pha lỏng và hệ số hoạt độ của cấu tử i được tính từ mô hình COSMO-SAC. Tính toán tính chất lý tưởng của pha khí được giả thuyết phân số mol pha hơi yi nhận được từ tỷ lệ áp suất riêng phần và áp suất tổng cộng. Hình 4. Sơ đồ tính toán cân bằng pha bằng mô hình COSMO Giản đồ pha được tính toán từ áp suất hơi của đơn chất 0ip ở nhiệt độ: 0 tot/i i i iy p x pγ= và với áp suất hơi riêng phần 0i i i ip p x γ= (40) Áp suất hơi ở nhiệt độ đã cho được tính theo phương trình Antoine. Hệ số phương trình nhận được từ 3 giá trị nhiệt độ và áp suất thực nghiệm tương ứng. Cấu trúc phân tử Tính toán hóa lượng tử Bề mặt phân tử lý tưởng: Phân bố mật độ điện tích trên bề mặt phân tử Cơ sở dữ liệu COSMO Tính toán nhiệt động thống kê. Giá trị sigma hỗn hợp Kết quả dữ liệu cân bằng pha của hỗn hợp C O SM O C O SM O -SA C 23 0ln( ) / ( )ip A B T C= − + (41) Ở đây A, B, C là hệ số phương trình, T là nhiệt độ ở áp suất 0ip tương ứng. Đánh giá sai số RMSE đối với kết quả dự đoán từ mô hình COSMO dựa vào công thức: 2 cal exp 1 (ln ln ) 2 n RMSE γ γ= −∑ (42) Ở đây n là số điểm dữ liệu, lnγcal và lnγexp là giá trị tính toán COSMO và giá trị thực nghiệm hay số liệu đã được công bố trong các tài liệu. II.4. Lý thuyết mô phỏng Monte Carlo[39,40,41,43,44,] 1. Cơ học thống kê mô phỏng toàn cục Monte Carlo Kỹ thuật mô phỏng Monte Carlo toàn cục GEMC (Gibbs Ensemble Monte Carlo) là một kỹ thuật hiệu quả đặc biệt để nghiên cứu cân bằng pha [44,45]. Trong kỹ thuật này việc xem xét một vùng lớn như một hộp ở các điều kiện không đổi N, V, T. Hộp được chia làm 2 vùng nhỏ I và II có thể tích là VI và VII có số hạt là NI và NII, và năng lượng EI và EII. Được biểu diễn ở Hình 5. Hình 5. Sơ đồ hệ mô tả mô phỏng toàn cục Monte Carlo GEMC Thể tích và số hạt tổng cộng của hệ được thỏa mãn: V = VI + VII và N = NI + NII Các vùng riêng biệt được xem xét như thể tích của các pha cân bằng của hệ vĩ mô; bề mặt phân cách pha giữa hai vùng không có ý nghĩa về mặt vật lý, chỉ là bề mặt tưởng tượng. Như vậy, các hạt đặt bên trong mỗi vùng nhỏ được xem xét sâu bên 24 trong các pha đặc trưng bởi các điều kiện của hệ nhỏ I và II. Hai hộp nhỏ này tạo thành và tồn tại cùng nhau và đặc trưng ở các điều kiện không đổi NVT. Trong quá trình thực hiện mô phỏng, có 3 loại dịch chuyển khác nhau: - Thay thế vị trí các hạt: mỗi hộp nhỏ được xem xét độc lập. Vì các hạt chuyển động trong các hộp như nhau, nên có thể chọn hộp I để xem xét tính chất này, trong không gian này trạng thái của hạt luôn thay đổi với xác suất exp[-βνI], ở đây νI là năng lượng cấu hình của hệ trong pha I. Các trạng thái này được tạo ra khi sử dụng giải thuật Metropolis: o Một hạt trong hộp I được chọn ngẫu nhiên o Hạt được thay thế vị trí ngẫu nhiên o Dịch chuyển hạt được chấp nhận với xác suất tối thiểu Imov(1,P )min ; ở đây ImovP là tỷ lệ xác suất giữa trạng thái cũ và trạng thái mới: I I new mov I old exp( ) exp( ) P βνβν −= − Hoặc I I mov exp( )P β ν= − Δ (43) Ở đây ΔνI là năng lượng trao đổi đối với sự chuyển dịch. Những sự dịch chuyển này làm cân bằng các hạt về nhiệt trong hai hộp. Một sự dịch chuyển bị loại bỏ thì sự dịch chuyển sẽ được thực hiện lại trong chuỗi các trạng thái Markov. - Thay đổi thể tích: quá trình này được bắt đầu với hộp I, như là phần của toàn cục không gian cân bằng đẳng nhiệt, đẳng áp (NIPT). Vì sự dịch chuyển của hộp I là tương quan với sự dịch chuyển của hộp II. Có thể xem hộp II là phần đẳng nhiệt, đẳng áp NIIPT. Sự dịch chuyển bao gồm: o Chọn sự thay đổi thể tích ngẫu nhiên ΔV đối với hộp I. o Đối với hộp I, tỷ lệ xác suất của trạng thái mới và cũ được xác định: (44) o Thể tích tổng cộng của hai hộp luôn được bảo toàn, sự thay đổi thể tích đối với hộp I, có nghĩa là thay đổi thể tích hộp II đồng thời -ΔV. Như vậy thay đổi tạo ra các trạng thái đẳng nhiệt, đẳng áp (NIIPT) ở hộp II: 25 o Giả sử áp suất P trong pha I và II là bằng nhau; đây là ngưỡng chuẩn cân bằng của hai pha I và II. o Xác suất thay đổi thể tích đối với hộp II là (45) o Tỷ lệ toàn thể xác suất đối với sự thay đổi thể tích (46) o Sự dịch chuyển thay đổi thể tích được chấp nhận với xác suất min(1,Pvol) o Nếu sự dịch chuyển không chấp nhận, trạng thái cũ được chấp nhận. - Trao đổi hạt: xét hộp I được đặc trưng (μVIT). Tạo ra một hạt trong hộp I. o Vị trí trong hộp một được chọn ngẫu nhiên o Đặt một hạt ở một điểm đã chọn ngẫu nhiên o Một sự dịch chuyển tạo một cấu hình mới đối với hộp I. Tỷ lệ xác suất đối với trạng thái mới và cũ là: (47) Ở đây λ, μ là bước sóng của hạt và hóa thế; hệ số hoạt độ được định nghĩa là z: o Vì tổng số hạt luôn được bảo toàn, việc tạo ra một hạt ở hộp I phải tương ứng loại bỏ một hạt ở hộp II. Hộp II là một mẫu biểu diễn toàn cục ở điều kiện không đổi (μVIT). Với μ trong mỗi hộp như nhau vì hai hệ cùng tồn tại. Để hoàn thiện sự dịch chuyển trao đổi, một nguyên tử trong hộp II được chọn ngẫu nhiên và được loại bỏ khỏi hộp II. o Tỷ lệ xác suất của các trạng thái mới và cũ của hộp I là 26 Hoặc (48) o Xác suất của tất cả các giai đoạn: (49) o Sự trao đổi được chấp nhận với xác suất min(1,Pex). Các dịch chuyển loại bỏ được xem xét lại trong chuỗn Markov. o Để thỏa mãn khả năng tráo đổi của hệ vi mô, việc tạo thành và loại bỏ trong hai hộp phải được thực hiện bằng xác suất bằng nhau. Điều này dẫn đến việc chọn hộp I và hộp II để cố gắng tạo ra xác suất trong suốt quá trình mô phỏng. o Nếu hộp II được chọn, các ký hiệu I và II được chuyển đổi trong phương trình (49) 2. Hàm thế tương tác phân tử [4,5,16] Trong mô phỏng hệ cân bằng pha, phương pháp nhiệt động thống kê được sử dụng và dựa trên giả thuyết tất cả các cấu hình có thể có của hệ thỏa mãn các điều kiện biên đã cho như năng lượng, thể tích và số hạt là luôn xảy ra là bằng nhau. Hệ toàn cục sẽ là ở hầu hết cấu hình có thể về phương diện thống kê. Nhiệt động học thống kê luôn nghiên cứu tính chất vi mô của hệ nhiệt động sử dụng lý thuyết xác suất. Nhiệt động thống kê đưa ra những giải thích mức độ phân tử ở mức độ định lượng nhiệt động như nhiệt, năng lượng hay entropi. Mục tiêu của nhiệt động thống kê là để hiểu biết và giải thích các tính chất của vật liệu vĩ mô có thể đo ở các số hạng tính chất của các hạt và sự tương tác giữa chúng. Để xem xét tính chất vật lý của hệ, đặc biệt các lực van der Waals, các kỹ thuật mô phỏng phân tử thực hiện các giải thuật để tạo ra bề mặt năng lượng đã biết như là các thế tương tác nguyên tử. Các giải thuật mong muốn đơn giản, nhưng thực tế 27 là phức tạp cần thiết để hoàn thành điều kiện trước đó. Một hàm tổng quát đối với một thế tương tác nguyên tử: (50) Ở đây Vi các đóng góp nguyên tử đơn; V2 và V3 là đóng góp tương tác giữa 2 và 3 hạt. Hàm thế tương tác phân tử nhận được hoặc từ tài liệu khi khớp cho các số liệu thực nghiệm hoặc từ các tính toán hóa học lượng tử hoặc có thể tính toán trực tiếp các số hạng khác nhau trong hàm thế. Đặc biệt khi thực hiện khớp một hàm thế với dữ liệu thực nghiệm là rất quan trọng. Tính toán thế năng là tất yếu để hiểu được tính chất hút và đẩy giữa các phân tử. Tương tác phân tử là kết quả của hai ảnh hưởng vùng gần (short-range) và vùng xa (long-range). Ảnh hưởng tĩnh điện, hút và đẩy là những mẫu cho sự tương tác vùng xa. Tương tác tĩnh điện nhận được từ phân bố điện tích giữa các phân tử. Ảnh hưởng có thể hoặc hút hoặc đẩy giữa các cặp. Tuy nhiên đóng góp quan trọng nhất là ảnh hưởng hút của những sự đẩy tăng lên từ sự dịch chuyển electron. Tương tác vùng gần (short-range) xảy ra ở khoảng cách tách biệt phân tử nhỏ, có sự xen phủ của các hàm sóng phân tử. Về lý thuyết có thể tính toán được tương tác phân tử theo nguyên lý thứ nhất. Tuy nhiên trong thực tế nguyên lý thứ nhất hoặc phương pháp ab initio giới hạn cho những hệ đơn giản. Hơn nữa ảnh hưởng của tương tác phân tử được biểu diễn bằng một vài kiểu thế tương tác phân tử. Hàm thế tương tác phân tử đã được phát triển và sử dụng rộng rãi để mô phỏng: - Thế tương tác phân tử Lennard-Jones 12-6 [39]: Thế này được sử dụng nhiều, thích hợp tốt cho sự tương tác giữa các phân tử. (51) Ở đây U(rij) là thế năng cấu hình; rij là khoảng cách giữa các hạt; εij và σij là độ sâu của thế và vị trí điểm giao nhau của thế. 28 - Thế tương tác phân tử Buckingham Exponential-6:[5] Nhiều công thức ở đó số hạng đẩy r-12 được thay thế bằng hàm số mũ thực tế hơn về lý thuyết, những tham số này bao gồm trong thế Buckingham: (52) Ở đây có 3 tham số hiệu chỉnh: mijr là bán kính năng lượng cực tiểu; εij là độ sâu của thế; αij điều khiển độ dốc của thế đẩy. - Thế tương tác phân tử hai tâm BFW[21] Đã có những hàm thế tương tác phân tử được phát triển chính xác được đưa ra cho một nguyên tử đặc biệt hay một nhóm các nguyên tử. Một tính chất trong tất cả các thế là sự đẩy mà được xử lý bằng số hạng mũ. Thế Barker-Fisher-Watts hai tâm (two-body): [ ]5 2 2 62 6 0 0 ( ) ( 1) exp (1 ) jii j i j C U r A x x x ε α δ + + = = ⎡ ⎤= − − −⎢ ⎥+⎣ ⎦∑ ∑ (53) Ở đây x = r/rm và các tham số khác nhận được bằng cách khớp thế cho dữ liệu thực nghiệm đối với các tâm phân tử, các hệ số virial và hệ số tương tác vùng xa (long-range). - Thế tương tác phân tử ab initio[5,7,40,41,42,43,46] Hầu hết các thế tương tác hai tâm (two-body) nhận được từ việc khớp thế cho số liệu thực nghiệm hai tâm (two-body). Một cách lựa chọn khác là nhận được một hàm thế bằng cách khớp cẩn thận hàm đã chọn cho dữ liệu bề mặt thế năng tương tác phân tử được tính toán từ phương pháp lượng tử ab initio. Các hàm thế ab initio đã được phát triển cho các tương tác nhiều vị trí trên phân tử: o Hàm thế ab initio tương tác 3 vị trí được phát triển bởi nhóm Naiker và Sum [43] sử dụng mô phỏng chất lỏng HCl: 3 3 2 min 10 1 1 0 ( ) {1 (1 exp[ ( )]) } ( ) 4 i jij ij ij ij ij ij ij i j ij q q U r r r f r r β α δ πε= = ⎡ ⎤= − − − − − +⎢ ⎥⎢ ⎥⎣ ⎦∑∑ (54) 29 Với 10 10 0 ( ) ( ) 1 ! ij ij k r ij ij ij ij k r f r e k δ δδ − = = − ∑ Ở đây βij tham số độ sâu của thế; αij tham số độ rộng của thế; minijr vị trí thế năng cực tiểu đối với tất cả tương tác giữa vị trí i và j trên phân tử a, b; f10(δijrij) hàm damping Tang-Toennies [54]. o Các hàm thế ab initio tương tác 5 vị trí được phát triển bởi nhóm Tat và Deiters [46] dùng mô phỏng các chất lỏng với phân tử nhiều tâm: 5 5 1 2 1 1 6,8,10 0 ( ) [ ( ) ( ) ] 4 ij ij ij r i jij n ij e ij a ijn i j n ij ij q qCU r D e f r f r r r α πε − = = = = + +∑∑ ∑ (55) 5 5 ( ) 2 1 2 1 1 6,8,10,12 0 ( ) ((1 ) 1) ( ) ( ) 4 ij ij ij ij r i jij n ij e ij b ijn i i n ij ij q qCU r D e f r f r r r α β πε − − = = = ⎡ ⎤⎛ ⎞= − − + +⎢ ⎥⎜ ⎟⎜ ⎟⎢ ⎥⎝ ⎠⎣ ⎦ ∑∑ ∑ (56) 25 5 2 1 22 1 1 3,4,5 0 ( ) ( ) ( ) 4 ij ij ij r r i jij n ij e ij a ijn i i n ij ij q qCU r D e f r f r r r α β πε − + = = = ⎡ ⎤= + +⎢ ⎥⎢ ⎥⎣ ⎦∑∑ ∑ (57) Với hàm damping 2( 2) 151( ) (1 )ij ij r ijf r e δ− − −= + , 2 ( ) 1 ij ijra ijf r e β−= − và 102 0 ( ) ( ) 1 ! ij ij k r ij ij b ij k r f r e k δ δ− = = − ∑ Ở đây ijeD , αij và βij là tham số hiệu chỉnh; qi, qj điện tích vị trí; rij khoảng cách giữa các vị trí; ijnC là các hệ số đẩy thuộc vùng xa (long-range); o Các hàm thế ab initio tương tác nguyên tử được phát triển bởi nhóm Eggenberger [16] đối với khí hiếm neon: 2 2 2 2 4 6 0 0 0 10 8 6 1 3 5 7 8 7 0 0 0 ( ) r r ra a a a a a ij r r rU r a e a e a e a a a a a a − − −⎛ ⎞ ⎛ ⎞ ⎛ ⎞− − −⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞= + + + + +⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠ ⎝ ⎠ (58) Ở đây a0 là bán kính Bohr và các tham số còn lại không có ý nghĩa vật lý. Đối với hệ đơn chất các tham số có thể nhận được từ việc khớp thế cho các số liệu thực nghiệm. Nhưng đối với các hệ hỗn hợp các cấu tử khác nhau, tham số thế được xác định bằng một quy tắc trộn tổ hợp Lorentz – Berthelot 1,ij ij i jkε ε ε= (59) 2, ( ) / 2ij ij i jkσ σ σ= + (60) Ở đây k1,ij và k2,ij là các hệ số trộn của quy tắc tổ hợp. 30 3. Kỹ thuật mô phỏng [39,45] Phương pháp Gibbs có thể được sử dụng cho những hệ nguyên tử hay phân tử. Đối với những hệ phân tử, những dịch chuyển kiểu NVT cũng sẽ bao gồm những thay đổi về hướng và cấu dạng. Những dịch chuyển kiểu μVT sẽ cũng bao gồm các hướng ngẫu nhiên của phân tử. Chúng ta có thể dễ dàng mở rộng phương pháp Gibbs cho hỗn hợp. Tất cả các dịch chuyển là như nhau trong trường hợp đơn chất, không kể những trao đổi hạt đối với các hộp NI và NII số các phân tử của các cấu tử được trao đổi. Trong hỗn hợp của A và B có các cấu tử A, B là I II I IIA A B BN , N , N , N với I II I IIA A B BN=N + N + N + N . Để chắc chắn khả năng chuyển về hệ vĩ mô, các bước tiếp theo phải được thực hiện theo những sự dịch chuyển trao đổi: - Chọn hoặc là hộp I hoặc là hộp II để tạo ra một xác suất bằng nhau - Chọn xác suất tùy ý cố định mà những cấu tử được trao đổi. Chọn A 20%, B 80% theo thời gian. Cách chọn khác dựa vào việc chọn lựa một nửa thời gian cho mỗi cấu tử. Không thể chọn một hạt ngẫu nhiên từ một trong các hộp để tráo đổi, vì khả năng chuyển đổi của hệ vi mô phải được thực hiện để tạo ra hay loại bỏ mỗi cấu tử một cách tách biệt. - Ngay khi một cấu tử được chọn, một cách ngẫu nhiên chọn một vị trí để tạo ra và chọn một cách ngẫu nhiên một trong những cấu tử đó để loại bỏ. Phương pháp mô phỏng có thể bao gồm áp suất không đổi đối với hỗn hợp. Điều này là hiệu quả cho việc chỉ định áp suất khi mô phỏng. Trong quá trình mô phỏng thể tích hệ tổng cộng không thay đổi, những thay đổi thể tích của hai hộp xảy ra độc lập với nhau. Các điều kiện chấp nhận là như nhau đối với mô phỏng NPT của hai vùng độc lập. Áp suất toàn cục Gibbs không đổi chỉ có ý nghĩa trong trường hợp của các hỗn hợp, vì đối với các cấu tử đơn số các biến số có thể đưa ra độc lập. Đối với mô phỏng Monte Carlo, điểm lưu ý đầu tiên là bắt đầu với cấu hình hộp tâm mặt với số hạt được đưa ra cho mỗi vùng và đặt ngẫu nhiên trong hộp. Các vòng lặp mô phỏng bao gồm sự thay thế vị trí của tất cả các hạt (sự dịch chuyển NVT), sự thay đổi thể tích và số hạt trao đổi. Các thay thế vị trí được chọn ngẫu nhiên từ phân bố đồng nhất. Thay thế vị trí lớn nhất được hiệu chỉnh trong cả hai vùng để cho một tỷ lệ chấp nhận trung bình 50%. Sau khi tất cả các hạt đã được thay thế vị trí, sự thay đổi thể tích được thực hiện. Những thay đổi thể tích được 31 chọn ngẫu nhiên từ phân bố đồng nhất và những thay thế lớn nhất được hiệu chỉnh để có tỷ lệ chấp nhận 50%. Những tính toán năng lượng và áp suất được thực hiện sau khi thay đổi thể tích. Những hiệu chỉnh vùng xa cho thế tương tác phân tử cần được thực hiện. Điều này thực hiện lấy tích phân thế tương tác từ một nửa độ dài của hộp (gọi là thế cut-off) và thiết lập g(r) = 1. Để thực hiện những trao đổi hạt trong quá trình mô phỏng nó được quyết định với một xác suất bằng nhau ngẫu nhiên mà hộp sẽ nhận một hạt và hộp kia sẽ loại bỏ một hạt. Một điểm được chọn ngẫu nhiên trong hộp nhận, một hạt được chọn ngẫu nhiên để loại bỏ trong hộp mất. Đối với những trao đổi đa cấu tử, các hộp nhận và hộp mất được chọn ngẫu nhiên và với một xác suất bằng nhau, thì sự giống nhau của hạt trao đổi được xác định ngẫu nhiên khi sử dụng một xác suất cố định cho mỗi cấu tử. Xác suất được thiết lập tốt nhất và thích hợp để thay thế vị trí cho tất cả các cấu tử. Hóa thế có thể dễ dàng tính toán theo giai đoạn mô phỏng và được sử dụng để kiểm tra quá trình mô phỏng. Để thực hiện điều này năng lượng tương tác được tính toán theo các bước trao đổi hạt tương đương với năng lượng hạt “ma” của phương pháp Widom. Đó là thế năng của một bước trao đổi thế của hạt i, νi,ex có thể được sử dụng theo công thức: βμi,r = -ln〈exp(-βνi,ex)〉 (61) Ở đây μi,r là hóa thế của cấu tử i; νi,ex là thế tương tác của các cấu tử i với phần còn lại của chất lỏng trong hộp nhận, bao gồm những hiệu chỉnh vùng xa. Bởi vậy hóa thế có thể được đánh giá trong cả 2 vùng. Điều này có thể sử dụng để kiểm tra sự cân bằng của hệ. 4. Tính toán cân bằng lỏng-hơi và tính chất nhiệt động [38,39] Từ kết quả mô phỏng Monte Carlo toàn cục, nhiệt độ tới hạn Tc/ K, tỷ trọng tới hạn ρc/ g.cm-3 và thể tích tới hạn Vc/ cm3.mol-1 của các đơn chất được tính bằng việc khớp phi tuyến cho các dữ liệu tỷ trọng cân bằng của các pha ở các nhiệt độ khác nhau dựa vào mối quan hệ: βρρ ρρρ )( )( 2 21 21 c cc TTB TTA −=− −+=− (62) 32 Ở đây ρl và ρv là tỷ trọng pha lỏng và pha hơi; β là hệ số mũ tới hạn (β = 0.325). A và B là các tham số hiệu chỉnh. Áp suất tới hạn Pc/ MPa được tính toán từ phương trình Antoine [39,]. Áp suất hơi và nhiệt hóa hơi được ΔvapH được tính bằng phương trình Clausius-Clapeyron: ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ +Δ−= 212 1 11ln TTR H P P vap (63) Trạng thái chuẩn P0 = 0.101 Mpa, quan hệ này được viết lại: R S TR H P P vapvap Δ+Δ−= 1ln 0 (64) Ở đây T1 và T2 là các nhiệt độ ở áp suất P1 và P2. Độ dốc và giá trị chặn lnP tương ứng với biến 1/ T là các hệ số ΔvapH và ΔvapS III. CÁC NHÓM HỢP CHẤT Trong đề tài này cấu trúc phân tử điện tử 3 chiều của các hợp chất được sử dụng trong nghiên cứu được nhận từ cơ sở dữ liệu của Viện tiêu chuẩn và kỹ thuật quốc gia (National Institute of Standards and Technology, NIST) của Mỹ [6,13]. Để thực hiện các mục tiêu nghiên cứu của đề tài, các hợp chất được phân thành 2 nhóm hợp chất chính hiện đang được chú ý nhiều trong các ngành công nghiệp kỹ thuật sản xuất hóa chất, kỹ thuật dược phẩm và xử lý môi trường. Các nhóm hợp chất chính: - Các đơn chất vô cơ - Nhóm hợp chất là chất lỏng hữu cơ - Nhóm dược chất hữu cơ Nhóm hợp chất được sử dụng nhiều trong thực tế, chiếm lượng lớn trong kỹ thuật sản xuất hóa chất cơ bản hiện nay; các đơn chất và hỗn hợp của các chất này dùng trong công nghiệp thực phẩm, công nghiệp hóa chất và làm dung môi trong công nghiệp. Việc nghiên cứu và hiểu biết về cân bằng lỏng-hơi hay lỏng-lỏng của các đơn chất và hỗn hợp khác nhau của chúng là vấn đề quan trọng đóng góp vào việc tạo ra hỗn hợp nhiên liệu lỏng và các ứng dụng làm dung môi công nghiệp. 33 Các hợp chất lỏng hữu cơ hiện nay được sử dụng nhiều trong các ngành công nghiệp như làm các chất tẩy rửa và cũng được làm các dung môi khác nhau trong công nghiệp sản xuất sơn. Các chất thải công nghiệp chứa các chất ion lỏng khi đi vào môi trường cũng làm ô nhiễm môi trường nước. Vì vậy hiểu biết về tính chất hóa lý của các hợp chất lỏng ion là việc quan trọng đóng góp vào việc giải quyết các vấn đề môi trường, và giúp các ngành công nghiệp tạo ra những sản phẩm công nghiệp quan trọng như các loại sơn nước, sơn bề mặt, và sơn chống thấm cho các hệ thống bê tông, hay chế tạo ra các loại sơn dẫn điện. Nhóm các hợp chất tự nhiên và đặc biệt là các dược chất có hoạt tính sinh học hiện nay đang là đối tượng được nghiên cứu tìm kiếm và thăm dò. 34 PHÖÔNG PHAÙP I. ĐIỀU KIỆN THỰC HIỆN I.1. Dữ liệu cấu trúc Để thực hiện đề tài này, các cơ sở dữ liệu cấu trúc phân tử được lấy từ cơ sở dữ liệu trực tuyến NIST [6,7] và của hội hóa học Mỹ [13]. I.2. Phần mềm Phần mềm được sử dụng trong nghiên cứu này: 1. Chương trình thương mại Phần mềm Sử dụng HyperChem Pro 8.03 Phục hồi cấu trúc phân tử 3D từ cơ sở dữ liệu DMol3 Tối ưu hóa cấu trúc, tính toán COSMO MS-EXCEL 2003 Chuyển đổi dữ liệu COSMO-SAC Tính toán hệ số hoạt độ Sigma-profile Fortran Tính toán giá trị bề mặt điện tích Origin 8.0 Xây dựng biểu đồ, đồ thị và giản đồ pha 2. Chương trình được viết ở mã nguồn C++ và Fortran Chương trình Sử dụng AutoFitting Thực hiện khớp các dữ liệu Virial 1.5 Tính toán Virial bằng thuật toán Monte Carlo Fortran 90 Tính toán VLE và LLE GEMC-NVT 2.0 Thực hiện mô phỏng Monte Carlo điều kiện NVT GEMC-NPT 2.0 Thực hiện mô phỏng Monte Carlo điều kiện NPT I.3. Sơ đồ thực hiện chung Đề tài được thực hiện theo sơ đồ chung sau: 35 Hình 6. Sơ đồ tính toán cân bằng pha theo các mô hình II. PHƯƠNG TRÌNH TRẠNG THÁI Để thực hiện các tính toán tính chất cân bằng pha của các hệ đơn chất và hỗn hợp các chất khác nhau, đề tài sử dụng hai phương trình trạng thái cơ bản là Peng- Robinson và phương trìn