计算化学公社
标题: 使用Gaussian+PySOC在TDDFT下计算旋轨耦合矩阵元 [打印本页]
作者Author: sobereva 时间: 2018-3-15 05:11
标题: 使用Gaussian+PySOC在TDDFT下计算旋轨耦合矩阵元
2018-Dec-29注：有不少人表示PySOC运行不了。只要读者使用本文的操作系统、严格按本文的操作步骤计算，至少本文的例子是绝对能运行成功的。读者应当先确保本文例子能运行正常，由此确认安装方式和自己的操作过程无误后再去跑你自己的体系。如果你当前没有本文用的系统，可以用VMware虚拟机装，半小时就能装好。
2019-Feb-10注：笔者后来又写了《使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱》（http://sobereva.com/462），用起来明显比用Gaussian+PySOC方便得多，还更快、更准确，且完全免费，因此在笔者来看PySOC没什么使用价值了，本文仅适合了解计算原理之用。读者不要再以任何形式问笔者PySOC的使用问题，笔者不予回复，笔者一律建议用ORCA算SOC矩阵元。
使用Gaussian+PySOC在TDDFT下计算旋轨耦合矩阵元
Using Gaussian+PySOC to calculate spin-orbit coupling matrix elements with TDDFT
文/Sobereva @北京科音 First release: 2018-Mar-15 Last update: 2019-Mar-30
0 前言
本文会简要介绍旋轨耦合矩阵元的基本知识和计算原理，然后介绍怎么使用Gaussian结合完全开源免费的PySOC程序在研究激发态最常用的TDDFT方法下计算各个单-三重态之间的旋轨耦合矩阵元。其实PySOC使用起来一点也不难，希望看过此文的Gaussian用户不会出现就为了计算旋轨耦合矩阵元而去买昂贵的ADF程序的情况。本文对于理论部分只是简要提及，便于读者能够理解计算结果和计算意义，系统的相对论量子化学知识介绍可以参看相关专著，诸如Dyall写的Introduction to Relativistic Quantum Chemistry就很好。限于篇幅和表达方式的限制，具体、细碎、系统的东西在文中不好呈现，相关内容会在北京科音(http://www.keinsci.com)的高级量子化学培训班里讲授。
1 旋轨耦合矩阵元的基础知识
薛定谔方程是标量方程，求解得到的波函数是“一分量”的。而考虑相对论效应时与薛定谔方程对应的叫做Dirac方程，解出来的波函数用含有四个分量的spinor（旋量）表示，体现出电子与正电子、alpha与beta自旋间的耦合。由于基于Dirac方程计算实际化学体系耗时较高（但也没一些文章和书籍里说得那么离谱），为了使求解更容易进行，不同的人提出了不同方法将Dirac方程进行变换，使四分量哈密顿中的电子和正电子部分脱耦合，从而得到描述电子部分的二分量方程，此时其alpha和beta自旋还是耦合在一起的。变换成二分量的方法非常多，诸如FW(Foldy-Wouthuysen)、X2C、ZORA(Zero-order regular approximation)、DK(Douglas-Kroll)、IOTC等。这些二分量哈密顿算符可以写为非相对论部分、标量相对论部分和旋轨耦合（H_SO）这三类算符的加和。因此，旋轨耦合本质上来自于相对论效应。
在不考虑相对论效应的情况下，电子波函数的空间部分和自旋部分是可以分离的。旋轨耦合体现的是电子自旋与电子轨道运动的相互作用，因此当考虑旋轨耦合算符后，电子波函数的空间和自旋部分就不能分离了。因此，当考虑旋轨耦合后，就没有诸如单重态、三重态这样的概念了。当不考虑旋轨耦合时，不同自旋多重度之间的态之间哈密顿矩阵元精确为0，换句话说就是没有相互作用。但是把旋轨耦合算符引入哈密顿之后，不同自旋多重度的态之间就可以发生相互作用，因此最终得到的本征态将是不同自旋多重度的态的线性组合。
所谓旋轨耦合矩阵元，就是指<i|H_SO|j>这种积分，它体现出i、j两个电子态之间的旋轨耦合作用的大小。在能够计算旋轨耦合矩阵元的一些程序和计算了旋轨耦合矩阵元的文章中，它往往也被称为旋轨耦合常数（但注意旋轨耦合常数这个词也有其它指代）。
计算旋轨耦合矩阵元有什么用？最常见的用处就是通过微扰方式计算磷光发射速率或磷光寿命。磷光发射是三重态激发态到单重态基态的跃迁，当不考虑旋轨耦合时，由于这两个态自旋部分正交，因此计算跃迁偶极矩时候结果精确为0，导致振子强度为0，故是跃迁禁阻的。而当以微扰的方式考虑旋轨耦合的时候，则单重态里面就会掺进去一定三重态成份，而三重态里面会掺进去一定单重态成份，因此再计算跃迁偶极矩的时候，由于会引入单-单重态和三-三重态的跃迁偶极矩，因此结果就不严格为0了。至于具体如何通过本文计算的旋轨耦合矩阵元计算磷光速率，笔者之后会另文介绍。此外，计算旋轨耦合矩阵元还被用于考察系间窜越速率（正比于相应两个态之间旋轨耦合矩阵元的模方），通过二阶微扰方法还可以通过旋轨耦合矩阵元计算多重态因为旋轨耦合效应导致的子态的能级分裂（这被称为ZFS，零场分裂）。
计算旋轨耦合矩阵元的时候H_SO的算符定义并不唯一，取决于用什么方法把Dirac方程变换为二分量。把Dirac方程对于分子体系最完整的哈密顿形式，即Dirac-Coulomb-Breit哈密顿做FW变换得到的是Breit-Pauli (BP)哈密顿，多数程序里和文献里都是取BP哈密顿中的旋轨耦合算符H_SO来计算旋轨耦合矩阵元，常用的ORCA也是如此，这样处理比较简单。但也有其它情况，比如Dalton程序也支持用DK哈密顿的旋轨耦合算符来计算，而ADF则是基于ZORA来计算。
很多程序都可以计算旋轨耦合矩阵元。大多数程序支持的是在MCSCF或多参考方法级别下计算，比如ORCA、GAMESS-US、Molpro、Molcas等。然而众所周知，这类方法不黑箱，需要使用者能够恰当定义活性空间，而且计算较耗时，用于诸如较大的配合物体系比较困难。然而，具有较强磷光发射能力的体系大多却正是含有过渡金属的配合物体系，这类体系主流的研究手段是TDDFT。因此，在TDDFT下计算旋轨耦合矩阵元这一点非常重要。可惜的是，能做到这点的程序很有限。常见的就三个，一个是贼贵的ADF，一个是虽然开源免费，但是使用门槛比ADF高得多的Dalton，还有一个是ORCA，在这方面做得很好，见《使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱》（http://sobereva.com/462）。
Gaussian在旋轨耦合计算方面是个显著软肋，只支持MCSCF下计算旋轨耦合矩阵元，而且还很不好用。好在，在J. Chem. Theory Comput., 13, 515−524 (2017)一文中，作者发布了一个名为PySOC的开源免费的小工具，可以基于Gaussian的TDDFT的输出文件在TDDFT级别下计算旋轨耦合矩阵元，解决了很多人长期发愁的一大问题，使用方法将在本文后面详细介绍。
2 旋轨耦合矩阵元的一些特点
旋轨耦合矩阵元是个复数，有虚部也有实部，有些程序（如ADF）会分别输出这两部分，有的程序（如PySOC）则只给出它的模，即实部和虚部的平方和开根号。
旋轨耦合算符可写为其x、y、z三个分量加和的形式，多数量化程序（如ORCA、Dalton等）会分别计算、分别输出；而有的程序（如PySOC）并不分别计算并分别输出，而是在计算的时候就合在一起了。三个分量有的对应实部有的对应虚部，程序在给出三个分量时经常将虚部用负数表示。
我们最关心的一类旋轨耦合矩阵元是单重态与三重态间的，PySOC计算的也正是这种矩阵元。之所以三重态叫三重态，是因为它包含三个子态(sublevel)，自旋磁量子数分别为+1、0和-1。单重态与三重态的不同子态之间的旋轨耦合矩阵元并不相同。有的程序（如PySOC）会分别输出单重态与每个子态的旋轨耦合矩阵元，也有的程序只给出总大小，相当于单重态与每个子态旋轨耦合矩阵元的平方和开根号。
由于旋轨耦合矩阵元数量级比较小，因此一般用cm^-1为单位表示，与常用单位的转换关系是：1Hartree=219474.6363cm^-1、1eV=8065.5447cm^-1。
3 关于旋轨耦合积分
确定了H_SO的形式，就可以计算两个电子态之间的旋轨耦合矩阵元了。通过推导，它可以写为基函数间的旋轨耦合算符对应的积分（旋轨耦合积分）的线性组合，组合系数涉及轨道展开系数和组态系数。
H_SO其中包含单电子旋轨耦合算符H_SO1和双电子旋轨耦合算符H_SO2，前者体现电子自旋与它绕着原子核的轨道运动的耦合，后者体现电子自旋与这个电子绕着其它电子的轨道运动的耦合。前者很好算，花不了什么时间，然而对于较大体系，计算后者耗时非常高，牵扯到双电子积分的导数。然而，H_SO1和H_SO2的贡献都在同一个数量级，若直接忽略后者甚至能对旋轨耦合矩阵元带来达到1/3程度的误差。为了降低双电子旋轨耦合积分的耗时，有人提出自旋轨道平均场(SOMF)的处理方法，修改了H_SO2算符的形式，使得双电子旋轨耦合以平均场的方式近似计算（类似于Hartree-Fock的思想），这使耗时大为降低，而带来的误差可忽略不计（对于不是特别重的原子来说）。这种自旋轨道平均场的做法在实际量化程序中一般是以称作Atomic-mean field integral (AMFI)的方式具体实现的，在计算时只保留SOMF的单中心项，使得耗时更低，尽管误差也会稍大一点。
有效原子电荷(Zeff)方法也是非常常用的近似计算旋轨耦合积分的方法。它的思想是，由于实际研究发现H_SO1和H_SO2部分有比较好的比例关系，因此，可以修改H_SO1项，来把H_SO2的效果等效地吸入H_SO1里面。这样仅需要花费计算H_SO1项的耗时，就可以同时体现出H_SO2的效果了。Zeff这种做法是目前最便宜的考虑旋轨耦合的做法，比AMFI便宜得多，但是误差比较明显，可能能达到百分之十几、二三十甚至更多（和体系、元素、参数有关。对于很重的过渡金属误差可能较大）。Zeff的还一个很大好处是可以结合赝势使用来节约耗时，同时等效考虑标量相对论效应。比如计算一个Pt的配合物，如果用Zeff，那么可以给Pt用相对论赝势，这样即便不用相对论哈密顿，Au的标量相对论效应也考虑了，而且在此基础上通过很便宜的Zeff方式就能把旋轨耦合效应也考虑了，真是美哉！反之，用AMFI的时候不仅考虑旋轨耦合这部分相对更昂贵，对Pt还得用全电子基组（算一个的耗时顶得上算一把轻原子），而且雪上加霜的是还得用相对论哈密顿考虑标量相对论效应。可见，以Zeff方式计算旋轨耦合积分来得到旋轨耦合矩阵元，对于大体系是很理想的做法。虽说精度差一些，但对于大体系来说，一般也没必要那么较真。
使用Zeff的一个关键是有效电荷。把H_SO2的效果吸入H_SO1就是通过把H_SO1里面原子核电荷改为有效电荷来实现的。有效电荷是前人拟合出来的，比如可以通过计算一批体系的零场分裂能，通过令计算值和实验值尽可能相符来得到。有效电荷的好坏直接影响Zeff方法的精度，而且，仅当要算的元素有已经拟合好的有效电荷的时候才能用Zeff方法处理。有效电荷的数值一定程度上受到当初拟合它用的计算级别的影响，但影响不至于特别大，只要自己用的基组和当初在拟合时相差不是特别大就可以用，比如当初是在6-31G*下拟合的，你结合cc-pVDZ、def2-SVP、6-311G*等使用也没问题。有些有效电荷是对于赝势基组来拟合的，你用的时候可以结合其它赝势基组用，但你用的赝势和拟合时候用的赝势所赝化的内核电子数必须相同，比如不能当初是对小核赝势基组拟合的有效电荷，你却结合大核赝势基组来用。
有不同的人提出不同的有效电荷，最常用的一套是Koseki在90年代搞的一套。在JPC, 96, 10768 (1992)中对6-31G*拟合了前三周期元素有效电荷，在JPC, 99, 12764 (1995)中对SBKJC赝势拟合了一直到碘的所有主族的有效电荷，在JPCA, 102, 10430 (1998)中对SBKJC赝势拟合了所有d、ds族过渡金属的有效电荷。拟合时候都是用MCSCF做的，结合如今常用的DFT也没什么问题。SBKJC赝势如今不怎么用，它对过渡金属是小核，对主族是大核，在Gaussian里直接写SDD、lanl2所用的赝势也是这种情况，因此可以把这套有效电荷搭配常用的SDD、lanl2赝势来用。
4 PySOC的特征和原理
PySOC可以在https://github.com/jzpathfinder/pysoc下载。程序虽然开头俩字母是py，强调程序里用了Python，但实际上大部分都是Fortran写的。此程序的主要用处是结合Gaussian计算TDDFT级别下的单重态-三重态间的旋轨耦合矩阵元，基态要求是闭壳层单重态。程序目前支持G09，笔者发现似乎没法正常处理G16的输出（网友提供了可以兼容G16的方法，见http://bbs.keinsci.com/thread-19813-1-1.html）。PySOC应当也可以在Windows下用，但本文都只涉及Linux的情况。
PySOC实际上还支持DFTB+程序，DFTB+是用来做DFTB方法的计算的。DFTB从精度和速度上类似于半经验版的DFT，将DFTB和TDDFT弄在一起就是TD-DFTB，耗时比一般的TDDFT低得多。根据PySOC原文的测试，基于TD-DFTB得到的旋轨耦合矩阵元比起TDDFT的相差很多，并没什么卵用，所以本文也不去谈怎么把PySOC和DFTB+结合使用。
PySOC的运作机制是这样的：首先用户用Gaussian09做TDDFT计算得到同等数目的单重态和三重态激发态，并保留rwf文件。然后启动PySOC.py，这个脚本就会调用Gaussian自带的rwfdump工具从rwf文件中提取各种接下来用的信息，比如组态系数、分子轨道展开系数，并且从Gaussian输出文件中提取一些信息，比如激发能、基组定义。然后PySOC自动调用附带的修改版MolSOC以Zeff方式计算旋轨耦合积分，这里MolSOC是一个可以独立运行的专门用来计算旋轨耦合积分的程序。接下来PySOC就按照其原文里的公式，将旋轨耦合积分、组态系数、分子轨道展开系数组合到一起，得到包含基态在内的各个单重态与各个三重态激发态间的旋轨耦合矩阵元。
原版的PySOC和原版的MolSOC都有一些不完美的地方，于是笔者进行了修改，下文一切操作的步骤一律都是相对于这里提供的修改版来说的。笔者的PySOC+MolSOC修改版下载地址：http://sobereva.com/soft/sob_PySOC_MolSOC.zip
笔者的修改版相对于原版在功能上主要有以下改变：
1 把SOC矩阵元输出文件的态的序号输出格式改为了两位（原版只留了一个字符的位置，因此态数如果是两位数就会被显示为星号，不便考察）
2 PySOC默认会计算激发态之间的跃迁偶极矩，但是这没什么用，而且笔者发现虽然对于TDA计算得到的跃迁偶极矩是正确的，但是对于TD计算，得到的跃迁偶极矩并不合理。于是笔者删除了计算这个的功能，PySOC计算大体系的总耗时由此可降低约40%。PS：如果真需要计算激发态间的跃迁偶极矩，可以用Multiwfn，结果可靠，速度也快。见《利用Multiwfn计算Gaussian输出的激发态之间的跃迁偶极矩》（http://sobereva.com/227）。
3 修改了MolSOC里定义有效电荷的sozeff.f文件。此程序自带的有效电荷是其作者自己拟合的前五周期主族的（用B3LYP结合DZVP全电子基组下拟合，见JCC,30,832 (2009)），以及极个别其它零碎的元素。显然只有这些有效电荷根本不够用，结合赝势算配合物都不行。笔者遂把上面提到的所有Koseki拟合的有效电荷都补充进去了。
5 PySOC的安装、配置
用PySOC最好系统别太老，否则一方面系统自带的Python版本太老会导致PySOC的Python脚本无法正常运行（要求Python版本>=2.7，可以用python --version命令查看本机python版本），另一方面由于PySOC和MolSOC预编译版是在较新的系统下编译的，在老系统上运行会因为找不到较新的glibc库而出错。笔者建议用CentOS/RHEL >7.0版。如果是较老的系统，也可以升级Python，并且自己编译PySOC和MolSOC。
注：很多本文的读者反映PySOC没运行成功，很大一部分是因为用的Python版本问题不合适所致的。笔者这里用的CentOS 7.2自带的Python是2.7.5版。如果你运行不成功，尝试用这个或相近版本的Python。
首先将修改版的PySOC+MolSOC解压，比如PySOC解压后的目录为/sob/pysoc，MolSOC解压后的目录为/sob/molsoc_modified。分别进入这两个目录，都运行chmod +x * -R命令，使这俩目录和子目录下的所有文件都有可执行权限。
把pysoc/bin/pysoc.py里的scrip_soc变量改为实际的PySOC可执行文件路径，比如scrip_soc = '/sob/pysoc/bin'。对于当前情况，/sob/pysoc/bin目录里面会看到有个soc_td，这正是笔者编译出来的PySOC中由Fortran写的代码对应的可执行文件。
修改/sob/pysoc/input_template/init.py，把脚本里的g09root改为g09所在目录（和安装Gaussian时的g09root环境变量指向的路径一致），把molsoc_path改为实际的修改版molsoc的路径，比如molsoc_path = '/sob/molsoc_modified/molsoc0.1/bin/molsoc0.1.exe'。
修改用户目录下的.bashrc文件，添加export PATH=$PATH:/sob/pysoc/bin/。之后重新进入终端，就可以在任意目录下直接输入pysoc.py来调用这个运算脚本了。
进入/sob/molsoc_modified/molsoc0.1/bin目录，会看到有三个可执行文件，对应不同有效电荷的情况，你当前的计算适合哪种情况，就把哪个可执行文件改名为molsoc0.1.exe。这三个文件对应的情况为：
molsoc0.1_1.exe：保留了作者自己对前五周期主族拟合的有效电荷（用于全电子计算，对DZVP拟合），加入了Koseki对d、ds族金属结合SBKJC小核赝势拟合的有效电荷
molsoc0.1_2.exe：保留了作者自己对前三周期拟合的有效电荷，对第四、五周期主族改为了Koseki对SBKJC大核赝势拟合的有效电荷，加入了对d、ds金属koseki对小核SBKJC赝势拟合的有效电荷。
molsoc0.1_3.exe：把前三周期元素改为了Koseki对6-31G*拟合的有效电荷，对第四、五周期主族改为了Koseki对SBKJC大核赝势拟合的有效电荷，加入了对d、ds金属koseki对小核SBKJC赝势拟合的有效电荷。
一般来说，推荐用molsoc0.1_3.exe，对应于前三周期都用全电子基组，而对之后的元素用SDD、lanl2等赝势的情况，这是最常见的情形。
再次提醒，使用PySOC在用赝势的时候，赝化的电子数一定要和SBKJC对应，缺乏基础知识的话仔细看《赝势的函数形式以及在量子化学程序中定义的方式》（http://sobereva.com/188）、《谈谈赝势基组的选用》（http://sobereva.com/373）、《在赝势下做波函数分析的一些说明》（http://sobereva.com/156）。比如你把molsoc0.1_3.exe改名为了molsoc0.1.exe以令PySOC能够调用它，若你要计算个含碘的体系，那对碘就不能用def2系列或者cc-pVnZ-PP系列了，因为它们对主族是小核，而Koseki拟合有效电荷时用的SBKJC对主族是大核。
如果你想自己编译PySOC也可以，最好用ifort编译器。在/sob/pysoc/src目录下，修改Makefile里的编译器和库的目录，然后直接运行make，就会产生soc_td，将之挪到/sob/pysoc/bin目录下就可以用了。
如果要自己编译MolSOC，就进入/sob/molsoc_modified/molsoc0.1/source源代码目录，里面会看到有sozeff_1.f、sozeff_2.f、sozeff_3.f三个文件，它们是用来设定有效电荷的，分别对应于上述molsoc0.1_1.exe、molsoc0.1_2.exe、molsoc0.1_3.exe的情况，打算用哪套有效电荷组合就把哪个.f文件改名为sozeff.f。之后运行make来编译，默认会调用ifort编译，并产生可执行文件/sob/molsoc_modified/molsoc0.1/source/molsoc0.1.exe。
6 PySOC的使用
用PySOC计算旋轨耦合矩阵元之前需先自行用Gaussian09运行TD(50-50,nstates=x)任务计算单重态和三重态各x个，x设多大up to you。由于程序只能处理笛卡尔型基函数，因此必须写上6d 10f关键词。虽然程序手册里没明确强制要求，但实际发现nosymm也得写，否则PySOC貌似没法正常解析输出文件。另外还得写GFInput让Gaussian把基组定义输出出来。输出文件名字必须为gaussian.log，并且必须保留下来rwf文件，名字为gaussian.rwf。以下是输入文件开头部分的例子：
%rwf=gaussian.rwf
# td(50-50,nstates=5) wB97X/TZVP 6D 10F nosymm GFInput
注：PySOC对态数、泛函没有任何限制，但是由于MolSOC只支持角动量<=f的基函数的旋轨耦合积分的计算，因此用的基组不能包含f以上角动量基函数（由于高角动量函数对于TDDFT计算以及旋轨耦合矩阵元的影响很小，这个要求无关大碍）。
G09算完后，确保gaussian.rwf和gaussian.log都在当前目录了，将init.py从/sob/pysoc/input_template里拷到当前目录下，然后用文本编辑器打开它。Gaussian做TD(50-50)计算时nstates设了几，就把此文件里n_s =和n_t =后面的序列从1一直写到多少。比如nstates设了5，init.py里就要设n_s = [1, 2, 3, 4, 5]和n_t = [1, 2, 3, 4, 5]。有的时候我们可能要算几十个态，手动从1写到最后相当麻烦，为了方便大家使用，我把从1到100的序号附到这个文件里了，并在开头加了#来注释掉。比如算了40个态，想计算S0~S40与T1~T40之间的旋轨耦合矩阵元，那就把从1到40的序号复制到n_s和n_t后面的方框里就行了。
最后，运行pysoc.py即可，程序就会根据当前目录下init.py的信息利用gaussian.rwf和gaussian.log文件进行运算了，屏幕上输出的都是一堆烂七八糟对用户没什么用的debug信息。PySOC的计算耗时和体系大小、计算的态数都有关系，小体系一两分钟就能算完，而大体系、上百个态的情况甚至需要算半个小时乃至更多。可惜PySOC没有被并行化，没法利用多核的优势。算的态数很多的时候，建议以诸如pysoc.py > pysoc.out方式运行来把输出信息重定向到某个文件里，要不然可能由于屏幕上输出信息量太大导致拖慢计算。
算完后当前目录下会看到一大堆文件，有的是PySOC调用rwfdump导出的数据，有的是MolSOC程序的输入输出文件，有的是PySOC运行产生的一些中间信息。其中soc_out.dat是最重要的输出文件，这就是各个态之间的旋轨耦合矩阵元。这是其中一条输出：
sum_soc, <S0|Hso|T1,1,0,-1> (cm-1): 115.94270 81.98387 0.01644 81.98387
代表S0和T1之间总SOC是115.94 cm-1，它由这行后面的<S0|Hso|T1_Ms=1>、<S0|Hso|T1_Ms=0>、<S0|Hso|T1_Ms=-1>三个值求平方和开根号得到，这里Ms指T1子态的自旋磁量子数。另外当前目录下还会产生ene_out.dat，这里面是各个态的激发能。
下面我们通过两个实例演示怎么使用PySOC。本文所有Gaussian计算用的都是G09 E.01 Linux版，系统是CentOS 7.3。计算前读者应当根据本文第5节的方式把PySOC运行环境进行了恰当配置、恰当修改了相关的.py文件。下面例子用的有效电荷对应的是上文提到的molsoc0.1_3.exe的情况（即用户已经把此文件改名为了molsoc0.1.exe）。更老的G09版本笔者未经测试，建议不要用老于D.01版的，尤其是早期G09程序的TD(50-50)选项有bug，PySOC肯定运行不正常。
7 实例1：计算H2CO的单-三重态间的旋轨耦合矩阵元
本例我们考察一个很简单的体系，甲醛。在PySOC原文的表3中，给出了这个体系的计算结果，如下所示
(, 下载次数 Times of downloads: 228)
这里我们要重现红框中的值。通过查看甲醛的分子轨道和电子跃迁信息，可以知道红框中的值对应的是<S1|H_SO|T2>。从表中的数据也可以看到，在B3LYP下，PySOC基于Zeff的结果（45cm-1），和更严格的在Dalton中通过直接计算单、双电子旋轨耦合积分得到的结果（54cm-1）还是比较相符的。虽然误差不能说非常小，但从实用性角度来说，准确度也够了。（注：PySOC原文11式是错的，文中给出的旋轨耦合矩阵元并没有乘以系数1/3）
首先我们在PBE0/TZVP级别下优化甲醛分子，输入文件如下。
# PBE1PBE/TZVP opt
test
0 1
C 0.00000000 0.00000000 -0.56221066
H 0.00000000 -0.92444767 -1.10110537
H -0.00000000 0.92444767 -1.10110537
O 0.00000000 0.00000000 0.69618930
然后把优化好的结构保存成新的Gaussian输入文件，内容如下，所用计算级别和原文表3一致。虽然我们只需要考察S1和T2，但正如《Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法》（http://sobereva.com/314）中所提到的，算的态数不应当正好卡着要研究的态。如果感兴趣的是第i态，那么计算到i+3态一般是比较稳妥的，因此本例我们总共计算5个单重态和5个三重态。
%rwf=gaussian.rwf
# td(50-50,nstates=5) B3LYP/TZVP 6D 10F nosymm GFInput
test
0 1
C 0.00000000 0.00000000 -0.52513500
H 0.00000000 0.93987900 -1.11261300
H 0.00000000 -0.93987900 -1.11261300
O 0.00000000 0.00000000 0.67200400
运行g09 < H2CO.gjf > gaussian.log，算完后当前目录下就有了gaussian.log和gaussian.rwf。然后把/sob/pysoc/input_template/init.py文件拷到当前目录，用文本编辑器打开它，确认n_s和n_t都被设为了[1, 2, 3, 4, 5]。然后，直接输入pysoc.py命令来运行之，由于体系很小而且考虑的态数也很少，眨眼间就算完了。屏幕上最后会提示有报错“forrtl: severe (153): allocatable array or pointer is not allocated”，这没有关系，这是因为笔者的修改版把PySOC计算跃迁偶极矩的部分给去掉了所致，这不影响旋轨耦合矩阵元的输出。
打开当前目录下产生的soc_out.dat，从中会看到这行
sum_soc, <S 1|Hso|T 2,1,0,-1> (cm-1): 49.81285 0.00000 49.81285 0.00000
即其中第一个数字49.81285 cm-1就是总的<S1|H_SO|T2>矩阵元。这个值和文献中的45 cm-1略有偏差，原因是原文计算时候并没有修改MolSOC里默认的有效电荷，即用的是MolSOC作者自己搞的有效电荷，而我们现在用的则是Koseki拟合的有效电荷。可见有效电荷的不同，对于结果是会产生一定影响的，但影响不至于很大。倘若我们运行pysoc.py之前是把molsoc0.1_1.exe改名为了molsoc0.1.exe，即计算H,C,O元素的时候用MolSOC作者自己搞的有效电荷，结果将为44.73483 cm-1，和PySOC原文里的精确一致。
8 实例2：计算Os配合物的单-三重态间的旋轨耦合矩阵元
有较强发磷光能力的分子一般都是含有靠后周期的d金属的配合物，因为这样的金属会引入显著的旋轨耦合效应，这类配合物体系是需要考察单-三重态旋轨耦合的主要场合之一。Gaussian+PySOC以Zeff方式计算这类体系是否靠谱？我们这里重现一下Phys. Chem. Chem. Phys., 16, 26184 (2014)文章里的旋轨耦合数据。此文研究了一系列含有Os、Au、Cu、Ag的配合物，通过Gaussian在B3LYP/6-31G*结合lanl2DZ做了几何优化，然后使用ADF在B3LYP结合TZP全电子基组并考虑ZORA标量相对论效应下算了最低10个单重态和三重态，之后计算了激发态间的旋轨耦合矩阵元。计算时候都以PCM模型考虑了二氯甲烷的溶剂效应。此文考察了不少体系，为了节约时间，我们就考察其中比较小的Os-4，结构如下
(, 下载次数 Times of downloads: 206)
本节例子涉及的文件都可以在此下载：
(, 下载次数 Times of downloads: 243)
首先我们用gview等程序构建Os-4结构，保存输入文件，优化此体系同时做振动分析确认无虚频。我们用的计算级别和原文一致，输入文件如下（其实Os用SDD会更好，参见《谈谈赝势基组的选用》http://sobereva.com/373）
# B3LYP/genecp opt freq scrf(solvent=CH2Cl2)
[...略]
C O F H
6-31G*
****
Os
lanl2DZ
****
Os
lanl2DZ
将优化后的结构保存成新的输入文件，做电子激发计算。我们用的泛函和原文一样也是B3LYP，也是算最低10个单重态和三重态。原文用的是ADF，这程序很非主流，其基组都是基于Slater函数的，因此文中用的TZP基组在Gaussian这样的基于高斯函数的程序中没有严格对应的。我们这里就用基于高斯函数的3-zeta基组中质量上乘的def2-TZVP来做计算（这对于当前计算来说其实有点浪费了），整体和TZP能达到同一个档次。由于def2-TZVP对Os是小核赝势基组，笔者的修改版MolSOC里又加入了Koseki对所有d族对小核赝势拟合的有效电荷，因此此时正好可以用PySOC来计算（如果在Gaussian里用更昂贵、更准确的DKH2标量相对论方法结合全电子基组来做此体系的TDDFT计算，PySOC则没法用，因为没有对MolSOC添加对过渡金属全电子基组拟合的有效电荷）。
Os-4的电子激发计算输入文件如下（Os-4_TD.gjf）
%rwf=gaussian.rwf
# b3lyp/def2TZVP scrf(solvent=CH2Cl2) TD(nstates=10,50-50) 6d 10f nosymm gfinput
[...略]
然后运行g09 < Os-4_TD.gjf > gaussian.log。由于用的基组较大，因此耗时还是不短的。算完后把init.py拷到当前目录，把其中的n_s和n_t后面都改为[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]，然后运行pysoc.py > pysoc.out，过一会儿就算完了。
文献中在表1中给出了各个Os配合物发生主要系间窜越的两个态之间的旋轨耦合矩阵元，对Os-4给出的是<S1|H_SO|T2>，数值为620 cm-1，在S0结构下它们的能级差（即垂直S1-T2能级差ΔE_ST）为0.11eV。
(, 下载次数 Times of downloads: 195)
我们看看我们算的。打开ene_out.dat，会看到S1和T2能量分别为3.9079eV和3.9945eV，能级差为0.09eV，和文献给出的值之间的绝对差异很小。然后打开soc_out.dat，会看到
sum_soc, <S 1|Hso|T 2,1,0,-1> (cm-1): 593.44440 419.62830 0.69267 419.62830
即我们算出来的<S1|H_SO|T2>为593cm-1，和文献给出的620cm-1之间的相对误差才不到5%。虽然相符这么好不免有些巧合，但至少也足够说明我们用的赝势+Zeff的廉价考虑标量相对论效应和旋轨耦合效应的组合对于配合物是可靠的。
笔者在更廉价的基组上也做了计算测试，配体用6-311G*，Os用SDD，这个档次对此体系耗时仅为def2-TZVP时的1/4，算出来的<S1|H_SO|T2>为539.9cm-1，ΔE_ST=0.093eV。可见单-三重态能级差用6-311G*结合SDD已经足够算准，旋轨耦合矩阵元对基组更敏感点，但当前级别至少也足够给出定性合理的结果。
9 总结
Gaussian+PySOC这种组合可以很容易地计算TDDFT级别下的旋轨耦合矩阵元，虽然PySOC计算本身也花时间，但耗时比TDDFT过程低一个数量级，即曰只要TDDFT算得动，单-三重态旋轨耦合矩阵元就可以轻易获得。虽然PySOC的结果基于比较糙的Zeff方法，但精度对于考察包括过渡金属配合物在内的较大体系来说够用了。
如果你不是非得用Gaussian不可，那么相对于本文的做法，用ORCA是好得多的选择，更快、更方便还更准确，看《使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱》（http://sobereva.com/462）。
免费的Dalton也可以算TDDFT下的单-三重态旋轨耦合矩阵元，但门槛稍高，而且算激发态间的旋轨耦合矩阵元耗时远高于用Gaussian+PySOC或ORCA。
作者Author: chiweijie 时间: 2018-3-15 08:14
好文
作者Author: 让你变成回忆 时间: 2018-3-15 08:37
我使用pysoc计算过一些体系的soc，但感觉普适性还有待考证。
如果pysoc在计算准确度上有保证的话，估计今后像ADF就会又少很大一批用户了。
作者Author: alwens 时间: 2018-3-15 09:40
千呼万唤始出来。
作者Author: superrice 时间: 2018-3-15 09:45
这个好，顶一下！
作者Author: 让你变成回忆 时间: 2018-3-15 14:46
sob老师您好，在PySOC的安装、配置部分的第三段
把pysoc/bin/soc.py里的scrip_soc变量改为实际的PySOC可执行文件路径
应该是pysoc/bin/pysoc.py里的scrip_soc，在soc.py中并没有找到scrip_soc
是否是笔误？
作者Author: wuy069 时间: 2018-3-15 15:34
如果能给出辐射速率常数，那就完全没必要用ADF程序了
作者Author: sobereva 时间: 2018-3-16 02:06
做这个事情用Dalton非常理想，完全不需要买ADF
作者Author: sobereva 时间: 2018-3-16 02:17
是笔误，已改
作者Author: abin 时间: 2018-3-17 12:58
本帖最后由 abin 于 2018-3-20 16:16 编辑
更正，之前测试，
放做了文件。
G09, G16, 对于同一个xyz文件， 采用同样的输入控制， 然后用pysoc处理结果，完全一致。
请放心使用。
特此更正。
如给你带来误解，敬请谅解。
===============
可以处理G16输出结果，
需要求改几处配置。
不过如同楼上提到的一样，
计算结果的普适性有待考证。
而且，我测试发现， 同一个结构， 分别使用G09， G16来处理， 最后得到的SOC数据偏差很大， 前者可能是0， 后者可能是好几百，
简单说， 前者可能认为， 比如T1-S0 SOC是0， 二后者可能会告诉你T1-S0 SOC是好几百。
G09 G16 用pysoc处理后的结果，
和ADF给出的结果， 也不同。
也就是说， 计算SOC， 三者的结果相互不一致，就是定性也不一致。
另，算速率常数啊， 看介绍资料MOMAP不是可以处理吗？
至于计算结果，用了才知道。
还有MOMAP的手册啊， 高明的厉害啊。
作者Author: sobereva 时间: 2018-3-17 13:21
以G09为准，程序也是针对G09进行的测试。自行修改PySOC代码肯定是可以兼容G16的。
如果非要用G16，结果相差很大的情况肯定是自己修改的PySOC没弄对，可以在源代码中加入一些调试语句试图弄清楚原因，特别是对比PySOC读进去的组态和轨道系数。只要G16和G09用的int=acc2e、积分格点相同，TDDFT本身给出的信息是精确相同的（起码对于真空下的计算来说），因此也注定不会导致PySOC运行正常时给出的结果有明显差异。
PySOC的普适性是没问题的，原文对有机分子有充分测试。我测试过的各种配合物体系结果也都合理。如果出现不合理、不准确则是Zeff方法的事，是Zeff的普适性问题，而不是PySOC的普适性问题，PySOC的实现是严格的。如果嫌Zeff精度还不够，或者想弄清楚和其它程序给出的SOC矩阵元的差异原因，用Dalton便可，Zeff、AMFI以及完整计算旋轨耦合积分都支持。
作者Author: abin 时间: 2018-3-17 14:31
我没仔细看。
看到Pysoc中调用Gauss， 貌似只是读取rwf。如果G09 G16 rwf读取之后的格式和内容没有明显变化， 结果应该不会差太多的。
或者也可以偷懒处理。 如果需要激发态的构型， 而且需要TD优化， 可以用G16优化， 然后读取构型， 用G09再搞一次单点的TD激发能计算， G16 TD+opt+freq毕竟快很多的。 个人愚见， 仅供参考。
作者Author: ketone 时间: 2018-3-20 21:13
好文！给作者个赞！
作者Author: yjr 时间: 2018-3-21 12:36
仿照计算第一个例子时，pysoc.py后出现报错，如下
File "/home/yangjr/sob_PySOC_MolSOC/pysoc/bin/soc.py", line 359
trans_dict = {trans_key[k]: ci_coeff[np+k] for k in range(ndim)}
^
SyntaxError: invalid syntax
forrtl: No such file or directory
请问什么原因呢？
作者Author: sobereva 时间: 2018-3-21 14:27
8成是系统的python版本太老
作者Author: yjr 时间: 2018-3-21 15:39
好的，谢谢您的及时回复，我去更新一下
作者Author: 哆啦哆啦 时间: 2018-3-21 16:58
运行第一个例子，出错信息如下：
/home/liu/software/molsoc_modified/molsoc0.1/bin/molsoc0.1.exe: /lib64/libc.so.6: version `GLIBC_2.14' not found (required by /home/liu/software/molsoc_modified/molsoc0.1/bin/molsoc0.1.exe)
Traceback (most recent call last):
File "/home/liu/software/pysoc/bin/soc.py", line 513, in <module>
hso[i] = read_file('soint', cont0, nt, 0)
File "/home/liu/software/pysoc/bin/soc.py", line 40, in read_file
with open(filein, 'r') as f:
IOError: [Errno 2] No such file or directory: 'soint'
/home/liu/software/pysoc/bin/soc_td: /lib64/libc.so.6: version `GLIBC_2.14' not found (required by /home/liu/software/pysoc/bin/soc_td)
请问什么原因？
作者Author: 让你变成回忆 时间: 2018-3-21 17:02
需要安装GLIBC库。我不清楚这个是否是python自带的。我使用的时候是单独安装了这个库。
作者Author: sobereva 时间: 2018-3-21 17:41
系统太老，GLIBC库版本太低
作者Author: YuhangYao 时间: 2018-5-7 21:58
本帖最后由 YuhangYao 于 2018-5-7 22:07 编辑
计算第一个例子没有问题，但是做我自己体系时出现报错：
Traceback (most recent call last):
File "/gpfs/share/home/1501110340/software/PySOC/pysoc/bin/soc.py", line 522, in <module>
s_matr = read_file('molsoc_overlap.dat', 'AO_overlap', nt, 0)
File "/gpfs/share/home/1501110340/software/PySOC/pysoc/bin/soc.py", line 40, in read_file
with open(filein, 'r') as f:
IOError: [Errno 2] No such file or directory: 'molsoc_overlap.dat'
forrtl: No such file or directory
forrtl: severe (29): file not found, unit 171, file /gpfs/share/home/1501110340/LuL1/PySOC/soc_td_input.dat
Image PC Routine Line Source
soc_td 000000000041BBB6 Unknown Unknown Unknown
soc_td 0000000000427AEB Unknown Unknown Unknown
soc_td 0000000000407F0B Unknown Unknown Unknown
soc_td 000000000040254E Unknown Unknown Unknown
soc_td 00000000004024FE Unknown Unknown Unknown
libc.so.6 00002AE445625B35 Unknown Unknown Unknown
soc_td 0000000000402409 Unknown Unknown Unknown
我看了一下molsoc.out的输出：SOZEFF: SOZEFF IS NOT AVAILABLE FOR ATOMIC NUMBER 71。我的体系中含有一个金属Lu，所以是Zeff方法计算SOCME不适用于稀土金属吗？
作者Author: sobereva 时间: 2018-5-8 03:50
只要有有效电荷参数，所有元素原则上都可以算。你得自己补充Zeff参数。
作者Author: YuhangYao 时间: 2018-5-8 21:08
谢谢！补充个Zeff参数后果然可以了
作者Author: zhangfuqiang 时间: 2018-5-24 16:37
这两天测试了一下，有个小问题，也不知做的对不对： soc.py不读Landl2dz赝势，在其111行填入ECP, 结果如下：['AO basis set','primitive gaussians','ECP']) ，发现全部都进去了。运行正常。
作者Author: abin 时间: 2018-6-8 13:46
本帖最后由 abin 于 2018-6-13 15:12 编辑
遇到一个奇怪的错误。
原始的算例：
%chk=ch2o0-g16
%rwf=ch2o0-g16
# td(50-50,nstates=5) B3LYP 6-31G(d) 6D 10F nosymm GFInput
Run with G16
0 1
C -0.131829 -0.000001 -0.000286
O 1.065288 0.000001 0.000090
H -0.718439 0.939705 0.000097
H -0.718441 -0.939705 0.000136
得到的结果可以使用pysoc处理。
使用类似的计算，测试了一个pt配合物：
%rwf=pt-demo-pysoc
# td(50-50,nstates=5) B3LYP GenECP 6D 10F nosymm GFInput
Run with G16
0 1XYZ
Pt 0
LanL2DZ
****
O N C H 0
6-31G(d)
****
Pt 0
LanL2DZ
结果，pysoc无法处理， 报错信息为：
QM_code gauss_tddft
open file gaussian.log
[[1, 2.2291], [3, 2.5906], [4, 2.6787], [5, 2.7284], [7, 2.892]]
2 2.4176
['ANG', 'Zeff', 'DIP']
Error when reading gaussian output
奇怪的是： 如果保留pysoc处理ch2o0-g16的时候，产生的一堆输出文件，比如以下文件，
molsoc.out *coeffs molsoc_* mo_* *overlap* ao_basis.dat dip_ao.dat ci_coeff.dat ene_out.dat init.pyc molsoc.inp MO_energy s_matr.dat soc_ao.dat soc_td_input.dat trans_dip_out.dat soint
然后，在同样的目录，再次处理pt的结果，居然就么有报错信息，会得到结果。
但是该结果，至少和Dalton, ADF的结果，相差甚大。
说明：
1 也测算过其他分子体系，没有问题；
2 以上测试，用原始的pysoc版本， 或者Sob提供版本，都是同样的问题。
奇怪的就是， 当前这个pt体系无法计算。
这个错误信息：
[[1, 2.2291], [3, 2.5906], [4, 2.6787], [5, 2.7284], [7, 2.892]]
2 2.4176
['ANG', 'Zeff', 'DIP']
Error when reading gaussian output
看上去，貌似没有提取到Singlet的能量信息。
有遇到类似错误信息的同仁吗？ 有什么处理方案吗？
更新： 使用Sob博文中的算例， Os配合物。
无论是我自己修改的pysoc版本，还是Sob编译后的pysoc版本，
都可以正常处理Os的计算结果。
但是换到我测试用的pt配合物，就是报错了。
作者Author: guoy14iccas 时间: 2018-6-12 09:04
本帖最后由 guoy14iccas 于 2018-6-12 09:06 编辑
AO_overlap 514R
/bin/sh: line 0: export: `gaussian.rwf': not a valid identifier
/bin/sh: line 0: export: `514R': not a valid identifier
XY_coeffs 635R
/bin/sh: line 0: export: `gaussian.rwf': not a valid identifier
/bin/sh: line 0: export: `635R': not a valid identifier
MOB_coeffs 526R
/bin/sh: line 0: export: `gaussian.rwf': not a valid identifier
/bin/sh: line 0: export: `526R': not a valid identifier
MO_energy 522R
/bin/sh: line 0: export: `gaussian.rwf': not a valid identifier
/bin/sh: line 0: export: `522R': not a valid identifier
Error when reading gaussian output
Traceback (most recent call last):
File "/home/yuanguo/program/pysoc-1.0/pysoc/bin/soc.py", line 491, in <module>
write_file(mo_ene, 'mo_ene.dat', '{} ')
NameError: name 'mo_ene' is not defined
forrtl: No such file or directory
forrtl: severe (29): file not found, unit 171, file /home/Co_Works/HH/opt-svp/SOC/Te/FPDI-S/soc_td_input.dat
Image PC Routine Line Source
soc_td 000000000048B27E Unknown Unknown Unknown
soc_td 0000000000489D16 Unknown Unknown Unknown
soc_td 0000000000449EF2 Unknown Unknown Unknown
soc_td 000000000041ECCB Unknown Unknown Unknown
soc_td 000000000041E232 Unknown Unknown Unknown
soc_td 00000000004287F8 Unknown Unknown Unknown
soc_td 000000000040755A Unknown Unknown Unknown
soc_td 00000000004022D9 Unknown Unknown Unknown
soc_td 000000000040227C Unknown Unknown Unknown
libc.so.6 00002B9FD00F14B0 Unknown Unknown Unknown
soc_td 0000000000402179 Unknown Unknown Unknown
我用的python 3版本还是报错，请老师帮忙看看为什么？
作者Author: sobereva 时间: 2018-6-12 10:15
用CentOS 7.x再试（或Python 2.7.x）
太新的可能也不行
作者Author: 哇哇吐 时间: 2018-6-29 10:02
您好，我也遇到了如此错误，请问您解决了么，我计算的是一个含碘的体系，一开始怀疑是不是基组的错误，改了之后还是不行。
QM_code gauss_tddft
open file gaussian.log
[[1, 2.8384], [2, 3.5714], [3, 3.6104], [4, 3.6178], [5, 3.6377], [6, 3.6631], [7, 3.6658], [8, 3.8167]]
9 4.0877
['ANG', 'Zeff', 'DIP']
Error when reading gaussian output
作者Author: abin 时间: 2018-7-1 21:56
没有。
我放弃使用这个代码了。
因为看不懂。
作者Author: sobereva 时间: 2018-7-2 03:45
含Pt体系没问题
至少这个配合物输入文件按照我文中的做法能正常处理（用与文中不同的系统、其它Gaussian版本另谈）
(, 下载次数 Times of downloads: 18)
作者Author: sobereva 时间: 2018-7-2 03:49
老有人问怎么PySOC没法运行，很大程度上要么没严格按博文的步骤走、没仔细理解博文，要么用的系统和文中的不同（会导致自带的python版本和PySOC兼容性有问题），或者Gaussian用的版本和文中的不同。想排除系统差异的问题，装个vmware虚拟机，用我文中的系统就完了。
作者Author: abin 时间: 2018-7-4 08:50
可能我表达有误。
我的意思是说， “我遇到有一个Pt的体系，pysoc计算搞错”， 但是同样是含有pt, 另一个类似的结构没有任何问题。
另，pysoc与高斯版本没啥必要的关系，只要能从rwf中提取数据，就可以计算。 当然了，前提是必须有“有效电荷”定义。
作者Author: abin 时间: 2018-7-10 09:52
测试了一下。
你提到的算例 1.gjf可以处理；同样，把坐标文件替换为其他的pt complex的坐标，也能计算；
但是类似的计算控制，这里提到的pt-demo-pysoc.com， 依旧是无法处理。
报错信息大概如下：
QM_code gauss_tddft
open file gaussian.log
[[1, 2.2291], [3, 2.5906], [4, 2.6787], [5, 2.7284], [7, 2.892]]
2 2.4176
['ANG', 'Zeff', 'DIP']
Error when reading gaussian output
Traceback (most recent call last):
File "/home/wang/composition/pysoc-master/pysoc/bin/soc-g16.py", line 493, in <module>
write_file(mo_ene, 'mo_ene.dat', '{} ')
NameError: name 'mo_ene' is not defined
forrtl: No such file or directory
作者Author: tiandikuoyuan 时间: 2018-7-23 19:34
Error when reading gaussian output
Traceback (most recent call last):
File "/home/pysoc/bin/soc.py", line 491, in <module>
write_file(mo_ene, 'mo_ene.dat', '{} ')
NameError: name 'mo_ene' is not defined
forrtl: No such file or directory
同样的错误，请问现在能解决吗？rwf文件太大会不会影响计算过程？
作者Author: 让你变成回忆 时间: 2018-8-18 14:08
本帖最后由 让你变成回忆 于 2018-8-18 14:09 编辑
sob老师您好，请教您一个问题：
我想把PySOC改来支持G16，现在正在做一些各种g09(version E0.1)和g16(version A03)的所产生的rwf的对比工作，例如PySOC中用到的各种参数等。
为了使g09和g16的输出结果完全一致，是否仅仅需要在g09中加上int=acc2e=12, int=ultrafine，这两个看您写的解决SCF不收敛的帖子里提到过，我不是很清楚是否还需要其他的一些设置。
谢谢sob老师~
作者Author: sobereva 时间: 2018-8-19 12:43
对于真空下算单点，g09加上这俩就和G16结果一样了
作者Author: 让你变成回忆 时间: 2018-8-19 12:53
行，我去测试一下。 谢谢sob老师~
作者Author: feishuangyu 时间: 2018-9-12 16:24
请问你怎么解决的，我的也是这个报错
作者Author: guoy14iccas 时间: 2018-9-13 08:12
放弃使用这个，换软件算了
作者Author: jiajihui 时间: 2018-12-28 17:12
老师您还，麻烦问下windows版下该如何计算，按要求安装后，在windows下无法运行，出现了，文件名，目录名的错误提示
作者Author: sobereva 时间: 2018-12-29 01:53
非要用Windows，用vmware装个Linux虚拟机跑。用win10的WSL或许也可以，我没试过
作者Author: pyscf 时间: 2019-1-8 22:01
老师：请问假如基态不是singlet 而是doublet的话 还能用tddft+pysoc吗？谢谢
作者Author: sobereva 时间: 2019-1-9 00:06
pysoc只支持闭壳层
作者Author: pyscf 时间: 2019-1-9 06:17
好的，谢谢老师
作者Author: xinren 时间: 2019-1-14 22:17
想要算单重态和三重态的旋轨耦合，看到guassian用CASSCF的方法，或者是结合其他的软件计算，想要问问有没有可能只用DFT和TDDFT计算即可
作者Author: sobereva 时间: 2019-1-15 02:30
这帖子里说的就是这个问题啊，仔细看了自然就明白了
作者Author: archilles 时间: 2019-5-7 09:24
我遇到了类似的错误，无法读取gaussian output。
我的问题是g09root路径选错了，多加了一个'/09'，把多余的字符删掉之后就可以运行了。
作者Author: hanmingxi 时间: 2019-10-28 09:32
sob老师，请问这里面算的ene_out.dat文件中各个态的激发能可以直接用吗？还是用“Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法”这篇博文里的方法更好？
作者Author: sobereva 时间: 2019-10-29 07:49
都一样，本来就是从Gaussian输出的信息里直接提取的
作者Author: hanmingxi 时间: 2019-10-29 10:02
但我比较了下两种方法算出的结果差别还是挺大的，查了下之前计算的关键词也没有错，是有什么问题吗我没有注意到，还是两种方法本来就不一样，这样的话用哪种比较好？
作者Author: gray 时间: 2019-10-29 12:43
感谢
作者Author: sobereva 时间: 2019-10-30 05:43
方法显然一样，pysoc自己根本不重新计算能量
注意是否读错、单位弄错没有
作者Author: hanmingxi 时间: 2019-10-30 09:17
单重态：
%chk=c3-B-single.chk
%nprocshared=16
%mem=8GB
# opt b3lyp/def2svp td(nstates=15) nosymm
c3-B-single
0 1
N 25.89646400 4.66987100 -0.86072900
O 30.48938300 1.81943200 6.96125700
N 30.68140600 -0.45920700 0.99261100
N 32.14381500 5.15659900 13.41839200
O 26.41582000 7.09946000 5.01982500
C 29.71621400 3.83233300 6.00395100
C 26.19829900 3.43423100 -1.45272600
O 30.41035800 8.07044900 8.00496100
C 28.73963400 4.35703000 5.14402300
H 28.30790700 3.72139400 4.36997300
C 25.70885300 3.43000800 -2.78589300
C 28.26623800 5.66957200 5.31425700
C 26.21762000 5.05272900 0.46110700
C 28.79506000 6.45459700 6.34763700
H 28.43649300 7.47642700 6.47979400
C 30.18226300 2.40183600 5.93505800
C 26.85835500 5.78806000 3.11067900
C 30.21141000 4.62791900 7.04968800
H 30.96198900 4.19414000 7.71158600
C 27.85069400 5.30738000 2.23923000
H 28.88888200 5.24714200 2.56973500
C 30.76015400 4.78174300 11.43753500
H 30.35627200 3.87882100 11.89809900
C 30.55082300 2.36571500 3.40880400
H 30.68115700 3.44955800 3.40184300
C 29.76486400 5.94409700 7.22322400
C 30.26655500 1.69891700 4.61240600
C 31.69432500 5.55956200 12.14030900
C 25.54481500 5.91950000 2.62361800
H 24.77951300 6.31332100 3.29516200
C 30.33151800 5.17167200 10.16914800
H 29.58862800 4.55944600 9.65537800
C 25.21786500 5.45719800 -1.80211400
C 27.53883800 4.95434600 0.92669000
H 28.32279600 4.61713000 0.24664900
C 30.27983300 -0.41521400 3.40738400
H 30.19583400 -1.50344800 3.40495700
C 30.69338200 1.65673200 2.21626200
H 30.90584000 2.18330700 1.28429800
C 25.08659400 4.72098600 -3.00991000
C 30.55043800 0.25986000 2.20341400
C 30.14970300 0.29759900 4.59504600
H 29.96940600 -0.22220900 5.53805100
C 30.80689000 6.35791600 9.58278800
C 25.22069100 5.54103700 1.32542000
H 24.19000900 5.60813600 0.97306200
C 27.12952600 6.25245100 4.51115600
C 26.97226800 1.18090000 -1.71253300
H 27.45467000 0.29019700 -1.30284800
C 31.71922200 7.14702100 10.30773300
H 32.07322300 8.07593000 9.85660600
C 31.50495900 -1.22353700 -1.01119200
C 31.76252500 -0.38121600 0.10306900
C 25.87337500 2.28665500 -3.58106400
H 25.50227600 2.27159800 -4.60894500
C 26.82221400 2.30760100 -0.90188800
H 27.17526300 2.29911400 0.13028200
C 24.75067800 6.77326200 -1.69604400
H 24.87316600 7.34978100 -0.77812600
C 32.66677200 3.89706700 13.74584000
C 29.73152200 -1.34024100 0.45805800
C 32.17316100 6.74868000 11.55977100
H 32.91177700 7.34821300 12.09444600
C 30.20745900 -1.83188000 -0.78715800
C 26.50871300 1.16933300 -3.04126600
H 26.64408200 0.27201200 -3.64940900
C 32.96249800 0.33223600 0.21603800
H 33.17261400 0.96044500 1.08281800
C 30.33569000 6.88070800 8.25951400
C 33.00786300 3.88343900 15.12471800
C 24.12710500 7.33597400 -2.81102100
H 23.75378000 8.36132300 -2.74947900
C 24.45603900 5.30951900 -4.11538100
H 24.34858700 4.75343100 -5.05009100
C 32.14374300 5.94929400 14.57537300
C 33.89416800 0.20798900 -0.81587700
H 34.83702900 0.75664300 -0.74914900
C 32.45771800 -1.33099400 -2.03437900
H 32.27167800 -1.97550700 -2.89708200
C 31.67820800 7.25697000 14.76346800
H 31.24810900 7.83278400 13.94301400
C 32.75437600 5.76472600 16.93111000
H 33.15884800 5.19115500 17.76884200
C 32.67136600 5.19090700 15.65441200
C 28.47351900 -1.70640400 0.95384400
H 28.09710100 -1.31393300 1.89947100
C 32.90368900 2.78625300 12.92601700
H 32.66537600 2.80129700 11.86173500
C 31.77109600 7.80459500 16.04404100
H 31.41485400 8.82403600 16.21245300
C 32.30795600 7.07199100 17.11876200
H 32.36777700 7.53118800 18.10827500
C 23.97442100 6.61358700 -4.00873600
H 23.48010000 7.08282900 -4.86257500
C 33.64667200 -0.61018900 -1.93365300
H 34.39666900 -0.68400600 -2.72452800
C 28.16931400 -3.09843400 -1.03320300
H 27.54405500 -3.79082600 -1.60156400
C 27.70485400 -2.59061300 0.19432200
H 26.72078100 -2.89050800 0.56305700
C 29.41622300 -2.71991500 -1.52981400
H 29.77131900 -3.10613500 -2.48846300
C 33.46622600 1.65110500 13.51198900
H 33.65711900 0.77258300 12.89049800
C 33.57180600 2.73032000 15.68892600
H 33.83770600 2.70770600 16.74873900
C 33.79394400 1.61687800 14.87993500
H 34.23135000 0.71168200 15.30759900
SOC：
%rwf=gaussian.rwf
%nprocshared=16
%mem=8GB
# td(50-50,nstates=15) b3lyp/def2svp 6D 10F nosymm GFInput
gaussian
0 1
N 25.89646400 4.66987100 -0.86072900
O 30.48938300 1.81943200 6.96125700
N 30.68140600 -0.45920700 0.99261100
N 32.14381500 5.15659900 13.41839200
O 26.41582000 7.09946000 5.01982500
C 29.71621400 3.83233300 6.00395100
C 26.19829900 3.43423100 -1.45272600
O 30.41035800 8.07044900 8.00496100
C 28.73963400 4.35703000 5.14402300
H 28.30790700 3.72139400 4.36997300
C 25.70885300 3.43000800 -2.78589300
C 28.26623800 5.66957200 5.31425700
C 26.21762000 5.05272900 0.46110700
C 28.79506000 6.45459700 6.34763700
H 28.43649300 7.47642700 6.47979400
C 30.18226300 2.40183600 5.93505800
C 26.85835500 5.78806000 3.11067900
C 30.21141000 4.62791900 7.04968800
H 30.96198900 4.19414000 7.71158600
C 27.85069400 5.30738000 2.23923000
H 28.88888200 5.24714200 2.56973500
C 30.76015400 4.78174300 11.43753500
H 30.35627200 3.87882100 11.89809900
C 30.55082300 2.36571500 3.40880400
H 30.68115700 3.44955800 3.40184300
C 29.76486400 5.94409700 7.22322400
C 30.26655500 1.69891700 4.61240600
C 31.69432500 5.55956200 12.14030900
C 25.54481500 5.91950000 2.62361800
H 24.77951300 6.31332100 3.29516200
C 30.33151800 5.17167200 10.16914800
H 29.58862800 4.55944600 9.65537800
C 25.21786500 5.45719800 -1.80211400
C 27.53883800 4.95434600 0.92669000
H 28.32279600 4.61713000 0.24664900
C 30.27983300 -0.41521400 3.40738400
H 30.19583400 -1.50344800 3.40495700
C 30.69338200 1.65673200 2.21626200
H 30.90584000 2.18330700 1.28429800
C 25.08659400 4.72098600 -3.00991000
C 30.55043800 0.25986000 2.20341400
C 30.14970300 0.29759900 4.59504600
H 29.96940600 -0.22220900 5.53805100
C 30.80689000 6.35791600 9.58278800
C 25.22069100 5.54103700 1.32542000
H 24.19000900 5.60813600 0.97306200
C 27.12952600 6.25245100 4.51115600
C 26.97226800 1.18090000 -1.71253300
H 27.45467000 0.29019700 -1.30284800
C 31.71922200 7.14702100 10.30773300
H 32.07322300 8.07593000 9.85660600
C 31.50495900 -1.22353700 -1.01119200
C 31.76252500 -0.38121600 0.10306900
C 25.87337500 2.28665500 -3.58106400
H 25.50227600 2.27159800 -4.60894500
C 26.82221400 2.30760100 -0.90188800
H 27.17526300 2.29911400 0.13028200
C 24.75067800 6.77326200 -1.69604400
H 24.87316600 7.34978100 -0.77812600
C 32.66677200 3.89706700 13.74584000
C 29.73152200 -1.34024100 0.45805800
C 32.17316100 6.74868000 11.55977100
H 32.91177700 7.34821300 12.09444600
C 30.20745900 -1.83188000 -0.78715800
C 26.50871300 1.16933300 -3.04126600
H 26.64408200 0.27201200 -3.64940900
C 32.96249800 0.33223600 0.21603800
H 33.17261400 0.96044500 1.08281800
C 30.33569000 6.88070800 8.25951400
C 33.00786300 3.88343900 15.12471800
C 24.12710500 7.33597400 -2.81102100
H 23.75378000 8.36132300 -2.74947900
C 24.45603900 5.30951900 -4.11538100
H 24.34858700 4.75343100 -5.05009100
C 32.14374300 5.94929400 14.57537300
C 33.89416800 0.20798900 -0.81587700
H 34.83702900 0.75664300 -0.74914900
C 32.45771800 -1.33099400 -2.03437900
H 32.27167800 -1.97550700 -2.89708200
C 31.67820800 7.25697000 14.76346800
H 31.24810900 7.83278400 13.94301400
C 32.75437600 5.76472600 16.93111000
H 33.15884800 5.19115500 17.76884200
C 32.67136600 5.19090700 15.65441200
C 28.47351900 -1.70640400 0.95384400
H 28.09710100 -1.31393300 1.89947100
C 32.90368900 2.78625300 12.92601700
H 32.66537600 2.80129700 11.86173500
C 31.77109600 7.80459500 16.04404100
H 31.41485400 8.82403600 16.21245300
C 32.30795600 7.07199100 17.11876200
H 32.36777700 7.53118800 18.10827500
C 23.97442100 6.61358700 -4.00873600
H 23.48010000 7.08282900 -4.86257500
C 33.64667200 -0.61018900 -1.93365300
H 34.39666900 -0.68400600 -2.72452800
C 28.16931400 -3.09843400 -1.03320300
H 27.54405500 -3.79082600 -1.60156400
C 27.70485400 -2.59061300 0.19432200
H 26.72078100 -2.89050800 0.56305700
C 29.41622300 -2.71991500 -1.52981400
H 29.77131900 -3.10613500 -2.48846300
C 33.46622600 1.65110500 13.51198900
H 33.65711900 0.77258300 12.89049800
C 33.57180600 2.73032000 15.68892600
H 33.83770600 2.70770600 16.74873900
C 33.79394400 1.61687800 14.87993500
H 34.23135000 0.71168200 15.30759900
这两个关键词都是一样的，但结果却是输出不一样，不知道问题出在哪里
作者Author: sobereva 时间: 2019-10-31 06:20
明显关键词不一样，怎么会是一样的
连任务类型、算的态、用的基函数都不同（一个用了6D 10F）
作者Author: hanmingxi 时间: 2019-10-31 11:49
嗯嗯，还想问您下这两个使用哪个得到的激发能比较好，是SOC的还是S1单重态的？
作者Author: hanmingxi 时间: 2019-10-31 16:24
老师，请问我在用Gaussian+PySOC计算时出现下面的报错，之后我又运行了下例子是正常的，这个要怎么修改？下面是我的.log文件
Len(soint) 0
Traceback (most recent call last):
File "/home/zyguo/bin/pysoc/bin/soc.py", line 522，in <module>
s_ matr = read file( 'molsoc overlap.dat，AO_ overlap，nt, 0)
File "/home/zyguo/bin/pysoc/bin/soc.py". line 40，in read_ file
with open(filein, 'r') as f:
IOError: [Errno 2] No such file or directory: ' molsoc overlap. dat'
forrtl: No such file or directory
forrtl: severe (29): file not found, unit 171， file /home/zyguo/bin/pysoc/ input template/s
oC_ td_ input.dat
作者Author: sobereva 时间: 2019-11-1 07:29
pysoc这玩意我已经弃了，赶紧改用ORCA，更快更方便，还可以做到更准确
使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱
http://sobereva.com/462（http://bbs.keinsci.com/thread-12114-1-1.html）
作者Author: sobereva 时间: 2019-11-1 07:30
你还没有理请逻辑关系
每个结构下都能算激发能，没有好不好，只能说哪个对当前研究更有意义。诸如你要研究荧光发射能，显然得用S1极小点结构下的
看此文补补基本知识
Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法
http://sobereva.com/314（http://bbs.keinsci.com/thread-2413-1-1.html）
作者Author: 陈丢丢 时间: 2019-11-28 18:22
sob老师好，我想请问下pysoc或者ORCA软件可以在DFT下计算旋轨耦合矩阵元？单线态和三线态的。如果可以是直接把输入文件的关键词换成DFT的吗？
作者Author: sobereva 时间: 2019-11-29 08:25
不能
作者Author: 陈丢丢 时间: 2019-11-29 16:42
谢谢老师。那我想算DFT下的SOC，可以用什么软件呢？
作者Author: sobereva 时间: 2019-11-30 03:08
用NWChem做SODFT在DFT计算中考虑旋轨耦合效应
http://sobereva.com/368（http://bbs.keinsci.com/thread-5573-1-1.html）
要么Dirac等支持二分量或四分量相对论计算的程序，旋轨耦合效应自然就体现进去了
作者Author: q1w2e3 时间: 2020-1-3 18:33
sob老师好，我用ubuntu虚拟器安装了python2.7.17和pysoc，按照您的教程修改后，运行pysoc.py命令时出现报错：
/bin/sh: 1: /sob/pysoc/bin/soc.py: not found
/bin/sh: 1: /sob/pysoc/bin/soc_td: not found
请问是什么原因呢？谢谢！
作者Author: sobereva 时间: 2020-1-3 18:39
我只对CentOS用户提供帮助
作者Author: 13104468322 时间: 2020-6-9 15:46
Error when reading gaussian output
Traceback (most recent call last):
File "/home/user039/soft/export/pysoc/pysoc/bin/soc.py", line 491, in <module>
write_file(mo_ene, 'mo_ene.dat', '{} ')
NameError: name 'mo_ene' is not defined
forrtl: No such file or directory
forrtl: severe (29): file not found, unit 171, file /home/user039/guo-ceshi/lingf/hbps/pysoc/wang/s1/soc_td_input.dat
Image PC Routine Line Source
soc_td 0000000000426956 Unknown Unknown Unknown
soc_td 000000000043288B Unknown Unknown Unknown
soc_td 0000000000407B24 init_prep_mp_read 26 init_prep.f90
soc_td 0000000000402D23 MAIN__ 14 soc_td.f90
soc_td 0000000000402CCE Unknown Unknown Unknown
libc.so.6 00000037B7A1ED5D Unknown Unknown Unknown
soc_td 0000000000402BD9 Unknown Unknown Unknown
sobereva老师，出现这种报错是怎么回事@sobereva
作者Author: sobereva 时间: 2020-6-9 16:03
要么操作不对，要么python版本不兼容
建议用ORCA
作者Author: Believe 时间: 2020-7-20 22:18
sob老师，您好！我想知道怎样通过Gaussian+PySOC计算的旋轨耦合矩阵元计算磷光速率？“通过本文计算的旋轨耦合矩阵元计算磷光速率，笔者之后会另文介绍。”谢谢您~
作者Author: sobereva 时间: 2020-7-21 00:33
需要自己写程序实现，公式可参考Chem. Rev. 2017, 117, 6500
我也写过这么个程序，但发现结果不如Dalton通过响应理论算的稳健，也就没有再写博文介绍
作者Author: Believe 时间: 2020-7-21 08:46
好的，谢谢sob老师~
作者Author: shuifangren 时间: 2021-10-15 16:45
sob老师，您好，您在本贴中提到“至于具体如何通过本文计算的旋轨耦合矩阵元计算磷光速率，笔者之后会另文介绍。” 我想用PySOC得到的SOC结果去计算磷光速率，但是没有找到您其他博文关于这部分内容的，所以想求助一下~
作者Author: sobereva 时间: 2021-10-16 03:14
用Dalton算。靠SOC矩阵元结合微扰公式，根据我的经验，明显不如Dalton基于响应函数算得准
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |