Post

利用高斯调度xTB完成量子化学计算

Posted
By Bane Dysta
7 min read
利用高斯调度xTB完成量子化学计算

最近一个过渡态找麻了，体系大算得慢，还偶尔会藏个无垒的，此时想到了sobereva老师的高斯external脚本，于是拿来用了一下，非常好用，特此记录。

接口文档

后续可能考虑把atb2弄进来，这样就能算激发态优化了，所以大概记一下接口的使用方法。

根据说明书，external脚本使用关键词external=’‘实现。默认情况下调用方式为：

1

$Gau_External layer InputFile OutputFile MsgFile FChkFile MatElFile
  • layer：指ONIOM的layer。真实系统为R,中间层为M，模型系统为S
  • InputFile：在scratch产生的数据交换文件，用来向外部程序传递信息
  • OutputFile：要求外部程序写入信息的文件名
  • MsgFile：消息文件，创建此文件时会将文件内容复制到高斯log内
  • FChkFile：检查点文件，看原文似乎意思是某些过link402的任务需要往这里写信息，而不是OutputFile，其他没看懂，用到再说
  • MatElFile：矩阵文件，可以在脚本后面加关键词要求Gaussian导出一些单电子积分相关的东西之类的

显然对于不耍杂技的情况，最重要的是InputFile和OutputFile，这两个文件直接承担了信息交换功能。

InputFile格式为：

1
2
3
4
5

#atoms  derivatives-requested  charge  spin
原子数   导数阶数                电荷    自旋
atomic#  x  y  z  MM-charge MM-atom_type 
元素符号    坐标    MM电荷     MM原子类型
...

OutputFile格式为：

1
2
3
4
5
6
7
8
9
10
11

Items                          Pseudo Code                              Line Format
energy, dipole-moment (xyz)    E, Dip(I), I=1,3                         4D20.12
能量和偶极矩
gradient on atom (xyz)         FX(J,I), J=1,3; I=1,NAtoms               3D20.12
梯度
polarizability                 Polar(I), I=1,6                          3D20.12
极化率
dipole derivatives             DDip(I), I=1,9*NAtoms                    3D20.12
偶极导数
force constants                FFX(I), I=1,(3*NAtoms*(3*NAtoms+1))/2    3D20.12
力常数(Hessian)

只要接口程序解析InputFile，调用外部程序计算，写入OutputFile，Gaussian即可利用外部程序的信息开始计算。

gau_xtb

原版gau_xtb接口有三个文件，分别为genxyz、extderi和xtb.sh，在运行时需要全部放在计算目录下。关键词：

1
2

%nproc=1
# opt(ts,calcfc,noeigen,nomicro) external='./xtb.sh'

注意：

  • 如果进行opt任务，则需要写入nomicro，避免Gaussian调用自己的MM优化器
  • %nproc需要设为1，因为Gaussian本体并不进行QC计算，只作为driver引导external程序进行计算

笔者在运行时遇到了问题：

1
2
3
4
5
6
7

Running external command "./xtb.sh R"
         input file       "/home/apprepo/g16/scratch/Gau-123180.EIn"
         output file      "/home/apprepo/g16/scratch/Gau-123180.EOu"
         message file     "/home/apprepo/g16/scratch/Gau-123180.EMs"
         fchk file        "/home/apprepo/g16/scratch/Gau-123180.EFC"
         mat. el file     "/home/apprepo/g16/scratch/Gau-123180.EUF"
Failed to open output file from external program.

几番周折，最终定位到的问题是脚本直接在scratch目录似乎创建不了Gau-123180.EOu，但如果在本地创建Gau-123180.EOu再cp过去则没有问题，很邪门。

1
2
3
4
5
6
7

# 报错Failed to open output file from external program
./extderi $3 $atoms $derivs

# 正常计算
./extderi test $atoms $derivs
cp test $3

怀疑是权限相关的问题，但没有进一步探究，因为感觉三个文件还是有点麻烦，于是合并成了一整个文件，后面也没再出现Failed to open output file from external program的问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

!Fixed xtb interface for Gaussian external calls
!Based on original xtb.sh + genxyz + extderi logic
!Written by Dr. Tian Lu at Beijing Kein Research Center for Natural Sciences (www.keinsci.com)

program xtb_interface_fixed
implicit real*8 (a-h,o-z)
parameter (b2a=0.529177249d0)

! Element symbols array (same as original genxyz)
character*2 ind2name(0:86)
data ind2name /"Bq","H ","He", &
"Li","Be","B ","C ","N ","O ","F ","Ne", &
"Na","Mg","Al","Si","P ","S ","Cl","Ar", &
"K ","Ca","Sc","Ti","V ","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr", &
"Rb","Sr","Y ","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I ","Xe", &
"Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu", &
"Hf","Ta","W ","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn"/

character*256 inputfile,outputfile,cmd,c80
real*8 polar(6)
real*8,allocatable :: ddip(:),hess(:,:)

! Get arguments - handle both 6-arg (Gaussian) and 2-arg (legacy) formats  
nargs = command_argument_count()
if (nargs.eq.6) then
    call get_command_argument(2,inputfile)  ! Second arg is input file
    call get_command_argument(3,outputfile) ! Third arg is output file  
else if (nargs.eq.2) then
    call get_command_argument(1,inputfile)
    call get_command_argument(2,outputfile)
else
    stop
end if

! Read input file (first line: natm derivs charge spin)
open(10,file=inputfile,status='old')
read(10,*) natm,ider,icharge,ispin
close(10)

! Output start message
write(*,*) ">>>>> xTB interface start calculation"

! Generate XYZ file (following original genxyz logic)
open(10,file=inputfile,status='old')  
open(11,file='mol.xyz',status='replace')
read(10,*) natm,ider,icharge,ispin
write(11,*) natm
write(11,*)
do i=1,natm
    read(10,*) idx,x,y,z  ! Read only first 4 fields, ignore MM data
    write(11,*) ind2name(idx),x*b2a,y*b2a,z*b2a
end do
close(10)
close(11)

! Clean up previous files (following original xtb.sh)
call system("rm -f charges energy xtbrestart gradient hessian xtbout")
call system("rm -f xtb_normalmodes g98_canmode.out g98.out wbo xtbhess.coord")

! Run xtb (following original xtb.sh logic)
iuhf = ispin - 1
if (ider.eq.2) then
    write(cmd,'(a,i0,a,i0,a)') "xtb mol.xyz --chrg ",icharge, &
     " --uhf ",iuhf," --hess --grad > xtbout"
else if (ider.eq.1) then  
    write(cmd,'(a,i0,a,i0,a)') "xtb mol.xyz --chrg ",icharge, &
     " --uhf ",iuhf," --grad > xtbout"
else
    write(cmd,'(a,i0,a,i0,a)') "xtb mol.xyz --chrg ",icharge, &
     " --uhf ",iuhf," > xtbout"
end if

call system(cmd)

! Extract results (following original extderi logic)
open(11,file=outputfile,status='replace')

! Extract energy from xtbout (exact format from extderi.f90)
open(10,file='xtbout',status='old')
do while(.true.)
    read(10,'(a)') c80
    if (index(c80,'total energy').ne.0) exit
end do  
close(10)
read(c80,'(30x,f24.12)') ene

! Output Gaussian-style energy line for banelib
write(*,'(" SCF Done:  E(xTB) = ",f15.8)') ene

write(11,'(4d20.12)') ene,0d0,0d0,0d0

! Extract gradients if needed
if (ider.ge.1) then
    open(10,file='gradient',status='old')
    read(10,*)
    read(10,*)
    do iatm=1,natm
        read(10,*)
    end do
    do iatm=1,natm
        read(10,*) fx,fy,fz
        write(11,'(3d20.12)') fx,fy,fz
    end do
    close(10)
end if

! Extract Hessian if needed (following original extderi logic)
if (ider.eq.2) then
    polar=0
    write(11,'(3d20.12)') polar
    allocate(ddip(9*natm))
    ddip=0
    write(11,'(3d20.12)') ddip
    allocate(hess(3*natm,3*natm))
    open(10,file='hessian',status='old')
    read(10,*)
    read(10,*) ((hess(i,j),j=1,3*natm),i=1,3*natm)
    write(11,'(3d20.12)') ((hess(i,j),j=1,i),i=1,3*natm)
    close(10)
end if

close(11)

! Final cleanup (following original xtb.sh)
call system("rm -f charges energy xtbrestart gradient hessian xtbout mol.xyz tmpxx vibspectrum")
call system("rm -f xtb_normalmodes g98_canmode.out g98.out wbo xtbhess.coord .tmpxtbmodef") 
call system("rm -f .engrad *.engrad xtbtopo.mol xtbhess.xyz")

! Output end message
write(*,*) "xTB interface done <<<<<"

end program

编译方法：

1

gfortran -O2 -fdefault-real-8 -o xtb_external xtb_external.f90

其他技巧

如果目标体系hessian非常耗时，但已经得知比较可靠的初猜，可以考虑利用xtb产生初始hessian给过渡态计算读取：

1
2
3
4
5
6
7
8
9
10

%nproc=1
%chk=xtbhess.chk
# freq external='./xtb.sh'

...

--link1--
%oldchk=xtbhess.chk
%chk=ts.chk
# opt=(TS,rcfc,noeigen) freq b3lyp/6-31g*
Quantum Chemistry
This post is licensed under CC BY 4.0 by the author.
Total hits!