量子化学中的主要近似 - 图文

量子化学中的主要近似

量子化学中的主要近似：

1. 单Slater行列式近似

2. 自洽场近似

3. MO-LACO近似

各种近似表示列表： 零级 一级 高级

1. 波函数 单slater行列式 多行列式

稳定分子 CI CASSCF

2. Hamilton量 HF MP2，MP4 QCISD CASPT2

RHF，UHF CCSD

3. DFT 局域密度近似（LDA） GGA（B3LYP，BLYP，PBE）

4. MO-LCAO STO-3G 6-31g系列 aug-cc-pVDZ系列

密度泛函中的单Slater行列式 密度泛函中的单Slater行列式

1．密度泛函的波函数是无相互作用的波函数，仅仅是为得到电子密度而引入

r(r)=|y(r)*y(r)|

2．绝大多数情况下只需要单Slater行列式 3．含过渡金属体系的计算，可引入非整数电子填充方式，即按照能量的指数衰减函数把电子填入轨道，比如HOMO填1.6，LUMO填0.4个电子等等 Hamilton量的近似

Hamilton量的近似

1．Schrodinger方程是多体作用的方程，其Hamilton算符H是多体相互作用的算符。

2． Hatree-Fock算符F是单电子算符。核与其它电子对它的作用都用一个等效势能来代替。

3．单电子的Schrodinger方程是可以计算的。

4．要使整个多电子体系在单电子“各自为政”情况下合理共存，要使用自洽场(SCF)方法 密度泛函理论

密度泛函理论(DFT)的近似

1 基本原理：体系的基态能量由密度唯一确定。

2 基本方程：Kohn-Sham方程

3 E[r(r)]=E动能[r(r)]+E静电[r(r)]+EXC[r(r)]

4 动能项和静电项都与HF方法一样，不同之处在于交换相关项EXC[r(r)]

交换相关泛函

1 局域密度近似(LDA)：EXC[r(r)]

2 梯度校正（GGA)： EXC[r (r), ? r (r) ]

3 梯度校正并未增加很多计算量，因此一般都使用到梯度校正

4 常用泛函有：B3LYP(杂化)，BLYP，PBE

5 所有泛函都包含几个经验参数，由小分子拟合得到。

6 从计算结果与实验结果的对比上确定对体系合适的交换相关泛函, 如几何结构

DFT与从头算比较

1 精度： DFT ~ (HF+MP2) 大多数情况下!

2 可靠性：DFT: 需要把结果与实验值（如几何结构）对照，先验证泛函的合理性和适用性，再进行进一步预测其它性质

（HF+MP2)：大多数情况下结果都可靠，甚至可以推倒实验的解释结果

3 计算量：DFT << (HF+MP2)

4 弱作用：主要还是用从头算。 关联能的计算

直接的电子相关处理比Hartree-Fock或DFT需要更大的单电子基组，以产生趋于收敛的结果，因此这样的计算可能相当昂贵。对于给定

的基组，关联能计算通常比HF计算更加昂贵，因此很多没有经验的人试图用小基组进行关联能计算。但这是完全不合理的，对于有意义的计算，至少应当使用加上几个极化函数的三-zeta基组（例如，cc-pVTZ）。

MOLPRO 对波函的近似和优化有很多不同的方法，例如M?ller-Plesset (MP)微扰理论，组态相互作用(CI)，或者耦合簇(CC)方法。密度泛函(DFT)方法也考虑了电子关联，虽然它是一种不如从头方法成系统，定义也不是很好的方法。

需要注意，HF近似，以及所有使用HF行列式作为零级近似的单参考方法，通常仅使用于平衡结构附近的计算。在大多数情况下，它们不能正确离解分子键，或者描述电子激发或（近）简并态。在这种情况下，使用多组态SCF波函(MCSCF)作为零级近似的多参考方法，可以提供合理的选择。完全活性空间SCF (CASSCF)是MCSCF的一种特殊变体。在MOLPRO中可以使用各种多参考电子相关方法，例如多参考微扰理论(MRPT,CASPT2)和多参考组态相互作用(MRCI)，以及它们的变体，例如多参考耦合对泛函(MR-ACPF)。 molpro和gaussian基组的异同

1 在molpro中基组的表达是basis=name

在molpro的网页中可以找到能用的基组。

其中比较常用的基组是dunning等人的相关一致基组：cc-pvdz（双zeta）（在molpro中可以表示为vdz），cc-pvtz（三zeta，vtz），cc-pvqz（vqz）等。 对于某些特性的计算，如电子亲和力，极化率，或者分子间势，还需要加上弥散函数，这些函数可以表示为：avdz，avtz，avqz等。molpro默认的是cc-pvdz。

另外一个常用的是pople等人的高斯基组，如6-31G**等，但是需要注意的是molpro默认使用球谐函数（5d，7f等），而gaussian使用笛卡儿函数（6d，10f等）。这会产生略为不同的能量。笛卡尔函数也能用于molpro，但是关键字cartesian必须在基组定义以前给出。

2 molpro也可以对不同的原子使用不同的基组： basis，C=vtz，O=avtz，H=vdz等

3 molpro也可以人工定义基组（指数和收缩因子）

总体上，gaussian可以使用的基组，molpro都可以使用，但是由于molpro的计算精度比gaussian高，所以molpro使用的基组一般是比较大的基组。

molpro中处理弱相互作用之一 SAPT

在molpro中处理弱相互作用目前我知道的有3种方法： 1 SAPT（只能限制在dimer体系） 2 BSSE （超分子方法） 3 LMP2 （局域电子相关方法） 下面分别一一介绍：

1 SAPT （没有BSSE误差）

在molpro中的输入文件：

SAPT的输出文件：

与ccsdt结果的比较：

molpro中处理弱相互作用之二 超分子方法--BSSE

2 超分子方法中考虑BSSE修正 这和gaussian中差不多，但是输入文件不同 bsse的介绍：

molpro中的BSSE的输入文件：

不同方法的BSSE修正的比较：

molpro中处理弱相互作用之三 局域电子相关方法-LMP2

3 LMP2

其实这个方法，我还没有认真读文献，他的精髓我还不是特别的懂，但是我知道，他的计算精度很高其中DF-SCS-LMP2可以和CCSDT的结果相媲美。 DF-LMP2也不错啦，还有什么F12（还不知道是什么东西呢）精度也很高。

这个可以用来处理大一些的体系，但是有一点很遗憾：在2002版本中并没有DF-SCS-LMP2。 参考文献见：JPCA 1998，102，5997-6003 （Martin Schutz，中间的u上面有两点哦..)(LMP2的方法介绍）

Physical Chemistry Chemical Physics 2006, 8, 4072-4078 JG Hill, JA Platts, HJ Werner （DF-SCS-LMP2方法介绍）

\

the reference structure should preferably have higher symmetry ， so that for symmetry reasons the adiabatic and diabatic states are identical.

20.) Geometry optimization for large molecules. Use DFT or df-lmp2 for larger molecules.

22.) Local Minima in mcscf/casscf orbital optimization.

It can happen that in MCSCF/CASSCF calculations local minima are found in the orbital optimization. Often it helps first to optimized the orbitals for a smaller active space and then increase the active space gradually. It also helps to first run a calculation with all CLOSED set to FROZEN, and then a subsequent calculation that starts from these orbitals.

23.) How can we get a trial version for three months?

This should now have been sent to everybody who is registered and here.

24.) Can you tell me more about installation and parallelization?

All aspects of installation, including for parallel execution, are described in the installation

guide. Running the program in parallel is described in the User Manual, section 2.2. One can either launch jobs with the supplied molpro script, or directly using standard local mpirun or mpiexec.

25.) What are the roles of the various output files that are produced?

Normally, the output from the job is written to the .xml file. This file consists of most parameters and results written in a well-defined standard format, capable of post-processing by the many XML tools and libraries that are available. Most of the output, however, particularly all those parts that give additional human-readable information, appear as comments in the .xml file. As the job is running, the .xml is filtered such that just these comment parts are written to the .out file. If the .xml or .out files already exist at the start of the job, the old versions are renamed to a file with an underscore and a number. Additionally the program creates a .log file during geometry optimizations and numerical frequencies; this contains the output for the individual geometry points, not normally of interest, but sometimes useful if the job fails. As well as this, other files can be created as the result of PUT or CUBE comman

molpro的对称性

molpro中对称性还是有一些复杂的，当然了，复杂是相对于我这样的比较搓的人来讲的，对于那些点群学的好的，可能就不是复杂，而是觉得的理所当然哈

molpro不能使用非阿贝尔群，换句话说就是只能使用阿贝尔群。Abelian point group需要看书哦。

molpro自动 recognizes and exploits molecular symmetry. 当然也可以根据需要把对称性降低或者取消对称性

molpro里边把不可约表示用1-8表示： 下面是D2h点群的对称性的表示：

下面是其他对称性的表示：

molpro中scf程序

见下面的示例：

alpha=108.3121 degree r=1.370848397 ang geometry={O; !define z-matrix O2,o,r; O3,o,r,O2,alpha} hf;occ,6,2,4,1;closed,6,1,4;wf,24,3,2 uccsd(t) optg freq,symm=auto zp(1)=zpe e(1)=energy （1）其中hf行就是scf过程

对于闭壳层体系，hf后边一般不需要其他的关键字

但是对开壳层体系，推荐使用wf，occ，closed或open卡，用来唯一定义波函数。 A wf卡定义波函数的电子数和总对称性 wf，elec，sym，spin 其中 elec 定义电子数

sym 波函数总的对称性（不可约表示符号）

spin 定义自旋对称性，spin=2*s（单重态=0，双重态=1，三重态=2，其实也就是但占据的数目）

B occ卡定义每个对称性的占据轨道数

occ，n1，n2，......,n8;

对于高对称性的情况，为避免收敛问题，这个卡总要包含在输入文件中。 ni是不可约表示i中的占据轨道的数量。轨道的总数必须等于（elec+spin）/2 C closed定义闭壳层的轨道 closed，n1，n2，......,n8;

这一项可以用于开壳层计算，指定每个对称性中闭壳层轨道的数量。这可以在不使用open卡的情况下，强制特定的态。

（2）freq,symm=auto

在数值计算Hessian的过程中，分子的对称性可能会降低。如果给定symm=auto，程序在每一个能量/梯度计算中将使用分子波函数可能的最高对称性，因此这个选项将使计算量最小。 如果使用symm=no，则在计算中不使用对称性，这是默认的。

轨道局域化

轨道局域化

轨道局域化是依据boys和pipek-Mezey判据进行计算 局域化是在每个对称类中进行

若是需要全部局域化，就不要使用对称性 局域化程序用locali命令调用：

locali，method ！method可以是boys可以是PM。 例： ***,h2o memory,4,m gthresh,energy=1.d-8 geometry={nosym;

H1,, 0.000000, 1.415022, 1.558516 O1,, 0.000000, 0.000000, 0.465140 H2,, 0.000000, -1.415022, 1.558516 }

basis=vdz

hf;

locali,pipek; ！scf过程使用局域化

locali,boys;

mp2;local;orbital,local(boys) ！METHOD，method（boys，PM），这难到是另外一种写法？ ！ Molecular orbitals read from record 2100.2 Type=RHF/LOCAL(BOYS)

e(1)=energy

ibaso(1)=0 ！不知道这个是干嘛的 prog(1)='LMP2/BOYS' mp2;local;orbital,local(pm)

e(2)=energy prog(2)='LMP2/PM' ibaso(2)=0

locali,boys;save,2101.2 ！Localized orbitals saved to record 2101.2 locali,pipek;save,2101.2

mp2;local,ibaso=1;orbital,2101.2,local(boys) ！ Molecular orbitals read from record 2101.2 Type=RHF/LOCAL(BOYS） e(3)=energy

prog(3)='LMP2/BOYS' ibaso(3)=1

mp2;local,ibaso=1;orbital,2101.2,local(pm) e(4)=energy prog(4)='LMP2/PM'

ibaso(4)=1 show,e

molpro中如何判断波函数是否稳定

在molpro中计算中，

（1）如何确定波函数是不是稳定呢（只是看单点能是不是单调变化的吗？）

（2）初始占据是不是正确？应该怎么确定呢？ 解答：

对于态平均mcscf计算，hf波函是否稳定并不重要，只要收敛就行，开壳层体系可以通过加减电子的办法变成闭壳层，作为mcscf初始波函。对于闭壳层的单参考计算，当homo-lumo gap > 1 eV的情况下，molpro的初始猜测通常都是合理的，波函数通常是稳定的。但是如果低于1 eV，或者开壳层体系，波函有可能不稳定。这类体系的初始占据没有好的判断方法，因为不同方法得到的基态和电子组态经常不相同。可以查文献的报道，如果计算量不大可以通过mcscf/caspt2/mrci/ccsd计算来确定。

波函不稳定有时候也包括结构不稳定，比如把水分子作为线型分子计算，波函就是不稳定的。通过频率计算发现有虚频，可以知道波函不稳定。

molpro中双体的相互作用能的计算 （dummy）cp校正

dimer相互作用的求算相互作用能的bsse（cp）方法的设置：

下面是个例子：h2o-ar memory,1,m

proc addit ！定义了一个进程 if (#calc.le.0) then calc=0; endif; calc=calc+1

energies(calc)=energy dipoles(calc)=dmz programs(calc)=program endproc ！进程结束

gthresh,energy=1e-10 basis=aug-cc-pVDZ

geometry={h1;o,h1,roh;h2,o,roh,h1,theta;ar,o,r,h1,180-theta/2,h2,180} roh=1 angstrom, theta=104 degree, r=4.5 angstrom

rhf;addit;mp2;cphf,1;addit;eall=energy; ！计算整个体系的能量，dimer的能量（cphf方法在2002中没有关键字，addit就是上面定义的进程

dummy,1,2,3 ！ 把h o设成了鬼原子

rhf;addit;mp2;cphf,1;addit;earghost=energy; ！ 求算AR的能量，在h2o的核电荷为0的情况下

dummy,4

rhf;addit;mp2;cphf,1;addit;eh2oghost=energy; ！ 同理求算h2o的能量

deghost=(earghost+eh2oghost-eall)*tocm ！这个是考虑了平衡修正后的体系的相互作用的大小

geometry={ar};rhf;mp2;cphf,1;ear=energy; ！计算单体ar的能量

geometry={o;h1,o,roh;h2,o,roh,h1,theta};rhf;addit;mp2;cphf,1;addit;eh2o=energy; ！计算单体h2o的能量

deghost=(earghost+eh2oghost-eall)*tocm ！这个是考虑了平衡修正后的体系的相互作用的大小

de=(ear+eh2o-eall)*tocm ！这个是没有考虑平衡修正后的体系的相互作用的大小 cph2o=(eh2o-eh2oghost)*tocm ！ 单体水的修正能 cpar=(ear-earghost)*tocm ！ 单体Ar的修正能 show,programs,energies,dipoles

进行平衡修正的结构优化计算

这是进行平衡修正的结构优化计算 以HF二聚体为例:

***,HF dimer mp2/CP optimization basis=avdz

gthresh,energy=1.d-8

maxit=20 !max number of iterations ! INITIAL VALUES OF GEOMETRY VARIABLES RFF= 5.31431160 R1= 1.75768738 R2 = 1.75298524 THETA1 = 7.03780227 THETA2 = 111.25930975 geometry={x;noorient f1 f2 f1 rff

h1 f1 r1 f2 theta1

h2 f2 r2 f1 theta2 h1 180.} !几何构型

do iter=1,maxit !optimization loop ！优化的最大循环次数 text, CALCULATION AT LARGE SEPARATION

rff_save=rff !save current rff distance ！把当前的参数给了rff_save

rff=1000 !dimer calculation at large separation ! 把双体的距离拉得特别远 text, HF1 ！下面开始计算HF1的能量 （没有所谓的cp校正）

dummy,f2,h2; !second hf is now dummy （把HF2设成dummy） hf;accu,16 !scf for first monomer mp2; !mp2 for first monomer

ehf1inf=energy !save mp2 energy in variable

dip1inf=sqrt(dmy^2+dmz^2) !total dipole moment of first monomer

forces; !compute mp2 gradient for first monomer （计算能量梯度） text, HF2 ！ 同理重复HF2 dummy,f1,h1; !first hf is now dummy hf;accu,16 !scf for second monomer mp2; !mp2 for second monomer ehf2inf=energy !save mp2 energy in variable

dip2inf=sqrt(dmy^2+dmz^2) !total dipole moment of second monomer forces; !compute mp2 gradient for second monomer

add,1 !add to previous gradient （当前的梯度根据给定的因子1添加到前一个上面） 为什么要add，1， 我不知道

rff=rff_save !reset HF - HF distance to current value （把两个单体的距离设到正常）

text, CPC calculation for HF1 MONOMER ！下面开始cp计算，先是HF1，然后HF2，为什么这里要用add，-1，我不知道。

dummy,f2,h2; !second hf is now dummy hf;accu,16 !scf for first monomer mp2; !mp2 for first monomer

ehf1=energy !save mp2 energy in variable

dip1=sqrt(dmy^2+dmz^2) !total dipole moment of first monomer forces; !compute mp2 gradient for first monomer

add,-1 !subtract from previous gradient （当前的梯度根据给定的因子-1添加到前一个上面）

text, CPC calculation for HF2 MONOMER dummy,f1,h1; !first hf is now dummy hf;accu,16 !scf for second monomer mp2; !mp2 for second monomer ehf2=energy !save mp2 energy in variable

dip2=sqrt(dmy^2+dmz^2) !total dipole moment of second monomer forces; !compute mp2 gradient for first monomer add,-1 !subtract from previous gradient

dummy !reset dummies

text, DIMER CALCULATION ！双体的计算 hf;accu,16 !scf for dimer mp2; !mp2 for dimer

edimer=energy !save mp2 energy in variable dipdim=sqrt(dmy^2+dmz^2)!total dipole moment of dimer

forces; !compute mp2 gradient for dimer

add,1 !add to previous gradient （这个地方是要加的） einf=ehf1inf+ehf2inf !total energy of unrelaxed momomers decpc=einf-ehf1-ehf2 !counter poise correction

de=edimer-ehf1-ehf2 !Interaction energy relative to unrelaxed monomers

etot=edimer+decpc !total cpc corrected energy ！在这以前都没有优化，就如前面的dummy计算平衡修正能差不多，但是计算了 forces (为什么呢）。 opt; !find next energy ！ 下面开始优化 show,optconv !show convergence parameter if(optconv.lt.0.0001) goto,lab: !test for convergence

enddo !end of optimzation loop lab:

text, compute optimized monomer energy rhf=r1 geometry={h1 F1,H1,rhf}

hf;accu,16 !scf for relaxed monomer mp2; !mp2 for relaxed monomer ehf=energy !save mp2 energy in variable

dipmon=sqrt(dmy^2+dmz^2)!total dipole moment of relaxed monomer optg !optimize HF structure

text, compute interaction energy including relaxation erelax = 2*ehf - einf !relaxation energy

derel=de-erelax !interaction energy relative to relaxed monomoers text, results

show,r1,r2,rhf,rff,theta1,theta2 !optimized geometry parameter

molpro中直接积分direct

除了CI和三重激发微扰（T）外，molpro中所有的方法都可以使用直接积分的方法。 当对小分子使用大基组时，direct可以通过需要很多CPU时间，但是减少磁盘空间的需求。对于大分子尺度，有150－200个以上的基函数的计算通常应当使用direct。 例如：

***, Formaldehyde memory,2,m basis=vdz

gthresh,energy=1.d-8 geomtyp=xyz geometry={x 4

FORMALDEHYDE

C 0.0000000000 0.0000000000 -0.5265526741 O 0.0000000000 0.0000000000 0.6555124750 H 0.0000000000 -0.9325664988 -1.1133424527 H 0.0000000000 0.9325664988 -1.1133424527

}

direct ！直接积分计算，一般若是它在特定程序的输入区之外，它的作用是全局的，但是它可以在特定的 程序之内，那么它只影响该程序，gdirect的作用总是全局的。 hf

mp2;local,save=6000.2,start=6000.2 ！ save，和start指令 optg

show,energy,dmz

molpro 中 file 和 restart卡

1 file 卡

*** file,1,h2o.int,new file,2,h2o.wfu,new geometry={

在这个里边，file卡从例子中可以看到它的用法，

其中最要的是后边那个new，当然它也可以用其他的如：unkown，old（如果文件存在则继续重新开始计算）等等代替。 2 restart的使用

*** file,1,h2o.int,new file,2,h2o.wfu,new geometry={ H1 O,H1,2; H2,O,2,H1,100} hf e(1)=energy method(1)='hf' mp2 e(2)=energy method(2)='mp2' ccsd(t) e(3)=energy method(3)='ccsd(t)' ***

file,2,h2o.wfu ！ 这样就重新开始了一个工作，读取h2o.wfu的信息 ci e(4)=energy method(4)='ci (restart,2)' *** file,1,h2o.int

file,2,h2o.wfu ！ 这样就重新开始了一个工作，读取h2o.int和h2o.wfu的信息 ccsd(t) e(5)=energy method(5)='ccsd(t) (restart,1,2)' show,e

其中最要的是后边那个new，当然它也可以用其他的如：unkown，old（如果文件存在则继续重新开始计算）等等代替。 2 restart的使用

*** file,1,h2o.int,new file,2,h2o.wfu,new geometry={ H1 O,H1,2; H2,O,2,H1,100} hf e(1)=energy method(1)='hf' mp2 e(2)=energy method(2)='mp2' ccsd(t) e(3)=energy method(3)='ccsd(t)' ***

file,2,h2o.wfu ！ 这样就重新开始了一个工作，读取h2o.wfu的信息 ci e(4)=energy method(4)='ci (restart,2)' *** file,1,h2o.int

file,2,h2o.wfu ！ 这样就重新开始了一个工作，读取h2o.int和h2o.wfu的信息 ccsd(t) e(5)=energy method(5)='ccsd(t) (restart,1,2)' show,e

本文来源：https://www.bwwdw.com/article/u7uv.html