Luận văn Theo dõi chuyển động của hạt nhân hydro trong quá trình đồng phân hóa vinylidene/ acetylene bằng laser xung cực ngắn

ron sinh ra một lưỡng cực. Lưỡng cực này dao động cùng với sự dao động của electron trong trường ngoài. Chính sự dao động của lưỡng cực này phát ra các bức xạ điều hòa mà ta gọi là các HHG. Quá trình phát xạ HHG được minh họa ở hình 1.5. Hình 1.5. Mô hình ba bước Lewenstein. Do sự lập lại chuẩn tuần hoàn của quá trình này trong trường laser nhiều chu kỳ, phổ lưỡng cực phát ra là gián đoạn, chứa bậc lẻ của tần số laser. Tần số dao động tức thời  của lưỡng cực liên hệ với động năng của electron theo biểu thức kE . Điều này có nghĩa là động năng của electron tại thời điểm va chạm chuyển thành năng lượng của photon phát ra, tương ứng với một sóng hài có tần số  . Những tần số khác nhau của các sóng hài tương ứng với các quỹ đạo khác nhau của electron khi trở về, làm cho động năng của chúng tại thời điểm va chạm là khác nhau. Về mặt toán học, phổ các sóng hài do một đơn nguyên tử hay phân tử phát ra chính là khai triển Fourier của gia tốc của lưỡng cực. Cường độ của phổ I và pha của sóng hài  tại tần số  có thể viết tách biệt ra như sau:     ,24  dI  (1.9)     ,arg  d (1.10) trong đó, moment lưỡng cực chuyển trạng thái trong vùng của phổ là            .exp rdrikrerkad g   (1.11) Ở đây bó sóng electron khi quan sát một phân tử được mở rộng trong sự chồng chất của các sóng phẳng như      rdxikexpkac   . Trong đó  k là số sóng (mometum) ứng với tần số điều hòa  và   ka là biên độ phức của nó. Giả sử như định hướng trục phân tử dọc theo trục x trong trường laser  tE , phân cực trên mặt phẳng x-y với một góc  hợp bởi trục định hướng của phân tử với vectơ phân cực của laser. Thành phần song song của moment lưỡng cực được viết dưới dạng:                  ....-*,iS-exp-sin-cos sincos 2/ st ** 0 //                 tatattdtd tdtd i diD yx yx (1.12) Một cách gần đúng khi bỏ qua sự suy giảm của trạng thái cơ bản, ta có thể xem     1-tat*a  , trong đó:       tA-,tpdtd st  và        -tA-,tpd-td st là các moment lưỡng cực của trạng thái cơ bản và trạng thái kế tiếp.       t -t '' st /dttA,tp là monent động lượng cổ điển tại những điểm đứng yên A là vectơ thế . Tương tự ta có thành phần vuông góc của moment lưỡng cực:                      ...-tata,tiS-exp-tE-tdsin-tdcos tdcos-tdsin 2/i ditD * styx 0 * y * x              (1.13) Những tác động cổ điển tại những điểm đứng yên lên electron truyền trong laser là:        ' t - 2 ' 2 -, , dtI tAtp tS t pst              , (1.14) với pI là thế ion hóa của phân tử. Photon năng lượng điều hòa có năng lượng cho bởi: ),(17.3max tUI PP  (1.15) với )(tUP là thế trọng động được tính tại thời điểm t photon phát ra. Hình dạng HHG phát ra từ acetylene khi sử dụng laser xung nhiều chu kỳ (2.1014W/cm2; 10fs) được đưa ra ở hình 1.6 như một ví dụ. Ta nhận thấy, đỉnh của những bậc điều hòa đầu tăng nhanh, theo sau đó là một vùng phẳng, cuối cùng là vị trí cut-off được xác định theo công thức .17.3 PP UI  (1.16) Hình 1.6. HHG thu được từ acetylene. Từ đó, ta có HHG điển hình gồm 3 vùng, như trong hình 1.7. Hình 1.7. HHG điển hình. Vùng thứ nhất gồm những bậc đầu tiên nơi mà có cường độ lớn hơn nhiều so với cường độ những thành phần phổ khác, vì hiệu suất chuyển đổi của chúng tỷ lệ với cường độ laser, mà cường độ này tăng lên tới bậc nhiễu loạn phi tuyến của quá trình; vùng này tuân theo lý thuyết nhiễu loạn bậc thấp (LOPT). Vùng thứ hai gọi là vùng phẳng (pleateau), đây là một dải rộng các đỉnh gần như có cùng chiều cao. Vùng thứ ba gọi là vùng cut-off, gồm một nhóm nhỏ các bậc điều hòa lớn với chiều cao giảm nhanh theo thứ tự bậc. 1.3. Chương trình tính toán LewMol 2.0 Đây là chương trình mô phỏng tính toán HHG do nguyên tử, phân tử phát ra khi tương tác với laser xung cực ngắn, cường độ mạnh. Mô hình này như sau: trước tiên trường điện của laser sẽ ion hóa phân tử làm cho điện tử ở vân đạo ngoài cùng (HOMO) bay ra vùng liên tục; sau đó trường điện tiếp tục gia tốc điện tử vừa được giải thoát trong bước thứ nhất; tiếp theo trường điện đổi chiều kéo điện tử quay lại tương tác với phân tử mẹ và phát ra sóng thứ cấp với tần số gấp nhiều lần tần số laser ban đầu. Vì sóng phát xạ ngay ở bước thứ ba khi phân tử va chạm với điện tử nên HHG mang thông tin về phân tử. Thông tin này được lấy trong thang thời gian cỡ femto giây (thời gian của phản ứng hóa học) cho nên có tầm quan trọng đặc biệt. Theo Lewenstein, laser tương tác chủ yếu với electron lớp ngoài cùng (tức HOMO). Vì vậy, để thu được sóng hài do nguyên tử, phân tử tương tác với laser, để đơn giản, thay vì thiết lập quá trình tương tác giữa laser với nguyên tử, phân tử; ta chỉ thiết lập những tính toán đối với HOMO của nguyên tử, phân tử. Ta thực hiện các bước sau : - Sử dụng Gaussian - một phần mềm được ứng dụng trong hóa lượng tử - tính toán các thông tin của phân tử như khoảng cách giữa các nguyên tử, các MO, đặc biệt là HOMO. Tất cả các thông tin này sẽ chứa trong file *.out. - Sau khi thu được các thông tin trên. Ta sử dụng get_wf_xie-modifyO2.f. Đây là một đoạn mã được thiết lập trên lập trình FORTRAN dùng để thu được các thông tin về hệ hàm cơ sở, các MO và đặc biệt là HOMO. - Cuối cùng ta sử dụng LewMol_2.2.f, đây cũng là một đoạn mã viết trên FORTRAN 7.0 sử dụng mô hình Lewenstein để tính toán HHG phát ra do nguyên tử, phân tử tương tác với laser cường độ mạnh. Source code này được viết bởi nhóm nghiên cứu đại học Kansas Hoa Kỳ và nhóm Đại học Sư phạm TP.HCM. Phần cơ sở lý thuyết của đoạn mã này đã được trình bày ở phần trên. Chương 2 Mô phỏng quá trình đồng phân hóa vinylidene/acetylene Nội dung chính của chương này là mô phỏng quá trình đồng phân hóa vinylidene/acetylene. Trước tiên, mô hình mẫu phân tử C2H2 được giải thích thông qua lý thuyết lai hóa các orbital trong phân tử. Sau đó, chúng tôi đi trình bày lý thuyết mặt thế năng (Potential Energy Surface – viết tắc là PES). Từ lý thuyết này, chúng tôi đi khảo sát tổng quát động học phân tử C2H2 thông qua kết quả PES tính được từ chương trình Gaussian. Các trạng thái đồng phân cũng như các trạng thái chuyển tiếp của quá trình đồng phân hóa vinylidene/acetylene trên đường phản ứng được chỉ rõ từ kết quả PES này. Sau đó, phương pháp mô phỏng động học phân tử được trình bày cùng phép gần đúng Born-Oppenheimer giúp chúng tôi có kết quả mô phỏng động lực học cho phân tử C2H2 bằng tiến trình BOMD trong chương trình Gaussian. Ở đây, các điều kiện ban đầu để quá trình đồng phân hóa vinylidene/acetylene diễn ra cũng được thu nhận. 2.1. Mô hình mẫu phân tử C2H2 2.1.1. Sự hình thành mối liên kết trong phân tử C2H2 Để giải thích sự định hướng và độ bền của các hợp chất hữu cơ, Pauling và Slater vào năm 1931 đã đưa ra khái niệm về sự lai hóa các orbital trong nguyên tử [25, 29]. Các orbital hóa trị ở các phân lớp khác nhau thường có năng lượng và hình dạng khác nhau; để tạo ra liên kết bền thì các orbital này cần phải đồng nhất giống nhau nhưng định hướng khác nhau trong không gian. Sự lai hóa của các orbital trong nguyên tử theo đó đã tạo nên sự tổ hợp một số orbital trong nguyên tử. Khi orbital 2s của nguyên tử C tổ hợp với một hoặc nhiều orbital 2p thì sẽ xảy ra ba trường hợp sau: Orbital 2s + 1 Orbital 2p → 2 Orbital lai hóa sp + 2 Orbital 2p còn lại. Orbital 2s + 2 Orbital 2p → 3 Orbital lai hóa sp2 + 1 Orbital 2p còn lại. Orbital 2s + 3 Orbital 2p → 4 Orbital lai hóa sp3. Trong đó, orbital lai hóa sẽ được dùng trong liên kết sigma với nguyên tử khác, còn các orbital còn lại được dùng cho liên kết pi. Cụ thể là orbital lai hóa sp thường được dùng để liên kết với hai hoặc một nhóm nguyên tử , orbital lai hóa sp2 thường liên kết với ba nguyên tử và orbital lai hóa sp3 thường liên kết với bốn hoặc nhóm nguyên tử. Lai hóa sp3 là sự tổ hợp một orbital s với ba orbital p của một nguyên tử tham gia liên kết tạo thành bốn orbital lai hóa sp3 định hướng từ tâm đến bốn đỉnh của hình tứ giác đều, các trục đối xứng của nó tạo với nhau một góc ≈ 109.5°. Ví dụ: trong hình 2.1, phân tử CH4 có orbital 2s lai hóa với ba orbital 2p tạo thành bốn orbital lai hóa sp3, bốn orbital lai hóa sp3 xen phủ với orbital 1s của nguyên tử hyđro tạo thành bốn liên kết sigma, lúc này, góc liên kết trong phân tử là 109°28". Hình 2.1. Cấu trúc orbital của CH4 và sự lai hóa sp3. Lai hóa sp2 với sơ đồ được đưa ra ở hình 2.2 là sự tổ hợp một orbital s với hai orbital p của một nguyên tử tham gia liên kết tạo thành ba orbital lai hóa sp2 nằm trong một mặt phẳng, định hướng từ tâm đến các đỉnh của tam giác đều với góc liên kết là 120°. Lai hóa sp2 được gặp trong các phân tử BF3, C2H4... Hình 2.2. Sơ đồ lai hóa sp2. Lai hóa sp là sự tổ hợp một orbital s với một orbital p của một nguyên tử tham gia liên kết tạo thành hai orbital lai hóa sp nằm thẳng hàng đối xứng với nhau hướng về hai phía và tạo với nhau một góc 180°. Lai hóa sp được gặp trong các phân tử BeH2, C2H2, BeCl2... Sơ đồ của sự lai hóa sp này được thể hiện ở hình 2.3. Hình 2.3. Sự lai hóa sp của orbital cacbon. Acetylene có công thức hóa học là C2H2, trong đó mỗi nguyên tử cacbon có hai liên kết. Như vậy, phân tử này có sự lai hóa sp: hai orbital lai hóa sp xen phủ với nhau tạo ra một liên kết sigma sp-sp giữa hai nguyên tử cacbon, phần còn lại của hai orbital lai hóa sp này xen phủ với orbital s của hai hydro để tạo ra hai liên kết sigma với hai nguyên tử hyđro. Hai orbital p (py và pz) không tham gia lai hóa còn lại xen phủ bên với nhau từng đôi một tạo ra hai liên kết pi như minh họa ở hình 2.4. Sự kết hợp của liên kết sigma sp-sp và hai liên kết pi p-p tạo ra mối liên kết ba giữa hai nguyên tử cacbon như ở hình 2.5. Hình 2.4. Liên kết pi tạo bởi hai orbital p. Hình 2.5. Sự hình thành mối liên kết của phân tử acetylene. 2.1.2. HOMO phân tử C2H2 Phần mềm Gaussian giúp ta có thể xây dựng cấu trúc hình học giả định của phân tử bằng cách tính toán một cấu trúc cân bằng của phân tử này. Ta cần chọn một hệ hàm cơ sở, chức năng tối ưu hóa (optimization) của Gaussian sẽ cho ta vị trí của tất cả các nguyên tử trong phân tử khi hệ đạt được trạng thái cân bằng (trạng thái năng lượng cực tiểu). Từ các dữ kiện thu được, ta tiếp tục dùng Gaussian để tìm các orbital của phân tử cần quan tâm. Chẳng hạn với cấu trúc giả định đưa vào trong file input như sau: X Y Z C 2.820513 1.324786 0.000000 C 4.015513 1.324786 0.000000 H 1.759513 1.324786 0.000000 H 5.076513 1.324786 0.000000 Chọn phương pháp tính toán và hệ hàm cơ sở DFT/ 6-311+g(2df,2pd) kết hợp với các từ khóa cần thiết như “opt” và “freq”, chương trình Gaussian sẽ xử lý số liệu đưa vào cho ta những thông tin hữu ích chứa trong file output, ví dụ cấu trúc hình học hợp lý của phân tử: Chiều dài liên kết C-C = 1.1834 Angstroms (thực nghiệm: 1.1979Angstroms [7]) Chiều dài liên kết C-H = 1.0558 Angstroms (thực nghiệm: 1.0575Angstroms [7]) X Y Z C 0.00000000 0.00000000 0.598061 C 0.00000000 0.00000000 -0.598061 H 0.00000000 0.00000000 1.659912 H 0.00000000 0.00000000 -1.659912 Ta nhận thấy cấu trúc thu được rất gần với số liệu thực nghiệm. Ta có thể sử dụng tiếp những số liệu này để tìm orbital phân tử MO, đặc biệt là HOMO như ở hình 2.6 và hình 2.7. Với HOMO thu được ta có thể mô phỏng quá trình tương tác giữa laser và phân tử để thu được sóng hài bậc cao HHG phát ra. Hình 2.6. HOMO vẽ trên các mặt Oxy, Oxz và Oyz của phân tử acetylene. HOMO này có năng lượng -0,414 Hartree. Hình 2.7. Các MO khác của phân tử acetylene. Bảng 2.1. Thông tin phân tử của acetylene và vinylidene. acetylene C2H2 vinylidene H2CC Khối lượng phân tử 26,038g/mol 26,038g/mol Thế ion hóa 11.4eV 11.47 eV Cấu trúc cân bằng X Y Z C 0.00000000 0.00000000 0.598061 C 0.00000000 0.00000000 -0.598061 H 0.00000000 0.00000000 1.659912 H 0.00000000 0.00000000 -1.659912 X Y Z C 0.000000 0.000000 -0.474329 C 0.000000 0.000000 0.815862 H 0.935848 0.000000 -1.026370 H -0.938128 0.000000 -1.022825 2.2. Mặt thế năng và các trạng thái đồng phân của phân tử C2H2 2.2.1. Mặt thế năng Năng lượng của phân tử là hàm phụ thuộc vào vị trí cũng như các tọa độ bên trong của các hạt nhân (như độ dài liên kết, góc liên kết và góc phẳng nhị diện). Khi hạt nhân chuyển động, các điện tử sẽ lập tức được sắp xếp lại; do đó, sự thay đổi năng lượng của hệ phân tử sẽ tương ứng với mỗi sự thay đổi của cấu trúc phân tử, sự thay đổi này có thể được quan sát rõ rệt thông qua mặt thế năng (PES). Như vậy, PES là hình vẽ trực quan mô tả hàm thế năng của phân tử theo cấu trúc phân tử. Số chiều của PES phụ thuộc vào số nguyên tử trong phân tử. Nếu phân tử có N nguyên tử thì nó sẽ có ba bậc tự do cho chuyển động tịnh tiến của toàn phân tử và có hai hoặc ba chuyển động quay nếu nó có cấu trúc tuyến tính hay không tuyến tính. Do đó, PES sẽ có (3N-5) hay (3N-6) chiều, dẫn đến sự phức tạp trong quá trình tính toán số học một PES lượng tử đầy đủ. Cụ thể, đối với phân tử hai nguyên tử, ta có đồ thị hai chiều giữa khoảng cách liên nguyên tử trên trục x, và giá trị năng lượng ở mỗi khoảng cách đó trên trục y. Lúc này PES đơn giản chỉ là một đường thế năng phụ thuộc vào độ dài liên kết như trong hình 2.8. Đối với phân tử lớn hơn, mặt sẽ có nhiều chiều hơn. Khi đó, mỗi điểm trên mặt này tương ứng với một cấu trúc hai chiều của phân tử với độ cao của mặt tại điểm đó cho ta năng lượng của phân tử. Hình 2.8. PES trong trường hợp đơn giản - phân tử hai nguyên tử. Hình 2.9. PES và các vùng đặc trưng. 2.2.2. Đặc điểm của mặt thế năng Hình 2.9 minh họa cho mặt thế năng cùng các đặc điểm từng vùng. Trên mặt này, ta sẽ chú ý đến một số điểm đặc biệt: điểm năng lượng cực tiểu với vùng trũng tương ứng với vị trí cân bằng (trạng thái bền) của phân tử; điểm năng lượng cực đại tương ứng với trạng thái chuyển tiếp tại vùng yên ngựa (saddle point). Năng lượng của các trạng thái này sẽ giúp ta có được năng lượng hoạt hóa của phản ứng. 2.2.2.1. Cấu trúc cân bằng Những phân tử có thể quan sát bằng thực nghiệm tương ứng với cực tiểu địa phương (local minimum) trên PES, tại đó ta có được cấu trúc cân bằng của phân tử. Để tìm ra vị trí này, trước tiên ta cần phải dự đoán cấu trúc của nó, sau đó sử dụng kỹ thuật tối ưu hóa (optimization) tìm được cực tiểu địa phương ở vùng lân cận. Tuy nhiên, vì quá trình tối ưu hóa thường định ra một vị trí tĩnh lân cận mà nó có thể là trạng thái cực tiểu hay chuyển tiếp, nên phép phân tích sự dao động thông qua các tần số nên được tiến hành trên cấu trúc đã tối ưu hóa để xác định rằng nó là cực tiểu địa phương. Ngoài ra vì cực tiểu địa phương không nhất định là cực tiểu toàn cục (global minimum) nên dự đoán chính xác chất ban đầu là cần thiết cho việc nhận được thông tin chính xác như độ bền tương đối, năng lượng phản ứng, độ lớn ngưỡng… 2.2.2.2. Cấu trúc trạng thái chuyển tiếp Cấu trúc trạng thái chuyển tiếp của một phản ứng hóa học tương ứng với điểm yên ngựa (saddle point) trên mặt thế năng. Như đã biết, việc tìm trạng thái chuyển tiếp cho một phản ứng hóa học là khá quan trọng, đặc biệt đối với các hệ phân tử lớn, vì vậy điểm yên ngựa có vai trò khá quan trọng trong việc giải thích cơ chế của phản ứng và động năng hóa học. Điểm yên ngựa cần tương ứng với đường chuyển động liên kết từ chất tham gia phản ứng đến sản phẩm. Một cách tốt hơn và chắc chắn hơn để xác định vị trí chuyển tiếp có chính xác hay không là thực hiện tính toán đường phản ứng. 2.2.2.3. Đường phản ứng Đối với một PES bất kỳ của phản ứng từ trạng thái A đến trạng thái B, ta có thể tìm được nhiều đường nối liên kết các trạng thái đầu và cuối. Tuy nhiên chỉ có đường ứng với sự thay đổi năng lượng tối thiểu từ chất phản ứng đến sản phẩm là đường mô tả cách dễ dàng nhất để xảy ra phản ứng, đường nối này được gọi là đường tọa độ phản ứng (reaction coordinate path hay Intrinsic reactant coordinate IRC) hay đường năng lượng cực tiểu (minimum energy path MEP). Theo đó, đường IRC là đường có độ dốc lớn nhất kết nối điểm yên ngựa với những vị trí cực tiểu địa phương gần nhất trên cả hai phía của dốc; đó là đường màu đỏ được chỉ ra trên Hình 2.9 kết nối cực tiểu của chất phản ứng hóa học đến cấu trúc chuyển tiếp A và cực tiểu của sản phẩm A. Đường IRC sẽ cung cấp cho ta thông tin về sự khác biệt rõ rệt giữa các mức năng lượng tương ứng với các trạng thái trong một quá trình phản ứng. Hình 2.10. Đường IRC điển hình. Như vậy, để nghiên cứu một phản ứng biến đổi của phân tử, ta cần nhận biết những vị trí cực tiểu của chất phản ứng, cũng như của sản phẩm tạo thành và điểm yên ngựa nối chúng với nhau. Đường phản ứng thường được tính toán để xác nhận điểm yên ngựa là điểm thật sự kết nối chất phản ứng và chất tạo thành. Vì phương trình Schrodinger chỉ có thể giải gần đúng cho hệ nhiều electron, độ chính xác của nó phụ thuộc vào bậc của phép gần đúng (phụ thuộc vào phương pháp và hệ hàm cơ sở), nên những vị trí sử dụng cho các điểm dừng (trạng thái cực tiểu hay chuyển tiếp) có thể hơi khác với cấu trúc thật của chúng. Do đó, cần có sự so sánh giữa các cấu trúc được dự đoán và những cấu trúc được quan sát từ thực nghiệm để đạt được độ chính xác cao trong phương pháp được dùng. 2.2.3. PES và các trạng thái đồng phân của phân tử C2H2 Đối với phân tử C2H2, PES sẽ phụ thuộc vào vị trí của cả bốn nguyên tử. Tuy nhiên, phân tử C2H2 có mối liên kết C-H rất linh động, đó là nguyên nhân dẫn đến các trạng thái bền năng lượng khác nhau của phân tử ứng với các vị trí khác nhau của hydro – gọi là các trạng thái đồng phân. Do đó, ở đây, để tiến hành mô phỏng động học phân tử, để đơn giản ta sẽ xét sự phụ thuộc của thế năng phân tử vào vị trí của hai nguyên tử hydro bằng cách cung cấp vận tốc cho một nguyên tử hydro. Khi nguyên tử hydro này nhận được năng lượng đủ lớn, nó sẽ chuyển động và chuyển phân tử sang dạng đồng phân khác. Trong quá trình nghiên cứu, để đặc trưng cho vị trí của nguyên tử hydro này, chúng tôi chọn hai tham số: R là khoảng cách từ nguyên tử hydro này đến khối tâm của liên kết C- C, A là góc hợp bởi vector R và trục phân tử. Điều này được thể hiện rõ trong hình 2.11. Khi đó, mặt thế năng PES thu được bằng cách sử dụng phần mềm Gaussian với phương pháp gần đúng DFT kết hợp với phép gần đúng Born-Oppenheimer và hai hệ hàm cơ sở B3LYP/6-311+G(2df,2pd) và PBE1PBE/aug-cc-pvtz. Khi này, ta bỏ qua năng lượng mất mát trong quá trình chuyển động của phân tử. Hình 2.11. Phân tử C2H2 ở trạng thái bất kỳ. * Sự khác biệt giữa hai hệ hàm cơ sở 6-311+G(2df,2pd) và aug-cc-pvtz: - 6-311+G(2df,2pd): là hệ hàm cơ sở Split-valence triple – zeta trong phép giải gần đúng HF, trong đó có tính thêm hiệu ứng momen và hiệu ứng khuyếch tán. Bộ hàm này bỏ qua vai trò của những điện tử bên trong trong việc tính toán. - aug-cc-pvtz: là một hệ hàm cơ sở được sử dụng phổ biến hiện nay (được phát triển bởi Dunning). Ngoài những khả năng như 6-311+G(2df,2pd), bộ hàm này còn tính thêm hiệu ứng tương quan, do đó kết quả tính được có độ chính xác cao hơn. Và cũng chính từ kết quả này, chúng tôi thấy chỉ khi tính đến hiệu ứng tương quan khi sử dụng hệ hàm cơ sở aug-cc-pvtz thì mới có thể thấy rõ sự cách biệt năng lượng giữa các trạng thái trong quá trình chuyển đồng phân C2H2 như kết quả PES thu được ở hình 2.12. Hình 2.12. PES ba chiều của C2H2. Để thấy rõ hơn sự cách biệt năng lượng giữa các đồng phân và các trạng thái chuyến tiếp của phân tử C2H2, chúng tôi vẽ mặt cắt thế năng của phân tử như ở hình 2.13. Hình 2.13. Mặt cắt thế năng của C2H2. Hình 2.14. Đường IRC của phân tử C2H2. Hình 2.14 là đường IRC của quá trình đồng phân hóa vinylidene/acetylene. Từ kết quả này ta thấy năng lượng của vinylidene cao hơn năng lượng của acetylene là 1.92eV. Điều này hoàn toàn phù hợp với kết quả thực nghiệm [7, 11] cho rằng sự chênh lệch năng lượng này là khoảng 2eV. Bảng 2.2 Thông số cấu trúc của các cấu hình từ vinylidene đến acetylene. Cấu hình RCC ( o A ) RCH ( o A )  (độ) R ( o A ) H (độ) Vinylidene 1.290 1.088 120.6 1.517 38.5 Chuyển tiếp 1 1.274 1.074 152.6 1.260 63.3 Cực tiểu địa phương 2 1.258 1.072 171.9 1.110 81.8 Chuyển tiếp 2 1.257 1.075 178.4 1.140 100.4 Acetylene 1.198 1.061 180.0 1.660 180.0 2.3. Gần đúng Born-Oppenheimer và tiến trình BOMD trong Gaussian 2.3.1. Phương pháp mô phỏng động học phân tử (MD) Phương pháp MD đầu tiên được đưa ra bởi Berni Alder và Tom Wainwright vào cuối những năm 1950 [4-5], trong quá trình nghiên cứu sự tương tác của những quả cầu rắn. Sự kiện này khi đó đã thu hút nhiều mối quan tâm của giới khoa học về mẫu chất lỏng. Sau đó, vào 1964, Rahman đưa ra sự mô phỏng đầu tiên về tương tác giữa các nguyên tử bằng thế năng liên tục cho argon lỏng [27]. Trong suốt những năm 1970, khi máy tính trở nên phổ biến, MD được phát triển cho những hệ phân tử phức tạp hơn. Năm 1974, Rahman và Stillinger tiến hành mô phỏng động lực học phân tử đầu tiên của hệ chất lỏng thật [30]. Sau đó là mô phỏng đầu tiên cho protein vào năm 1976 bằng cách sử dụng hàm năng lượng thực nghiệm. Từ đó đến nay, sự mô phỏng này được ứng dụng rộng rãi không chỉ trong vật lý, mà còn thường được dùng trong khoa học vật liệu, hóa sinh, vật lý sinh học, trong nghiên cứu protein và các phân tử sinh học. MD là một trong những công cụ chủ yếu trong nghiên cứu lý thuyết phân tử. Đây là phương pháp mô phỏng bằng máy tính dựa trên cơ học thống kê, trong đó các nguyên tử và phân tử tương tác với nhau, sự phụ thuộc thời gian của hệ phân tử được tính toán bằng cách giải số học phương trình chuyển động Newton. Từ đó, ta có các thông tin ở cấp độ vi mô trong thang thời gian nguyên tử như vị trí, vận tốc, sự dao động của các nguyên tử và sự thay đổi cấu trúc của phân tử. Những thông tin ở cấp độ vi mô này có thể được đưa về cấp độ vĩ mô dưới dạng áp suất, năng lượng, nhiệt dung thông qua cơ học thống kê, trong đó có nguyên lý Ergodic. Nguyên lý này nói rằng: khi hệ ở trạng thái cân bằng, giá trị trung bình trên tập hợp của một đại lượng vật lý của một hệ tại một thời điểm nào đó trùng với giá trị trung bình của đại lượng này tính theo thời gian của một hệ duy nhất. Mô phỏng MD trong nhiều khía cạnh rất giống với một thí nghiệm thực sự. Khi ta tiến hành một thí nghiệm thật, ta chuẩn bị sẵn mẫu vật cần nghiên cứu, kết nối mẫu vật này với các công cụ đo (nhiệt kế, máy đo độ nhớt…), đo các tính chất quan tâm trong suốt một khoảng thời gian xác định. Nếu phép đo của chúng ta có sai số thống kê thì thời gian càng dài, phép đo càng chính xác. Trong MD, chúng ta cũng có phương pháp tương tự. Đầu tiên ta chọn hệ mẫu chứa N hạt, những điều kiện ban đầu, giải phương trình chuyển động Newton cho hệ cho đến khi tính chất của hệ không thay đổi theo thời gian (cân bằng hệ). Sau sự cân bằng đó, chúng ta tiến hành phép đo thực sự. Giới hạn chủ yếu của phương pháp là yêu cầu về một nguồn máy tính mạnh. Kích cỡ mô phỏng (số hạt), bước thời gian và khoảng thời gian cần được xác định để việc tính toán có thể hoàn tất trong một khoảng thời gian hợp lý. Tuy nhiên, để tạo ra kết luận có giá trị thống kê, khoảng thời gian phải hợp với động lực học của quá trình tự nhiên. Thông thường, hệ phân tử và chuyển động của những hạt cấu thành nên hệ được mô tả chính xác bởi cơ học lượng tử. Tuy nhiên, ở đây chúng ta sẽ dùng cơ học cổ điển Newton để mô tả chuyển động của hạt nhân trong MD bằng cách sử dụng các phép gần đúng cho sự mô tả lượng tử. Hàm sóng phân tử sau khi được giải có thể tách riêng thành hai phần ứng với phần hạt nhân và phần điện tử theo Born-Oppenheimer, do đó chuyển động của chúng cũng có thể được tách ra và giải riêng lẻ: ta xem hạt nhân là hạt cổ điển, còn các biến số electron có thể được tách rời hoặc là được coi như các bậc tự do trong phương trình Schrodinger. Theo đó, MD được phân thành hai nhóm: 1/ MD cổ điển: lực trên hạt nhân được tính từ hàm thế năng. 2/ MD Ab initio: lực trên hạt nhân nhận được từ các tính toán cấu trúc điện tử. 2.3.2. MD cổ điển MD dựa trên phương trình chuyển động Newton. Nếu biết lực trên mỗi nguyên tử , ta sẽ xác định được gia tốc mỗi nguyên tử trong hệ. Tích phân của phương trình chuyển động lúc này mang lại quỹ đạo mô tả vị trí, vận tốc và gia tốc của các hạt theo thời gian. Từ quỹ đạo này, các giá trị trung bình của các đại lượng sẽ có thể được xác định, đồng thời trạng thái của hệ cũng có thể được dự đoán tại mọi thời điểm tương lai hoặc quá khứ. Lực có thể được viết dưới dạng gradient của thế năng phụ thuộc vào vị trí:   ,...,,2,1,...,, 21 NIRRRV R RMF N I III       (2.1) với III RMF ,, là lự