计算化学公社

标题: 计算有机分子的LUMO HOMO的疑问 [打印本页]

作者
Author: exity 时间: 2017-4-17 10:23
标题: 计算有机分子的LUMO HOMO的疑问
有请大神解惑，我试着尽量把问题描述的清楚一点。我想对一个有机分子进行计算后得到各个原子对总的LUMO HOMO的贡献，然后通过比较这些值来判断氧化剂进攻改分子时可能发生断键的位点，主要参考"Rapid Removal of Tetrabromobisphenol A by Ozonation in Water: Oxidation Products, Reaction Pathways and Toxicity Assessment"这篇文献。
文献里说可以通过计算LUMO HOMO的平方和之类的值判断断键位点并给出了下图：
(, 下载次数 Times of downloads: 99)

(, 下载次数 Times of downloads: 103)

于是我按照文献里方法用如下的关键词对一个三嗪类化合物进行了计算

# opt freq b3lyp/6-311g(d,p) scf=qc pop=all scrf=(solvent=water)

计算过程很顺利，也没有出错。

结果算完了以后，用“Molecular Orbital Coefficients"为关键词在输出文件里搜索，什么也找不到，我自己也不知道应该怎么去找LUMO HOMO的值，烦请大神解惑。是我的关键词不对，还是关键词参数不对？

此外，我有一系列的化合物，如下图
(, 下载次数 Times of downloads: 84)
我在用同一种氧化剂对这些物质进行降解研究，得到的结论是都在氨基的右侧发生断键后生成苯胺和对应的烷烃，但是母体的降解速率各不相同。

我分别对这些分子进行计算，得到每个分子中的氮原子和相邻的碳链上的碳原子对各个分子的LUMO HOMO的贡献，能否通过比较这些值来对降解速率进行说明？

我是计算新手，上面肯定有一些蠢问题，还请大神不要笑我！谢谢！

PS 输入输出文件见附件

作者
Author: sobereva 时间: 2017-4-17 18:29
Multiwfn是计算轨道成份最方便最强大的工具，这部分功能的操作和介绍看：
谈谈轨道成份的计算方法
http://sobereva.com/131
Multiwfn入门贴：
Multiwfn入门tips
http://sobereva.com/167

另外关于反应位点预测的文章看http://bbs.keinsci.com/forum.php?mod=viewthread&tid=384，特别是《亲电取代反应中活性位点预测方法的比较》（http://www.whxb.pku.edu.cn/CN/abstract/abstract28694.shtml）

作者
Author: exity 时间: 2017-4-17 19:41

sobereva 发表于 2017-4-17 18:29
Multiwfn是计算轨道成份最方便最强大的工具，这部分功能的操作和介绍看：
谈谈轨道成份的计算方法
http:/ ...

谢谢Sob,我也正在研究你编写的这个软件，命令行的输入方式很燃，像黑客帝国。不过，关于我这个问题，您能不能就事论事的点评一下呢？我真的很希望能搞清楚这个事情。

作者
Author: 冰释之川 时间: 2017-4-17 20:26

exity 发表于 2017-4-17 19:41
谢谢Sob,我也正在研究你编写的这个软件，命令行的输入方式很燃，像黑客帝国。不过，关于我这个问题，您能 ...

Sob老师回答得很清楚了：
（1）算各个原子对总的LUMO HOMO的贡献参考：
谈谈轨道成份的计算方法
http://sobereva.com/131

（2）不会用Multiwfn的参见：
Multiwfn入门tips
http://sobereva.com/167

（3）研究反应位点还有其他更好的方法，你可以参考使用其它方法：
反应位点预测的文章看http://bbs.keinsci.com/forum.php?mod=viewthread&tid=384，特别是《亲电取代反应中活性位点预测方法的比较》（http://www.whxb.pku.edu.cn/CN/abstract/abstract28694.shtml）

作者
Author: exity 时间: 2017-4-17 20:34

冰释之川发表于 2017-4-17 20:26
Sob老师回答得很清楚了：
（1）算各个原子对总的LUMO HOMO的贡献参考：
谈谈轨道成份的计算方法

好吧，也就是说我用opt freq b3lyp/6-311g(d,p) scf=qc pop=all scrf=(solvent=water)这样的关键词算出来的fch文件也能通过那个软件算出各自的LUMO HOMO恭贡献？

作者
Author: sobereva 时间: 2017-4-17 23:02

exity 发表于 2017-4-17 20:34
好吧，也就是说我用opt freq b3lyp/6-311g(d,p) scf=qc pop=all scrf=(solvent=water)这样的关键词算出来 ...

是

根本没有pop=all的关键词，pop=full写了也毫无意义。没事不要用scf=qc。

另见http://bbs.keinsci.com/forum.php?mod=viewthread&tid=334

作者
Author: exity 时间: 2017-4-17 23:30

sobereva 发表于 2017-4-17 23:02
是

根本没有pop=all的关键词，pop=full写了也毫无意义。没事不要用scf=qc。

谢谢sob！以后不会用好吧做句首发语词了，一定注意。

作者
Author: exity 时间: 2017-4-18 12:17

sobereva 发表于 2017-4-17 23:02
是

根本没有pop=all的关键词，pop=full写了也毫无意义。没事不要用scf=qc。

SOB老师，再次打扰您一下！按照您给的指示，我进行了计算。得到我的有机分子FCHK文件后用Multiwfn打开了，用Hirshfeld方法进行轨道分析，那么我又不清楚那个是LUMO 那个是 HOMO 轨道，我就用GVIEW打开FCHK文件，在右键 edit 里 MOs中，找到这样的图 (, 下载次数 Times of downloads: 88)
是不是通过这里就知道第38是LUMO 而 37是HOMO，然后回到Multiwfn中，选择对应的原子编号，然后轨道输入 37,38 输出的就是对应的LUMO HOMO贡献？
比如，我想知道氮原子(8号)对LUMO HOMO的贡献，我就输入，8 enter, 37,38 enter 得到这样的输出

Input atom index, e.g. 4
8
Input orbital range e.g.  20,25
37,38
Orb# Type    Ene(a.u.)    Occ Composition Population
37 Alpha&Beta    -0.2625 2.000    16.089%    0.321773
38 Alpha&Beta    -0.0329 0.000    3.206%    0.000000
Population of this atom in these orbitals: 0.321773

是不是Ene就是我要的能量，而Composition就是对轨道贡献的百分比？

同时Ene对应的就是这个文献插图的中的最上面结构中的蓝色和黑色小数字，不知道我这样理解对不对？

麻烦SOB老师点拨一二！
(, 下载次数 Times of downloads: 89)

作者
Author: sobereva 时间: 2017-4-18 20:21

exity 发表于 2017-4-18 12:17
SOB老师，再次打扰您一下！按照您给的指示，我进行了计算。得到我的有机分子FCHK文件后用Multiwfn打开了 ...

轨道成份算得对

ene是轨道能量，跟你图上黑字和蓝字完全没关系。具体什么含义在图注中应当会写

哪个是LUMO哪个是HOMO进入主功能0的时候直接在屏幕上提示了，根本用不着gview。看
使用Multiwfn观看分子轨道
http://sobereva.com/269

作者
Author: exity 时间: 2017-4-19 09:35

sobereva 发表于 2017-4-18 20:21
轨道成份算得对

ene是轨道能量，跟你图上黑字和蓝字完全没关系。具体什么含义在图注中应当会写

SOB老师，谢谢您，我这里有一个建议，不知道对不对！
对于我这样的超级没入门者来说，有些基本概念还是欠缺不少，但是迫于进度压力，必须进行一些量化计算。
我发现有一个地方可能对初学者（超级菜鸟）有误导：例如，你的“http://sobereva.com/269”例子中，对LUMO. HOMO轨道的判断如下

图形窗口上方中，点击Orbital info.--show all，在文本窗口中就会一目了然地输出所有轨道的信息。例如：
Orb: 19 Ene(a.u./eV): -0.340271    -9.2592 Occ:  2.000000 Type: A+B
Orb: 20 Ene(a.u./eV): -0.246291    -6.7019 Occ:  2.000000 Type: A+B
Orb: 21 Ene(a.u./eV): -0.246291    -6.7019 Occ:  2.000000 Type: A+B
Orb: 22 Ene(a.u./eV):    0.003620    0.0985 Occ:  0.000000 Type: A+B
Orb: 23 Ene(a.u./eV):    0.003620    0.0985 Occ:  0.000000 Type: A+B
显示了a.u.和eV为单位的所有轨道能量、占据数和自旋类型。此例A+B就是Alpha+Beta，即闭壳层轨道。从这里也可以立刻看出21和22号分别是HOMO和LUMO，20是HOMO-1，23是LUMO+1。
刚好21,22的能量从负值过渡到正值，我开始以为就是通过这样的正负变化来判断LUMO HOMO的（大误！），刚好在我自己的例子中

Orb: 55 Ene(a.u./eV): -0.280158    -7.6235 Occ:  2.000000 Type: A+B
Orb: 56 Ene(a.u./eV): -0.251841    -6.8529 Occ:  2.000000 Type: A+B
Orb: 57 Ene(a.u./eV): -0.250066    -6.8046 Occ:  2.000000 Type: A+B
Orb: 58 Ene(a.u./eV): -0.031105    -0.8464 Occ:  0.000000 Type: A+B
Orb: 59 Ene(a.u./eV): -0.010169    -0.2767 Occ:  0.000000 Type: A+B
Orb: 60 Ene(a.u./eV):    0.016059    0.4370 Occ:  0.000000 Type: A+B
Orb: 61 Ene(a.u./eV):    0.043075    1.1721 Occ:  0.000000 Type: A+B
Orb: 62 Ene(a.u./eV):    0.049899    1.3578 Occ:  0.000000 Type: A+B

正确的LUMO HOMO应该是58,57，可是我最开始的时候以为是59,60（通过正负变化判断的，其实这是错误的判断方法），和GVIEW里的结果是不一样的。然后突然我似乎顿悟了一样，是不是通过2.00000-0.00000的变化来判断LUMO HOMO，这里的2和0就是指占据轨道的电子数，2就是有电子占据的所以57是HOMO, 58就是LUMO。这样就比较没有误会了，SOB老师可以添加一句话，解释一下2,0的意义，就会对我这样的超级菜鸟帮助特别大，不知道我说的对不对，请批评。

另外，现在基本上我已经知道了怎么看各个原子对LUMO HOMO的贡献了，那么如果我的分子结构比较大，例如有60多个原子，我想知道其中的40个原子对LUMO HOMO的贡献，岂不是要输入40*2次（原子序数和轨道数），然后输出的结果也要分别的黏贴复制到一起才能再继续下面的统计分析的工作，有没有简便一点的方法，烦请不吝赐教！

作者
Author: sobereva 时间: 2017-4-19 10:00

exity 发表于 2017-4-19 09:35
SOB老师，谢谢您，我这里有一个建议，不知道对不对！
对于我这样的超级没入门者来说，有些基本概念还是 ...

没有任何误导，懂得HOMO、LUMO含义的人都不会犯这种错误

HOMO、LUMO序号是多少，屏幕上都直接输出了，都根本不用看占据数，没法更贴心了

Hirshfeld轨道成份分析界面里很清楚有一个定义片段的功能，定义完片段，看轨道成份的时候一下子就输出片段的贡献了

作者
Author: exity 时间: 2017-4-19 11:17

sobereva 发表于 2017-4-19 10:00
没有任何误导，懂得HOMO、LUMO含义的人都不会犯这种错误

HOMO、LUMO序号是多少，屏幕上都直接输出了， ...

明白了，谢谢您！你真是一本真人版的量化百科全书！

作者
Author: exity 时间: 2017-4-19 11:21

sobereva 发表于 2017-4-19 10:00
没有任何误导，懂得HOMO、LUMO含义的人都不会犯这种错误

HOMO、LUMO序号是多少，屏幕上都直接输出了， ...

Orbital 57 is HOMO, energy: -0.250066 a.u. -6.804646 eV
Orbital 58 is LUMO, energy: -0.031105 a.u. -0.846420 eV
LUMO/HOMO gap: 0.218961 a.u. 5.958227 eV 574.881600 kJ/mol

原来在这个地方，非常感激，真的是不能再贴心了。\笑中带泪

然后我用VMD输出LUMO HOMO图，发现可以实现原子编号（也是参考了你的教程），不过不能加原子名称，比如我要把C, N, S这样的原子符号标上去，请问应该怎么实现？

作者
Author: exity 时间: 2017-4-19 13:52

sobereva 发表于 2017-4-19 10:00
没有任何误导，懂得HOMO、LUMO含义的人都不会犯这种错误

HOMO、LUMO序号是多少，屏幕上都直接输出了， ...

SOB老师，我还有一个问题想请教。

应该怎么计算有机分子中各个原子在LUMO HOMO下的原子轨道？

比如下图这样
(, 下载次数 Times of downloads: 98)
这个计算结果就把各个分子轨道下原子轨道的贡献给出了，我想要得到这样的细节数据应该怎么输入关键词。

我现在用的关键词如下，得到OUT文件里没有上图这样的结果。

%nprocshared=8
%chk=DIA.chk
# opt freq b3lyp/6-311g(d,p) scrf=(solvent=water) scf=qc

作者
Author: lao7 时间: 2017-4-19 14:43

exity 发表于 2017-4-17 20:34
好吧，也就是说我用opt freq b3lyp/6-311g(d,p) scf=qc pop=all scrf=(solvent=water)这样的关键词算出来 ...

咱们一样弱，跟不上拍，呵呵！

作者
Author: exity 时间: 2017-4-19 14:45

lao7 发表于 2017-4-19 14:43
咱们一样弱，跟不上拍，呵呵！

不怕慢，就怕站

作者
Author: lao7 时间: 2017-4-19 14:47

exity 发表于 2017-4-19 14:45
不怕慢，就怕站

嗯，有时候担心大神们没有耐心。你说的这个功能，Multiwfn一个软件就能搞定。但我现在对这个H、L分析问题还是一知半解。追你帖子，好好看看！

作者
Author: exity 时间: 2017-4-19 14:50

lao7 发表于 2017-4-19 14:47
嗯，有时候担心大神们没有耐心。你说的这个功能，Multiwfn一个软件就能搞定。但我现在对这个H、L分析问题 ...

我研究出了一些，有一条路已经走通了，结果发现目的地不是我要去的地方，这才是悲剧呢。你现在卡到哪一步了？我们可以讨论一下。

作者
Author: lao7 时间: 2017-4-19 14:54

exity 发表于 2017-4-19 14:50
我研究出了一些，有一条路已经走通了，结果发现目的地不是我要去的地方，这才是悲剧呢。你现在卡到哪一步 ...

我的比较滞后，还在计算，计算还没有结束，但是下一步分析轨道占据情况。你现在面临的问题，估计就是我要面临的。我要做的是电荷分解分析，但同样用到L、M的结果。我现在连别人计算的结果都没有完全搞明白如何解读。

作者
Author: exity 时间: 2017-4-19 14:59

lao7 发表于 2017-4-19 14:54
我的比较滞后，还在计算，计算还没有结束，但是下一步分析轨道占据情况。你现在面临的问题，估计就是我要 ...

好的， keep in touch

作者
Author: sobereva 时间: 2017-4-19 21:49

exity 发表于 2017-4-19 13:52
SOB老师，我还有一个问题想请教。

应该怎么计算有机分子中各个原子在LUMO HOMO下的原子轨道？

如果你要轨道系数，pop=full直接就输出了

作者
Author: sobereva 时间: 2017-4-19 21:51

exity 发表于 2017-4-19 11:21
Orbital 57 is HOMO, energy: -0.250066 a.u. -6.804646 eV
Orbital 58 is LUMO, energy: - ...

把我提供的tcl脚本里的%1i改为%1i%a

作者
Author: beidiwu 时间: 2017-4-20 03:25
sob大神不但学问深，还有极好的耐心，令人尊敬！

作者
Author: exity 时间: 2017-4-23 20:03

sobereva 发表于 2017-4-19 21:49
如果你要轨道系数，pop=full直接就输出了

SOB老师，请问计算静电势时必须使用弥散函数吗？
比如我用b3lyp/6-311g(d,p)或b3lyp/6-311g**(d,p)会对我最终结果有影响吗？

作者
Author: sobereva 时间: 2017-4-24 05:56

exity 发表于 2017-4-23 20:03
SOB老师，请问计算静电势时必须使用弥散函数吗？
比如我用b3lyp/6-311g(d,p)或b3lyp/6-311g**(d,p)会对 ...

6-311g**(d,p)没有这种写法
给氢加弥散没用，6-311+G**比6-311G**对静电势的描述会略有改善，但还不如用def2-TZVP

作者
Author: exity 时间: 2017-4-25 16:55
本帖最后由 exity 于 2017-4-25 17:00 编辑

sobereva 发表于 2017-4-24 05:56
6-311g**(d,p)没有这种写法
给氢加弥散没用，6-311+G**比6-311G**对静电势的描述会略有改善，但还不如用 ...

谢谢您，SOB大神。我还有4个问题，请您指导一下。
Q1：在你的文章“使用multiwfn的定量分子表面分析功能预测反应位点，分析分子间相互作用”中第3节分子表面的静电势分析实例，中有提到“通常B3LYP/6-31G**级别的波函数已经能很好地重现分子表面的静电势分布了”，这里您推荐的是def2-TZVP，请问我应该怎么取舍？
Q2：计算羟基自由基的时候，在GV里画了一个水分子后删除了一个氢，然后在输出文件中把电荷改为0，自旋多重度改为2，请问这样设置对不对？
Q3：用# opt freq b3lyp/6-311g(d,p) scrf=(solvent=water) pop=full这样的关键词计算后，根据您的文章输出了静电势的彩色图，发现我的分子中氮原子附近的静电势是正的，可是在Murray, J.S. and Politzer, P. (2002) Encyclopedia of Computational Chemistry, John Wiley & Sons, Ltd.这篇综述中，不是说Vmin会出现在O,N,F,S...等原子附近吗？请问这要怎么理解文献和实际计算结构的差异？
(, 下载次数 Times of downloads: 51) (, 下载次数 Times of downloads: 56)
Q4：请问能否实现这样一种计算：固定两个分子之间的距离，例如5A左右，以此情况模拟分子在靠近后在各自的电场影响下本来的ESP发生了什么变化？
问题比较多，麻烦您了。

最后有一个小问题啊，真心没有冒犯您的意思，我知道你英语一定非常非常好的，可您不觉得SOB这个名字缩写有点怪怪的吗？还是和“I am a bit*h”有异曲同工之妙？这个问题可以不用理会。/抱拳 /抱拳

作者
Author: sobereva 时间: 2017-4-25 17:33

exity 发表于 2017-4-25 16:55
谢谢您，SOB大神。我还有4个问题，请您指导一下。
Q1：在你的文章“使用multiwfn的定量分子表面分析功能 ...

1 def2-TZVP定量精度更高，6-31G**足够保证定性正确
2 对
3 不能一概而论。跟具体所处的化学环境有关。你图中5N、7N都是正常情况，而8N、9N这种成键方式极为特殊
4 用Multiwfn对二聚体的静电势和单体的静电势求差，如同作密度差图。

完全不怪。能错误理解的人脑回路过于特殊

作者
Author: exity 时间: 2017-4-26 06:12

sobereva 发表于 2017-4-25 17:33
1 def2-TZVP定量精度更高，6-31G**足够保证定性正确
2 对
3 不能一概而论。跟具体所处的化学环境有关。 ...

未暴走初号机的意思？

作者
Author: exity 时间: 2017-4-26 06:23
本帖最后由 exity 于 2017-4-26 06:25 编辑

sobereva 发表于 2017-4-25 17:33
1 def2-TZVP定量精度更高，6-31G**足够保证定性正确
2 对
3 不能一概而论。跟具体所处的化学环境有关。 ...

老师啊，8N 9N就是碳上面接的氨基，这个特殊在什么地方呢？能不能细说一下。事实上，我还有几种化合物和这个结构非常接近，就是这个8N9N上还接着一个碳链（是一种除草剂，阿特拉津），我算出的那个ESP也是这样的，下面两个氮是正的。可我实际实验结果羟基自由基（亲电试剂），进攻的就是这两个氮上的碳链，后发生断键形成上个贴子中这样的产物，请问这算是计算结果和实际结果冲突了吗？是需要改进或换一种算法还是ESP这个解释工具只能弃之不用？另外我用fukui函数倒是能解释实验结果，这时计算结果和实验实际又符合了。很苦恼啊。

作者
Author: exity 时间: 2017-4-28 09:02

sobereva 发表于 2017-4-25 17:33
1 def2-TZVP定量精度更高，6-31G**足够保证定性正确
2 对
3 不能一概而论。跟具体所处的化学环境有关。 ...

大侠，我错了，我不该调侃你的，不该在这样一个严肃论坛里用戏谑的话，请您高抬贵手帮帮我吧。保证以后不再犯了。

作者
Author: sobereva 时间: 2017-4-30 19:57

exity 发表于 2017-4-26 06:23
老师啊，8N 9N就是碳上面接的氨基，这个特殊在什么地方呢？能不能细说一下。事实上，我还有几种化合物和 ...

你把氢取消显示但又没提，别人还以为那是悬着的单个氮，此时显然电子结构特殊
氮周围静电势为正显然是氨基上的带正电的氢导致的

作者
Author: exity 时间: 2017-5-2 14:30

sobereva 发表于 2017-4-30 19:57
你把氢取消显示但又没提，别人还以为那是悬着的单个氮，此时显然电子结构特殊
氮周围静电势为正显然是氨 ...

谢老师，不计前嫌，悉心解答。

作者
Author: jjf_sxnu 时间: 2017-5-6 11:30
受教了，回去我也试试。

欢迎光临计算化学公社 (http://bbs.keinsci.com/)