计算化学公社
标题: 使用sobMECP程序结合Gaussian程序搜索极小能量交叉点 [打印本页]
作者Author: sobereva 时间: 2015-3-18 07:47
标题: 使用sobMECP程序结合Gaussian程序搜索极小能量交叉点
使用sobMECP程序结合Gaussian程序搜索极小能量交叉点
Using sobMECP program combined with Gaussian program to search for minimum energy crossing points
文/Sobereva@北京科音 First release: 2015-Mar-18 Last update: 2022-May-27
本文主要介绍怎么用sob修改版的Harvey的MECP程序（称sobMECP）结合Gaussian搜索两种不同自旋多重度的态之间的极小能量交叉点(MECP)。第一节先简单介绍一下基本理论和算法方面的知识，第二节介绍一下sobMECP程序的特点和使用流程，第三节给出一些具体实例。
1 原理
1.1 MECP的相关理论知识
我们先看看什么情况下两个态的势能面才会交叉。设两个态波函数分别为|1>和|2>，体系哈密顿算符为H，则这样的两态模型构成2*2哈密顿矩阵，其四个矩阵元包括两个对角元H11=<1|H|1>、H22=<2|H|2>，以及两个等价的耦合矩阵元H12=<1|H|2>=H21=<2|H|1>。让这两个态的能量E1、E2相同需要同时满足两个条件：(1) H11=H22 (2) H12=0。
这些哈密顿矩阵元的数值都是体系几何结构的函数。对于两个自旋多重度不同的态，在不考虑旋轨耦合情况下，由于波函数自旋部分不同，导致耦合矩阵元H12总为0，因此只需要满足H11=H22就够了。对于非线性分子，坐标是3N-6维的，我们可以对坐标进行某种变换，使得H11-H22只依赖于某一维坐标，只要这个坐标恰为某个值时就能满足H11=H22，而不管其它3N-7个坐标的数值是如何的。这也就是说，不考虑旋轨耦合下两个自旋多重度不同的态的交叉并不是一个特定的点（对应某个具体的结构），而是在3N-7维的超曲面中交叉的。换句话说，能让两个态能量相同的结构有无限个，自然不可能都考察，我们研究的只是其中最有意义的一个点，这一般就是极小能量交叉点(MECP, Minimum energy crossing point)，它是指两个态能量相同的3N-7维超曲面中能量最低的那个结构。但也不一定MECP就只有一个，正如同几何优化时往往可以得到多个能量极小点，两个态之间也可以有多个MECP，能找到哪个取决于MECP搜索时的初猜结构是否离它比较近。
一般搜索MECP的时候是不考虑旋轨耦合的。如果在此之后在计算能量时将旋轨耦合算符引入哈密顿，由于此时两个态之间耦合矩阵元H12=H21不再为0，两个态就会发生混合，产生A、B两个新的态，波函数表示为c1*|1>+c2*|2>形式，且能量相应地出现小幅分裂（即E_A≠E_B）。这种情形类似于自旋多重度相同的两个态之间的避免交叉极小点（Avoided crossing minimum），因为都是H12=H21≠0.
系间窜越(Intersystem crossing)是重要的光化学过程，是指电子激发后，由于不同自旋多重度的态的势能面之间存在交叉，导致体系经历这样的结构时以非辐射方式改变自旋多重度。因此，寻找MECP的结构对光化学研究十分重要。
1.2 MECP结构的搜索方法
搜索MECP结构涉及两个问题：（1）理论方法与基组（2）优化算法。
搜索MECP可以用的理论方法有两类，一类是做态平均CASSCF计算（或MS-CASPT2、MRCI等进一步考虑动态相关的），同时计算多个根，并搜索指定的两个自旋多重度不同的根之间的MECP。另一类是用普通的理论方法，如DFT、MP2、CCSD，分别对两个不同自旋多重度的态做单点计算，并由此搜索这两个态的MECP。一般来说基于DFT搜索MECP是首选，便宜、方便，结果也不错，可以用在颇大的体系。泛函、基组的选择和常规DFT计算并无区别。如果对结果要得很精确，可以再考虑CCSD或CASSCF。
（PS：这里说的都是搜索自旋多重度不同的态之间的MECP。搜索自旋多重度相同的态之间的MECP或避免交叉极小点只能用CASSCF，但如果能量分裂较大用TDDFT等单参考态方法也凑合，这些与本文无关就不多说了。）
Robb提出的搜索MECP的算法其实很简单，和普通几何优化算法思路颇为类似，这也是最常用的MECP搜索方法。定义两个有效梯度，f和g。f是(E1-E2)^2对核坐标求导，具体写为E1-E2乘以两个态的梯度矢量之差。沿着f移动坐标会令两个态能量差减小。g是E1对核坐标求导，并且投影掉了它在f方向上的分量以使f和g正交。沿着g移动会令第一个态的能量降低。若将f和g考虑到一起，按照f+g作为体系有效受力进行优化，那么在优化过程中就会一边降低两个态的能量差一边降低体系总能量，显然最后得到的就是MECP结构。优化过程用的具体数值算法和一般几何优化一样是基于牛顿法或赝牛顿法，只不过Hessian矩阵也相应地和普通几何优化不同，详细公式见Theor. Chem. Acc., 99, 95(1998)。
MECP搜索一般同时使用位移、受力、两个态能量差作为收敛判断标准，三者都非常小的时候才算收敛。但如果把能量差必须接近0这条收敛标准去掉，那么上面介绍的MECP搜索算法也可以直接用于搜索避免交叉极小点。
搜索MECP用的初猜就取平衡结构或过渡态结构就行，离预期的MECP结构越近越好，但并不需要很精确。上述MECP搜索算法还是比较稳健的。如果搜索不成功可以再考虑调整初猜结构。
2 MECP程序的使用
Harvey的MECP程序基于Gaussian，只能用于Linux平台，搜索可以在HF、DFT、MP2等级别下进行，运行过程中会生成Gaussian输入文件、调用Gaussian进行计算，并且读取受力，对结构按照搜索算法进行位移。优化算法基于L-BFGS赝牛顿法。MECP原版可以在这里下载：http://pan.baidu.com/s/1o682NB8。使用时应引用Harvey对苯基阳离子搜索MECP的文章Theor. Chem. Acc., 99, 95(1998)。值得一提的是MECP搜索算法并非是Harvey提出来的，MECP程序里用的，也即上一节介绍的算法是Robb提出的标准的MECP搜索方法。之所以Harvey的MECP程序知名度较高，是因为它支持的是最常用的Gaussian，而与此同时Gaussian自身仅支持用CASSCF搜索MECP（众所周知Gaussian的CASSCF不给力，而且CASSCF本身又相对复杂、耗时），这才导致了MECP程序在文献中用得很多。
笔者发现这MECP程序用起来颇不方便，每处理一个新体系要改好多地方，用法上也有不少别扭的地方。笔者遂对此程序进行了很多修改，使之变得方便好用得多，修改版称为sobMECP，下载地址：http://sobereva.com/soft/sobMECP.zip。下文讨论全都是对于sobMECP而言的。如果你使用了sobMECP发表了文章，请在文中以这种格式引用：Tian Lu, sobMECP program, http://sobereva.com/286 (accessed month day, year)
2.1 sobMECP程序的基本使用流程
将sobMECP解压到任意一处，进入其中，然后按下列步骤使用之。
1 设定好Input_Header_A和Input_Header_B文件，内容分别是两个态的分子坐标之前的部分。%chk必须得写，chk文件具体怎么命名无所谓。force和guess(read)这两个关键词是必须写的，其它关键词，包括方法和基组、辅助SCF收敛的选项等，根据实际情况去写。必须用#n。
如果需要用gen或genECP从分子坐标末尾读取基组/赝势定义，则需要在Input_Tail文件中以常规方式写明基组/赝势定义，比如
C H 0
cc-pVTZ
****
Cu 0
SDD
****
[空行]
Cu 0
SDD
[空行]
[空行]
如果不需要自定义基组/赝势，则Input_Tail文件留空即可。
2 在geom里填好体系初始结构，注意末尾空一行，原子名必须用元素的序号代替。例如：
6 -0.17110831 0.45776462 0.00000000
8 -0.01984499 -0.74914326 0.00000000
1 -1.10138660 1.00356716 0.00000000
1 0.76276200 0.95429984 0.54191098
1 0.76276200 0.95429984 -0.54191098
[空行]
3 运行./prepare.sh。这将会在当前目录下编译产生MECP.x可执行文件、运行进度文件ProgFile、主脚本runMECP.sh，并且清空当前目录下的零碎文件和JOBS目录。
4 运行./runfirst.sh，这将会在初始结构下对两个态进行计算，产生Job0_A.gjf、Job0_B.gjf、Job0_A.log、Job0_B.log以及相应的chk文件。
5 运行./runMECP.sh，即开始搜索MECP，直到全部标准收敛为止，收敛信息会不断输出到屏幕上以便监控。随着搜索的进行，geom文件会不断更新为当前结构，ab_initio文件内容会不断更新为当前两个态的能量和梯度，ProgFile里记录的进度信息也会不断更新。
运行过程中产生的各种细节信息都会输出到当前目录下ReportFile文件中。包括每一步的结构、两个态的能量、收敛情况、有效梯度（即前文提到的f+g，但这里f额外乘了个刻度因子）、差值梯度（两个态受力之差）、平行梯度（前文的g）。
运行过程中产生的每一步的Gaussian输入输出文件都在JOBS目录下。
MECP搜索过程中产生的每一步的结构都会记录到traj.xyz下作为轨迹文件，可以用VMD程序打开来方便地观看搜索过程中的结构变化，第一帧对应于初始结构。
如果运行中途中断，应重新依次执行prepare.sh、runfirst.sh、runMECP.sh，这会基于geom文件中储存的最后一步结构重新进行搜索。
2.2 sobMECP程序的相关细节
默认情况prepare.sh是调用gfortran编译器来编译MECP程序。如果机子里没有gfortran，可以改用ifort编译，做法是在prepare.sh里在gfortran前头写上#将之注释掉，而把ifort那行开头的#去掉。
如果要调整收敛判断阈值，应修改temp/MECP.f的第300行。总共有5个标准。TDE是能量变化，TDXMax和TDXRMS是位移最大值和方均根，TGMax和TGRMS是受力最大值和方均根。默认情况下，受力和位移判断标准比Gaussian09几何优化默认判断标准松将近一倍。
如果要调整搜索步数上限，应修改temp/runMECP.sh第30行的数字。默认是100步。
有sobMECP用户反映把优化过程的置信半径减小五倍，往往令收敛更容易、降低震荡倾向。我未充分考证这一点，读者可以试试。做法是在MECP.f里搜索STPMX = 0.1d0，将数值改小。
进行如上修改后需再次运行./prepare.sh方可生效。
extract_energy是Linux下awk工具的脚本文件，用于从Gaussian输出文件中提取能量。提取不同理论方法输出的能量需要用不同的脚本。temp目录下自带了四种：
extract_energy_SCF：提取HF/DFT/半经验的
extract_energy_MP2：提取MP2的
extract_energy_CIS：提取CIS的
extract_energy_TD：提取TDHF/TDDFT的（对于G16用户，此文件里的TD-KS必须手动改为TD-DFT）
当前计算在哪个级别下进行，就把上述哪个文件拷到上一级目录下改名为extract_energy。程序包里直接带的extract_energy对应的是extract_energy_SCF。
如果第一个态和第二个态用的理论方法不同，比如第1个态用TDDFT计算S1，第二个态用UKS计算T1，那么应同时准备extract_energy和extract_energy2，前者会被用于提取第一个态的能量，后者会被 用于提取第二个态的能量。
如果用的Gaussian不是09版，需把sub_script中的g09替换为相应命令。笔者只测试了Gaussian 09 D.01，可以完全正常运行，其它版本G09以及G03、G16应该也可以正常使用。
3 MECP搜索实例
下面给出一些sobMECP程序使用例子，涉及到的文件都在sobMECP的test文件夹里。
3.1 CH2（卡宾）
这一节寻找CH2的单-三重态间的MECP。
在sobMECP文件夹中编辑Input_Header_A文件成为如下内容，对应单重态计算
%mem=6GB
%nproc=4
%chk=singlet.chk
#n B3LYP/6-311G** force guess(read)
First State
0 1
然后，类似地将Input_Header_B文件写为下面这样，对应三重态计算
%mem=6GB
%nproc=4
%chk=triplet.chk
#n B3LYP/6-311G** force guess(read)
Second State
0 3
注意上面两个文件末尾无空行，后同。
然后编辑geom文件，将下面的初始结构写进去
6 0.00000000 0.00000000 0.13397933
1 0.00000000 -0.92611695 -0.40193800
1 -0.00000000 0.92611695 -0.40193800
[空行]
然后用chmod +x *将sobMECP下的文件都加上可执行权限。确保Gaussian 09可以在命令行下用g09直接正常调用。之后依次运行
./prepare.sh
./runfirst.sh
./runMECP.sh (运行这步之前可能还得执行一次chmod +x *，因为MECP.x是新编译出来的)
收敛情况随着计算不断在屏幕上输出，仅用了5步就收敛了，屏幕上显示5个YES。此时geom文件中的结构就是MECP结构了。从ab_initio文件中可见，两个态的能量分别是-39.1443410840和-39.1443771072 a.u.，确实基本一致。搜索过程的详细信息在ReportFile里可看到。
将traj.xyz拖到VMD程序的主窗口里可以看到搜索过程结构的变化，可见一开始H-C-H角度是119.9度，而在最后一帧，即MECP结构下，角度减小为100.1度。
3.2 FeO+
这一节寻找FeO+阳离子的四-六重态间的MECP。对O用6-311G*，对Fe用SDD赝势。对应四重态的Input_Header_A文件应当为
%chk=A.chk
#n B3LYP/genecp force guess(read)
First State
1 4
对应六重态的Input_Header_B文件应为
%chk=B.chk
#n B3LYP/genecp force guess(read)
Second State
1 6
还要写一个Input_Tail文件，用来设定要从分子坐标后面读取的内容，这里就是自定义基组和赝势信息：
O 0
6-311G*
****
Fe 0
SDD
****
[空行]
Fe 0
SDD
[空行]
然后在geom里写入初始坐标
26 0.000000 0.000000 0.000000
8 0.000000 0.000000 0.670000
[空行]
然后也是依次运行prepare.sh、runfirst.sh、runMECP.sh。
初始结构中Fe-O长度为0.67埃，从输出结果可见，最后Fe-O长度变为了1.322埃。对于这样只有一个几何变量的体系，显然也可以自行对两个态做势能面扫描，然后拟合成曲线来确定交叉点。
注：在2.0埃附近实际上还有个交叉点，但是在附近区域，当前理论方法下两个态的势能曲线几乎完全平行，此时MECP搜索算法完全失效，因此没法让初猜Fe-O长度在2.0附近来寻找这个点。这个交叉点只能通过自行考察势能曲线获得。
3.3 C6H5+
苯基阳离子C6H5+正是使用MECP需要引的Theor. Chem. Acc., 99, 95(1998)这篇文章中研究的体系，文中使用不同方法考察了它的单-三重态MECP结构。类似前面的例子，也是先写Input_Header_A、Input_Header_B、geom文件，然后依次运行prepare.sh、runfirst.sh、runMECP.sh。涉及的文件在sobMECP/test/C6H5+中可找到，这里就不累述了。
这里用的初始结构就是把苯去掉一个氢而已，在B3LYP/6-31G**下经过十几步就找到了想要的MECP。体系依然是C2v构型，但是环发生了一定变形。
sobMECP的test目录下还有两个单-三重态MECP搜索例子，其中H3CO+就是原版MECP程序中自带的一个例子，Pt_coord是个含铂配合物，算是个较大体系的例子。
3.4 激发态单重态与三重态的交叉
前面考察的单-三重态的MECP，单重态都是直接用DFT算的基态。但是激发态单重态与三重态的交叉也很重要。比如S1-T1的交叉对于TADF（热活化延时荧光）的研究是关键，和反系间窜越效率问题密切相关。
test目录下有几个这样的例子，C6H5+_A2-3B1、C6H5+_B1-3A2、Pyrrole_A2-3B2，横杠前是单重态的电子态，横杠后是三重态电子态。这些任务计算第一单重态激发态用的是TDDFT，算三重态用的和之前一样是UKS。计算的时候要把相应目录下extract_energy和extract_energy2都拷到sobMECP目录下，前者提取第一个态的能量（TDDFT给出的），后者提取第二个态的能量（DFT算的）。注意像这些高对称结构，不同初始结构下T1电子态的不可约表示经常不同。找S1与三重态的MECP时由于不断读取上一步的初猜，所以跟踪的是初始结构下的T1态（这个态到了MECP结构下，未必是能量最低的三重态了）。而TDDFT计算时，S1的电子态的不可约表示则可能发生改变。
虽然原理上也可以用TDDFT算三重态能量，但之所以这里用UKS来算，一方面是更省时间，另一方面是结果往往更合理。
4 关于MECP处的自由能
有人问怎么计算MECP的自由能。一般没必要算它的自由能，绝大部分涉及到MECP的文章中也只讨论电子能量。但如果非要说怎么算的话，应当是在MECP处计算两个自旋态的能量简并的3N-7维空间里的3N-7个振动模式的频率，这是有意义的，因为在这3N-7维空间中MECP相当于极小点。有了振动频率，就可以用常规方式计算热力学量了，比如可以把振动频率等信息提供给Shermo来计算，见《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）。然而，这3N-7个频率据我所知没有简单现成的工具可以计算。
网上有文章称在MECP处用Gaussian里的freq=projected关键词可以得到MECP的自由能（此关键词本来目的是对IRC上的点计算垂直于IRC方向的频率），这明显不具备普适性。例如用上面的C6H5+的例子，在MECP处对单重态和三重态分别做freq=projected时，自由能分别为-231.188990和-231.186629 Hartree，可见存在显著差异，即结果依赖于自旋多重度的选取，结果不唯一。因此对于非特殊情况，靠freq=projected关键词是没法得到MECP的自由能的。其问题在于，从MECP位置可以对两个自旋态用IRC关键词跑downhill路径，见《谈谈Gaussian产生downhill路径的功能》（http://sobereva.com/571），对于MECP来说freq=projected算的是垂直于downhill路径的3N-7维空间的振动频率。而MECP处两个自旋态的downhill路径不仅方向往往不同，而且也往往都不垂直于唯一决定H11-H22的方向，即垂直于downhill路径的子空间不等同于两个态能量简并的3N-7维空间，因此靠freq=projected一般是没法正确计算MECP的自由能的。
作者Author: nkallwar 时间: 2015-3-20 23:59
前几天在QQ群里问过关于sob老师一个自旋态转变的问题，继续小白提问时间哈：
讨论的是一个 L-Cr(IV)-bridge--Cr(II)-L' 的混价金属配合物， 首先计算得到了七重态和三重态的结构和相应能量，二者的结构非常的接近，而能量（E(ZPE))也只相差了0.74kcal/mol，吉布斯自由能差别更小，只有0.03kcal/mol, 这两种自旋态几乎就是简并的。
期间努力寻找过五重态的结构，未遂，各种scf不收敛，sob老师认为此种结构的五重态很难得到，我也认为可能是五重态的电子结构太差，所以很难收敛。 三重态可以用 Cr（II，αααα）--Cr（IV，ββ）这种反铁磁性耦合来说明，七重态则是六个d电子均是自旋平行。
配合物可以认为是七重态和三重态两种状态的混合，两种状态的占有符合boltzman分布，当温度改变时，两种状态的分布概率会相应的改变，因此，会有一部分 三重态转变为七重态，或者反过来。 这里面要讨论的就是 三--七重态之间的转变是否要经历一个五重态的阶段。
先是看了一篇硕士论文里关于自旋翻转的必要条件是 |ΔS|=1，否则旋轨耦合矩阵元H12=0，无法完成自旋翻转。
而sob老师提出“如果通过势能面交叉来改变自旋多重度，和|ΔS|这个没关系。两个态没有耦合（H12=0）时满足H11=H22时能面就能够相交，由此导致自旋翻转，也无需经过五重态过程”。
昨天看到sob老师这个关于MECP的帖子，里面的例子倒都是|ΔS|=1，不知道像 三重态-七重态这种过程能不能直接用sob-MECP程序直接研究？
另外，关于自旋翻转的原理，又请教了一个做计算方法研究的老师，他给的回复是这样的：
”coupling matrix element一般可以用来衡量跃迁可能性的大小。这是一般性的结论，但势能面交叉的情况比较复杂。
不同的自旋态必须通过旋轨耦合作用才能翻转，因为其他的相互作用，包括vibronic coupling都不含自旋，是没有办法翻转自旋的。所谓的旋轨耦合矩阵元是衡量领头的跃迁几率（一阶微扰理论），如果这个领头项为零，那么就必须通过其它自旋态作为中介来达成跃迁（数学上对应高阶微扰理论），几率就会相对比较小。
但当势能面交叉的地方，原子核的运动非常迅速，这时候计算的这些电子态包括这个势能面本身是“不真实”的。换句话说，真实的分子波函数是几个电子态的组合。在这种情况下，即便两个自旋态需要多次翻转才能联系在一起，仍然有可能有显著的耦合。“
不是太明白“，这时候计算的这些电子态包括这个势能面本身是“不真实”的。换句话说，真实的分子波函数是几个电子态的组合”这段话的意思，不过结论好像还是要把 三重态和七重态连接起来，需要另一个自旋态（五重态）的阶段？
我考虑一下，如果没有电子转移导致Cr的价态变化的话， 那Cr（II）-Cr（IV）五重态可能是 (αβ，α，α） （αα），在八面体场中d4 有一种低自旋的基态是 3T1g，不过这种电子组态似乎不像是 （α，α，α，α）-（α，α）------>(α，α，α，α）（β，β）的中间态呢？
所以，如果在找不到五重态的情况下，还能找到 三重态--七重态的 转换机制吗？
作者Author: sobereva 时间: 2015-3-21 03:16
自旋翻转可以说有两种形式：
(1)两个态势能面相距较远，比如磷光发射，这种情况适用于Fermi golden规则。当不考虑旋轨耦合时H12就为0，由公式可见就不会发生跃迁，即翻转不了。
(2)透热过程，这对应于势能面交叉，即本文讨论的。此时适用于Landau-Zener方程（也可以通过Newton-X等程序进行激发态动力学模拟）。当H12为0时会完全交叉上，当耦合越强，H12越大，分裂也越大，发生透热过程的几率也会越低。
这两种方式都会导致自旋翻转，受H12的影响不同，(1)是H12越大越快，而(2)是H12越小越快。当势能面相距较远时(1)主导，当相距很近时(2)是主导。当能够发生(2)的时候，其速率是远远要大于(1)的。内转换速率比荧光快得多和这个道理是类似的。
所以说，一定要分清楚情况。H12=0不会发生(1)，但在交叉时(2)可能非常快。
|ΔS|=1这个规则我没印象。就算对于(1)的情况是必须的，和(2)也完全无关，(2)绝对不受这个条件支配，ΔS>1的两个态之间也照样可以交叉。我不觉得有必要非得经过五重态。
红字那段话我不做评论。
作者Author: lanxingren 时间: 2015-3-21 08:59
sob老师，文中提到“如果要调整搜索步数上限，应修改temp/runMECP.sh第30行的数字。默认是80步。”，在runMECP.sh中的第30行是“while ($num < 50)”这个是不是指搜索步数上限是50步。
作者Author: sobereva 时间: 2015-3-21 21:56
恩，应该是50步。帖子里已改。
当初向网盘里上传文件的时候文件里设的还是50步，后来用的时候我改成80步了，所以帖子和上传的文件里这个数不对应。
作者Author: studyba 时间: 2015-4-23 15:35
sob老师,请问一下，我运行到 [quote]把prepare.sh里开头的natm改成实际原子数。然后运行./prepare.sh。 这一步时，并没有在当前目录下编译产生MECP.x可执行文件、运行进度文件ProgFile、主脚本runMECP.sh。 显示natm=32:Command not found. 请问这是什么原因呢？
作者Author: sobereva 时间: 2015-4-23 15:58
你目前是bash环境否？
作者Author: studyba 时间: 2015-4-23 16:52
是的。 刚刚解决了这个问题，可以完全运行了。但是一直都不收敛。geom 里输入的是优化后的低自旋态的过渡态坐标，正确吗？
作者Author: sobereva 时间: 2015-4-23 18:46
并不要求优化过，更不需要是过渡态结构。建议先跑跑自带的例子。
作者Author: studyba 时间: 2015-4-24 14:07
好的！谢谢老师。
作者Author: studyba 时间: 2015-4-24 14:08
好的！谢谢老师
作者Author: studyba 时间: 2015-4-27 17:42
我运行例子了，算出来了。我参照例子改了该改的地方、为什么就是错误呢？一直找不到原因，请教sob老师
Problem with the ab initio Job.
cat:AddtoReportfile:No such file or directory
作者Author: sobereva 时间: 2015-4-27 17:49
光从这提示不好说，必须得反复尝试、改脚本之类。
作者Author: yaojihule 时间: 2015-6-19 15:27
sob老师 ，Landau-Zener公式中，我看Harvey 的文献，E=Eh-Emecp代表体系在跃迁方向上的动能，Emecp是最低能量交叉点处的相对能量，Eh指什么的能量？谢谢老师。
作者Author: sobereva 时间: 2015-6-19 21:39
这是哪篇文献的？我想不太起来了，我看看文献
作者Author: yaojihule 时间: 2015-6-23 09:56
本帖最后由 yaojihule 于 2015-6-23 09:59 编辑
嗯嗯 老师刚看到帖子，谢谢 我把文献发给你 在QQ上
作者Author: yaojihule 时间: 2015-6-23 10:01
谢谢老师
作者Author: 玥末夏初 时间: 2015-12-8 09:43
您好！用您的MECP计算，运行runMECP.sh时出错
[user2@localhost sobMECP1]$ bash ./runMECP.sh
./runMECP.sh: line 56: syntax error: unexpected end of file
试了dos2unix没有用，想问问您有没有什么解决办法，谢谢！
作者Author: sobereva 时间: 2015-12-8 16:03
运行时不要写上bash前缀
作者Author: 玥末夏初 时间: 2015-12-8 16:25
不写bash前缀，出错
[user2@localhost sobMECP1]$ runMECP.sh
-bash: runMECP.sh: command not found
作者Author: 玥末夏初 时间: 2015-12-8 16:53
我知道我哪里出错了，谢谢您了！
作者Author: 小范范1989 时间: 2016-2-15 17:32
sob老师好，我按照您说的步骤，我采用test中的例子想试试，出现下面的错误。
fanjz@node00:~/sobMECP> ls
extract_energy geom Input_Header_B prepare.sh sub_script test
extract_gradient Input_Header_A Input_Tail runfirst.sh temp
fanjz@node00:~/sobMECP> chmod +x *
fanjz@node00:~/sobMECP> ./prepare.sh
: command not found2:
fortcom: Severe: No such file or directory
... file is 'MECP.f'
compilation aborted for MECP.f (code 1)
: command not found7:
: command not found9:
: command not found13:
: command not found16:
说明一下：我已经修改了ifort，还有原子数目。我怀疑是不是我高斯0的问题，因为我们采用的是我们学院自己搭建的计算平台，高斯09好几个版本都安装在我们的集群上，平日提交任务都是使用脚本的形式，我把我采用的高斯09D01的脚本给你看看，谢谢老师。
#!/bin/sh
#PBS -N report
#PBS -l nodes=1:ppn=10
#PBS -l walltime=168:00:00
#PBS -q GAUSSIAN
#PBS -j oe
# Envrionment setting for Gaussian09
#
export g09root=/home/software/g09D01
PATH=$g09root/g09:$PATH
export PATH
source $g09root/g09/bsd/g09.profile
#
# End of envrionment setting for Gaussian09
mkdir -p /tmp/$PBS_JOBID
cp -r $PBS_O_WORKDIR/* /tmp/$PBS_JOBID
cd /tmp/$PBS_JOBID
trap "cp -r /tmp/$PBS_JOBID/* $PBS_O_WORKDIR;rm -rf /tmp/$PBS_JOBID" SIGTERM
g09 report.gjf
cp -r /tmp/$PBS_JOBID/* $PBS_O_WORKDIR
rm -rf /tmp/$PBS_JOBID
作者Author: sobereva 时间: 2016-2-15 17:37
牵扯到脚本提交的事我不清楚
报错提示肯定是.sh脚本牵扯到的命令的可执行文件没找到
你可以prepare.sh在本机上运行，runfirst.sh在集群上通过提交运行
作者Author: 小范范1989 时间: 2016-2-15 17:55
谢谢sob老师。
作者Author: healthyworld 时间: 2016-5-11 18:24
好东西
作者Author: blueyangliu 时间: 2016-10-15 20:12
Sob老师，我在自己的单机上运行算例没有问题。在集群上提交计算时， 最后一个runMECP.sh，出问题了，不知道如何在PBS文件中把runMECP.sh这个加入进去，Sob老师怎么弄呢？
作者Author: sobereva 时间: 2016-10-15 20:53
不知道，我对PBS不熟
作者Author: blueyangliu 时间: 2016-10-15 21:18
那只能在单机运行了吗
作者Author: sobereva 时间: 2016-10-15 23:58
总有办法通过提交方式运行，琢磨琢磨就有办法
作者Author: blueyangliu 时间: 2016-10-17 08:14
谢谢！
作者Author: 风飞 时间: 2016-12-16 11:11
老师：您好！如果在虚拟机中装linux系统，运行这个软件，是否有影响，因为怕虚拟机装linux兼容性不好，谢谢！
作者Author: sobereva 时间: 2016-12-16 11:32
毫无影响，完全正常
作者Author: 风飞 时间: 2016-12-16 21:11
谢谢！
作者Author: HERO 时间: 2016-12-24 20:06
本帖最后由 HERO 于 2016-12-24 20:41 编辑
sob老师，在运用sobMECP时遇到几个问题，1. 如果我想要寻找S1态和T1态的交叉点，input_Hearder_A输入S1态的优化构型，input_Hearder_B输入T1态的优化构型，那么geom初始构型怎样输入？ 是否是根据经验可能的构型。是不是程序通过A和B的构型寻找交叉点并输出到geom中？。在C6H5+_A2-3B1例子中，输入文件A和B的构型与geom_final是一样的,而与geom_init不一样，为何？2. 如程序正常运行结束，ab_initio中两个态的能量相差比较大，是不是意味着没有交叉点。 另外，怎样对结果进行分析？ 得到了一个交叉点的构型，是否能够给出ISC的比例等信息？ 3. 如果寻找的是S1和S0的交叉点，是否能够说明荧光淬灭或减弱的问题。 谢谢！
作者Author: sobereva 时间: 2016-12-25 06:15
input_Header里面输入的是计算用的关键词，而不是几何结构。几何结构是两个态共享的，在geom里定义。
你哪看到“A和B的构型与geom_final是一样的”？计算用的初始构型就是geom_init，改名成geom就能运算了。
程序只找个极小交叉点结构而已，不要以为一下子就能把什么ISC的比例、荧光减弱之类一大堆信息都得到，当前只是提供个初步信息而已，额外的都得自己再折腾，用诸如landau-zener之类公式、或者做surface hopping动力学之类，都比这复杂得多。
作者Author: HERO 时间: 2016-12-25 18:34
谢谢sob 学习中。
作者Author: 杨小狗 时间: 2017-2-20 08:44
请问您是如何解决的？
作者Author: ter20 时间: 2017-2-26 15:02
求问，sobMECP能做63原子，286电子这么巨大的体系吗？
作者Author: sobereva 时间: 2017-2-26 15:31
这个和sobMECP没关系，看机子以及计算级别
原理上，只要做几何优化能做得动，用sobMECP就没问题
作者Author: ter20 时间: 2017-2-27 08:11
OK，我先试试，谢谢Sob!
作者Author: ter20 时间: 2017-3-6 20:16
一个小问题，就是我在运行 ./runMECP.sh & 命令后，收敛信息会不断输出到屏幕上，这时如果我登出当前服务器远程窗口（exit），我发现计算任务就会终止，不知道如何解决这个问题，求助Sob老师。
PS: 我的ssh远程登录很恶心，一段时间不操作连接就断了，我试了反空闲设置，但是没用（我用的mac，不知道是不是设置得有问题，以前用windows也是这个毛病）。所以我没法一直保持屏幕输出状态
作者Author: sobereva 时间: 2017-3-6 20:28
nohup ./runMECP &
然后用exit登出服务器，就会一直运行，输出信息都在nohup.out里。
作者Author: ter20 时间: 2017-3-6 20:29
好的，我试试，非常感谢！
作者Author: SongChao 时间: 2018-4-14 11:50
本帖最后由 SongChao 于 2018-4-14 11:52 编辑
Hello:
Sorry for not installing Chinese input so I have to write in English.
sobMECP works perfectly with C6H5+, CH2 and H3CO+ examples. But when I reproduce MECP of CF3C(O)N and CF3NCO which published in a angewandte paper. prepare.sh and runfirst.sh runs normally, without any error comparing to those examples. But I get "fortran error" after running runMECP.sh immediately:
- songchao@linux-jygu:~/Documents/ResearchInProgress/attempt_MECP/CF3CON/MECP> ./runMECP.sh
- ProgFile Exists - OK
- First Input OK
- At line 230 of file MECP.f (unit = 8, file = 'ProgFile')
- Fortran runtime error: Bad integer for item 1 in list input
- Error termination. Backtrace:
- #0 0x2b89d0a03607 in ???
- #1 0x2b89d0a04115 in ???
- #2 0x2b89d0a04869 in ???
- #3 0x2b89d0abec10 in ???
- #4 0x2b89d0ac2b93 in ???
- #5 0x404a7c in ???
- #6 0x4052cb in ???
- #7 0x4054b8 in ???
- #8 0x2b89d1481724 in ???
- #9 0x400bc8 in ???
- at ../sysdeps/x86_64/start.S:118
- #10 0xffffffffffffffff in ???
- Problem with Fortran program
- Step Number 0 -- MECP not yet converged
- An error has occurred, possibly in the Gaussian Job
It seems the error is about MECP.f file in temp folder. I didn't find out how to correct it to make calculation perform normally.Is that a bug or something I didn't do correctly?
作者Author: SongChao 时间: 2018-4-14 15:12
I use atom symbol in the "geom" file, which is atomic number in the example. It solved.
作者Author: 风影月瞳 时间: 2018-5-1 21:23
sob老师，我运行CH2和FeO+的完美进行，但是跑激发态S1和T1的交叉时候总出现
ProgFile Exists - OK
First Input OK
Problem with the ab initio Job.
cat: AddtoReportFile: No such file or directory
Step Number 0 -- MECP not yet converged
，去看ab_initio里看到S1的能量没提取出来，试了g09和g16两个版本都是同样的错误，提取S1能量的代码是这样的
{
if ($3 == "E(TD-HF/TD-KS)") {
print $5
}
}
，对于代码我基本不懂，求解？
作者Author: sobereva 时间: 2018-5-1 21:56
去网上找点资料看看awk命令的语法就懂了
作者Author: 风影月瞳 时间: 2018-5-2 10:21
谢谢老师，已经搞定，我的高斯版本里TD能量项的关键词是TD-HF/TD-DFT，改一下就正常了。
作者Author: sobereva 时间: 2018-5-2 14:32
G16净瞎改...
作者Author: ter20 时间: 2018-8-17 09:28
使用了sobMECP，请问文献该如何引用？是不是引用Tian Lu, Feiwu Chen, Multiwfn: A multifunctional wavefunction analyzer, J. Comput. Chem., 33, 580-592 (2012)就可以了，谢谢！
作者Author: sobereva 时间: 2018-8-17 10:55
Tian Lu, sobMECP program, http://sobereva.com/286 (accessed 月 日, 年)
等以后有空，打算彻底重写一遍sobMECP代码，支持更多程序、功能和选项，作为个独立程序发表
作者Author: ter20 时间: 2018-8-17 11:06
谢谢Sob！
作者Author: ter20 时间: 2018-8-17 12:09
有个问题请教Sob，以卡宾CH2的MECP计算为例，单重态计算时为何不考虑用unrestricted DFT来描述开壳层的单重态双自由基的情况 (singlet diradical)？谢谢！
作者Author: sobereva 时间: 2018-8-17 21:54
是应该这样考虑，更得当
作者Author: ter20 时间: 2018-8-19 07:22
本帖最后由 ter20 于 2018-8-19 07:24 编辑
同时引用了MECP，sobMECP和Multiwfn，Sob麻烦看一下是否合适，谢谢 (请忽略具体文献格式)！
The crossing point (CP) was obtained with the minimum energy crossing point (MECP) method [1] as implemented in the sobMECP module [2] in the Multiwfn 3.3.9 program package[3]
[1] Harvey, Theor Chem Acc (1998) 99:95-99
[2] Tian Lu, sobMECP program, http://sobereva.com/286, (accessed on 12 July 2017)
[3] Tian Lu, Feiwu Chen, Multiwfn: A multifunctional wavefunction analyzer, J. Comput. Chem. (2012), 33, 580-592
作者Author: sobereva 时间: 2018-8-19 11:34
sobMECP并不是Multiwfn的组成部分，不需要引（这么引的话怕是有人会说洒家以非正当方式赚引用）
1、2都没问题，说成sobMECP program就行了，不用module这个词
作者Author: ter20 时间: 2018-8-19 13:46
好的，明白了，期待sobMECP早日发表，这样全球用户都可以使用了，非常有用的程序包！
作者Author: ter20 时间: 2018-8-29 13:13
就之前那个open-shell singlet diradical的问题，我测试了我自己的体系，发现用默认sobMECP程序 （即triplet/closed-shell singlet）得到的singlet态，做波函数稳定性测试，出现了RHF->UHF instability。我手动改了Input_Head_A的关键字为guess(mix,always)，得到了另一组结果，比默认程序得到结构的能量低。所以我建议大家在用sobMECP时可以根据自己体系的特征，确定是否需要考虑用开壳层单重态，纯属个人观点，希望Sob指正。
作者Author: loovfnd 时间: 2018-8-30 08:30
这里的MECI可以用软件求得吗？
http://bbs.keinsci.com/thread-10688-1-1.html
作者Author: sobereva 时间: 2018-8-30 15:29
那篇文中是一个特例，直接通过一维势能面扫描就能得到
考虑多维的时候才需要考虑用sobMECP
作者Author: 西柚西柚啊 时间: 2019-1-18 15:10
老师您好，我在测试ch2的例子时候，进行到./runfirst.sh，出现了这样的问题呢，请问应该如何解决呢
作者Author: sobereva 时间: 2019-1-18 19:25
没正确安装G09
作者Author: 西柚西柚啊 时间: 2019-1-18 19:51
老师打扰了，我在集群上如何调用g09，请问有什么方法或者途径吗？谢谢老师
作者Author: sobereva 时间: 2019-1-18 19:58
问管理员
你得能直接用g09命令启动g09才行
作者Author: 陈丢丢 时间: 2019-11-28 23:26
老师好，我用sobMECP软件算出了单线态和三线态之间的MECP，我现在想算MECP处单线态和三线态之间的梯度差，即两个态在MECP处的斜率差。我可以怎么计算呢？
作者Author: sobereva 时间: 2019-11-29 08:23
Gaussian里对这个结构用force关键词计算受力，分别算单单重态和三重态
作者Author: 陈丢丢 时间: 2019-11-29 21:47
谢谢老师。我今天加force关键词算了MECP处的不同重态，想要计算出2个势能面之间斜率的差值，是可以将两个out文件里面的那个数值进行相减呢？每个势能面在MECP处的斜率就是“Cartesian Forces: Max 0.009240955 RMS 0.001622638”里面的Max这个值吗？
作者Author: sobereva 时间: 2019-11-30 03:05
看你具体怎么定义斜率。势能面是3N-6维的，不同方向斜率都不一样。
作者Author: 陈丢丢 时间: 2019-11-30 09:02
本帖最后由 陈丢丢 于 2019-11-30 14:06 编辑
谢谢老师！我想的是要在MECP处沿着该势能面的切线的斜率，是可以看哪个数值呢？
作者Author: sobereva 时间: 2019-12-2 02:13
你先搞清楚交叉点附近的势能面结构
作者Author: 陈丢丢 时间: 2019-12-2 09:02
老师的意思是交叉点附近的结构也需要去算吗？
作者Author: sobereva 时间: 2019-12-3 07:34
本身是个高维曲面，切线的意义也不明确。说梯度方向那倒行
作者Author: 332544875 时间: 2020-9-9 19:16
请问Sob老师，这个sobMECP程序支持Gaussian的ONIOM模式吗？谢谢
作者Author: 函数与激情 时间: 2020-9-9 20:56
支持，需要自己修改下脚本
作者Author: 332544875 时间: 2020-9-10 08:59
谢谢你的回复，请问是需要修改脚本的什么内容？
作者Author: 函数与激情 时间: 2020-9-10 15:55
sob老师提供的MECP.x其实是个optimizer, 可以自己自定义接口其他软件的能量和梯度。你需要修改runMECP.sh使其产生出来的geom文件是oniom格式的，和QM不同有第二列，活动和固定的原子（0还是-1），以及最后一列分层（H M L），这步应该是最关键的，其余的修改HEAD.gjf比较容易了。我之前修改过并测试ONIOM中是可以运行的，也优化出了晶体中交叉结构，但是例子已经找不到了。希望可以帮助到你。
作者Author: 332544875 时间: 2020-9-10 16:58
谢谢老师，ONIOM基本计算的设置我是了解的。我在疑惑读取的ONIOM计算的能量还是读取TD-HF/TD-DFT和SCF done关键词后的能量吗？因为我看到ONIOM计算后输出文件里有多种能量，如high system：model energy，extrapolated energy，SCF done，TD-DFT等，具体每个代表什么含义我比较疑惑。
作者Author: 函数与激情 时间: 2020-9-10 19:40
能量读取的是“ONIOM: extrapolated energy =”后面的部分，这是ONIOM体系的总能量
作者Author: 332544875 时间: 2020-9-11 19:24
老师你好，我在使用sobMECP搜索计算ONIOM体系时发现生成的geom文件里冻结原子参数（0 -1）和高低层参数（H L）没有生成，请问老师你是怎么处理和解决的？谢谢
作者Author: 函数与激情 时间: 2020-9-12 09:35
把冻结和高低层分别整成两个单独的文件 每次运行runMECP.sh时将原来的geom与这两个文件合成一个新的geom文件
作者Author: Freeman 时间: 2021-2-8 11:23
社长您好。
我研究的交叉点的单重态部分是个自旋极化单重态，就要在Input_Header里写上stable=opt关键词，但是stable和force关键词是不相容的；如果把stable=opt换成guess=mix nosym um062x，那么guess=mix又和guess(read)不相容了。
看过sobMECP的几个sh脚本之后，我想到的解决方法是，先单独运行一下带有stable guess=mix的gjf文件，看看单靠guess=mix能不能使波函数稳定。如果是，那就在prepare.sh之前在Input_Header里加上guess=mix nosym um062x来产生出稳定波函数的chk文件（因为prepare.sh尚不涉及到guess(read)），然后在runMECP.sh之前把Input_Header里的guess=mix删掉。
但是这么做的问题是，guess=mix不比stable=opt稳妥，有时可能无法得到稳定波函数。所以想请问一下，有没有更好的办法。
作者Author: sobereva 时间: 2021-2-8 21:27
runfirst就相当于自动去掉guess(read)并且跑当前输入文件，用来产生初始的chk文件而已
你自己先直接手动用stable=opt算出个对称破缺单重态的chk文件，就没必要执行runfirst了，直接runMECP即可，这样计算的第一步就会通过guess(read)读你自己产生的chk
作者Author: shuifangren 时间: 2021-2-24 10:27
老师，我想求助一下，在执行第三步./runMECP.sh之后，输出的traj.xyz结构一直没有更新（和初始结构一样），Reportfile中没有内容，这是什么原因呢？
作者Author: sobereva 时间: 2021-2-24 12:08
看屏幕上的提示，弄清楚是正在算着还是怎么回事
在跑自己的任务前先把例子跑一遍
作者Author: shuifangren 时间: 2021-2-25 08:49
好的，老师，问题解决了，谢谢您~
作者Author: sui 时间: 2021-3-16 14:46
请问Sob老师，我运行了测试文件没有问题，算自己的文件时提示
ProgFile Exists - OK
First Input OK
Problem with the ab initio Job.
cat: AddtoReportFile: No such file or directory
Step Number 0 -- MECP not yet converged
Problem with the ab initio Job.
是哪里有问题呀？
作者Author: sobereva 时间: 2021-3-17 06:06
看中途产生的Gaussian输出文件判断原因
作者Author: sui 时间: 2021-3-18 16:24
谢谢老师
作者Author: Freeman 时间: 2022-1-20 22:12
社长您好。我最近在用bagel优化mecp，可是它的优化算法太不给力了，我就想把高斯作为跳板，间接用sobmecp来优化。现在一个问题是，高斯调用外部程序计算受力，不会生成带有波函数的chk文件，但是sobmecp要求gjf文件里要写guess(read)。这该如何是好？可以不写这个关键词吗？
作者Author: sobereva 时间: 2022-1-20 22:37
不是非得有guess(read)，这只不过是让sobMECP运行像普通的Gaussian优化一样，每一步读取上一步收敛的波函数，使得波函数的状态能延续下去，并且减少每一步达到SCF收敛的耗时。如果你不在乎这个，完全可以去掉
作者Author: Freeman 时间: 2022-1-24 14:51
社长您好。可以正常运行了，以下是运行中遇到的问题。
优化到某个点就开始振荡。这个点的S0/S1 gap已经很小了，我就以它为起点，再开一个优化任务，步长调成0.000500。但是优化完全是朝着反方向进行的。
- (base) [yzhangnn@login-0 sobMECP]$ grep State ReportFile
- Energy of First State: -913.3140316190
- Energy of Second State: -913.3140222100
- Energy of First State: -913.3141242030
- Energy of Second State: -913.3139641900
- Energy of First State: -913.3142225020
- Energy of Second State: -913.3138938610
- Energy of First State: -913.3142930820
- Energy of Second State: -913.3138421970
- Energy of First State: -913.3143875760
- Energy of Second State: -913.3137692700
- Energy of First State: -913.3144722690
- Energy of Second State: -913.3137019970
- Energy of First State: -913.3144550870
- Energy of Second State: -913.3137119090
- Energy of First State: -913.3145141190
- Energy of Second State: -913.3136749680
- Energy of First State: -913.3145143560
- Energy of Second State: -913.3136808920
- Energy of First State: -913.3145192510
- Energy of Second State: -913.3136830340
- Energy of First State: -913.3145278840
- Energy of Second State: -913.3136832340
请问有什么参数可以用来帮助它收敛吗？至少不要让它越优化，gap越大。
作者Author: sobereva 时间: 2022-1-24 18:33
可能是初始Hessian不好所致
可以从更早一些的结构作为起始结构试试
作者Author: lianghuashou 时间: 2022-3-8 19:26
sob老师，我刚刚安装测试了sobMECP，我的编译器是ifort，编译之后，用软件自带的例子做测试，用了Pt_coord 和 FeO的例子， prepare.sh 和 runfirst.sh 都可以正常运行，但是运行runMECP.sh的时候，报错信息为@: Expression Syntax. 此外没有其他的报错信息。检查了runMECP的文件，没有编程基础，不知道哪里有语法错误
还请sob老师帮忙看下
作者Author: sobereva 时间: 2022-3-9 04:52
我不知道你怎么运行的
这是csh脚本，当前运行环境必须有csh
作者Author: 好兄弟 时间: 2022-3-25 12:22
老师我之前运行sobMECP时都能正常运行，但我没注意sobMECP在后台运行时又运行了一个，导致之后运行时就出现了如下报错，自己尝试解决未果，希望老师能帮我排查下原因，教授我解决方案
(, 下载次数 Times of downloads: 38)
(, 下载次数 Times of downloads: 1)
(, 下载次数 Times of downloads: 1)
作者Author: sobereva 时间: 2022-3-25 15:43
重算
作者Author: 好兄弟 时间: 2022-3-25 15:57
我重算了依然不行，，，
作者Author: sobereva 时间: 2022-3-25 16:02
log文件看不出问题
之前的文件彻底清理干净，先把帖子里的例子跑出来再说
作者Author: 好兄弟 时间: 2022-4-10 10:51
老师您好，在输出五个yes得到mecp下结构后，是否还需要通过查看某些数值等方式去验证是否计算的是正确的mecp下的结构？
作者Author: sobereva 时间: 2022-4-11 00:26
达到收敛标准，而且肉眼看上去这结构就是你想要的，就够了。
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |