计算化学公社

标题: 在Gaussian中对不同能量区间分批计算激发态 [打印本页]

作者
Author:
sobereva    时间: 2016-9-23 22:14
标题: 在Gaussian中对不同能量区间分批计算激发态
在Gaussian中对不同能量区间分批计算激发态
Calculating excited states in batches for different energy ranges in Gaussian

文/Sobereva @北京科音   2016-Sep-23


1 激发态的分批计算

从Gaussian09的后期版本开始,包括D.01、E.01,对CIS、TD、TDA激发态计算增加了设定能量下限关键词DEmin。比如TD(nstates=50,DEmin=5329),则程序只会计算出激发能高于5329/1000=5.329eV的50个激发态。

DEmin最大的实际用处就是补算激发态。比如对一个大分子,想计算UV-Vis光谱,一开始设了nstates=50,结果发现态数算少了,模拟出的光谱的波长下限不够,比如只计算到了300nm的地方,但还想要240~300nm的区段,可能再加估摸50个态才够。这时候不必重新算100个态,更省时的做法是先打开之前的输出文件看看最高的激发能,比如最高的是4.1213eV,就写一个新的输入文件,里面用TD(nstates=50,DEmin=4122) guess=read。这里4122即4.122eV,是比之前算的最高激发能高一丝的值。guess=read读取之前的chk用来省下重新做一遍基态SCF的时间。

有人认为一次算一部分激发态,总耗时比起一次算一大批激发态更省时间。最近一个西方的Multiwfn用户就这么干,他算一个很大的体系的UV-Vis光谱,要算500个态,他觉得一次算这么多吃不消,于是分成5批,以上述接力方式每次算100个态。至于这种做法是否真能节约时间,我做了一下测试,机子是Intel双路16核2.6G,64GB机子。

第一个测试体系比较小,是含30个原子的共轭体系,TD-PBE1PBE/6-311G*。一次性算150个态耗时为17m37s。而以接力方式每次算50个态,耗时为:
第一次50个态:6min54
第二次50个态(用了guess=read):7min39
第三次50个态(用了guess=read):8min48
可见总耗时反倒还高于一次性算150个态。值得注意的是,算越高能量区间耗时越多,因为算DEmin以上的态的时候实际上并非完全忽略掉DEmin下面的态,只不过TD迭代求解过程利用了一些技巧,虽然不用管下面的那些态的死活从而节省了迭代到收敛的时间,但不代表能够完全忽略掉它们。

第二个测试体系是胡萝卜素,含有96个原子,TD-M062X/6-31G*一次性算200个态耗时为2h51m54s,而以接力方式每次算100个态,耗时为:
第一次100个态:1h7m23s
第二次100个态(用了guess=read):1h26m12s
可见总耗时确实比一次算200个态低,但只低了一点。

限于时间和条件,笔者没有做更多测试,上述测试看到接力方式计算可能导致耗时比一次性计算更高也可能更低,这必定和体系、接力方式有关系。笔者不排除像那个Multiwfn用户一样对巨大体系分5批算100态确实比一次性算500个态耗时更低的可能。

这种接力方式得到的所有的态的能量,和一次性计算得到的相对比,通过上述测试体系并未发现有任何差异,所以用DEmin得到的结果是安全的。


2 用Multiwfn对分批计算的情况绘制总光谱图

以上述接力方式计算激发态由于会有多个输出文件,没法直接用gview绘制成完整的UV-Vis光谱图(或ECD图),而通过十分灵活的Multiwfn则可以解决这个问题。没用过Multiwfn画光谱图的读者先看看《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD和VCD光谱图》(http://sobereva.com/224)。

我们就拿前面第一个测试体系为例。首先载入第一次算50个态的Gaussian输出文件,进主功能11,选UV-Vis,然后选-2 Export transition data to plain text file将激发能和振子强度都导出到当前目录下transinfo.txt中。然后将此文件改名为1.txt。此时里面内容如下
    50     2
       3.097200       0.884800       0.666667
       3.944900       0.000000       0.666667
       4.220900       0.200600       0.666667
       4.381800       0.051000       0.666667
       4.451200       0.000100       0.666667
       4.455700       0.151000       0.666667
...[略]
含义在Multiwfn手册3.13.2节介绍了。50代表此文件里记录了50个态,后面的2代表此文件里对每个态同时记录了激发能、振子强度、FWHM。

之后再载入第二次算50个态的Gaussian输出文件,还是如上步骤得到transinfo.txt,改名为2.txt。之后再对第三次算50个态的输出文件也这么干,得到3.txt。

之后,把2.txt和3.txt的第一行都去掉,然后把2.txt和3.txt的内容按顺序复制到1.txt的末尾。再把1.txt第一行里的50改成当前文件里的态数,即150。然后保存文件。此时1.txt就含有所有激发态信息了,可以用Multiwfn绘制150态合在一起的光谱图了。

启动Multiwfn,载入1.txt,然后按照一般的绘制UV-Vis光谱图的操作,看到的光谱就和一次性算150个态的时候完全相同了。
作者
Author:
我本是个娃娃    时间: 2016-9-23 22:19
还有这么强大的功能,坐看接下来的讨论
作者
Author:
乐嘻嘻嘻    时间: 2020-6-25 21:37
感谢老师的这篇博文,我找了个体系测试了一下,如图先计算了6个激发态,后来在计算了3个激发态的基础上补3个激发态,但是结果如图,发现振子强度和波长有区别,不知道这种情况应该怎么办,麻烦老师指点一下。谢谢老师。
作者
Author:
liyuanhe211    时间: 2020-6-26 03:45
乐嘻嘻嘻 发表于 2020-6-25 21:37
感谢老师的这篇博文,我找了个体系测试了一下,如图先计算了6个激发态,后来在计算了3个激发态的基础上补3 ...

有啥区别
作者
Author:
乐嘻嘻嘻    时间: 2020-6-26 08:57
liyuanhe211 发表于 2020-6-26 03:45
有啥区别

没有第五激发态,直接第六了…
作者
Author:
liyuanhe211    时间: 2020-6-26 10:32
乐嘻嘻嘻 发表于 2020-6-26 08:57
没有第五激发态,直接第六了…

考虑到第五和第六激发态在能量上几乎简并,先挑到了第六个也并不奇怪。但这不是“振子强度、波长等有区别”
另外为啥不直接算456了,而是要把3再算一遍
作者
Author:
乐嘻嘻嘻    时间: 2020-6-26 20:35
liyuanhe211 发表于 2020-6-26 10:32
考虑到第五和第六激发态在能量上几乎简并,先挑到了第六个也并不奇怪。但这不是“振子强度、波长等有区别 ...

谢谢您,拿了个体系试试,没在意是345还是456~
作者
Author:
Mr.zhen    时间: 2024-3-4 15:48
本帖最后由 Mr.zhen 于 2024-3-5 09:42 编辑

请问老师,我尝试后发现设置了DEmin的输入文件计算出来的结果与没有设置的完全相同,关键词是#p TD(nstates=50,DEmin=3243) guess=read b3lyp/6-31g(d,p),输入文件与前50个激发能已上传,不知是哪里出了问题,希望您能指点一二
作者
Author:
sobereva    时间: 2024-3-5 02:16
Mr.zhen 发表于 2024-3-4 15:48
请问老师,我尝试后发现设置了DEmin的输入文件计算出来的结果与没有设置的完全相同,关键词是#p TD(nstates ...

上传压缩后的输出文件
作者
Author:
Mr.zhen    时间: 2024-3-5 09:33
本帖最后由 Mr.zhen 于 2024-3-5 09:43 编辑
sobereva 发表于 2024-3-5 02:16
上传压缩后的输出文件

已上传前后两次计算的输出文件
作者
Author:
sobereva    时间: 2024-3-5 11:46
Mr.zhen 发表于 2024-3-5 09:33
已上传前后两次计算的输出文件

明显已经生效了
最低激发能高于你设的值

  1. Excitation energies and oscillator strengths:

  2. Excited State   1:  66.002-A      3.2463 eV  381.92 nm  f=0.0036  <S**2>=1088.824
  3.     306B -> 319B      -0.10866
  4.     307B -> 318B       0.16128
  5.     307B -> 321B      -0.16611
  6.     307B -> 323B       0.10875
  7.     307B -> 324B       0.13982
  8.     307B -> 329B       0.10818
复制代码





欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3