Scarlet Land

0%

电子积分库Libcint的用法

发表于 分类于

libcint是Sun Qiming写的电子积分库，有Fortran接口。我想以后使用它写点小玩意。现在先尝试搞懂它的用法，并分享自己的一点经验。

当然最重要的是看libcint自己的手册（和GitHub仓库里的例子），然后可以参考这篇文章：Libcint电子积分库使用教程

计算积分前要先将每一个primitive Gaussian函数（直接从基组文件里读取的）归一化，然后将重叠矩阵的对角元归一化。

调用CINTgto_norm函数，它接受primitive Gaussian函数的总角动量angl（笛卡尔型Gaussian函数的总角动量$l=l_x+l_y+l_z$）和指数项expnt（$\alpha$）。乘到primitive Gaussian函数对应的收缩系数cntr_coeff即可。

1
cntr_coeff = cntr_coeff * CINTgto_norm(angl, expnt)

我写了一个简单的测试（比较随意，有些地方是写死的），计算了水分子的def2-SVP基组的电子积分。

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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
module mod_info
    use iso_fortran_env, only: r8 => real64
    implicit none

    ! molecule: water
    ! basis set: def2-SVP

    integer, parameter :: natm = 3
    integer, parameter :: nshl = 12
    integer, parameter :: nprm = 22
    integer, parameter :: nbas = 24
  
    integer, dimension(nshl), parameter :: cntr_odr = &
    [5,1,1,3,1,1,3,1,1,3,1,1]
  
    integer, dimension(nshl), parameter :: angl = &
    [0,0,0,1,1,2,0,0,1,0,0,1]
  
    integer, dimension(nshl), parameter :: shl_belong_to_atom = &
    [1,1,1,1,1,1,2,2,2,3,3,3]
  
    integer, dimension(nshl), parameter :: sh_indx = &
    [1,2,3,4,7,10,15,16,17,20,21,22]
  
    real(kind=r8), dimension(nprm), parameter :: expnt = &
    [ 2266.1767785, 340.87010191, 77.363135167, &
      21.479644940, 6.6589433124, 0.80975975668, &
      0.25530772234, 17.721504317, 3.8635505440, &
      1.0480920883, 0.27641544411, 1.2000000, &
      13.0107010, 1.9622572, 0.44453796, &
      0.12194962, 0.8000000, 13.0107010, &
      1.9622572, 0.44453796, 0.12194962, &
      0.8000000 ]
  
    real(kind=r8), dimension(nprm), parameter :: coeff = &
    [ -0.53431809926E-02, -0.39890039230E-01, -0.17853911985, &
      -0.46427684959, -0.44309745172, 1.0000000, &
      1.0000000, 0.43394573193E-01, 0.23094120765, &
      0.51375311064, 1.0000000, 1.0000000, &
      0.19682158E-01, 0.13796524, 0.47831935, &
      1.0000000, 1.0000000, 0.19682158E-01, &
      0.13796524, 0.47831935, 1.0000000, &
      1.0000000 ]
    
    real(kind=r8), dimension(3,natm), parameter :: geom = reshape(&
    [ 0.00000000, -0.00000000, -0.11085125, &
      0.00000000, -0.78383672, 0.44340501, &
      0.00000000, 0.78383672, 0.44340501 ], [3,natm])
  
    integer, dimension(natm), parameter :: charge = [8,1,1]

end module mod_info

module mod_libcint
    use mod_info
    implicit none
  
    integer, dimension(:,:), allocatable :: atm
    integer, dimension(:,:), allocatable :: bas
    real(kind=r8), dimension(:), allocatable :: env
  
    integer, parameter :: CHARGE_OF  = 1
    integer, parameter :: PTR_COORD  = 2
    integer, parameter :: NUC_MOD_OF = 3
    integer, parameter :: PTR_ZETA   = 4
    integer, parameter :: ATM_SLOTS  = 6
    integer, parameter :: ATOM_OF    = 1
    integer, parameter :: ANG_OF     = 2
    integer, parameter :: NPRIM_OF   = 3
    integer, parameter :: NCTR_OF    = 4
    integer, parameter :: KAPPA_OF   = 5
    integer, parameter :: PTR_EXP    = 6
    integer, parameter :: PTR_COEFF  = 7
    integer, parameter :: BAS_SLOTS  = 8
    integer, parameter :: PTR_ENV_START = 20
  
contains
  
    subroutine set_libcint_input()
        integer :: off, prim_off
        integer :: iatom, ishl, iprim, ioff
        real(kind=r8), external :: CINTgto_norm
        
        allocate(atm(ATM_SLOTS,natm))
        allocate(bas(BAS_SLOTS,nshl))
        allocate(env(PTR_ENV_START+3*natm+2*nprm))
  
        off = PTR_ENV_START
        do iatom = 1, natm
            atm(CHARGE_OF,iatom) = charge(iatom)
            atm(PTR_COORD,iatom) = off
            atm(NUC_MOD_OF,iatom) = 1
            env(off+1:) = geom(:,iatom)
            off = off + 3
        end do
  
        prim_off = 1
        do ishl = 1, nshl
            bas(ATOM_OF,ishl) = shl_belong_to_atom(ishl) - 1
            bas(ANG_OF,ishl) = angl(ishl)
            bas(NPRIM_OF,ishl) = cntr_odr(ishl)
            bas(NCTR_OF,ishl) = 1
  
            bas(PTR_EXP,ishl) = off
            ioff = 0
            do iprim = prim_off, prim_off + cntr_odr(ishl)-1
                env(off+1+ioff) = expnt(iprim)
                ioff = ioff + 1
            end do
  
            off = off + cntr_odr(ishl)
  
            bas(PTR_COEFF,ishl) = off
            ioff = 0
            do iprim = prim_off, prim_off + cntr_odr(ishl)-1
                env(off+1+ioff) = coeff(iprim) * CINTgto_norm(angl(ishl), expnt(iprim))
                ioff = ioff + 1
            end do
  
            off = off + cntr_odr(ishl)
  
            prim_off = prim_off + cntr_odr(ishl)
        end do
  
    end subroutine set_libcint_input
  
    subroutine normalize()
        integer :: ishl, iprim
        integer :: di, dj
        real(kind=r8), dimension(:,:), allocatable :: buf1e
        integer, dimension(2) :: shls
        integer, external :: CINTcgto_spheric
  
        do ishl = 1, nshl
            shls(1) = ishl - 1
            shls(2) = ishl - 1
            di = CINTcgto_spheric(ishl-1, bas)
            dj = CINTcgto_spheric(ishl-1, bas)
            allocate(buf1e(di,dj))
            call cint1e_ovlp_sph(buf1e,shls,atm,natm,bas,nshl,env)
            do iprim = 1, cntr_odr(ishl)
                env(bas(PTR_COEFF,ishl)+iprim) = env(bas(PTR_COEFF,ishl)+iprim) / sqrt(buf1e(1,1))
            end do
            deallocate(buf1e)
        end do
  
    end subroutine normalize

end module mod_libcint

program main
    use mod_info
    use mod_libcint
    implicit none
    integer :: i, j, k, l, di, dj, dk, dl, x, y, z, w
    real(kind=r8), dimension(nbas,nbas) :: S
    real(kind=r8), dimension(nbas,nbas) :: T
    real(kind=r8), dimension(nbas,nbas) :: V
    real(kind=r8), dimension(nbas,nbas,nbas,nbas) :: eri
    real(kind=r8), dimension(:,:), allocatable :: buf1e
    real(kind=r8), dimension(:,:,:,:), allocatable :: buf2e
    integer, dimension(4) :: shls
    integer, external :: CINTcgto_spheric
  
    call set_libcint_input()
    call normalize()
  
    open(11, file="int.txt")
    write(11,*) atm(CHARGE_OF,:) 
    write(11,*) atm(PTR_COORD,:) 
    write(11,*) atm(NUC_MOD_OF,:)
    
    do i = 1, nshl
        write(11,"(*(I4))") bas(ATOM_OF,i), bas(ANG_OF,i), bas(NPRIM_OF,i), &
                            bas(NCTR_OF,i), bas(PTR_EXP,i), bas(PTR_COEFF,i)
    end do
  
    do i = 1, size(env)
        write(11,*) i, env(i)
    end do
    
    do i = 1, nshl
        shls(1) = i - 1
        di = CINTcgto_spheric(i-1, bas)
        do j = 1, nshl
            shls(2) = j - 1
            dj = CINTcgto_spheric(j-1, bas)
            allocate(buf1e(di,dj))
            x = sh_indx(i); y = sh_indx(j)
            call cint1e_ovlp_sph(buf1e,shls,atm,natm,bas,nshl,env)
            S(x:,y:) = buf1e(:,:)
            call cint1e_kin_sph(buf1e,shls,atm,natm,bas,nshl,env)
            T(x:,y:) = buf1e(:,:)
            call cint1e_nuc_sph(buf1e,shls,atm,natm,bas,nshl,env)
            V(x:,y:) = buf1e(:,:)
            deallocate(buf1e)
        end do
    end do
  
    do i = 1, nshl
      shls(1) = i - 1
      di = CINTcgto_spheric(i-1, bas)
      do j = 1, nshl
          shls(2) = j - 1
          dj = CINTcgto_spheric(j-1, bas)
          do k = 1, nshl
              shls(3) = k - 1
              dk = CINTcgto_spheric(k-1, bas)
              do l = 1, nshl
                  shls(4) = l - 1
                  dl = CINTcgto_spheric(l-1, bas)
                  allocate(buf2e(di,dj,dk,dl))
                  x = sh_indx(i); y = sh_indx(j)
                  z = sh_indx(k); w = sh_indx(l)
                  call cint2e_sph(buf2e,shls,atm,natm,bas,nshl,env,0_8)
                  eri(x:,y:,z:,w:) = buf2e(:,:,:,:)
                  deallocate(buf2e)
              end do
          end do
      end do
    end do
  
    do i = 1, nbas
        write(11,"(*(F15.10))") S(i,:)
    end do
    write(11,*)
    do i = 1, nbas
        write(11,"(*(F15.10))") T(i,:)
    end do
    write(11,*)
    do i = 1, nbas
        write(11,"(*(F15.10))") V(i,:)
    end do
    write(11,*)
    do i = 1, nbas
        do j = 1, nbas
            do k = 1, nbas
                write(11,"(*(F15.10))") eri(i,j,k,:)
            end do
        end do
    end do
    close(11)
    
end program main

我们可以用Psi4的电子积分做对比。有些矩阵元的顺序不一样，这个很正常（2p壳层的3个基函数的顺序不同）。

Psi4获得电子积分的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
import psi4

mol = psi4.geometry('''
    O  0.00000000  -0.00000000 -0.11085125
    H  0.00000000  -0.78383672  0.44340501
    H  0.00000000   0.78383672  0.44340501
    unit = au
''')
basis_set = "def2-SVP"

wfn = psi4.core.Wavefunction.build(mol, basis_set)
mints = psi4.core.MintsHelper(wfn.basisset())

# get integrals from psi4 module
Enuc = mol.nuclear_repulsion_energy()
S = np.asarray(mints.ao_overlap())
V = np.asarray(mints.ao_potential())
T = np.asarray(mints.ao_kinetic())
I = np.asarray(mints.ao_eri())