電子密度差Δρ深入計算一例：HCl

我們知道，波（狀態(tài)）函數和幾率，是作為量子力學第一個基本假設（公理）提出來的。
在量子化學，ρ=Ψ*Ψ，它天才地溝通了電子的波粒二象。然而通常比較忽視作為粒子性的ρ，這主要是因為計算得到的“凈電荷”以及其他電荷都有人為處理的失誤，與經驗、實驗相距甚遠，使人們對電荷失去信心，從而不理會電荷，這也累及到ρ。

其實ρ從Ψ的平方而來，它也與Ψ一樣，是從量子化學得到的最原始、原生態(tài)數據，它本身是有來歷、可信的。
這里的初步工作，是想從ρ出發(fā)，從新的視角和方法，認識電子在形成化學鍵時的變化。

    我在上一篇帖子《電子密度ρ和密度差Δρ》說，在GsGrid中，對成鍵前后的幾個ρ進行三維空間格點相減，從而得到Δρ的cube新文件。再次進入GsGrid，就會自動對Δρ的cube進行三維空間格點體積微元中的電子的正值部分和負值部分分別進行加和，這樣一來，Δρ不僅是一幅直觀的圖象，而且有了電子凈增（正值）、凈減（負值）的具體數值，又精確定量了。

這樣，對H2的Δρ進行計算，得到Δρ的電荷在鍵中間的（下稱“中移”）凈增加值為0.251229，在二H電荷凈減少值共為0.250936。
對二聚水電子轉移Δρ進行計算，轉移量為：水-1孤對電子處（凈減少）-0.39139，水-2的H—O處電子凈增加為0.39128。

    然而，象上面單純的電子向鍵中間移動、單純的電子在分子之間轉移，畢竟是少數，大量的情形是電負性不同的原子之間成鍵，同時伴隨著電子的中移和轉移，它們在Δρ是重疊的，要想方設法理清這種重疊，判明電子中移和轉移的確切量。
    這需要算很多分子來進行分析考察，現(xiàn)在先拿一個HCl的計算例，說明這種計算的基本思路、過程。

關于Cl的C03計算：
Cl由于有一個單電子，計算時可用ROHF或者HF方法。
用ROHF方法算，單電子能級（-0.19956）在兩個孤對電子簡并能級之上，比H的軌道能級（-0.49982）還高出很多。
用HF方法算，總能量比ROHF低，單電子能級（-0.57891）在兩個孤對電子簡并能級之下，也在H的軌道能級之下，故本計算取了HF方法計算。

選擇Δρ進行計算考察，因為它是原子或基團形成分子前后的電子密度差，反映了成鍵前后電子的凈變化。
選擇HCl這類分子，因為它形成σ鍵的MO單一，可以認為是Cl的單電子軌道與H的軌道衍生而來，Δρ只需要用相關的軌道，情況單一明確，而其它軌道都不會有電子得失上的變化。

圖一便是這樣作出來的Δρ的一個有代表性的截面。我們可以清楚地看到，成鍵后電子發(fā)生了什么樣的凈變化�？偟目磥恚�電子大量地由H轉移的Cl上了，H、Cl之間也有中移電子。



為了分析Δρ所展示的電荷，從圖一看，從左到右，依次是正（增加）、負（減少）、正(中移)、負（轉移與中移疊加）四個部分，首先就要將這四個部分分別算出相應的電荷量來。

在空間找到的電子量，本是指這樣一個積分的結果：∫∫∫ρ(r) dxdydz ，其中r是坐標矢量。在gsgrid中計算時，由于cube的格點已是離散化的，就成了∑[x]∑[y]∑[z]ρ(x,y,z) dxdydz，令dxdydz=dv，即dv是空間微元，也就成了∑[x]∑[y]∑[z]ρ(x,y,z)dv。所以格點文件中每個位置的ρ(x,y,z)乘以空間微元dv的積的加和就應當歸一或歸N，即1或者N個電子。

從cube文件的構成，或者更直接由GsGrid為每個格點配齊三維坐標計算的實際運作來看，它是從空間格點的起始點開始，每變一個x的步長，y、z就輪回一遍，三維坐標格點也正是這樣排列的。因此我們只要以x坐標為準，在x軸上選一個值a，對三維坐標格點的從起始點開始一直到x坐標值為a，對ρ進行加和，再乘以空間微元dv，就得到這一段空間的電荷量。于是從圖一我們找到正、負、正、負之間的x軸的分界點，分段求和，就可以計算這四段的電荷量。
當然為了劃分更加嚴格分明，實際的具體操作是將GsGrid的三維坐標格點文件在Excel中打開，并將這個ρ值列進一步分為正值、負值各一列。并對正值、負值各作出相應的截面圖，見圖二、圖三。在圖二、圖三上分別找到兩個正值和兩個負值的分界點。然后對兩個正值分別加和，兩個負值分別加和，再乘以空間微元dv，得到各自的電荷量。由于計算都是采用的原子單位，ρ是bolr^3中出現(xiàn)的電子量，空間微元dv也是以bolr^3為單位的，GsGrid會自動輸出dv的大小。


由上面所講可以看出，由于三維坐標格點文件只是以x坐標為基準在變化的，相應的三維坐標格點也正是這樣排列的，并且無論分子的實際坐標取向如何。如果分子取向不是現(xiàn)在的x軸方向，而是y軸方向，我們將完全不能由劃分x軸的區(qū)間來劃分計算空間，而且又不能去劃分y軸。
這在作Δρ的截面圖時，本沒有什么問題，只需要作圖時交換一下坐標，圖形就會旋轉90度，按我們的習慣來展示圖形。但這樣作完全不能改變三維坐標格點以x坐標為準變化的排列。

而且實際計算還真的出了這個問題！
當我們用G03計算Cl原子時，我們需要用的那個單電子波函數，它競是個Py，F(xiàn)原子的單電子波函數也是個Py，O原子有兩個單電子，競是個Py、Pz！
在G03中，我們目前還找不到什么指令或方法讓單個原子的波函數改變取向。如果我們按照Cl原子單電子波函數Py取向安排HCl分子，它將是y軸方向的，上面的劃分x軸的區(qū)間的方法就完全行不通了。在Excel中作改變將會太復雜。

幸好multiwfn 1.1新增加了這個功能，將波函數Py取向改為Px！具體作法是：

1、在用G03計算Cl原子時，增加關鍵詞output=wfn輸出波函數，在輸入的分子信息之后空一行，寫上輸出的波函數名稱和路徑，如：D:\Cl.wfn。
此時不再用關鍵詞cube計算格點文件，因為這時得到的單電子波函數是Py。

運行multiwfn，讀入Cl.wfn，按程序使用說明書操作，將選定的波函數Py取向改成Px。
在功能選擇菜單上選擇選改波函數的功能6，再選次級功能2，輸入分子軌道號碼，比如Cl單占據軌道就是7。這時出現(xiàn)該波函數的展開系數都在Py上，并且可以看到基函數Py對應的序號和應交換的相鄰Px的序號，反復選擇交換功能7，寫下要對應交換的Py序號和Px序號，用“，”號分隔。直到所有交換進行完畢，選擇功能8，自動將改好的文件存為new.wfn。視基組大小，只須改變幾個基函數即成。

再次運行multiwfn，讀入new.wfn，按程序使用說明書操作，就可計算Cl單占據軌道的電子密度三維格點文件（Cl.cub），此時取向已經是Px，它的格式與G03的cube文件完全一樣。

2、同樣對HCl、H按作Δρ的要求進行G03計算，也用output=wfn得到波函數，同樣進入multiwfn做成電子密度三維格點文件（HCl.cub、H.cub）。

3、運行GsGrid，將電子密度三維格點文件HCl.cub分別減去Cl.cub、H.cub，就得到HClΔρ的三維格點文件HCl-DD.cub。

4、運行GsGrid，讀入HCl-DD.cub，取xy截面，得到HCl-DD-xy.txt，它可以在sigmaplot中直接打開，作出等值線圖形。可以觀察HCl的Δρ。
如果Δρ圖形正確，將HCl-DD-xy.txt在Excel中打開，將Δρ值的列分成正值列和負值列。一方面可以用它們作一個Δρ的虛實線圖作為插圖。另一方面可以分別作出正值線圖（圖二）和負值線圖（圖三），找到計算ρ值和的分界點的值。HCl的Δρ正值線的分界點值確定為x=0.4，負值線分界點值確定為x=0.89。

5、運行GsGrid，讀入HCl-DD.cub，選擇功能1，就為每個格點配齊三維坐標了，它是一個文本文件如記為HCl-DD-xyz.txt。
同時GsGrid會自動對給出dv的值為：0.0014577，計算Δρ的正值加和和負值加和值，以及乘以dv后的正值（電子凈增加總量）、負值（電子凈減少總量）。我們還需要將這兩個總量進一步分開。

6、將HCl-DD-xyz.txt在Excel中打開，如果是80×80×80的格點數，在Excel達到512000行！前三列為xyz坐標，第四列為Δρ的值。進一步將第四列分為第五列（正值列）和第六列（負值列）。

點選正值列第一個值，看著x軸的值下拉滾動條到x=0.4的最后一行，點擊選中，點擊菜單上的∑進行這一段Δρ值的加和，再乘以dv，得到電子量為：0.4960
點選正值列x=0.4的最后一行，拉滾動條到此列末尾，點擊選中，點擊菜單上的∑進行這一段Δρ值的加和，再乘以dv，得到電子量為：0.0229
再，點選負值列第一個值，看著x軸的值下拉滾動條到x=0.89的最后一行，點擊選中，點擊菜單上的∑進行這一段Δρ值的加和，再乘以dv，得到電子量為：-0.2170
點選負值列x=0.89的最后一行，拉滾動條到此列末尾，點擊選中，點擊菜單上的∑進行這一段Δρ值的加和，再乘以dv，得到電子量為：-0.2960
這樣，圖一所呈現(xiàn)的從左到右正、負、正、負的四個值都分別求出了。

由此，初步計算HCl的轉移電子為：0.2845。
用同樣的方法，初步計算HF的轉移電子為：0.3664。
推算轉移電子的方法這里沒有講，目前還沒有確定下來。

目前還需要算更多的分子來考察這方面的問題，要討論的問題將會很多，可說的話自然也很多，且留在以后的帖子再說吧。

[ Last edited by yjcmwgk on 2010-6-16 at 21:25 ] 返回小木蟲查看更多