24小时热门版块排行榜

返回列表

【奖励】本帖被评价12次，作者zhou2009增加金币 11.1 个

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

zhou2009

版主 (著名写手)

[资源] 【zhou2009个人文集】电子密度差Δρ深入计算一例：HCl

电子密度差Δρ深入计算一例：HCl

我们知道，波（状态）函数和几率，是作为量子力学第一个基本假设（公理）提出来的。
在量子化学，ρ=Ψ*Ψ，它天才地沟通了电子的波粒二象。然而通常比较忽视作为粒子性的ρ，这主要是因为计算得到的“净电荷”以及其他电荷都有人为处理的失误，与经验、实验相距甚远，使人们对电荷失去信心，从而不理会电荷，这也累及到ρ。

其实ρ从Ψ的平方而来，它也与Ψ一样，是从量子化学得到的最原始、原生态数据，它本身是有来历、可信的。
这里的初步工作，是想从ρ出发，从新的视角和方法，认识电子在形成化学键时的变化。

我在上一篇帖子《电子密度ρ和密度差Δρ》说，在GsGrid中，对成键前后的几个ρ进行三维空间格点相减，从而得到Δρ的cube新文件。再次进入GsGrid，就会自动对Δρ的cube进行三维空间格点体积微元中的电子的正值部分和负值部分分别进行加和，这样一来，Δρ不仅是一幅直观的图象，而且有了电子净增（正值）、净减（负值）的具体数值，又精确定量了。

这样，对H2的Δρ进行计算，得到Δρ的电荷在键中间的（下称“中移”）净增加值为0.251229，在二H电荷净减少值共为0.250936。
对二聚水电子转移Δρ进行计算，转移量为：水-1孤对电子处（净减少）-0.39139，水-2的H—O处电子净增加为0.39128。

然而，象上面单纯的电子向键中间移动、单纯的电子在分子之间转移，毕竟是少数，大量的情形是电负性不同的原子之间成键，同时伴随着电子的中移和转移，它们在Δρ是重叠的，要想方设法理清这种重叠，判明电子中移和转移的确切量。
这需要算很多分子来进行分析考察，现在先拿一个HCl的计算例，说明这种计算的基本思路、过程。

关于Cl的C03计算：
Cl由于有一个单电子，计算时可用ROHF或者HF方法。
用ROHF方法算，单电子能级（-0.19956）在两个孤对电子简并能级之上，比H的轨道能级（-0.49982）还高出很多。
用HF方法算，总能量比ROHF低，单电子能级（-0.57891）在两个孤对电子简并能级之下，也在H的轨道能级之下，故本计算取了HF方法计算。

选择Δρ进行计算考察，因为它是原子或基团形成分子前后的电子密度差，反映了成键前后电子的净变化。
选择HCl这类分子，因为它形成σ键的MO单一，可以认为是Cl的单电子轨道与H的轨道衍生而来，Δρ只需要用相关的轨道，情况单一明确，而其它轨道都不会有电子得失上的变化。

图一便是这样作出来的Δρ的一个有代表性的截面。我们可以清楚地看到，成键后电子发生了什么样的净变化。总的看来，电子大量地由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中出现的电子量，空间微元dv也是以bolr^3为单位的，GsGrid会自动输出dv的大小。

由上面所讲可以看出，由于三维坐标格点文件只是以x坐标为基准在变化的，相应的三维坐标格点也正是这样排列的，并且无论分子的实际坐标取向如何。如果分子取向不是现在的x轴方向，而是y轴方向，我们将完全不能由划分x轴的区间来划分计算空间，而且又不能去划分y轴。
这在作Δρ的截面图时，本没有什么问题，只需要作图时交换一下坐标，图形就会旋转90度，按我们的习惯来展示图形。但这样作完全不能改变三维坐标格点以x坐标为准变化的排列。

而且实际计算还真的出了这个问题！
当我们用G03计算Cl原子时，我们需要用的那个单电子波函数，它竞是个Py，F原子的单电子波函数也是个Py，O原子有两个单电子，竞是个Py、Pz！
在G03中，我们目前还找不到什么指令或方法让单个原子的波函数改变取向。如果我们按照Cl原子单电子波函数Py取向安排HCl分子，它将是y轴方向的，上面的划分x轴的区间的方法就完全行不通了。在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。这时出现该波函数的展开系数都在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
这样，图一所呈现的从左到右正、负、正、负的四个值都分别求出了。

由此，初步计算HCl的转移电子为：0.2845。
用同样的方法，初步计算HF的转移电子为：0.3664。
推算转移电子的方法这里没有讲，目前还没有确定下来。

目前还需要算更多的分子来考察这方面的问题，要讨论的问题将会很多，可说的话自然也很多，且留在以后的帖子再说吧。

[ Last edited by yjcmwgk on 2010-6-16 at 21:25 ]

回复此楼