科学网 › 标签 › gaussian

标签: gaussian

相关帖子	版块	作者	回复/查看	最后发表

没有相关内容

相关日志

MaterialStudio & Gaussian中构建反应过渡态输入文件的方法: autodock 2020-4-15 21:00; 0. 化合物结构的构建，化合物结构的构建可先在ChemDraw中绘制平面结构并保存为cdx格式文件，保存后cdx文件用Chem3D打开，在Chem3D中使用Calculations——MM2——Minimize Energy对结构进行初步优化（这一步骤中Chem3D优化的分子结构的键长、键角），初步优化后的结构另存为mol2格式，在gaussian view中打开，另存为XX.com/XX.gjf文件，添加基组/泛函等关键词后再提交到工作站进行计算。正常计算结束后得到的log/out文件在gaussian view中打开。另存为mol/mol2/pdb格式的文件。产物文件和反应物文件的制作方法均相同。反应物或产物的结构也可以直接在MaterialStudio中制作，但是MS中制作复杂化合物的结构稍显麻烦。 1. 在materialstudio中打开优化后的反应物结构，复制反应物结构另存为 cis-opt.xsd 2. 同时打开优化后的产物结构，复制反应物结构另存为trans-opt.xsd，这两个结构的原子数应保持一致 3. 点击Tools-Reaction Preview，进行原子匹配。完成匹配后会生成cis-opt-trans-opt.xtd文件，在这一步骤中生成的过渡态结构数可根据需要自行调整。若使用TS命令进行过渡态搜寻，将生成的中间结构作为初始猜测结构，取出其坐标，添加相应的泛函/基组及其他关键词即可进行计算。 4. 选择软件栏中的gaussian程序进行计算 5. 在Task中选择TS Search，Task-More中选择QST2或者QST3，调整其他相应的基组/泛函参数，在Job Control中修改文件名字，保存，MS软件会直接生成相应的输入文件。 6. 点击Files保存输入文件，此时，MS中会自动生成XX.ginf格式的输入文件，修改文件格式为com或gjf。即可提交至Linux系统进行计算。（在MaterialStudio中制作输入文件时能根据关键词自动生成相应的文件，可大量减少格式错误。）采用QST2计算时的文件格式及部分关键词 %chk=cis-opt-trans-opt.chk %Mem=2048MB #p B3LYP/6-31G* SCF=(MaxCycle=65) OPT=(QST2) GEOM(PrintInputOrient) cis-opt-trans-opt TS Search 0 1 C 1.551400 -0.496700 0.071200 C 0.730300 0.557400 -0.095900 C -0.730100 0.557400 0.095900 C -1.551600 -0.496500 -0.071100 H 1.185600 -1.464800 0.404500 H 2.619400 -0.416300 -0.106900 H 1.165100 1.516300 -0.379200 H -1.164500 1.516500 0.379100 H -1.186600 -1.464700 -0.405000 H -2.619500 -0.415700 0.107400 Product 0 1 C -1.857400 0.108900 0.000000 C -0.610000 -0.401900 -0.000100 C 0.610000 0.401900 -0.000200 C 1.857400 -0.108900 0.000000 H -2.034200 1.182300 0.000000 H -2.736300 -0.527700 0.000300 H -0.474300 -1.484000 0.000300 H 0.474300 1.484000 0.000500 H 2.736300 0.527700 0.000200 H 2.034200 -1.182300 0.000100 采用QST3时部分关键词和输入文件格式 %chk=cis-opt-trans-opt.chk %Mem=2048MB #p B3LYP/6-31G* SCF=(MaxCycle=65) OPT=(QST3) GEOM(PrintInputOrient) cis-opt-trans-opt TS Search 0 1 C 1.551400 -0.496700 0.071200 C 0.730300 0.557400 -0.095900 C -0.730100 0.557400 0.095900 C -1.551600 -0.496500 -0.071100 H 1.185600 -1.464800 0.404500 H 2.619400 -0.416300 -0.106900 H 1.165100 1.516300 -0.379200 H -1.164500 1.516500 0.379100 H -1.186600 -1.464700 -0.405000 H -2.619500 -0.415700 0.107400 Product 0 1 C -1.857400 0.108900 0.000000 C -0.610000 -0.401900 -0.000100 C 0.610000 0.401900 -0.000200 C 1.857400 -0.108900 0.000000 H -2.034200 1.182300 0.000000 H -2.736300 -0.527700 0.000300 H -0.474300 -1.484000 0.000300 H 0.474300 1.484000 0.000500 H 2.736300 0.527700 0.000200 H 2.034200 -1.182300 0.000100 TS Location 0 1 C 1.105650 -0.742722 1.098915 C 0.856364 0.166043 0.139717 C -0.459354 0.705527 -0.211200 C -1.390922 0.162899 -1.013965 H 0.289277 -1.320363 1.525674 H 2.112332 -1.092088 1.305003 H 1.671463 0.507138 -0.497108 H -0.651278 1.687622 0.218517 H -1.142237 -0.373709 -1.923725 H -2.391570 -0.037623 -0.641341 *本文为本人在学习使用MaterialStudioGaussian过程中的一些笔记，其中计算中采用的关键词(基组/泛函等）不一定合理，请各位在使用过程中根据自己体系需要选择合理的关键词。; 个人分类: 科研笔记|8486 次阅读|0 个评论

北京超算上gaussian怎么跑: richor 2018-10-6 12:31; 首先要 source 一下 g09.pbs: #g09.pbs exportHOME=/home/fang_ap/soft exportPATH=$PATH:$HOME/g09:$HOME/g09/bsd exportGAUSS_EXEDIR=$HOME/g09 exportLD_LIBRARY_PATH=$LD_LIBRARY_PATH:$HOME/g09 exportGAUSS_SCRDIR=/tmp g09.profile 2018年12月11日开始，GAUSS_SCRDIR这句需要加，不然不能跑。或者.bashrc中加入，则这句不需要。然后 bsub g09.lsf 不加，会提交到cpuII队列。其中 g09.lsf ： #BSUB-Jhuang-BN #BSUB-W10080 #BSUB-n20 #BSUB-Rspan #BSUB-qc_fangap #BSUB-o%J.out #BSUB-e%J.err g09./*.gjf 中间的文件在跑完后会自动删掉！还可以不用lsf脚本的方式： bsub-Jzn_trp1.s-W10080-n20-Rspan -qc_fangap-o%J.out-e%J.errg09./zn_w1.gjf 但是有的 gjf 会莫名崩溃，而在其他（比如组内）机器上不会。比如 cu_w7_s1.gjf 在北京超算上的结果： beijing.log 组内机器的结果： hp63.log 对比发现，北京超算上是A.02版本，而自己装的是D.01版本。按该帖，旧版本的 Gaussian 有些已知或未知的 Bug 。 http://bbs.keinsci.com/thread-4829-1-1.html; 个人分类: gaussian|0 个评论

用gaussian进行几何优化: richor 2018-8-31 10:50; 【优化参数】 opt=redundant 与 cartesian 的区别是什么？对 gaussian 量化有影响： https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.540120510 the set of internal coordinates chosen to carry out the optimization can have a significant effect on the rate of convergence ; this is especially true for cyclic systems where variables are likely to be highly coupled . 很全的 manual: http://www.lct.jussieu.fr/manuels/Gaussian03/g_ur/k_opt.htm Note: No partial optimization or freezing of variables can be done with purely Cartesian optimizations ，所以 freeze 的话，用默认的是可以的。 By default , Gaussian performs the optimization in redundant internal coordinates. This is a change from previous versions of the program (Gaussian98). There has been substantial controversy in recent years concerning the optimal coordinate system for optimizations. See references therein. This one is good. https://onlinelibrary.wiley.com/doi/epdf/10.1002/jcc.540120510 For minimizations, redundant internal coordinates provide substantial improvements in optimization efficiency over Cartesian and nonredundant internal coordinates, especially for flexible and polycyclic systems. C. Peng, P. Y. Ayala, H. B. Schlegel, M. J. Frisch, Using redundant internal coordinates to optimize geometries and transition states. J. Comp. Chem. 1996, 17, 49 - 56. or pdf 【优化算法】 Berny algorithm: http://www.pitt.edu/~jordan/chem3430/optimization.pdf The berny algorithm in G03 is similar to EF （ Eigenvalue following ）基于 Hessian 矩阵的。; 个人分类: gaussian|1 次阅读|0 个评论

gaussian计算电荷布居: richor 2018-6-12 15:47; pop 关键字。比较完全的文档： https://wenku.baidu.com/view/be050ff09e3143323968937f.html amber力场使用的是resp方法。 resp电荷计算： http://signe.teokem.lu.se/~ulf/Methods/resp.html 需要amber tools里的antechamebr，建议直接安装AmberTools. net charge 设置： -nc 选项！; 个人分类: gaussian|4397 次阅读|0 个评论

pdb file的格式: richor 2018-5-29 19:14; 详见： http://www.wwpdb.org/documentation/file-format-content/format33/sect9.html 其中 ATOM 部分： gaussview 有时候处理 pdb 文件，有时候会有问题： https://www.researchgate.net/post/Why_some_pdb_files_are_not_accepted_by_GaussView But most of the atoms from a single molecule are accepted, only some of the atoms are not accepted. Should i do any modification in that particular atom.. Yes! For example, “Zn” should begin with 13 and written as Zn instead of ZN.; 个人分类: 分子模拟|2 次阅读|0 个评论

如何计算gaussian的自旋多重度: 热度 1 richor 2018-5-23 14:21; http://bbs.keinsci.com/thread-4836-1-1.html 体系的自旋多重度 = 单电子数 +1 ，或者是自旋多重度 =2s+1 。而对于一个含有过渡金属配合物则不能按照这样的规则来计算。按 sobereva 的教导：总有些菜鸟、弱智文献说“体系的自旋多重度 = 单电子数 +1 ”，这根本不对！碰见双自由基、反铁磁性耦合体系就瞎了。正确定义是“自旋多重度 =alpha 电子数 -beta 电子数 +1 ”。更详细的解释： https://www.researchgate.net/post/What_are_the_alpha_orbitals_and_beta_orbitals_and_their_importance_in_DFT_calculations 其实就是未配对电子数。鉴于大部分体系不存在未配对电子，因此一般设为1.; 个人分类: gaussian|21417 次阅读|11 个评论

gaussian的一些错误及解决办法: richor 2018-4-25 16:06; 1. Error ： New curvilinear step failed. 解决办法： Try to restart optimization from a slightly different input geometry. http://blog.sina.com.cn/s/blog_8be57be70100xabz.html 2. Terminating because of bad ONIOM partitioning 一般是原子坐标后面多了内容。 3. opt过程中， Bend failed for angle错误。原因是出现了180度角。 http://muchong.com/html/201405/7393235.html 解决办法，将优化结果重新做成gjf, 继续优化。; 个人分类: 个人|0 个评论

Gaussian高斯里计算resp电荷方法: CanonLife 2017-10-10 15:01; 原理：高斯根据计算所产生的每个原子周围的静电势进行拟合，在高斯输出文件 log 文件里有 ESP 电荷，但是 ESP 电荷做的不够好，如对于 -CH3 基团，一般认为三个 H 所带电荷相同，在 ESP 电荷里略有不同， RESP 做出来后三个 H 所带电荷相同。需注意的是，对于同一个原子， RESP 电荷与 ESP 电荷应相差不大，如同一个 H 原子，高斯计算出 ESP 电荷是 0.05 ，做出的该 H 的 RESP 电荷如果是 0.5 ，那就可能哪里出问题了，总之， ESP 电荷可作为参考。【计算方法】：第一步：用高斯计算 esp 电荷算 resp 的结构应该是优化之后的结构 . 基于此 , 有两种办法 : 第 1 种 : 优化和 resp 计算分开计算，先优化结构，再做 resp 计算。优化可以用 B3LYP 或者别的泛函，但拟合 RESP 电荷用的静电势应当在 HF 下计算，原因是此级别高估偶极矩，被认为可以等效体现溶剂效应。第 2 种 : 优化和 resp 同时做 , 这时泛函必须选 HF, 原因就是 resp 这一步必须采用 HF. 做 resp 的关键词有 g03 和 g09 两个版本： g03 版 : iop(6/33=2,6/42=10,6/41=10) pop=mk 每一个 iop 概念，查手册。例子： #p opt hf/6-31+g(d) iop(6/33=2,6/42=10,6/41=10) pop=mk g09 版： iop(6/33=2,6/42=6,6/50=1) pop=mk 例子： #p opt hf/6-31+g(d) iop(6/33=2,6/42=6,6/50=1) pop=mk 算完之后，可以在 log 文件里搜索【 ESP 】关键词，可以看到 esp 电荷，如下：第二步：做 RESP 电荷此处方法，用到 antechamber 程序，所以需先加载好 antechamber 路径，因为 antechamber 在 amber 软件包里有，所以可以单独加载 amber 路径，也可以单独加载 antechamber 路径，加载 amber 软件，例： export PATH=$PATH:/share/apps/amber9.exe export AMBERHOME=/share/apps/amber9 单独加载 antechamber 程序路径，例： export PATH=/opt/antechamber-1.27/exe:$PATH export ACHOME=/opt/antechamber-1.27/ 到此，可以采用 antechamber 进行拟合 resp 电荷。拟合 resp 电荷有两种办法：第一种： antechamber -i LIG.log -fi gout -o LIG.ac -fo ac -nc -1 -pf y respgen -i LIG.ac -o test1 -f resp1 respgen -i LIG.ac -o test2 -f resp2 espgen -i LIG.log -o LIG.esp resp -O -i test1 -e LIG.esp -o LIG.resp1out -t qout_stage1 resp -O -i test2 -e LIG.esp -o LIG.resp2out -q qout_stage1 -t qout_stage2 解释： -i 指定输入文件 -fi 指定输入文件类型 -o 指定输出文件 -fo 指定输出文件类型 -nc 指定体系所带总电量，此处为-1（如果没有-nc这一命令，则默认体系呈中性） -pf y 不产生中间文件（此命令可以不加，只是这样会产生一系列中间文件）做好的resp电荷存储在 qout_stage2 文件中, 该文件长这样：此文件中，只有原子电荷信息。同时存有resp电荷信息的是 LIG.resp2out 文件，如下： q(opt) 一列即所需resp电荷，与qout_stage2保持一致。该拟合相当于拟合了两次，第一次拟合得到LIG.resp1out文件，该文件中q(opt)一列电荷与qout_stage1文件中保持一致。第二次拟合得到LIG.resp2out文件，在LIG.resp2out中q(init)一列即LIG.resp1out中q(opt)一列。第二种： antechamber -i LIG.log -fi gout -o LIG.mol2-fo mol2 -c resp -nc -1 -pf y 所产生的 resp 电荷存储在 mol2 文件中：上图中最后一列即 resp 电荷。解释： -c resp 指定拟合原子电荷的方法，此处为采用 resp 方法。在上述拟合 resp 电荷过程中，可能会出现如下一些信息，这都不是错误信息，属于正常： Info: Bond types are assigned for valencestate 11 with penalty of 1 Info: the number of the ESP exceeds the MAXESP(20000),extendthe size and reallocate the memory automatically 第二种 info, 是因为在高斯 com 文件中设置 iop 时，采用格点数太大，出现此类信息为正常，高斯官方回复： It is just for information. The program reallocates memory itself. 那么，什么时候需要自己算 resp 电荷呢？在自己的以往工作中遇到需要算 resp 电荷的情况是：在用 amber 跑生物体系动力学时，最初建力场时，关于配体有机小分子的电荷，如果所研究的电荷转移能量转移等有涉及配体分子，则此时配体分子的电荷需要跑 resp 电荷。当然可以采用 bcc 电荷，但是一般认为 resp 比 bcc 更准确。采用 bcc 电荷方法： antechamber -i LIG.pdb -fi pdb -o LIG.mol2-fo mol2 -c bcc -nc -1 -pf y 解释： -c bcc 指定拟合原子电荷的方法，此处为采用 bcc 方法。完！; 个人分类: 高斯|24241 次阅读|0 个评论

def2-tzvp基组加弥散的方法，弱相互作用计算: yanrding 2016-8-24 11:21; 量化大神sobereva的帖子（http://bbs.keinsci.com/forum.php?mod=viewthreadtid=3487extra=）指出： Truhlar等人提出一种通用的给出原本不含弥散函数的基组以最低限度的弥散函数的策略，这样的基组以ma-开头，含义是minimal augmentation，详见J. Chem. Theory Comput., 7, 3027、Theor. Chem. Acc., 128, 295。也就是将原先基组中指数最小的s和p的指数除以3作为弥散函数的指数，但不对氢加弥散函数。这种处理后的def2基组可以直接从这里获得Gaussian格式的定义： http://comp.chem.umn.edu/basissets/basis.cgi 。如ma-def2-TZVP在里面被简写为ma-TZVP。对于大多数需要弥散函数的问题，实际上像这样只给重原子加一层s和p弥散就已经解决绝大部分问题了，因此如果想给def2加弥散但又不想花费过高代价的话，这种ma-方式加弥散是比较理想的。以下，针对使用ma-def2tzvp基组，提供一种具体的操作方式（另外的方式比如使用genecp当然也可以）。 1. 首先得到def2tzvp基组的具体形式（可以在网站下载，或者用高斯的def2tzvp + gfinput来输出）。把s和p的指数系数最小的两项系数除以3，作为弥散函数。（这个也可以到 http://comp.chem.umn.edu/basissets/basis.cgi 下载，但是可能有些元素不全）。 2. 编辑一个基组文件 ma-def2tzvp.gbs，放在$g09root/basis下面，内容包含所研究的所有原子的弥散函数信息（这个文件可以慢慢积累，以后还可以用）。包含C和O的文件的ma-def2tzvp.gbs文件如下（注意元素符号前面的-号，这样即使体系不含该元素，程序也不会报错退出）： -C 0 S 1 1.00 0.000000000000 0.3172148001D-01 0.1000000000D+01 P 1 1.00 0.000000000000 0.3352274557D-01 0.1000000000D+01 **** -O 0 S 1 1.00 0.000000000000 0.6168178787D-01 0.1000000000D+01 P 1 1.00 0.000000000000 0.5826140423D-01 0.1000000000D+01 **** 3. 使用的时候，大致是这个样子（考虑到弥散函数常常用于弱相互作用的计算，这里把B3LYP-D3(BJ)的形式也写进去了）： # b3lyp/def2tzvp extrabasis empiricaldispersion=GD3BJ opt freq test 0 1 C 0.00000000 0.42348200 0.00000000 O -1.12068900 -0.21083400 0.00000000 。。。 @/home/soft/g09/basis/ma-def2tzvp.gbs /N 上面最后一行如果不加 /N ，那么会把ma-def2tzvp.gbs 文件里面的所有内容写到log文件里面，看上去不好看。关于DFT-D3的使用，看了sob的帖子， http://sobereva.com/210 ，总结了一下：弱相互作用修正：empiricaldispersion=GD3BJ or GD3? 两者差别不大。Gaussian中有时对一个泛函只有一种，此时如果必须用另一种，则需要手工设置参数。以下为具体情况： 1. 只能用GD3的情况： (1) 对于明尼苏达系列泛函（M06...），用BJ会出错。 (2) 低版本（G09D及以下）用BJ时的freq计算会有问题，一些体系可能因此得到小的虚频。 2. 只能用GD3BJ的情况：有些泛函只能用GD3BJ，否则可能会出错。比如双杂化泛函DSDPBEP86，使用GD3BJ没问题，使用GD3会有下面的错误： R6DS8: Unable to choose the S8 parameter。以上问题随着高斯版本的升级可能会得到修正。以上问题之外，在GD3或GD3BJ都可以的情况下，一般推荐使用BJ。; 个人分类: Gaussian|21964 次阅读|0 个评论

进入三维超分辨：GB-STED: 热度 3 xipeng1 2016-1-18 09:36; 进入三维超分辨： GB-STED 席鹏 2016-01-14 前面的博文里面我讲过了 STED 这一诺奖技术。该技术的核心是，用一束光制造一个空心光环 ---- 老外叫做 ” 甜甜圈 donut“ （老外咋那么多吃货呢）。有了这个 donut 光，将它套在原来的高斯 PSF 上，利用受激辐射作为擦除机制，即可实现把荧光的 PSF 缩小，也就提高了分辨率。问题是，要想实现空心光环，特别是要越小越好的空心光环，绝非易事。目前，通常的做法是用一个 0-2 π 的位相调制器，结合圆偏振光进行。这样的话： (1) 对于光学上来说， 2 π 的位相调制就回到了 0 ，因此是一个螺旋式上升的连续结构； (2) 这一位相调制无论从哪个角度剖开，都是一个 0- π 的台阶。由于从透镜前端到焦点刚好是一个傅里叶变换：我们可以看到，当 u=0,v=0 时，中心就是对于所有这些复振幅的积分或求和。由于 , , 这两者刚好相互抵消，因此中间的极大值变成了极小值，形成了非常美丽的面包圈结构。图 1 STED 通常 0-2 π 采用涡旋位相光栅实现 donut 。这一看似美丽的数学解有一个脆弱的地方：对于生物成像，样品会带来新的位相调制。由于生物细胞中，所有细胞器的折射率分布变化大和散射特性不同，导致它们对每一个角度的位相都会产生影响。这个扰动一旦变大， donut 中心就不再为零， STED 将会擦除 PSF 的中心地带，导致分辨率不但不能有效提升，甚至可能会变差。一句话， donut 不零， STED 不灵。这也是为什么大家看到很多漂亮的 STED 的成像结果，但都是二维图像的原因。而另一方面，人们研究了一些抗干扰的光束，如 Bessel 光， Airy 光等，并将其应用在光片成像、 OCT 等领域。 Bessel 光的核心是，通过一个 axicon 对高斯光束进行调制，可以实现一个针状的光分布，且这一针状光斑经过样品时，不会受到样品折射率变化的干扰。图 2 Bessel 光束调制原理。中心区域通过干涉增强实现了抗干扰的 Bessel 光。由于 STED 最怕干扰的是中空 donut ，在此条件下，我们把 axicon 和涡旋位相波片结合起来，实现了一个中空的竹子状光束。这一光束同样具有 self healing 的效果。图 3 Gaussian-Bessel STED 焦点分布示意图。接下来的工作就豁然开朗了：我们在琼脂样品上测试分辨率，发现在 155 微米的深度，我们能够达到和表面同样的分辨率。我们尝试了折射率失配的 PDMS 样品，发现能够达到 100 微米深度的超分辨。最后，我们制作了一个类脑白质的仿体，发现能够达到 100 微米的穿透深度，实现超分辨。在过去，很多超分辨的工作受限于样品散射等因素，被局限在二维世界；能研究的样品只有一些细胞，和非常浅层的脑成像。我们希望通过这一技术，实现深层的组织超分辨成像。相关成果被 Laser Photonics Reviews 作为 2016 年 1 月的封面文章发表。北京大学物理学院施可彬课题组的于文韬同学为本文第一作者，施可彬研究员和席鹏研究员为共同通讯作者。该工作得到了 973 国家重点基础研究发展计划和国家自然科学基金委的支持。 Wentao Yu et al, Super-resolution deep imaging with hollow Bessel beamSTED microscopy, Laser Photonics Reviews, 10, 147 – 152, 2016. 相关链接： http://onlinelibrary.wiley.com/doi/10.1002/lpor.201500151/abstract; 个人分类: STED|11111 次阅读|3 个评论

Matlab：Gaussian Model: lixujeremy 2015-9-15 14:27; Summary 开贴记录数据拟合之高斯函数。已知一组数据， x 表示数值， y 表示频率，如 Fig. 1 所示，图形非常接近高斯分布。 Fig. 1 Method （ 1 ）拟合高斯模型，手动方法在 Command Window ，键入 cftool ，出现 Curve Fitting Tool 。如 Fig. 2 所示，配置 x 、 y 及 Gaussian 。 Fig. 2 配置完成后，软件自动得到返回参数， R-square 达 0.9814 ，非常高，右侧图形也显示拟合高斯曲线与原始情况较为吻合。（试试改变图中 Number of terms ？）此时得到的高斯函数： y=1946*exp(-((x-442.9)/82.13)^2) （ 2 ）代码方法见 Codefit.m ，与手动方法返回结果一致。如 Fig.3 所示， Fitted Curve 是为高斯曲线。 Over 。 Fig. 3 gauss.rar Further Discussion 尽管该数据通过高斯拟合方法可以使得被解释的方差高达 98% ，但是该高斯函数的二阶导数等于 0 的解仍不能得到实数，更准确的说通过拟合逼近这一曲线仍不能得到该曲线的拐点。参考 Blocking A Non-Stationary Signal Using Wavelets ，将这一曲线视为非平稳信号，经连续小波变化（ Mexican hat ）处理，直接得到这一信号的二阶导数，找到二阶导数为 0 的序列。 Fig. 4 Fig. 4 从左至右的 2 图表示 1~200 尺度下的二阶导数，显然在 Scale 大于 75 之后，黑线仅与两条 0 等值线相交， 3 、 4 图表示 Scale=100/150 的情况，对应拐点是 469 和 383 。 gaussblock.rar References Normal (Gaussian) distrib ution, from Wikipedia, the free encyclopedia . Gaussian function, from Wikipedia, the free encyclopedia . Gaussian Models, MathWork . 利用 Matlab 解方程 .; 个人分类: Matlab|5391 次阅读|0 个评论

[转载]高斯优化过渡态的经典总结: bestener 2014-3-8 12:30; 一般地，优化所得驻点的性质（极小点还是过渡态）要靠频率来确定；而对过渡态，要确定反应路径（即到底是哪个反应的过渡态）必需要做IRC了，不然靠不住的（往往用QST找到的过渡态并不一定就是连接输入反应物和产物的过渡态）。在我们用 QST2 或 QST3 来优化过渡态时，需输入反应物和产物，实际上反应物和产物的输入顺序是没有关系的。就是说，先输反应物后输产物和先输产物后输反应物得到的是同样的过渡态。这也好理解， QST2 里对过渡态的初始猜测实际上是程序自动将输入的反应物和产物的各变量取个平均，所以输入顺序是没有关系的。对 QST3 和 TS ，需人为指定过渡态的初始猜测。上面说的反应物和产物的输入顺序没有关系，有个前提条件，就是反应物和产物的自旋多重度一致。对 QST3 ，因为是人为指定过渡态的初始猜测，所以没有影响；而对 QST2 ，过渡态的自旋多重度默认和后面输入的一致。如果反应物和产物的自旋多重度一致，那随便先输哪个都没关系；而如果反应物和产物的自旋多重度不一致，这时该怎么办？？？当反应物基态多重度为 1 ，产物基态多重度为 3 时，到底将过渡态的多重度定为 1 还是 3 ？还是两个都试，哪个能量低取哪个？假如取 1 ，在下来做 IRC 验证时， Reverse 还好办， Forward 却是按多重度为 1 做的，这样怎么能和多重度为 3 的产物连接起来？要将其关联起来，需有一个为激发态，这样才能在自旋多重度上保持一致。可这种处理对吗？还有更好的思想吗？另外，我觉得 IRC 难用得很，很难控制，也不知自己钻进了哪条死胡同！！一直对此很迷茫，希望各位大侠援助啊！对于这点确实很迷茫，我觉得好像得用 cas 解决，有看过类似的文献，上面有用 cas 得出过反应物基态经激发态再回到基态得出产物的反映路径。我没有试过，对于 cas 我很是头疼，理论和实践都一无所知。作 IRC 能说明哪些问题？ irc 做完，就说明我们完成了这个反应路径地计算，在 SUMMARY OF REACTION PATH FOLLOWING: ，我们可以看到分子过渡态的键长键角与能量随着反应坐标的变化而变化，如果你将反应坐标与能量作图，就可以得到一条过渡态曲线，由于 irc 计算是结构沿着反应路径的方向，在每个点进行优化的，所以如果你找的过渡态是正确且优化是成功的话， TS 确实是连接两个 minimum 的。一般 IRC 的两个终点并不恰好就是反应物和生成物，除非不断加大 Maxpoints 。然而此时 IRC 又往往会出错，或者并不能达到指定的 Maxpoints ，例如，我想正反两个方向各算 20 各点， IRC 计算正常完成，但结果中并没有 41 个构型。如果这时作能量 - 反应坐标曲线，如何确定反应物和生成物的反应坐标？如果你的过渡态确实做得漂亮，那么就一定可以成功地由过渡态找到与之对应的反应物、生成物。可这只是理论上的东东，现实中，由于反应的势能面实在太复杂，你千辛万苦所得到的过渡态是很难天衣无缝地连接反应物、生成物的。不断加大 Maxpoints 简直是 Kill time 。一般地，我们可以将 IRC 最后得到的两个结构的所有数据分别与你做 QST 时所用的两个结构进行对比，只要结构数据的差别、能量差别比较小，就可以基本认定（因为反应的势能面太复杂）你 IRC 是成功的了。说得好！取点越多越好，但我发现往往作不到这一点，所以我也是将 IRC 最后得到的两个结构的数据分别生产物态和产物态的两个结构进行对比。而且我发现，作 IRC 时频率计算给出的讯息很重要，呵呵！过渡态寻找小结刚刚做了一段时间的过渡态，期间碰到了许多的困难。寻找过渡态不是一件容易的事（对于我和大多数刚涉及量化的人来说），因此我希望通过写这个经验小结能对大家有些帮助。 1. 首先遇到的问题是，用哪种方法来寻找过渡态？ GAUSSIAN 提供的方法是 QSTN 和 TSN 方法。两种方法各有优点和缺点。 QSTN 方法特别 QST3 方法要求输入反应物，过渡态的猜测结构，产物这三者的结构。特别麻烦。但很管用，一般不会出现不收敛的情况。对于 TSN （对应关键词为 OPT=TS ）方法，只要求输入过渡态的初始结构，但这个初始结构非常的关键，如果结构不好，则很容易出现不收敛的情况。所以我建议，如果是刚开始做过度态的话，用 QSTN 方法是好的选择，等有了 “ 感觉 ” 之后，再用 TSN 方法。 2. 怎么解决经常出现的错误？在找过度态的时候，经常碰到的一些问题就是不收敛（ 1 ） . ，有一个错误的本征值（错误信息为： there is a wrong sign eigenvalue in hessian matrix.....) （ 2 ），和 LINK9999 错误导致退出。（ 3 ）对于不收敛的情况，可以分为两类，比如提示信息里的 CONVERGENCE FAILER 提醒收敛到了 10 （ -5 ），而此时你设定的 SCF 循环次数也仅仅是 64 步，那么完全有希望通过加大 SCF 循环次数来达到收敛的目的。倘若只收敛到 10 （ -3 ）或 10 （ -2 ），此时加大循环次数可能就没用了。结果还是 CONVERGE FAILER 。此时可采用 SCF=QC ，来达到强制收敛的目的。因为 SCF=QC （ LINK508 ）的计算量比默认的 L502 要大，所以不到万不得以就不用它了。出现第二个错误可以直接用关键词 OPT=NOEIGEN 来实现。 LINK9999 出错是因为已经走完了默认的步数，但还未完成。系统会自动跳出。出现这种情况大多数就是因为优化步数和 SCF 步数超过了默认值。可用 OPT （ MAXCYCLE=100 ）和 SCF （ MAXCYCLE=300 ）来改错。 3. 怎么样控制过渡态的优化，使得过渡态不至于收敛到其他的分子结构中去？我用 GAUSS VIEW 可以解决这个问题，当刚开始运行 GAUSSIAN 时，你用 GVIEW 去打开输出文件时，你可以看到你的过渡态的初始输入结构，当一个循环过后（从上一个 LINK502 到下一个 LINK502 ），你再打开输出文件，你就可以清晰地看到优化一步后分子的构型，这样就可以随时监控过度态分子的结构，倘若已经有收敛到其他分子构型的趋势时，你就可以把它给 KILL 了，而不至于需要等全部工作结束后，打开输出文件才知道已经不是想要的过渡态了。如果收敛到其它的构型上去，可以考虑缩小 OPT 的步长 .iop(1/8=2 或 3 ）即可。 4. 还需要加其它的关键词吗？建议在 OPT 中加入 CALCFC 。这样可以加大找到过渡态的几率。本人深有体会！先写这么多了，难免有错误和不恰当之处，还希望大家来指点和补充！我找过渡态的经验是：找过渡态的关键点，是如何根据反应体系的立体结构和反应过程化学键的断裂和形成的轨道的空间形状，提出一个合理初始的过渡态的结构，再进行优化。提出的初始过渡态结构越接近真实结构，就越容易找到，否则花一年时间也可能找不到一个过渡态。而找过渡态的 QSTN 方法、 TSN 方法的差别不是太大。这几个帖子都很好，我也找了十几个简单的过渡态，不是直接找到的，而是先用目标过渡态的四元环结构做，然后对四元环结构逐步 “ 修饰 ” 直到那个四元环就是我需要的四元环为止，这个方法来得一般比较简单，但也有无法 “ 修饰 ” 到目标过渡态的情况，其他方法都是一样的。 10 （ -5 ）表示十的 -5 次方。每个循环都会列出 SCF 计算的结果，和收敛值啊。可以找的到的。 iop(1/8)=2 相当于在找过度态的时候，以 2A （ A 为基本单位长度）为单位来寻找过渡态。就象有一百个抽屉里只放了一个苹果，当然是一个一个抽屉打开，找到苹果的把握最大了。此时意味着 iop(1/8)=1 。有虚频说明结构处于不稳定状态。过渡态的初始结构猜测应该根据已有的知识和经验啊，哪种构型或构象最可能是过渡态，注意有时候猜测可能与计算结果完全相反。反应坐标势能面是研究化学反应历程的基础。对一个势能面来说，我们比较关心其关键点。势能面上的关键点是指势能的极值点，包括反应物 (R) 、产物 (P) 、中间体 (I) 与过渡态 (TS) 或鞍点（ saddle point ）等。连接反应物、过渡态和产物的反应途径，是一条能量最低的路径（ minimum energy path ），称为反应坐标（ reaction coordinate ）。各关键点是通过反应坐标联结起来的，过渡态是其上面的一级鞍点。＊不同类型的反应坐标我们常说的反应坐标、最小能量路径等用语，严格地说，是有一定区别的。（ 1 ）最速下降路径：也称最小陡降路径。从过渡态出发，体系会沿着最陡的斜坡向反应物深谷和产物深谷移动，移动路径与等能量面正交，方向与势能梯度相反。这条路径的任何地方，其垂直方向有一极小点，就是垂直切割最速下降路径的交点。（ 2 ）内禀反应坐标：内禀反应坐标就是连接反应物、过渡态和产物的最速下降路径。（ 3 ）最小能量路径（ MEP ）：通常选定一个反应坐标，使能量相对于其它坐标为最小。这样的反应路径，常常会引起弯曲甚至不连续，特别地与最速下降路径偏离甚远，而后者是真正的最小能量路径。举例来说， Relaxed Scan 的结果就是一条 MEP 。（ 4 ）经典轨迹：上述路径在势能面上对应着一组原子无限缓慢地移动（零动量），而轨迹是通过在包括核最初动量的合适初始条件下求解反应物的运动方程得到的。动力学效应，例如离心力，将使分子沿着不同于最速下降的路径移动。随着分子象台球一样从势能山脊上弹起又滑下山坡，轨迹可能十分繁长。 irc 输入和结果 %Chk=irc #p hf/sto-3g irc(forward,calcfc) guess=read geom=allcheck ****************************************** %nproc=2 %Chk=irc2 #p hf/sto-3g irc(reverse,calcfc) guess=read geom=allcheck scf(maxcyc=200) (results:) SUMMARY OF REACTION PATH FOLLOWING: 与 scan 看结果一样是看 OPti 上的那个 orient 。如果是单独作 irc(forward) 和 irc(reverse) ，取出最后 opt 结果上面的那个 orient 来看 irc 到底是个什么样的过程。第一个是过渡态的 orient 了。关于虚频的一个简单的理解首先，什么是频率。中学的时候我们学过简谐振动，对应的回复力是 f=-kx, 对应的能量曲线，是一个开口向上的二次函数 E=kx^2/2. 这样的振动，对应的 x=0 的点是能量极小值点（简单情况下也就是最小值点）。这时的振动频率我们也会求： ω=2π sqrt(k/m) 。显然它是一个正的频率，也就是通常意义下的振动频率。那么，一维情况下，如果能量曲线是一个开口向下的二次曲线呢？首先，从能量上看，这是个不稳定的点，中学的物理书上称为 “ 不稳平衡 ” 。用现在的观点看，就是这一点导数是零（受力为 0 ），且是能量极大值。如果套用上面的公式， “ 回复力 ”f=-k'x （实际上已经不是回复，而是让 x 越来越远了），这里 k' 是个负数， ω=2π sqrt(k'/m) 显然就是一个虚数了，即所谓的虚频。 Gaussian 里面给出一个负的频率，就是对应这个虚频的。实际情况下，分子的能量是一个高维的势能面，构型优化的时候，有时得到了极小值点，这样这个点的任意方向上，都可以近似为开口向上的二次函数，这样这里对应的振动频率就都是正的。对于极大值点，在每个方向都是开口向下的二次函数，那么频率就会都是负的 —— 当然一般优化很少会遇到这样的情况。对于频率有正有负的情况，说明找到的点在某些方向上是极大值，有些方向上是极小值。如果要得到稳定的能量最低构型，显然需要通过微调分子的构型，消去所有的虚频。如何微调？要看虚频的振动方向。想象着虚频对应的就是开口向下的二次函数，显然，把分子坐标按照振动的方向移动一点点，分子应该就可以顺着势能面找到新的稳定点，但是也不能太小。而所谓的过渡态，则是连接反应物和产物之间的最低能量路径上的能量极大值。好比山谷中的 A ， B 两点，它们之间的一个小土丘，就是过渡态，从 A 到 B 的反应，需要越过的是这个小土丘，而不是两边的高山。这样，过渡态就是在一个方向上是极大值，而在其它方向上都是极小值的点。因此，过渡态只有一个虚频。频率分析只能在势能面的稳定点进行，这样，频率分析就必须在已经优化好的结构上进行。频率分析的另外一个用处是判断稳定点的本质。稳定点表述的是在势能面上力为零的点，它即可能是极小值，也可能是鞍点。极小值在势能面的各个方向都是极小的。而鞍点则是在某些方向上是极小的，但在某一个方向上是极大的，因为鞍点是连接两个极小值的点。有关鞍点的信息 : 1. 负的频率 2. 频率相应简正振动的模式当一个结构产生负的振动频率时，可以表明在该振动方向可能存在着能量更低的结构。判断所得鞍点是不是需要的鞍点的方法，就是察看它的简正振动模式，分析是不是可以导向所需要的产物或反应物。进一步的，更好的办法是通过 IRC 计算来判断反应物，产物与得到的鞍点是否有关系。向大家推荐一本书，翻译成中文的（原文我没有见过），可能是《化学反应中的电子》。在寻找过渡态计算中 opt=ts, 需要计算频率，如果有一个虚频就说明优化所得的构型是个过渡态。有时在做部分限制优化 opt=minimal 时，也能得到一个虚频（一般发生在 frozen 部分的化学键上）。我想问一下，虚频的数值大小本身具有什么物理意义。后一种情况是否也可以认定为过渡态。我的理解：虚频就是在这个振动方向上力常数是负的过渡态之所以有一个虚频，是因为它是一个鞍点，鞍点上有一个负的梯度的方向分子沿这个方向振动时将转化为反应物或产物。就象从山上掉下来，受到的力可以认为是负的。但并不是有一个虚频就是过渡态，还要在这个虚频的振动方向上分别指向反应物和产物 Gaussian 程序频率计算中的几个问题在 Gaussian 计算中，为了确定优化得到的几何结构是势能面上的局域极小点还是鞍点，或者要得到相关的热力学性质，经常需要对优化后的几何结构进行振动分析。这里我们将讨论几个频率计算中常见的一些问题。希望能对初学 Gaussian 的人有所帮助。首先，原则上说，振动频率分析只对稳定结构有意义。这里所说的稳定结构包括是势能面上的局域极小点和鞍点。如下图 1 所示是一维自由度上的势能面， A 和 B 处在势能面的局域极小点，而处在势能面的鞍点上。他们在都处在平衡位置（原子核受力为零），不同的是， A 和 B 来说离开平衡位置会受到指向平衡位置处的力，而 C 离开平衡位置会受到远离平衡位置的力。因此 A 和 B 处在稳定平衡点， C 处在不稳定平衡点。实际上，一个分子可以有很多的自由度，如果在所有自由度上分子都处在稳定平衡，就是稳定的分子。频率分析得结果是所有频率都是正的，表明这是一个局域的极小点。如果分子只在一个自由度上处于不稳定平衡位置，其他自由度上都处在稳定平衡位置，说明该结构是一阶鞍点。分子在稳定自由度方向上的振动才是真实的振动，在不稳定自由度方向上的实际上是不会有振动的。不过我们可以对不稳定方向上的运动也按振动来做数学处理，会的到负的振动频率，我们称它为虚频。虚频的出现表明该结构为鞍点。图 1 势能面上的局域极小点和鞍点第二， Gaussian 计算中，频率的计算一定要在和分子结构优化相同的方法，基组下进行，否则计算的结果是没有意义的。我们知道，任何理论水平下的计算，都是在一定的近似下进行的，不同的理论水平的近似程度是不同的。在一种理论水平 A 下优化的稳定结构 Geom_A 会和另一种理论水平 B 下优化的稳定结构 Geom_B 有差别，也就是说 Geom_A 不会是理论水平 B 下的稳定结构。根据前面我们所讨论的，在理论水平 B 下对一个不稳定的结构进行频率分析是没有意义的。图 2 示意说明了不同理论水平下稳定点结构的不同。图 2 不同理论水平下优化的稳定结构是不同的第三，频率计算中可以考虑同位素效应（ Freq=ReadIsotopes ）。在波恩 - 奥本海默近似下，对于同一种元素采用不同的同位素对几何优化和电子结构计算没有影响，频率计算所需的力常数矩阵（ Hessian 矩阵）也不会变化，变化的只是约化质量。容易理解，重的同位素会导致低的振动频率。实际上，原子序数大的元素的同位素效应非常不明显，一般只需考虑 H 原子的同位素效应。第四，各种方法计算的频率和实验结果之间存在系统误差，需要乘以一个约化因子来进行校正（ Scale=f ）。一般来说，理论计算的频率值会比实验结果大。下面是一些理论水平下的约化因子。注意频率和零点能的约化因子是可以不同的。更多水平下的约化因子需要查文献获得。方法：约化因子约化因子（频率）（ ZPE ） HF/3-21G 0.9085 0.9409 HF/6-31G(d) 0.8929 0.9135 MP2(Full)/6-31G(d) 0.9427 0.9646 MP2(FC)/6-31G(d) 0.9434 0.9676 BLYP/6-31G(d) 0.9940 1.0119 B3LYP/6-31G(d) 0.9613 0.9804 SVWN/6-31G(d) 0.9833 1.0079 第五， Gaussian 的频率计算有时会遇到下面的警告： Warning -- explicit consideration of 3 degrees of freedom as vibrations may cause significant error 这时一般有两种可能：一种可能是，优化的几何结构不够精确，还没有达到稳定点。对于这种情况，需要考虑用 OPT=tight 或 OPT=VeryTight ，结合 Int=fine 或者 Int=VeryFine 进行更加精确的优化。第二种可能是，在计算的振动模式中，有部分的低频模式对应于内转动模式，在热力学分析中应该按自由转子或受阻转子模型处理，如果按谐振子模型处理会造成较大的误差。采用受阻内转子模型，可以用 Freq=Hindered-Rotor 计算或者自己手动计算这些振动模式对热力学性质的贡献。第六，对于计算出的力常数较小的振动模式，其势能面形状可能不满足谐振子模型的要求，力常数函数需要考虑非谐性的成分，这时可以用 Freq=Anharmonic 来进行非谐性处理。第七，在频率输出前，会有这样的输出 : Low frequencies --- -0.0008 0.0003 0.0013 40.6275 59.3808 66.4408 Low frequencies --- 1799.1892 3809.4604 3943.3536 上面行的六个模式实际对应于平动（前三个）和转动（后三个），理论上，这六个数都应该是零。如果是对过渡态进行的频率计算，第一行会先输出虚频，然后才给出平动转动的数值。一般来说，平动的三个数都是接近于零的。如果转动的三个数不接近于零（一般解析频率 10 个波数以内，数值频率 50 个波数以内），说明需要进行 OPT=Tight 或 OPT=VeryTight 计算。第二行的几个实频率应该和随后输出中的相应频率进行对比，如果基本一样，说明频率计算结果很好；如果差别较大，说明这些频率受到的平动和转动的污染较大。第八，频率计算后程序会进行一步优化计算，正常会得到四个判据都是 YES 的结果，这说明优化的结构很好。但是有时也会遇到四个判据不全是 YES 的情况，这时需要仔细分析，一般有三种情况： 1 、 Force 和 RMS Force 都是收敛的， Displacement 和 RMS Displacement 虽然是 NO ，但是都比较接近收敛判据； 2 、如果 Force 和 RMS Force 都是收敛的， Displacement 和 RMS Displacement 不收敛且数值远大于收敛判据； 3 、 Force 和 RMS Force 不收敛。对于第一种情况，我们可以不管它，仍然可以认为计算结果是可靠的。对于第二和第三种情况，一般说明优化过程中估算的 Hessian 是不准确的，优化的结构可能还没有达到稳定点，需要重新优化。如果反复优化仍然无法解决这个问题，建议在优化过程中使用 OPT=CalcAll 。第九，优化一个局域极小点，收敛后从力常数本征值看应该是局域极小点（没有负的本征值），但是频率计算中会出现虚频和负的本征值。这种情况下虚频一般是由转动模式造成的，说明分子中两个基团之间的相互位置不是很合适，需要绕转动的键相对转动到一个合适的位置，重新优化。绕某个键转动两个基团，有时可以很方便地用修改二面角的方法实现： OPT=Modredundant 结合分子描述后输入： * m n * value 。其中 m ， n 为连接两个基团的键的顶端原子。有时，用 OPT=CalcAll 也可以解决这个问题。第十， Gaussian 程序中，频率计算是和温度、压力无关的。因为 Gaussian 所考虑的是原子核在一定核构型和电子运动状态构成的势场内振动，一般来说，温度和压力的变化不会影响分子的结构和电子运动状态，所以振动的势能函数是与温度、压力无关的。 Freq=Isotopes 关键词要求输入的温度和压力影响的只是平动，转动和振动的平衡分布和它们对热力学性质的贡献。分子激发态计算的理论方法一 1. 激发态计算方法的分类按照化学模型，目前计算激发态的方法可以分为： ____DFT (TDDFT, MC-DFT ， DFT+CI...) Chemical____| Model |____Wave function____|~~~~Ab initio (MR-CI, CIS...) |____Semiempirical (Zindo/1 ， INDO/S ， CNDO/S...) 这与基态计算方法的分类是一样的，从中看不出激发态计算的特点。另一种分类是按照获得激发能（或激发态能量）的方式，可以分为两类： i) 电子组态方法。选择不同的电子组态的线性组合进行迭代，得到激发态的绝对能量（有些程序经过简单的后期处理也能给出相对能量）。分为多参考 (MR) 方法和单参考 (SR) 方法两种。常见的多参考方法有： MCSCF( 有 CASSCF ， RASSCF 等多种形式 ) ， MR-CISD ， CASPT2 等。标准的 MR 的计算过程一般是三步： HF— 〉 MCSCF （常用 CASSCF ） — 〉 MRCISD 或 CASPT2...R 是 MR 的特例，不需要中间的 MCSCF 步骤，计算中只需要对 MR 程序模块指定一个参考组态即可，计算量大大减少，但是精度也会下降。对于一些无法用 SR 很好描述的体系（如镧系化合物），计算是失败的，这时就必须用 MR 方法。 ii) 线性响应理论 (LRT) 。基态用普通的理论方法求得，然后通过求解特殊的线性响应方程，获得激发态相对于基态的垂直激发能（相对能量）。有的方法求解方程需要做迭代，有的则不需要。这一方法要求作为参考态的基态必须被很好地描述。常见的有 TDDFT ， CIS ， CCSD-LRT 等，共同点是基态都使用 SR 模型。基态用 MR 模型的有 MCLR (MCSCF-LRT) ， MR-AQCC-LRT 等方法，用于计算那些无法用 SR 模型很好描述基态的体系。 2. 一些方法的特点及软件 i) Zindo/1 ， INDO/S ， CNDO/S 都属于半经验方法，可用于研究原子数在 100 左右的大分子体系。软件： Zindo/1 Gaussian, Zindo 。 Zindo 结合 ZINDO-MN 模块可以用 Zindo/1 方法计算溶剂对激发态的影响。 INDO/S ， CNDO/S 商业 MOPAC 的 MOS-F 模块。免费的 MOPAC 6.0/7.0 无此功能 ii) CIS 基态能量用 Hartree-Fock 方法得到，激发态精度也是 HF 级别的，得到半定量的结果。适于研究原子数在 50 左右的中等分子体系。此外还有精度更高一些的 CIS(D) ，计算量当然更大。 Gaussian 03 的 CIS 还能考虑溶剂环境对激发态影响。程序化的这些方法一般只能计算单重态和三重态。扩展 CIS (XCIS) 方法还能计算双重态和四重态。软件： CIS Gaussian, Q-Chem ， CaChe ，商业的 MOPAC ， TURBOMOLE ，最新版 MWChemCIS(D) Q-Chem ， TURBOMOLEXCIS Q-Chem iii)TD-DFT 和 TD-HF 目前比较流行的激发态计算方法，适用于十几个原子构成的中小分子，小分子计算也能得到令人满意的结果。 TD-DFT 目前的缺陷是无法研究多电子激发和 Rydberg 态。 TD-HF 早期也叫 RPA ，是 TD-DFT 的特殊形式，精度较差。 Gaussian 03 中 TDHF/TDDFT 的激发态计算还能考虑溶剂环境影响。程序化的 TD-DFT 一般只提供对单重态和三重态计算的正式支持。最新版 NMChem 的 TD-DFT 除了计算闭壳层分子的单重态和三重态，还能计算开壳层分子的双重态和四重态。软件： TD-HF Gaussian, ADF, NWChem, Q-Chem ， DeMon ， GAMESS-US ， ACES II ， TURBOMOLE TD-DFT Gaussian, ADF, NWChem, Q-Chem 。 DeMon ， TURBOMOLE ， CADPAC ， BDF 。一些计算晶体的程序也可以做分子的 TD-DFT 计算，但是由于使用平面波基组，结果不太好。 iv) CC-LRT CC-LRT 包括 CCS-LRT ( 等价于 CIS) ， CC2-LRT ( 对 CCSD-LRT 的近似 ) ， CCSD- LRT ， CC3-LRT ( 对 CCSDT-LRT 的近似 ) 等不同激发级别的形式。最常用的是 CCSD-LRT ，有很多种形式，如 EOM-CCSD ， VOOD-CCSD 等。适用于中小分子和小分子的计算，结果与 TD-DFT 相当，可能还略好一些。软件： CC-LRT DALTON 。 CCSD-LRT 的其他形式可在 ACES II ， Q-Chem ， PSI ， MOLPRO 中找到。 TURBOMOLE 可做 CC2-LRT 计算。 v) SAC-CI 。精度略高于 TD-DFT 的方法，适于研究十几个原子构成的有机分子，还可以包含多电子激发的影响。软件： SAC-CI Gaussian 03 vi) DFT+CI 混合方法。结合了 DFT 和 CI 的优点，适用于中小分子激发态计算。软件： DFT+CI TURBOMOLE vii) MCSCF ( 有 CASSCF ， RASSCF 等 ) 。做（电子组态方法）多参考计算的基础。适用于中小分子激发态计算。常用电子组态形式的 CASSCF 。线性响应理论形式 (MCLR) 用的较少。软件： MCSCF GAMESS-US, GAMESS-UK, Gaussian, DALTON, MOLPRO, MOLCAS, COLUMBUS ， NWChem MCLR DALTON ， GAMESS-UK viii)MR-CI 及其变体 MR-CISD 是计算（几个原子构成的）小分子激发态的标准方法。由于这一方法存在着大小一致性问题，因而又发展了很多变体。比较简单的是对 MR-CISD 做事后修正 (+Q) ，主要有 Davidson 修正和 Pople 修正两种。复杂一些的有 MR-ACPF 、 MR-AQCC 、 MR-CEPA 等，其中表现最好的是 MR-AQCC ，以及后来发展的虚轨道修正 MR-AQCC (MR-AQCC-V) 。有的程序还允许对 MR-CI 做激发级别的修正，得到 MR-CIS ， MR-CISD(T) 等。软件： MR-CISD GAMESS-US, GAMESS-UK, DALTON, MOLPRO, MOLCAS, DIRAC ， MOLFDIR ， COLUMBUS 变体： MOLPRO, MOLCAS ， COLUMBUS 。 COLUMBUS 还能做 MR-AQCC-LRT 计算，此方法综合了 MR-AQCC 与线性响应理论的特点。 MOLPRO 对 MR-CI 及其变体使用了内部收缩，速度更快。 ix) MR-MPn 多参考微扰理论。目前主要有二阶微扰和三阶微扰。在程序作了优化的前提下，它们比相同级别的 MR-CI 计算速度快（例如在 MOLPRO 中， CASPT2 比 MR-CISD 快），但是精度略低。目前应用的难题是入侵态问题。有些程序通过与 MR-CI 混合可以部分地解决这一问题。软件： CASPT2 MOLPRO, MOLCAS ， COLUMBUS ， GAMESS-UK 。 GAMESS-US 的 MC-QDPT 方法和 Gaussian 的 CASSCF+MP2 也属于 MR-MP2 。 CASPT3 MOLPRO ， GAMESS-UK x) MR-CC 和 FCI MR-CI 的理论极限是完全组态相互作用 (FCI) ，常用的 MR-CISD 是把 FCI 在 CI 空间的展开截断到二级的形式。鉴于耦合簇 CC 能获得比 CI 更好的结果，而且不会有大小一致性的问题，如果把 FCI 展开到 CC 空间并在一定的级别截断，就得到 MR-CC 。由于 FCI 计算量相当大，所以目前也只能计算较轻原子构成的双原子分子。从理论上讲， MR-CCSD 是小分子激发态计算最实用、最精确的方法。但目前还有入侵态的问题没有解决。公开发布的程序中目前还没有此方法。 xi) MC-DFT ， CAS-DFT ：多参考的 DFT ，尚在研究中。 MOLCAS 声称将在 6.0 版包含此方法。 xii) MR-Gn 多参考的 Gaussian-n 理论，可想而知是相当昂贵的方法。目前只有 GAMESS- US 小组的一篇文章研究过这一理论。将来能否用于激发态计算并获得高精度结果还未知。 3. 方法的选择构成分子的原子数： 2.....10.......20........50..........100.......... |_____________________| |______________| 半经验 CIS |_____________________| MCSCF( 注 1) |______________| TD-DFT( 注 2), CCSD-LRT |_____| MR-CI 及其变体， MR-MPn( 注 2) | FCI 可见，越往下精度越高，适用分子越小。注 1 ： CASSCF 的计算量和活性空间大小 2A 和活性电子数 E 有关。 C(2A,E)=(2A)!/ 。 C(2A,E) 越大则计算量越大。 RASSCF 由于对组态进行了筛选，能适当减少计算量。对于小分子， MCSCF 一般只作为高级 MR 计算（ MR-CI ， MR-MPn ）的中间步骤。注 2 ：对于含有重元素的小分子激发态计算，还要考虑自旋轨道耦合 (SO) 造成的电子态分裂。例如 AtBr 的 3Pi 会分裂为 0- ， 0+ ， 1 ， 2 四个态。能够在激发态计算中包含 SO 的程序有 MR-CISD 及其变体 + SO: COLUMBUS ， MOLPRO ， GAMESS-US MR-MPn + SO: MOLPRO ， GAMESS-US MCSCF + SO: COLUMBUS ， MOLPRO ， GAMESS-US, MOLCAS ， TD-DFT + SO: ADF 2002 ，尚在测试中，目前只能对原子获得较好的激发能完全相对论 MR-CISD: DIRAC, MOLFDIR 完全相对论 TD-DFT: BDF ，尚在测试中如果是用 Gaussian 或者 GAMESS-US ，我的理解是： 1. 先用 AM1 进行几何优化。 2. 用优化的结构进行 Hartree-Fock 计算。 3. 用第二步的分子轨道进行 CIS 计算（ Gaussian 使用 Guess=Read 读入第二步的 checkpoint 文件； GAMESS 使用 guess=moread ，并把第二步的 MO 复制到 $VEC 中）。 Gaussian 的 CIS 自动给出激发能和谐振强度， GAMESS 只给出激发能，计算谐振强度需要另外的计算步骤。 4.Gaussian 中的 Zindo 和 CIS 使用方法相同，速度更快，但适用的原子有限，如果计算量不是太大的话，还可以用 TD ，结果更好。另外 Gaussian 的 CIS 只能算单态、三态，其它的像双态、四重态、五重态可能也会给出结果，但是分析很麻烦。用 GAMESS 进行 CIS 计算可以定义任意多重度，比较直观。 http://emuch.net/html/200805/829151.html http://blog.sina.com.cn/s/blog_763deefd010174v6.html; 个人分类: 研究兴趣|23830 次阅读|0 个评论

Gaussian fit: hustfliee 2013-11-21 19:39; x_data = '; y_data = 3*exp(-(x_data-2).^2/4) + 0.5*rand(size(x_data))-0.25; cf1 = fit(x_data,y_data,'gauss1') % natural fit fff= fittype('a1*exp((-(x-b1).^2)./c1)','problem','c1') c1 = 2.5 cf2 = fit(x_data,y_data,fff,'problem',c1) plot(x_data,y_data,'bo') hold all plot(cf1,'r') plot(cf2,'g') legend('Data points','Gaussian fit','Gaussian fit with fixed c1'); 2286 次阅读|0 个评论

[转载]GAUSSIAN处理激发态: DonarF1 2013-1-23 09:41; GAUSSIAN处理激发态 http://scuchj.blog.sohu.com/112134965.html 分类：资源 2009-03-12 20:21 §3-1.输入文件中的电子态多重度输入文件中的多重度，指的是具有这种多重度的最低能量的电子态，不一定就是基态，不过一般输入文件都用基态。激发态的多重度是在输入文件的route section行控制的，如：singlets，triplets。虽然在GAUSSIAN使用帮助中没有明说，但是从给出的例子中可以体会出来。另外，在HyperChem的帮助文件中有类似的说明，虽然是针对HyperChem的，但我想对GAUSSIAN一样适用。 §3-2.GAUSSIAN中研究激发态的方法 GAUSSIAN可以用CIS，CASSCF和ZINDO方法计算激发态。GAUSSIAN中的CIS实际是对分子激发态的零级组态相关处理。有些从头计算程序（如GAMESS）包含计算激发态的一级组态相关（FOCI）和二级组态相关（SOCI），我不清楚是不是和CIS相对应。从结果来看，FOCI和 SOCI比CIS略有改善。ZINDO方法适用范围窄，只能用于周期表前48种原子构成的分子，这里不谈。98版还新增加了含时（Time Depend）方法处理激发态（包括含时Hatree-Fock(TDHF)和含时密度泛函(TDDFT)）。我们一般计算的激发态称为Low-lying Excited States，也就是由原子价电子构成的分子电子激发态。更高的电子态是Rydberg态。从头算法和其他半经验方法对空轨道，特别是较高空轨道的计算结果不好，因而得到的Rydberg态的结果缺乏明确的物理意义。 §3-3.含时方法（TD）输入方式帮助文件对TDDFT这么重要的方法说得真是太简单了，让人去参考CIS的输入方法，连个例子都没有。CIS的输入方法是： # CIS(Triplets,NStates=10)/3-21G Test 那就用B3LYP泛函试试看吧： # TDB3LYP(Triplets,NStates=10)/3-21G Test 但是不行。其实正确的输入方法是： # B3LYP/3-21G TD(Triplets,NStates=10)Test 这个问题浪费了我很长时间。为了表示它的重要，单独分为一节。 §3-4.CIS，TDHF，TDDFT方法激发能比较一般来说，CIS和TDHF方法结果的准确度是比较差的，TDDFT的结果是最好的。而各种TDDFT相比，B3LYP，B3P86，MPW1PW91是最好的（依具体的分子和使用的基组而不同）。可以参考文献： J. Chem. Phys., 111(7), 1999, 2889 Chem. Phys. Lett., 297(1-2), 1998, 60 J. Chem. Phys., 109(23), 1998, 10180 这些文章都是计算C，H，O组成的有机物在基态平衡位置的垂直激发能。我对含有重元素In的分子激发态进行计算，发现如果选好基组，很多TDDFT的结果和实验值的差距甚至小于500个波数，而CIS的结果则可能要差5000多个波数。所以现在使用TDDFT方法的比较多。另外还有人用TDDFT方法计算 DNA这样的大分子。 §3-5.GAUSSIAN处理激发态的BUG 1。用94版进行CIS计算，50-50选项不支持LanL2DZ和3-21G*基组。98版则不会出现这个问题。这个问题可能还跟计算的分子有关。 2。同时计算单重激发态和三重激发态（也就是使用50-50这个选项），有时候会出现致命错误：某些电子态在某些位置的能量有不规则起伏，而单独计算单重态或三重态（使用Singlets或者Triplets选项）则不会出现这个问题。我在计算中曾经遇见过两次。下面这个例子是从网上看到的，我想这位德国的大哥也遇见了同样的问题。他对同一个任务执行了三种计算： 1) #Bp86/6-31G* td=(nstates=5,50-50) scf=direct density 2) #Bp86/6-31G* td=(triplets,nstates=1) scf=direct density --link1-- #Bp86/6-31G* td=(nstates=10) scf=direct guess=read geom=check --link1-- #Bp86/6-31G* td=(triplets,nstates=10) scf=direct guess=read geom=check 3) #Bp86/6-31G* td=(nstates=10) scf=direct #Bp86/6-31G* td=(triplets,nstates=10) scf=direct density 第一种计算产生10个态（5S和5T），后两种产生20个。关键字Density对结果没有影响。结果是： 1) 2) 3) Triplet-B1U 3.1077 3.1077 3.1077 Triplet-B3G 3.9568 3.9568 3.9568 Triplet-AG 4.1612 4.1612 4.1612 Triplet-B2U 4.1769 4.1769 4.1769 Triplet-B3G 4.2101 4.2101 4.2101 Singlet-B3G 4.3609 4.3609 4.3609 Singlet-B1U 4.3788 4.4864 4.4864 Singlet-B2U 4.4864 4.6409 4.6409 Triplet-AG 4.5752 4.7144 4.7144 Singlet-B3G 4.6409 -- 5.0010 通过比较发现，对前六个态，三种处理的结果相同，但是后面的态，使用了50-50的方法就和另外两种方法的结果不一样了。我个人认为，尽量不要使用50-50这个选项。推荐单重态、三重态分开计算，虽谈多花一倍时间，但是结果保险。 §3-6.画激发态的势能曲线最笨的方法莫过于算完一个点后，更改输入文件的坐标，计算下一个点了。这样费时费力。在优化基态几何构型的时候，可以使用Scan关键字进行扫描，实际上 Scan对激发态一样适用。这样开上一晚上计算机，到早上所有激发态的势能曲线数据就全有了。例如计算从R=0.6埃到1.2埃氢分子的10个单重激发态： # RCIS(NStates=10)/6-311G Scan H2 Excited States 0,1 H1 H2 H1 r r 0.6 60 0.01 需要注意三点： 1。结果给的是垂直激发能，必须换算成相对于基态平衡位置的能量。 2。扫面范围不要离基态平衡位置太远，一般从0.75Re到2Re就可以了，太远了会报错。 3。某些理论和某些基组的搭配可能不支持Scan，这时候就只好手工控制了。最后是数据处理，这一步是最枯燥的了。数据处理需要用到Excel、Oringin、Matlab。这些软件用不着钻研很深，知道几个最基本的命令就足够用了。 §3-7.高对称性的分子使用CIS方法默认的CIS计算求出最低的三个激发态。以它为例，GAUSSIAN在初始的猜测中选择电子态数目两倍的向量（也就是6个）进行迭代，直到有三个收敛。对于大多数分子体系来说，这种默认的向量数目足够了。但是对于高对称性的分子，则需要特殊的方案，否则不会得到比较高的收敛态。方法是增加NStates的数量，以增加初始向量的数量。推荐NStates值为最大阿贝尔点群操作的数量，也就是输出文件在Standardorientation之前的 Symmetry部分（在输入文件中使用#控制打印输出就可以得到，#T不行）。这是一个计算苯分子的例子： Stoichiometry C6H6 Framework group D6H Deg. of freedom 2 Full point group D2H NOp 8 Largest concise Abelian subgroup D2 NOp 4 Standard orientation ------------------------------ 因此在route section中使用CIS(NStates=8)。（四）求分子的光谱常数 §4-1.求离解能首先必须注意，原子化能和离解能不是同一个概念。帮助文件和Exploring Chemistry with Electronic Structure Methods一书中的解释并不确切。计算原子化能和离解能是不同的两个过程。以三原子分子GaCl2为例，计算离解能的过程是： GaCl2(X2A1)----Cl(2P)+GaCl(X1Σ+) 过程中的能量变化称为离解能。很明显，离解能还可以细分为第一步离解能、第二步离解能......，上面的是计算第一步离解能，计算第二步离解能的过程是： GaCl(X1Σ+)----Cl(2P)+Ga(2P) 而计算原子化能的过程是： GaCl2(X2A1)----Cl(2P)+2Ga(2P) 可见还是有区别的。不过对于双原子分子来说，这两个概念是相同的。第二个需要注意的是，离解能一般用符号D0或者De表示（前者不包括零点能，后者包括），但是在分子光谱中，还有一组D0和De。为避免混淆，这里需要作特别说明。第一组D0和De，也就是离解能，数值比较大，数量级一般在0.1到1eV。这个参数多出现在真空-紫外、振动光谱中。有一个公式： ωe^2 ωeχe = ------- （^是指数符号） 4De 其中的De就是离解能。另一组D0和De称为离心畸变常数，数值很小，作为转动光谱中的微扰项。如： h^3 De = ----------------------------------- 128(π^6)(μ^3)(c^3)(re^6)(ωe^2) 其中的De就是离心畸变常数。 GAUSSIAN中有个求PH2分子的原子化能的例子，可以作为求离解能的参考。也就是文件e7-01。结果是原子化能 = E(P)(Hartree)+2E(H)(Hartree)-E(PH2)(Hartree)-ZPE(Hartree) = (-341.25930)+(-0.50027)*2-(-342.50942)-(0.01322)(Hartree) 其中的零点能ZPE经过了矫正 = 148.3 Kcal/mol 和实验值144.7 Kcal/mol的差距相当小了。 §4-2.求电离能电离能缩写为I.P.。在GAUSSIAN中有两种求I.P.的方法。法一：见文件e7-03，这个例子求PH2分子的电离能。结果为： E ZPE PH2+ -342.14416 0.01347 （ZPE经过矫正） PH2 -342.50942 0.01322 （这是上面求原子化能的结果） I.P. = (E(PH2+)+ZPE(PH2+))-(E(PH2)+ZPE(PH2)) = (-342.14416+0.01347)-(-342.50942+0.01322)(Hartree) = 0.36551(Hartree) = 9.95(eV) 法二：OVGF方法。帮助文件里有详细的使用说明。 §4-3.求激发态的光谱学参数Te，ωe，Re 有两种方法。 1。CIS可以用结构优化的方法。Starl的文章里说得很详细了，这里不再细讲。仅仅需要补充的是，如何强制保留设置后的波函：在使用Guess=Alter设置波函之后，程序总会回到设置前的波函。这个时候需要使用SCF关键字中的Symm选项，强制程序保留设置以后的波函。使用方式是在route section部分添加SCF=Symm。 2。势能曲线扫描法。 TDDFT方法不能计算激发态的波函，所以不能用来研究激发态的平衡构型和分子属性。上面的结构优化方法用在TDDFT中就不适合了。但是可以用势能曲线扫描的方法。但是，多原子分子能量是势能曲面，控制多个输入坐标变量，扫描整个曲面是困难的，所以这个方法用于双原子分子。如何扫描激发态的势能曲线见前面的§3-6节。有了势能曲线以后，就可以找到势能最低点Re（为了精确，这一区域的点要取密一些），Re位置的能量相对于基态平衡位置能量的差值就是电子态高度Te。势阱的深度就是离解能De。这些都简单，关键是求振动参数ωe。求ωe也有两种方法。一种是用Morse势函数拟合的方法。Morse势函数形式为： V(R) = Te+De{1-EXP }^2 （对基态Te是0） Te、De、Re是已知的，拟合出了α后，用公式 α = ωe ((2(π^2)c*μ)/(h*De))^0.5 （高斯单位制）就可以得到ωe。这种方法太麻烦，我们不这样做。实际只要把扫描出来的Morse势函数对R求二阶导数，在Re处的导数值f2和ωe的关系有： 4(π^2)(ωe^2)c*μ f2 = -------------------- h*(10^18) 其中c、h、μ是国际单位制，ωe的单位是波数cm^-1。另外在扫描的势函数中，De的单位是波数cm^-1，R和Re的单位是埃。否则上面的f2需要乘上个系数。求f2用Origin软件就可以做。注意要取Re位置的f2。有了f2，就可以算出ωe，省去了计算De和拟合α的步骤，避免引入更大的误差。 §4-4.如何获得光谱项的对称性如果用GAUSSIAN计算比较轻的分子，一般都能给出正确的激发态光谱项对称性，但是对含有重原子的分子，可能得不到对称性。这个时候就需要根据轨道跃迁的的信息获得对称性。下面是用CIS方法LanL2DZ基组计算InCl分子的例子，输出文件给出InCl的电子轨道信息： Occupied (SG) (SG) (PI) (PI) (SG) Virtual (PI) (PI) (SG) (PI) (PI) (SG) (PI) (PI) (SG) (SG) (SG) 这里结果只给出前4个单重激发态： CIS wavefunction symmetry could not be determined. Excited State 1: Singlet-?Sym 4.9226 eV 251.87 nm f=0.3243 5 - 6 .43587 5 - 7 .54909 This state for optimization and/or second-order correction. Copying the Cisingles density for this state as the 1-particle RhoCI density. CIS wavefunction symmetry could not be determined. Excited State 2: Singlet-?Sym 4.9226 eV 251.87 nm f=0.3243 5 - 6 .54909 5 - 7 -.43587 Excited State 3: Singlet-SG 6.7609 eV 183.38 nm f=0.2830 3 - 6 -.28038 4 - 7 -.28038 5 - 8 .57818 CIS wavefunction symmetry could not be determined. Excited State 4: Singlet-?Sym 6.8241 eV 181.68 nm f=0.0000 3 - 7 .49690 4 - 6 -.49690 由于In原子是重原子，所以大多数电子态并没有给出对称性（第三个激发态除外）。我们发现第一激发态是从第五个轨道跃迁到第六、七个轨道，而第六、七两个轨道是简并的（PI、DELT轨道都是二重简并的）。也就是说第一激发态的电子组态是： InCl: (内层电子)(1σ)2 (2σ)2 (1π)4 (3σ)1 (2π)1 这个组态耦合出来的单重态光谱项有一种（当然还有一个三重的）：1Π，所以就是1Π态。再往下看，第二激发态和第一激发态的电子组态、能量完全相同，说明这是一个简并态，而实际上除了Σ态以外的态都是二度简并的。所以第一激发态和第二激发态实际都是InCl的第一激发态1Π态。第三个激发态给出对称性了。如果没给，方法和上面的一样，先找到它的电子组态，它是由两个组态叠加而成的： InCl: (内层电子)(1σ)2 (2σ)2 (1π)3 (3σ)2 (2π)1 InCl: (内层电子)(1σ)2 (2σ)2 (1π)4 (3σ)1 (4σ)1 第一个组态可以耦合出的单重态光谱项有：1Δ，1Σ+，1Σ-。因为这不是一个简并态，所以只能是1Σ+和1Σ-其中之一。这就需要进一步知道+、-对称性。方法是看跃迁系数： 3 - 6 -.28038 4 - 7 -.28038 由于正负号相同，所以是1Σ+态。那么第四激发态必然是1Σ-了。另外还可以看第二个组态耦合出来的光谱项，只有1Σ+，因此第三个激发态必然是1Σ+。两个结果一致。以上是自旋限制的计算。非限制计算（每个轨道再分成自旋Alpha和自旋Beta）分析类似。 *************************************************************************** Gaussian 98的补充 1. 谐振强度f 谐振强度f出现在计算激发态的结果中，如： Excited State 3: Singlet-?Sym 5.0236 eV 246.80 nm f=0.3428 5 - 7 0.68918 在Gaussian 98早期版本中，TD方法计算激发态的谐振强度f是错误的（其它计算激发态的方法没有这个问题）。从A7版开始做了修正。这点要注意。谐振强度能干什么呢？可以预测电子激发态跃迁的寿命。用公式 t=1.499/(f*E^2) 这里的E是电子态之间的能量间隔，单位是波数，f是电子态谐振强度，t是电子态寿命，单位是秒。默认的下能级是基态。如果计算各激发态之间的f，在 Route Section中加入关键字AllTransitionDensities，就得到各电子态之间的偶极矩m(x)，m(y)，m(z)。求解f用公式： f=(2/3)*E* 因为Gaussian给出的m(x)，m(y)，m(z)是au单位，所以激发态之间的能量间隔E要换算成au单位，才能代入上面的公式。 2. 用TDDFT做几何结构优化一般认为，因为TDDFT不能计算激发态的波函，所以不能用来研究激发态的平衡构型和分子属性。所以用Gaussian 98提供的TDDFT方法做研究激发态的研究工作，只能得到垂直激发能，对于计算平衡键长Re、键角和振动频率ωe无能为力。不过，最近的报道说，TDDFT做激发态几何优化在理论上证明可行。已经有消息说，Gaussian公司准备在新版本的Gaussian中加入TDDFT几何优化了。可以参考文献： Van Caillie C, Amos R. D., Geometric derivatives of density functional theory excitation energies using gradient-correct functionals, Chem. Phys. Lett., 317(1-2), 2000, p.159-164 和 Hsu, Hirata, Head-Gorden, J. Phys. Chem., A105, 2001, 451-458 关于线性多烯烃的TDDFT激发能和谐振强度。; 7532 次阅读|0 个评论

[转载]gaussian 出错信息和可能的解决办法(转）: jlcuit 2012-11-3 16:03; from http://blog.sina.com.cn/s/blog_8be57be70100xabz.html Gaussian calculations can fail with various error messages. Some error messages from .out and .log files - and possible solutions - have been compiled here to facilitate problem solving. These are divided into: Syntax and similar errors Memory and similar errors a title="Convergence problems: " href="file:///E:/–°–‡¤¹/Trouble%20shooting%20Gaussian%20calculations%20€”%20Norwegian%20High%20Performance%20Computing.htm#convergence-problems" Convergence problems a title="Errors in solvent calculations: " href="file:///E:/–°–‡¤¹/Trouble%20shooting%20Gaussian%20calculations%20€”%20Norwegian%20High%20Performance%20Computing.htm#errors-in-solvent-calculations" Errors in solvent calculations a title="ERROR MESSAGES IN LOGFILES " href="file:///E:/–°–‡¤¹/Trouble%20shooting%20Gaussian%20calculations%20€”%20Norwegian%20High%20Performance%20Computing.htm#error-messages-in-logfiles" Errors in log files ERROR MESSAGES IN OUTPUT FILES Syntax and similar errors:End of file in ZSymb. Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l101.exe Solution: The blank line after the coordinate section in the .inp file is missing. Unrecognized layer "X". Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l101.exe Solution: Error due to syntax error(s) in coordinate section (check carefully). If error is "^M", it is caused by DOS end-of-line characters (e.g. if coordinates were written under Windows). Remove ^M from line ends using e.g. emacs . To process .inp files from command line, use sed -i 's/^M//' File.inp (Important: command does not work if ^ M is written as characters - generate ^M on command line using ctrl-V ctrl-M ). QPERR --- A SYNTAX ERROR WAS DETECTED IN THE INPUT LINE. Solution: Check .inp carefully for syntax errors in keywords RdChkP: Unable to locate IRWF=0 Number= 522. Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l401.exe or FileIO operation on non-existent file. Error termination in NtrErr: NtrErr Called from FileIO. Solution: Operation on .chk file was specified (e.g. geom=check, opt=restart), but .chk was not found. Check that: %chk= was specifed in .inp .chk has the same name as .inp .chk is in the same directory as .inp run script transports .chk to temporary folder upon job start. Run scripts downloaded here should do this. The combination of multiplicity N and M electrons is impossible. Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l301.exe Solution: Either the charge or the multiplicity of the molecule was not specified correctly in .inp. Memory and similar errors:Out-of-memory error in routine RdGeom-1 (IEnd= 1200001 MxCore= 2500) Use %mem=N MW to provide the minimum amount of memory required to complete this step Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l101.exe or Not enough memory to run CalDSu, short by 1000000 words. Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l401.exe or allocation failure: Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l1502.exe Solution: Specify more memory in .inp ( %mem =Nmb). Possibly, also increase pvmem value in run script. Especially solvent calculations can exhibit allocation failures and explicit amounts of memory should be specified. galloc: could not allocate memory. Solution: The %mem value in .inp is higher than pvmem value in run script. Increase pvmem or decrease %mem . Probably out of disk space. Write error in NtrExt1 Solution: /scratch space is most likely full. Delete old files in temporary folder. Convergence problems:Density matrix is not changing but DIIS error= 1.32D-06 CofLast= 1.18D-02. The SCF is confused. Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/linda-exe/l502.exel Solution: Problem with DIIS. Turn it off completely, e.g. using SCF =qc, or partly by using SCF =(maxconventionalcycles=N,xqc), where N is the number of steps DIIS should be used (see SCF keyword ). Convergence criterion not met. SCF Done: E(RHF) = NNNNNNN A.U. after 129 cycles Convergence failure -- run terminated. Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/linda-exe/l502.exe Solution: One SCF cycle has a default of maximum 128 steps, and this was exceeded without convergence achieved. Possible solution: In the route section of input file, specify SCF=(MaxCycle=N), where N is the number of steps per SCF cycles. Alternatively, turn of DIIS (e.g. by SCF=qc) (see SCF keyword ). Problem with the distance matrix. Error termination via Lnk1e in /pkg/gaussian/g03/l202.exe Solution: Try to restart optimization from a different input geometry. New curvilinear step not converged. Error imposing constraints Error termination via Lnk1e in /pkg/gaussian/g03/l103.exe Solution: Problem with constrained coordinates (e.g. in OPT=modredun calculation). Try to restart optimization from a slightly different input geometry. Optimization stopped. -- Number of steps exceeded, NStep= N Error termination request processed by link 9999. Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l9999.exe Solution: Maximum number of optimization steps is twice the number of variables to be optimized. Try increasing the value by specifying OPT=( MaxCycle =N) in .inp file, where N is the number of optimization steps (see OPT keyword ). Alternatively, try to start optimization from different geometry. Errors in solvent calculations:AdVTs1: ISph= 2543 is engulfed by JSph= 2544 but Ae( 2543) is not yet zero! Error termination via Lnk1e in /global/apps/gaussian/g03.e01/g03/l301.exe Solution: Problem is related to building of the cavity in solvent calculations. One possible solution is to change the cavity model (default in g03 is UAO, can be changed by adding RADII keyword in section below coordinates in the .inp file, e.g. RADII=UFF, see SCRF keyword ). Hydrogen X has 2 bounds. Keep it explicit at all point on the potential energy surface to get meaningful results. Solution: In UAO cavity model, spheres are placed on groups of atoms, with hydrogens assigned to the heavy atom, they are bound to. If assignment fails (e.g. because heavy atom-H bond is elongated), cavity building fails. Possible solutions: a) use cavity model that also assigns spheres to hydrogens (e.g. RADII=UFF) or b) Assign a sphere explicity on problematic H atom (use SPHEREONH= N, see SCRF keyword ) ERROR MESSAGES IN LOGFILES= PBS: job killed: walltime N exceeded limit M signal number 15 received. Solution: Job did not finish within specified wall time. Retrieve .out and .chk files from temporary folder /global/work/$USER/$JOB (or $PBS_JOBID) and restart calculation if possible (using e.g. opt =restart or scf =restart). cp: cannot stat $JOB.inp: No such file or directory Solution: The .inp file is not in the directory from where the job was submitted (or its name was misspelled during submission. If error reads: cp: cannot stat $JOB. inp.inp , the .inp file was submitted with extension). ntsnet: unable to schedule the minimum N workers Solution: The value of %NprocLinda=N in the .inp file is higher than the number of nodes asked for during submission. Make sure these values match. Connection refused died without ever signing in Sign in timed out after 0 worker connections. Did not reach minimum (N), shutting down Solution: Error appears if you run parallel calculations but did not add this file to your $HOME directory: .tsnet.config containing only the line: Tsnet.Node.lindarsharg: ssh; 个人分类: gaussian|22325 次阅读|0 个评论

在Linux系统中安装Gaussian程序: NanoZhang 2012-3-3 17:28; 计算机是64位的AMD处理器，系统是64位的openSUSE，准备安装64位的Gaussian 03。EM64T是Intel的64位技术，AMD的64位技术是AMD64，它们统称为x86-64。我手头有g03e01-em64t.tar.gz，试图找“g03e01-amd64.tar.gz”却没有找到，于是只好装em64t了，后来证明一切正常，说明它们兼容。我在Linux系统中建有user1、user2等用户，目录分别为/home/user1、/home/user2等。 1. 把g03e01-em64t.tar.gz上传至/home目录下，执行“tar -zxvf g03e01-em64t.tar.gz”解压。产生/home/g03目录。 2. 建立文件夹/home/tmp做为Gaussian 03计算的scratch文件夹。 3. 设置环境变量。编辑/root/.bashrc文件，在其末尾添加如下内容：　export g03root=/home 　export GAUSS_SCRDIR=/home/tmp 　source $g03root/g03/bsd/g03.profile 安装完成，可以用root用户登录做计算了。不过user1、user2等用户登录后却不能做计算，因为各个用户目录间没有相互操作的权限，现在我就使用chmod和chgrp赋予他们权限。 4. 执行“chgrp -R users g03”和“chgrp -R users tmp”，让/home/g03、/home/tmp也加入users组。 5. 执行“chmod -R 750 g03”和“chmod -R 770 tmp”，750即rwxr-x---，770即rwxrwx---，使得同组可操作。 6. 把/root目录下的.bashrc文件拷贝到每个用户目录下，这样各个用户登录时就会自动加载环境变量。大功告成！从现在开始使用Linux下的Gaussian了，不再菜鸟！; 14955 次阅读|0 个评论

gaussian 简单计算小得: wangshu 2012-1-31 22:50; 下面的小知识点很简单，只是用得少，老是记不住，就写出来列于此, ，好记性不如烂笔头嘛 ============== 1) 0 1 这里的0是电荷，1是多重度，即2S+1=2*0+1=1 若一中性分子的电荷，多重度分别不为0，1，则可能是多输了H，或电子，可以用右键viewlabel来找，或菜单editatom list中找。 2) geom=connectivity 是指明确原子之前连接方式，如果在原子位置后标注连接方式则必须用这个命令；要删全删，要留全留。 3）opt指优化结构，如果算已有结构，如晶体的定点能，则不需此命令。 4）即使有晶体结构，也要优化结构，一是作类比较，二是验证此基组适用的可行性。 ============ 数据举例，来源，CHEM. COMMUN. , 2003, 2116–2118，CCDC：211977 %chk=CB.chk %mem=2500MB %nprocl=1 %nprocs=8 # b3lyp/3-21G geom=connectivity opt CB 0 2 O -1.76400000 11.13560000 4.01040000 O -1.00400000 13.23060000 4.01130000 C 3.24620000 12.74610000 5.82070000 C 3.31200000 14.36970000 5.86070000 C 4.62480000 14.92760000 6.26510000 C 5.67340000 14.16090000 6.72830000 C 5.61390000 12.69150000 6.68790000 C 4.50420000 12.04830000 6.18680000 C 4.67770000 10.93430000 5.28310000 C 3.66800000 11.01310000 4.25400000 C 2.85330000 12.17760000 4.50200000 C 2.42990000 12.93000000 3.42910000 C 2.49170000 14.40400000 3.46690000 C 2.98100000 15.05470000 4.58220000 C 3.89880000 16.15440000 4.40660000 C 4.90220000 16.07580000 5.43260000 C 6.20620000 16.47750000 5.16550000 C 7.30280000 15.68150000 5.67130000 C 7.03490000 14.54540000 6.41880000 C 7.81270000 13.34730000 6.19520000 C 6.93810000 12.20560000 6.34520000 C 7.10580000 11.09580000 5.54110000 C 5.94430000 10.44010000 4.98210000 C 6.24180000 10.01680000 3.63740000 C 5.28070000 10.08840000 2.64130000 C 3.95890000 10.59560000 2.95880000 C 3.48160000 11.36600000 1.84110000 C 2.74470000 12.52000000 2.07740000 C 3.00320000 13.71060000 1.28860000 C 2.84950000 14.85820000 2.15190000 C 3.68870000 15.96390000 1.98290000 C 4.22500000 16.63050000 3.13960000 C 5.58480000 17.03350000 2.84470000 C 6.55840000 16.95630000 3.84120000 C 7.87640000 16.45110000 3.51990000 C 8.33390000 15.66330000 4.65310000 C 9.08140000 14.52020000 4.44260000 C 8.82170000 13.32760000 5.23010000 C 8.97820000 12.17640000 4.37280000 C 8.13700000 11.07830000 4.51980000 C 7.60210000 10.41370000 3.34810000 C 7.93950000 10.86570000 2.07200000 C 6.93140000 10.93790000 1.04000000 C 5.62270000 10.55680000 1.31180000 C 4.51630000 11.35180000 0.82050000 C 4.75940000 12.49080000 0.07470000 C 3.98680000 13.69090000 0.30710000 C 4.86470000 14.83290000 0.14460000 C 4.72190000 15.94390000 0.95880000 C 5.89140000 16.61070000 1.50020000 C 7.16190000 16.12730000 1.18790000 C 8.17790000 16.04970000 2.22530000 C 8.94950000 14.84860000 2.00490000 C 9.39400000 14.10150000 3.09610000 C 9.32420000 12.64820000 3.04980000 C 8.81900000 12.00850000 1.92780000 C 8.35480000 12.79150000 0.79700000 C 8.41560000 14.17500000 0.83360000 C 7.31430000 14.96860000 0.33470000 C 6.18380000 14.34060000 -0.17310000 C 6.12490000 12.89200000 -0.22510000 C 7.18130000 12.13560000 0.24900000 C 2.36230000 13.56340000 6.73990000 C 2.67400000 13.59320000 8.21560000 C 2.75730000 12.43390000 8.97120000 C 2.98080000 12.51050000 10.33750000 C 3.14370000 13.72670000 10.95180000 C 3.06950000 14.89230000 10.20790000 C 2.82790000 14.82210000 8.84660000 C 0.87330000 13.58330000 6.39020000 C 0.19050000 12.22970000 6.60460000 C -1.27100000 12.21200000 6.10350000 C -1.38760000 12.10640000 4.61500000 C -0.96090000 13.17800000 2.55440000 H 2.61900000 11.58570000 8.53800000 H 3.04770000 11.72860000 10.85410000 H 3.33290000 13.82080000 11.88330000 H 3.20620000 15.76940000 10.62210000 H 2.84250000 15.67120000 8.31090000 H 0.42290000 14.29440000 7.00660000 H 0.74840000 13.89850000 5.43050000 H 0.24080000 12.02210000 7.60170000 H 0.72430000 11.45640000 6.12160000 H -1.75880000 13.06610000 6.37680000 H -1.67030000 11.38940000 6.45750000 H -1.87310000 13.08920000 2.23900000 H -0.58960000 14.07390000 2.25580000 H -0.34920000 12.41000000 2.23220000 ?s 2.97920000 8.99550000 8.13160000 ?s 2.52550000 7.15780000 0.26900000 ?s 2.52550000 7.15780000 0.26900000 ?s 0.07150000 8.21900000 -0.47260000 ?s 0.07150000 8.21900000 -0.47260000 1 73 2.0 2 73 1.0 74 1.0 3 4 1.0 8 1.0 11 1.0 63 1.0 4 5 1.0 14 1.0 63 1.0 5 6 1.0 16 1.0 6 7 1.0 19 1.0 7 8 1.0 21 1.0 8 9 1.0 9 10 1.0 23 1.0 10 11 1.0 26 1.0 11 12 1.0 12 13 1.0 28 1.0 13 14 1.0 30 1.0 14 15 1.0 15 16 1.0 32 1.0 16 17 1.0 17 18 1.0 34 1.0 18 19 1.0 36 1.0 19 20 1.0 20 21 1.0 38 1.0 21 22 1.0 22 23 1.0 40 1.0 23 24 1.0 24 25 1.0 41 1.0 25 26 1.0 44 1.0 26 27 1.0 27 28 1.0 45 1.0 28 29 1.0 29 30 1.0 47 1.0 30 31 1.0 31 32 1.0 49 1.0 32 33 1.0 33 34 1.0 50 1.0 34 35 1.0 35 36 1.0 52 1.0 36 37 1.0 37 38 1.0 54 1.0 38 39 1.0 39 40 1.0 55 1.0 40 41 1.0 41 42 1.0 42 43 1.0 56 1.0 43 44 1.0 62 1.0 44 45 1.0 45 46 1.0 46 47 1.0 61 1.0 47 48 1.0 48 49 1.0 60 1.0 49 50 1.0 50 51 1.0 51 52 1.0 59 1.0 52 53 1.0 53 54 1.0 58 1.0 54 55 1.0 55 56 1.0 56 57 1.0 57 58 1.0 62 1.0 58 59 1.0 59 60 1.0 60 61 1.0 61 62 1.0 62 63 64 1.0 70 1.0 64 65 1.0 69 1.0 65 66 1.0 75 1.0 66 67 1.0 76 1.0 67 68 1.0 77 1.0 68 69 1.0 78 1.0 69 79 1.0 70 71 1.0 80 1.0 81 1.0 71 72 1.0 82 1.0 83 1.0 72 73 1.0 84 1.0 85 1.0 73 74 86 1.0 87 1.0 88 1.0 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93; 个人分类: 化学|4815 次阅读|0 个评论

更多...

帐号		自动登录	找回密码
密码			注册

关闭 安全验证

标签: gaussian

相关帖子

相关日志

关闭安全验证