高斯fch文件与wfn波函数文件的介绍及转换方法 | 宜武汇-ag真人国际厅网站

了解高斯的格式化检查点文件fch和波函数文件wfn的格式和每项的意义对编写第三方量化辅助程序是很重要的,前者蕴含的信息丰富,后者更为简明而且通用。本文将解读wfn文件的格式和读、写代码,对fch中的一些要点进行介绍,介绍如何写出将fch转化为wfn文件的程序fch2wfn,以帮助大家认识这两种文件并利用fch和wfn写出相关程序。

1. wfn(wavefunction)文件的格式和意义

wfn文件最早作为bader的aimpac程序的输入文件出现,aimpac是历史最久的aim分析软件,后来其它的aim分析软件如multiwfn (

)、aim2000、aimall、morphy、xaim等等也都需要wfn文件。wfn文件包含的是所有高斯函数的信息和波函数向它们的展开系数,结构简洁,计算空间中波函数值、电子密度等属性很方便。gaussian、gamess等许多主流量化程序都可以输出wfn文件,例如gaussian中只需要在route section加上output=wfn,并在分子结构部分末尾空一行写上wfn文件输出路径即可,如c:\divokej_bill.wfn。

这里以高斯中hf/3-21g计算的ch4分子输出的wfn文件为例对格式进行介绍

第1行

title card required

wfn中第一行没有任何意义,只是文件注释

第2行

gaussian              5 mol orbitals     27 primitives        5 nuclei

开头的gaussian说明文件中记录的函数类型是高斯类型,一般wfn中都是高斯类型函数。

5 mol orbitals代表文件中的分子轨道数,wfn中只记录占据轨道,虚轨道不记录,因为aim分析的对象是电子密度,而虚轨道对计算电子密度没有贡献。

27 primitives代表原始高斯函数的总数目,可独立变分的基函数是由它们收缩成的。

5 nuclei是原子数

第3-7行

这些是原子的信息。例如其中的第一行,c是原子符号,后面的1是原子序号。centre后面的也是原子序号,与前面的序号是重复的,所以取其一即可。接下来是原子坐标,单位为波尔,最后是原子的核电荷数。使用赝势时核电荷数不等于原子序数,而等于被明确表达的价层电子数。

第8、9行

centre assignments    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  2  2  2  3  3

centre assignments    3  4  4  4  5  5  5

每项代表每个高斯函数所在的原子序号。由于一共有27个高斯函数,这里也就有27项。

第10、11行

type assignments      1  1  1  1  1  2  2  3  3  4  4  1  2  3  4  1  1  1  1  1

type assignments      1  1  1  1  1  1  1

这是每个高斯函数的类型。wfn文件中一般都是用笛卡尔型高斯函数,这里的数字与类型的对应关系是

1/2/3/4/5/6/7/8/9/10=s/x/y/z/xx/yy/zz/xy/xz/yz

11/12/13/14/15/16/17/18/19/20=xxx/yyy/zzz/xxy/xxz/yyz/xyy/xzz/yzz/xyz

生成wfn文件时应注意应使用6d/10f型高斯函数而不要用球谐型。wfn文件对于g及更高角动量函数没有统一的序号定义,而且由于涉及这样的函数机会不多,所以多数aim软件还不支持这些函数,如果体系中有这样的函数,高斯也无法输出wfn文件。

第12-17行

这是每个高斯函数的指数,即n*x^ix*y^iy*z^iz*exp(-α*r^2)中的α。显然也是27项。

第18-52行

这部分是分子轨道信息,这里只列出第3条以说明。

mo    3代表是第3条轨道,也就是它在高斯中的轨道序号。对于闭壳层,这个序号和此轨道在本文件中的序号一致。

mo 0.0没有任何意义,只是个标签。

no =    2.0000000是占据数,这里就是代表有两个电子占据,若是非限制性波函数即为1。若是限制性开壳层计算,高斯会误把单占据轨道占据数写为2,需手动修改为1。

orb. energy =   -0.548298是轨道能,单位为hartree。

接下来就是分子轨道向这27个高斯函数的展开系数,一定要注意,每个展开系数已经包含了高斯函数的归一化系数,也就是分子轨道向没有归一化的高斯函数(即x^ix*y^iy*z^iz*exp(-α*r^2)形式)展开的系数。

这里的轨道并非一定是分子轨道,如果用后hf方法,并且加了density=current表明用后hf密度,则输出的每个轨道就是自然轨道,其占据数一般不是整数,输出的轨道数目与可独立变分的基函数一致。因为后hf密度的自然轨道没有是否占据之分,只有占据多少之分,每条轨道都对电子密度有贡献,所以所有的自然轨道都会列出。如果是开壳层计算,还应当加上pop=noab,这样就会记录自然自旋轨道。

轨道记录顺序是能量由低往高。注意对于自旋轨道,wfn中先按这个次序记录完所有alpha轨道再记录beta轨道。由于第一个beta轨道的能量肯定比最后的alpha轨道的能量低,所以通过检验轨道能量与前一个轨道能量的关系就能判断beta轨道从哪里开始记录。如果是自然自旋轨道,则要判断它与前一个轨道的占据数关系。

wfn中轨道的轨道序号,记录的前后顺序是随意的。高斯函数的顺序也是随意的,只要在类型、所在中心、指数、轨道展开系数段落中的位置对应就行,不会影响结果。

第53行

end data

其中end data代表标准wfn文件要求的全部信息的说明已结束,在此之后允许附加一些信息。

第54行及之后

 the  hf energy =    -39.976405845082 the virial(-v/t)=   2.00033480

这些是附加信息,是随意的。一般量化程序输出的wfn都会按照这个格式输出能量和维里值。the  hf energy是体系hf或dft总能量,如果是后hf方法,这里仍然是hf能量。有的程序还会附上其它一些信息。

2.读、写wfn文件的fortran代码解析

首先根据wfn所拥有的数据,定义一些变量,以及两种自定义类型用于储存原子和高斯函数的数据。

type atomtype

character*2 name !元素名称

real*8 x,y,z,charge !原子坐标和电荷

end type

type basistype

integer center,functype !高斯函数所在原子序号,函数类型

real*8 exp !指数

end type

character*8 mode !用来记录函数类型,一般wfn里都是高斯函数,所以可以不定义。

integer nmo,nprims,ncenter !轨道、高斯函数、原子的数目

real*8 totenergy,virialratio !总能量、维里值

type(atomtype),allocatable :: a(:) !储存原子信息

type(basistype),allocatable :: b(:) !储存基函数信息

real*8,allocatable :: moocc(:),moene(:) !轨道占据数和能量

real*8,allocatable :: co(:,:) !系数矩阵,co(i,j)代表第i个分子轨道向第j个高斯函数的展开系数,已包括归一化系数。

还定义两个变量,并不是直接从wfn文件中读取,需要用代码来判断,在写一些相关程序时会用得到它们。

integer,allocatable :: motype(:) !轨道类型,0代表双占据的空间轨道,1是alpha自旋轨道,2是beta自旋轨道

integer wfntype !用来说明波函数类型,0/1/2代表r/u/rohf波函数,3/4代表是闭壳层/开壳层后hf波函数,分别代表储存的是空间自然轨道和自旋自然轨道。

下面是读取wfn文件的子程序,是作为主程序的内部函数使用的,能共享前面在主程序中定义的变量,主程序中已使用了隐式声明implicit real*8(a-h,o-z)。

了解了wfn读入方法,输出wfn也就会了,以下是输出子程序。

3. fch(formatted check point)文件简介

任何平台下高斯执行任务生成的二进制的chk文件都可以用formchk生成fch/fchk文件,fch文件是ascii文件,是人可读的。它独立于平台,可用于在不同平台上交换chk文件,也可方便地从中提取计算过程中有用的信息。fch文件在gaussian03与09中内容大体不变,但是记录项目的顺序有所改变。

fch中的信息远比wfn丰富,包含了获得wfn文件所需要的一切信息,里面也有不少冗余的信息。本文主要目的是介绍fch向wfn的转换方法,如何写fch2wfn程序,下面只谈一些相关要点,一些条目靠fch自带的解释就可以明白就不谈了。

第一行是标题,是任意的,可以自行写一些附加信息在此告诉fch2wfn程序对它做哪些特殊处理。

第二行是任务类型、方法、基组,格式是a10,a30,a30。看第11列是r、u还是ro就能判断波函数类型。

每一类数组型数据前面都有一行标签,读取格式都是a40,3x,a1,5x,i12。比如

primitive exponents                        r   n=          18

18代表下面有18个数据。r代表下面的是实数型数据。若是i、c、l则分别代表是整型、字符型、逻辑型数据。实数型数据读取格式是5(1pe16.8),整型是6i12。

fch中与wfn显著不同的一点是基函数的储存方法。wfn中是按高斯函数顺序储存的,而fch中是一层一层地储存的。shell types代表了每个壳层的类型,0=s, 1=p, -1=sp, 2=6d, -2=5d, 3=10f, -3=7f, 4=15g, -4=9g。number of primitives per shell代表的是每一壳层内的基函数的收缩度。

fch中每个壳层里笛卡尔型基函数的类型是按如下顺序排列s,x,y,z, xx,yy,zz,xy,xz,yz, xxx,yyy,zzz,xyy,xxy,xxz,xzz,yzz,yyz,xyz,对比wfn中的函数类型编号,会发现f轨道的顺序不一致,这在转换时需要特别注意。如前所述,wfn文件里不能记录g及更高角动量函数,但是在fch文件中可以,g函数的顺序为zzzz,yzzz,yyzz,yyyz,yyyy,xzzz,xyzz,xyyz,xyyy,xxzz,xxyz,xxyy,xxxz,xxxy,xxxx。

coordinates of each shell其实是多余的,因为已经有shell to atom map(壳层所属原子)和current cartesian coordinates(每个原子xyz坐标),可以直接推出。

total scf density就是scf密度矩阵,用了后hf方法并用density=current还会有后hf的密度矩阵被记录,如total mp2 density。由于是对称的,所以只记录了下三角部分(含对角元),顺序是(1,1)(2,1)(2,2)(3,1)(3,2)(3,3)(4,1)…,

限制性方法计算的轨道能量和展开系数记录在alpha orbital energies/alpha mo coefficients里(尽管双占据称为alpha轨道并不妥),如果是非限制性的计算还会多出beta的。用density=current并不会使fch文件像wfn一样使记录的轨道成为自然轨道,它仍然是scf的系数。如果想把自然轨道系数写入到alpha mo coefficients当中,则先用后hf方法带着density=current关键字跑一次,然后再运行guess=(save,only,naturalorbitals) chkbasis即可。此时fch里轨道能量部分就不再是能量,而是占据数。

加了pop=savenbo使高斯调用nbo 3.1进行nbo分析后,会把nbo轨道保存到fch的分子轨道系数部分。如果用的scf密度,则分子轨道能量部分保存的是nbo轨道能量,读入下文的fch2wfn程序前应把第一行改为savenboene;如果用的是后hf密度,这部分保存的则是nbo轨道占据数,读入fch2wfn前应把第一行改为savenboocc。fch中倒数的几个nbo轨道的能量或占据数误为1000,应根据nbo分析输出的结果手动修改。

primitive exponents、contraction coefficients与p(s=p) contraction coefficients的关系:

这里举hf/sto-3g算ch4的fch的这部分作为例子,这三类只列出第一行来说明,最开头的就是c的s和sp壳层。

primitive exponents                        r   n=          18

  7.16168373e 01  1.30450963e 01  3.53051216e 00  2.94124936e 00  6.83483096e-01

contraction coefficients                   r   n=          18

  1.54328967e-01  5.35328142e-01  4.44634542e-01 -9.99672292e-02  3.99512826e-01

p(s=p) contraction coefficients            r   n=          18

  0.00000000e 00  0.00000000e 00  0.00000000e 00  1.55916275e-01  6.07683719e-01

下面是一些linux版高斯自带的sto3g.gbs文件,

-c

s    3 1.00

 0.7161683735d 02  0.1543289673d 00

 0.1304509632d 02  0.5353281423d 00

 0.3530512160d 01  0.4446345422d 00

sp   3 1.00

 0.2941249355d 01 -0.9996722919d-01  0.1559162750d 00  注:第一列是高斯函数指数,后两列分别是s和p型的收缩系数

 0.6834830964d 00  0.3995128261d 00  0.6076837186d 00

 0.2222899159d 00  0.7001154689d 00  0.3919573931d 00

将fch的内容竖着看,而将sto3g.gbs的内容横着看,可见fch中与sto3g.gbs的内容是对应的。对于s壳层,p(s=p) contraction coefficients对应的项都是0,到了sp壳层就不为0了,代表的就是p的收缩系数。如果用的是例如dunning相关一致性基组,由于s和p壳层是独立的,没有sp壳层,则这部分就不会出现。

4.读取fch文件的fortran代码

为了完成fch->wfn的转换,不需要读入所有fch的数据,只需要读入其中有用的数据并经过转换,把将前面outwfn子程序要用的变量全都填好就可以了。

读取fch文件的代码中要用到在fch中定位标签的子程序,例如call loclabel(10,’alpha mo coefficients’)就会使当前读取位置处在alpha mo coefficients这行开头,如果没找到要找的内容,则全局变量noentryinfch=1,否则为0。loclabel子程序代码如下:

下面是读入fch文件的子程序readfch:

5.由fch的系数矩阵生成wfn的系数矩阵

后面的代码要用到计算高斯函数归一化系数的函数,如下所示。函数需要输入高斯函数类型代表的数字和指数。其中ft()代表阶乘函数。

从fch的系数矩阵amocoeff和bmocoeff转化到输出wfn文件要用的系数矩阵co过程略微复杂,故单独说明。下面的图对3-21g计算的ch4的高斯函数循环过程进行了说明,每一排球是一个壳层,每个球是一个基函数,球中每条线是一个高斯函数,箭头指明了高斯函数循环次序。要把fch中通过层层循环描述的基函数的系数转换到wfn中的每个高斯函数的系数,就需要先把每个高斯函数都循环一遍,循环过程中赋予每个高斯函数系数、所在中心、指数和类型。

下面的代码可以进行这个工作。最外层循环每一壳层(桔色箭头),在壳层中循环每个基函数(粉色箭头),在基函数中循环每个高斯函数(绿色箭头),在每个高斯函数中循环每个分子轨道。

将上述内容组合在一起,就成为了fch2wfn程序,完整的代码及编译好的程序见

,在cvf6.5、ifort编译通过。可以运行后输入fch文件名和输出的wfn文件名,也可以直接用命令行模式,如fch2wfn c:\sob.fch ..\saint.wfn。再次提醒,含自然轨道的fch开头一行应改为saveno。含nbo的fch若用了后hf密度要把开头改为savenboocc,若用了scf密度要把开头改为savenboene,以便fch2wfn正确处理。

转换的结果和直接用高斯输出的wfn文件几乎完全一致,一般只有末位由于数值精度问题相差1或2。但有时会发现含自然轨道信息的fch转化出来的wfn中有几条轨道系数和高斯输出的wfn不一致,但实际上都是正确的,只是自然轨道取向不同而已。

aimall中的模块aimqb也能将fch向wfn的转换,但aimqb转换出的wfn虽然内容正确,与fch2wfn一致,但格式与标准格式不符,可能造成一些aim程序无法读取。而且高斯函数的储存顺序与一般规则有异(当然这完全不影响结果),比如3-21g,一般sp壳层的高斯函数的类型储存顺序是x x y y z z | x y z,其中|代表分裂价基的分隔,而aimqb转换出来的顺序是x y z x y z x y z。fch2wfn的转换不仅结果正确,也完整保留了wfn文件的习俗。fch2wfn的功能也已经整合进了multiwfn中,在multiwfn中读入fch后选功能6,再选保存wfn文件即可。

原文链接:http://sobereva.com/55

网络摘文,本文作者:15h,如若转载,请注明出处:https://www.15cov.cn/2023/08/27/高斯fch文件与wfn波函数文件的介绍及转换方法/

发表评论

邮箱地址不会被公开。 必填项已用*标注

网站地图