LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
Code Releases Activity

numgrid

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-21 16:31:39 +01:00
parent 97d5958add
commit 941ab79dc1
33 changed files with 623 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.CO
View File

@ -2,4 +2,4 @@
2 7 7 0 0
# Znuc x y z
C 0. 0. 0.
O 0. 0. 2.134
O 0. 0. 2.8

2
examples/molecule.F2
View File

@ -2,4 +2,4 @@
2 9 9 0 0
# Znuc x y z
F 0. 0. 0.
F 0. 0. 2
F 0. 0. 2.640

2
examples/molecule.H2
View File

@ -2,4 +2,4 @@
2 1 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.402
H 0. 0. 2.3

2
examples/molecule.LiF
View File

@ -2,4 +2,4 @@
2 6 6 0 0
# Znuc x y z
Li 0. 0. 0.
F 0. 0. 2.974
F 0. 0. 2.965

121
input/basis
View File

@ -1,39 +1,98 @@
1 10
S 8
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
S 8
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
1 15
S 9
1 6601.0000000 0.0001170
2 989.7000000 0.0009110
3 225.7000000 0.0047280
4 64.2900000 0.0191970
5 21.1800000 0.0630470
6 7.7240000 0.1632080
7 3.0030000 0.3148270
8 1.2120000 0.3939360
9 0.4930000 0.1969180
S 9
1 6601.0000000 -0.0000180
2 989.7000000 -0.0001420
3 225.7000000 -0.0007410
4 64.2900000 -0.0030200
5 21.1800000 -0.0101230
6 7.7240000 -0.0270940
7 3.0030000 -0.0573590
8 1.2120000 -0.0938950
9 0.4930000 -0.1210910
S 1
1 2.8360000 1.0000000
1 0.0951500 1.0000000
S 1
1 0.3782000 1.0000000
1 0.0479100 1.0000000
S 1
1 0.0222000 1.0000000
P 3
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
1 6.2500000 0.0033880
2 1.3700000 0.0193160
3 0.3672000 0.0791040
P 1
1 1.1430000 1.0000000
1 0.1192000 1.0000000
P 1
1 0.3300000 1.0000000
1 0.0447400 1.0000000
P 1
1 0.0179500 1.0000000
D 1
1 4.0140000 1.0000000
1 0.3440000 1.0000000
D 1
1 1.0960000 1.0000000
1 0.1530000 1.0000000
D 1
1 0.0680000 1.0000000
F 1
1 2.5440000 1.0000000
1 0.2460000 1.0000000
F 1
1 0.1292000 1.0000000
G 1
1 0.2380000 1.0000000
2 15
S 9
1 74530.0000000 0.0000950
2 11170.0000000 0.0007380
3 2543.0000000 0.0038580
4 721.0000000 0.0159260
5 235.9000000 0.0542890
6 85.6000000 0.1495130
7 33.5500000 0.3082520
8 13.9300000 0.3948530
9 5.9150000 0.2110310
S 9
1 74530.0000000 -0.0000220
2 11170.0000000 -0.0001720
3 2543.0000000 -0.0008910
4 721.0000000 -0.0037480
5 235.9000000 -0.0128620
6 85.6000000 -0.0380610
7 33.5500000 -0.0862390
8 13.9300000 -0.1558650
9 5.9150000 -0.1109140
S 1
1 1.8430000 1.0000000
S 1
1 0.7124000 1.0000000
S 1
1 0.2637000 1.0000000
P 3
1 80.3900000 0.0063470
2 18.6300000 0.0442040
3 5.6940000 0.1685140
P 1
1 1.9530000 1.0000000
P 1
1 0.6702000 1.0000000
P 1
1 0.2166000 1.0000000
D 1
1 5.0140000 1.0000000
D 1
1 1.7250000 1.0000000
D 1
1 0.5860000 1.0000000
F 1
1 3.5620000 1.0000000
F 1
1 1.1480000 1.0000000
G 1
1 2.3760000 1.0000000

4
input/methods
View File

@ -7,10 +7,10 @@
# CIS RPA RPAx ppRPA ADC
F F F F F
# G0F2 evGF2 G0F3 evGF3
T F F F
F F F F
# G0W0 evGW qsGW
F F F
# G0T0 evGT qsGT
F F F
T F F
# MCMP2
F

5
input/molecule
View File

@ -1,4 +1,5 @@
# nAt nEla nElb nCore nRyd
1 5 5 0 0
2 6 6 0 0
# Znuc x y z
Ne 0.0 0.0 0.0
Li 0. 0. 0.
F 0. 0. 2.965

5
input/molecule.xyz
View File

@ -1,3 +1,4 @@
1
2
Ne 0.0000000000 0.0000000000 0.0000000000
Li 0.0000000000 0.0000000000 0.0000000000
F 0.0000000000 0.0000000000 1.5690105433

4
input/options
View File

@ -7,10 +7,10 @@
# CIS/TDHF/BSE: singlet triplet
T T
# GF: maxSCF thresh DIIS n_diis lin renorm
256 0.00001 T 5 F 3
256 0.00001 T 5 T 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F T F F F T 0.000
# ACFDT: AC Kx XBS
T F F
T F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

121
input/weight
View File

@ -1,39 +1,98 @@
1 10
S 8
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
S 8
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
1 15
S 9
1 6601.0000000 0.0001170
2 989.7000000 0.0009110
3 225.7000000 0.0047280
4 64.2900000 0.0191970
5 21.1800000 0.0630470
6 7.7240000 0.1632080
7 3.0030000 0.3148270
8 1.2120000 0.3939360
9 0.4930000 0.1969180
S 9
1 6601.0000000 -0.0000180
2 989.7000000 -0.0001420
3 225.7000000 -0.0007410
4 64.2900000 -0.0030200
5 21.1800000 -0.0101230
6 7.7240000 -0.0270940
7 3.0030000 -0.0573590
8 1.2120000 -0.0938950
9 0.4930000 -0.1210910
S 1
1 2.8360000 1.0000000
1 0.0951500 1.0000000
S 1
1 0.3782000 1.0000000
1 0.0479100 1.0000000
S 1
1 0.0222000 1.0000000
P 3
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
1 6.2500000 0.0033880
2 1.3700000 0.0193160
3 0.3672000 0.0791040
P 1
1 1.1430000 1.0000000
1 0.1192000 1.0000000
P 1
1 0.3300000 1.0000000
1 0.0447400 1.0000000
P 1
1 0.0179500 1.0000000
D 1
1 4.0140000 1.0000000
1 0.3440000 1.0000000
D 1
1 1.0960000 1.0000000
1 0.1530000 1.0000000
D 1
1 0.0680000 1.0000000
F 1
1 2.5440000 1.0000000
1 0.2460000 1.0000000
F 1
1 0.1292000 1.0000000
G 1
1 0.2380000 1.0000000
2 15
S 9
1 74530.0000000 0.0000950
2 11170.0000000 0.0007380
3 2543.0000000 0.0038580
4 721.0000000 0.0159260
5 235.9000000 0.0542890
6 85.6000000 0.1495130
7 33.5500000 0.3082520
8 13.9300000 0.3948530
9 5.9150000 0.2110310
S 9
1 74530.0000000 -0.0000220
2 11170.0000000 -0.0001720
3 2543.0000000 -0.0008910
4 721.0000000 -0.0037480
5 235.9000000 -0.0128620
6 85.6000000 -0.0380610
7 33.5500000 -0.0862390
8 13.9300000 -0.1558650
9 5.9150000 -0.1109140
S 1
1 1.8430000 1.0000000
S 1
1 0.7124000 1.0000000
S 1
1 0.2637000 1.0000000
P 3
1 80.3900000 0.0063470
2 18.6300000 0.0442040
3 5.6940000 0.1685140
P 1
1 1.9530000 1.0000000
P 1
1 0.6702000 1.0000000
P 1
1 0.2166000 1.0000000
D 1
1 5.0140000 1.0000000
D 1
1 1.7250000 1.0000000
D 1
1 0.5860000 1.0000000
F 1
1 3.5620000 1.0000000
F 1
1 1.1480000 1.0000000
G 1
1 2.3760000 1.0000000

2
scan_H2.sh
View File

@ -1,7 +1,7 @@
#! /bin/bash
MOL="H2"
BASIS="cc-pvqz"
BASIS="cc-pvdz"
R_START=1.0
R_END=2.4
DR=0.1

8
scan_LiF.sh
View File

@ -1,10 +1,10 @@
#! /bin/bash
MOL="LiF"
BASIS="cc-pvqz"
R_START=2.965
R_END=2.965
DR=0.001
BASIS="cc-pvdz"
R_START=2.5
R_END=3.6
DR=0.1
for R in $(seq $R_START $DR $R_END)
do

2
scan_LiH.sh
View File

@ -1,7 +1,7 @@
#! /bin/bash
MOL="LiH"
BASIS="cc-pvqz"
BASIS="cc-pvdz"
R_START=2.5
R_END=3.6
DR=0.1

31
src/QuAcK/G0F2.f90
View File

@ -21,7 +21,7 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
double precision :: eps
double precision :: VV
double precision,allocatable :: eGF2(:)
double precision,allocatable :: Bpp(:,:)
double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:)
integer :: i,j,a,b,p
@ -36,13 +36,19 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
! Memory allocation
allocate(Bpp(nBas,2),Z(nBas),eGF2(nBas))
allocate(Sig(nBas),Z(nBas),eGF2(nBas))
if(linearize) then
! Frequency-dependent second-order contribution
write(*,*) '*** Quasiparticle equation will be linearized ***'
write(*,*)
Bpp(:,:) = 0d0
Z(:) = 0d0
end if
! Frequency-dependent second-order contribution
Sig(:) = 0d0
Z(:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
@ -51,8 +57,8 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
eps = e0(p) + e0(a) - e0(i) - e0(j)
VV = (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)
Bpp(p,1) = Bpp(p,1) + VV/eps
Z(p) = Z(p) + VV/eps**2
Sig(p) = Sig(p) + VV/eps
Z(p) = Z(p) + VV/eps**2
end do
end do
@ -66,8 +72,8 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
eps = e0(p) + e0(i) - e0(a) - e0(b)
VV = (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)
Bpp(p,2) = Bpp(p,2) + VV/eps
Z(p) = Z(p) + VV/eps**2
Sig(p) = Sig(p) + VV/eps
Z(p) = Z(p) + VV/eps**2
end do
end do
@ -78,15 +84,16 @@ subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
if(linearize) then
eGF2(:) = e0(:) + Z(:)*(Bpp(:,1) + Bpp(:,2))
eGF2(:) = e0(:) + Z(:)*Sig(:)
else
eGF2(:) = e0(:) + Bpp(:,1) + Bpp(:,2)
eGF2(:) = e0(:) + Sig(:)
end if
! Print results
call print_G0F2(nBas,nO,e0,eGF2,Z)
call print_G0F2(nBas,nO,e0,Sig,eGF2,Z)
end subroutine G0F2

10
src/QuAcK/G0T0.f90
View File

@ -71,7 +71,7 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Compute linear response
call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR, &
call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR, &
nOOs,nVVs,eHF(:),ERI(:,:,:,:), &
Omega1s(:),X1s(:,:),Y1s(:,:), &
Omega2s(:),X2s(:,:),Y2s(:,:), &
@ -96,7 +96,7 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Compute linear response
call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR, &
call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR, &
nOOt,nVVt,eHF(:),ERI(:,:,:,:), &
Omega1t(:),X1t(:,:),Y1t(:,:), &
Omega2t(:),X2t(:,:),Y2t(:,:), &
@ -113,6 +113,9 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
X1t(:,:),Y1t(:,:),rho1t(:,:,:), &
X2t(:,:),Y2t(:,:),rho2t(:,:,:))
rho2s(:,:,:) = 0d0
rho2t(:,:,:) = 0d0
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
@ -133,7 +136,8 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Solve the quasi-particle equation
!----------------------------------------------
eG0T0(:) = eHF(:) + Z(:)*SigT(:)
eG0T0(:) = eHF(:) + SigT(:)
! eG0T0(:) = eHF(:) + Z(:)*SigT(:)
!----------------------------------------------
! Dump results

5
src/QuAcK/QuAcK.f90
View File

@ -481,7 +481,7 @@ program QuAcK
if(doppRPA) then
call cpu_time(start_ppRPA)
call ppRPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF)
call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_ppRPA)
t_ppRPA = end_ppRPA - start_ppRPA
@ -636,7 +636,8 @@ program QuAcK
call cpu_time(start_G0T0)
call G0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
! call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0

14
src/QuAcK/evGF2.f90
View File

@ -28,7 +28,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
double precision :: rcond
double precision,allocatable :: eGF2(:)
double precision,allocatable :: eOld(:)
double precision,allocatable :: Bpp(:,:)
double precision,allocatable :: Sig(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: e_diis(:,:)
@ -45,7 +45,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
! Memory allocation
allocate(Bpp(nBas,2),Z(nBas),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
allocate(Sig(nBas),Z(nBas),eGF2(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
@ -65,7 +65,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
! Frequency-dependent second-order contribution
Bpp(:,:) = 0d0
Sig(:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
@ -74,7 +74,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
eps = eGF2(p) + e0(a) - e0(i) - e0(j)
Bpp(p,1) = Bpp(p,1) &
Sig(p) = Sig(p) &
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps
end do
@ -89,7 +89,7 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
eps = eGF2(p) + e0(i) - e0(a) - e0(b)
Bpp(p,2) = Bpp(p,2) &
Sig(p) = Sig(p) &
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps
end do
@ -133,11 +133,11 @@ subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
if(linearize) then
eGF2(:) = e0(:) + Z(:)*(Bpp(:,1) + Bpp(:,2))
eGF2(:) = e0(:) + Z(:)*Sig(:)
else
eGF2(:) = e0(:) + Bpp(:,1) + Bpp(:,2)
eGF2(:) = e0(:) + Sig(:)
end if

6
src/QuAcK/excitation_density_Tmatrix.f90
View File

@ -49,6 +49,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
! + ERI(p,i,c,d)*X1(cd,ab)
+ (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
end do
end do
@ -60,6 +61,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
! + ERI(p,i,k,l)*Y1(kl,ab)
+ (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
end do
end do
@ -116,6 +118,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
cd = cd + 1
rho1(p,i,ab) = rho1(p,i,ab) &
+ (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
print*,rho1(p,i,ab),ERI(p,i,c,d),X1(cd,ab)
end do
end do
@ -125,6 +129,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
kl = kl + 1
rho1(p,i,ab) = rho1(p,i,ab) &
+ (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
print*,rho1(p,i,ab),ERI(p,i,k,l),Y1(kl,ab)
end do
end do
@ -140,6 +145,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
end do
end do

1
src/QuAcK/linear_response_C_pp.f90
View File

@ -25,6 +25,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
! Define the chemical potential
eF = e(nO) + e(nO+1)
eF = 0d0
! Build C matrix for the singlet manifold

1
src/QuAcK/linear_response_D_pp.f90
View File

@ -25,6 +25,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
! Define the chemical potential
eF = e(nO) + e(nO+1)
eF = 0d0
! Build the D matrix for the singlet manifold

7
src/QuAcK/linear_response_pp.f90
View File

@ -1,4 +1,4 @@
subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA)
! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013)
@ -8,6 +8,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
! Input variables
logical,intent(in) :: ortho_eigvec
logical,intent(in) :: BSE
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR
integer,intent(in) :: nOO
@ -127,11 +128,13 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
! Split the various quantities in p-p and h-h parts
call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
call sort_ppRPA(ortho_eigvec,nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
! Compute the RPA correlation energy
! EcppRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
! print*,+sum(Omega1(:)),- trace_matrix(nVV,C(:,:))
! print*,-sum(Omega2(:)),- trace_matrix(nOO,D(:,:))
EcppRPA = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
! EcppRPA = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))

11
src/QuAcK/ppRPA.f90
View File

@ -1,4 +1,4 @@
subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e)
! Perform pp-RPA calculation
@ -7,7 +7,6 @@ subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc
! Input variables
logical,intent(in) :: doACFDT
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas
@ -62,7 +61,7 @@ subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc
allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), &
Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
call linear_response_pp(ispin,.false.,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
Ec_ppRPA(ispin))
@ -91,9 +90,9 @@ subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc
Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin))
call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
call linear_response_pp(ispin,.false.,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, &
Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), &
Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), &
Ec_ppRPA(ispin))
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin))

21
src/QuAcK/print_G0F2.f90
View File

@ -1,4 +1,4 @@
subroutine print_G0F2(nBas,nO,eHF,eGF2,Z)
subroutine print_G0F2(nBas,nO,eHF,Sig,eGF2,Z)
! Print one-electron energies and other stuff for G0F2
@ -8,6 +8,7 @@ subroutine print_G0F2(nBas,nO,eHF,eGF2,Z)
integer,intent(in) :: nBas
integer,intent(in) :: nO
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: Sig(nBas)
double precision,intent(in) :: eGF2(nBas)
double precision,intent(in) :: Z(nBas)
@ -24,23 +25,23 @@ subroutine print_G0F2(nBas,nO,eHF,eGF2,Z)
! Dump results
write(*,*)'--------------------------------------------------------'
write(*,*)'--------------------------------------------------------------------------'
write(*,*)' Frequency-dependent G0F2 calculation'
write(*,*)'--------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','e_G0F2 (eV)','|','Z','|'
write(*,*)'--------------------------------------------------------'
write(*,*)'--------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Sigma (eV)','|','e_G0F2 (eV)','|','Z','|'
write(*,*)'--------------------------------------------------------------------------'
do p=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',eGF2(p)*HaToeV,'|',Z(p),'|'
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',eGF2(p)*HaToeV,'|',Z(p),'|'
enddo
write(*,*)'--------------------------------------------------------'
write(*,*)'--------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'G0F2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0F2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0F2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'--------------------------------------------------------'
write(*,*)'--------------------------------------------------------------------------'
write(*,*)
end subroutine print_G0F2

11
src/QuAcK/print_G0T0.f90
View File

@ -21,7 +21,16 @@ subroutine print_G0T0(nBas,nO,e,ENuc,ERHF,SigT,Z,eGW,EcRPA)
HOMO = nO
LUMO = HOMO + 1
Gap = eGW(LUMO)-eGW(HOMO)
if(nBas >= LUMO) then
Gap = eGW(LUMO) - eGW(HOMO)
else
Gap = 0d0
end if
! Dump results

43
src/QuAcK/print_evGF2.f90
View File

@ -1,44 +1,53 @@
subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,eGF2)
subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,Sig,Z,eGF2)
! Print one-electron energies and other stuff for GF2
! Print one-electron energies and other stuff for G0F2
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO,nSCF
double precision,intent(in) :: Conv,eHF(nBas),eGF2(nBas)
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nSCF
double precision,intent(in) :: Conv
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: Sig(nBas)
double precision,intent(in) :: eGF2(nBas)
double precision,intent(in) :: Z(nBas)
integer :: x,HOMO,LUMO
integer :: p
integer :: HOMO
integer :: LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGF2(LUMO)-eGF2(HOMO)
Gap = eGF2(LUMO) - eGF2(HOMO)
! Dump results
write(*,*)'-------------------------------------------'
write(*,*)'--------------------------------------------------------------------------'
write(*,*)' Frequency-dependent evGF2 calculation'
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','e_evGF2 (eV)','|'
write(*,*)'-------------------------------------------'
write(*,*)'--------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Sigma (eV)','|','e_evGF2 (eV)','|','Z','|'
write(*,*)'--------------------------------------------------------------------------'
do x=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',eHF(x)*HaToeV,'|',eGF2(x)*HaToeV,'|'
do p=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',eGF2(p)*HaToeV,'|',Z(p),'|'
enddo
write(*,*)'-------------------------------------------'
write(*,*)'-------------------------------------------------------------'
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------'
write(*,*)'--------------------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'evGF2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------'
write(*,*)'--------------------------------------------------------------------------'
write(*,*)
end subroutine print_evGF2

6
src/QuAcK/print_evGF3.f90
View File

@ -8,7 +8,7 @@ subroutine print_evGF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
integer,intent(in) :: nBas,nO,nSCF
double precision,intent(in) :: Conv,eHF(nBas),eGF3(nBas),Z(nBas)
integer :: x,HOMO,LUMO
integer :: p,HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
@ -26,9 +26,9 @@ subroutine print_evGF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
'|','#','|','e_HF (eV)','|','Z','|','e_evGF3 (eV)','|'
write(*,*)'-------------------------------------------------------------'
do x=1,nBas
do p=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',eHF(x)*HaToeV,'|',Z(x),'|',eGF3(x)*HaToeV,'|'
'|',p,'|',eHF(p)*HaToeV,'|',Z(p),'|',eGF3(p)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------'

8
src/QuAcK/renormalization_factor_Tmatrix.f90
View File

@ -42,7 +42,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do i=nC+1,nO
do cd=1,nVVs
eps = e(p) + e(i) - Omega1s(cd)
Z(p) = Z(p) + 1d0*(rho1s(p,i,cd)/eps)**2
Z(p) = Z(p) + (rho1s(p,i,cd)/eps)**2
enddo
enddo
enddo
@ -53,7 +53,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do a=1,nV-nR
do kl=1,nOOs
eps = e(p) + e(nO+a) - Omega2s(kl)
Z(p) = Z(p) + 1d0*(rho2s(p,a,kl)/eps)**2
Z(p) = Z(p) + (rho2s(p,a,kl)/eps)**2
enddo
enddo
enddo
@ -68,7 +68,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do i=nC+1,nO
do cd=1,nVVt
eps = e(p) + e(i) - Omega1t(cd)
Z(p) = Z(p) + 1d0*(rho1t(p,i,cd)/eps)**2
Z(p) = Z(p) + (rho1t(p,i,cd)/eps)**2
enddo
enddo
enddo
@ -79,7 +79,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do a=1,nV-nR
do kl=1,nOOt
eps = e(p) + e(nO+a) - Omega2t(kl)
Z(p) = Z(p) + 1d0*(rho2t(p,a,kl)/eps)**2
Z(p) = Z(p) + (rho2t(p,a,kl)/eps)**2
enddo
enddo
enddo

8
src/QuAcK/self_energy_Tmatrix_diag.f90
View File

@ -46,7 +46,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do i=nC+1,nO
do cd=1,nVVs
eps = e(p) + e(i) - Omega1s(cd)
SigT(p) = SigT(p) + 1d0*rho1s(p,i,cd)**2/eps
SigT(p) = SigT(p) + rho1s(p,i,cd)**2/eps
enddo
enddo
enddo
@ -57,7 +57,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do a=1,nV-nR
do kl=1,nOOs
eps = e(p) + e(nO+a) - Omega2s(kl)
SigT(p) = SigT(p) + 1d0*rho2s(p,a,kl)**2/eps
SigT(p) = SigT(p) + rho2s(p,a,kl)**2/eps
enddo
enddo
enddo
@ -72,7 +72,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do i=nC+1,nO
do cd=1,nVVt
eps = e(p) + e(i) - Omega1t(cd)
SigT(p) = SigT(p) + 1d0*rho1t(p,i,cd)**2/eps
SigT(p) = SigT(p) + rho1t(p,i,cd)*rho1t(p,i,cd)/eps
enddo
enddo
enddo
@ -83,7 +83,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do a=1,nV-nR
do kl=1,nOOt
eps = e(p) + e(nO+a) - Omega2t(kl)
SigT(p) = SigT(p) + 1d0*rho2t(p,a,kl)**2/eps
SigT(p) = SigT(p) + rho2t(p,a,kl)*rho2t(p,a,kl)/eps
enddo
enddo
enddo

11
src/QuAcK/soG0T0.f90
View File

@ -79,8 +79,8 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Compute linear response
call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV, &
seHF,sERI,Omega1,X1,Y1,Omega2,X2,Y2, &
call linear_response_pp(ispin,.true.,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV, &
seHF,sERI,Omega1,X1,Y1,Omega2,X2,Y2, &
EcRPA)
call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:))
@ -92,6 +92,10 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
X1(:,:),Y1(:,:),rho1(:,:,:), &
X2(:,:),Y2(:,:),rho2(:,:,:))
rho2(:,:,:) = 0d0
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
@ -112,7 +116,8 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Solve the quasi-particle equation
!----------------------------------------------
eG0T0(:) = seHF(:) + Z(:)*SigT(:)
eG0T0(:) = seHF(:) + SigT(:)
! eG0T0(:) = seHF(:) + Z(:)*SigT(:)
!----------------------------------------------
! Dump results

31
src/QuAcK/sort_ppRPA.f90
View File

@ -1,4 +1,4 @@
subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! Compute the metric matrix for pp-RPA
@ -7,6 +7,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! Input variables
logical,intent(in) :: ortho_eigvec
integer,intent(in) :: nOO
integer,intent(in) :: nVV
double precision,intent(in) :: Omega(nOO+nVV)
@ -91,24 +92,28 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
write(*,*) 'pp-RPA positive excitation energies'
call matout(nVV,1,Omega1(:))
write(*,*)
! write(*,*) 'pp-RPA positive excitation energies'
! call matout(nVV,1,Omega1(:))
! write(*,*)
write(*,*) 'pp-RPA negative excitation energies'
call matout(nOO,1,Omega2(:))
write(*,*)
! write(*,*) 'pp-RPA negative excitation energies'
! call matout(nOO,1,Omega2(:))
! write(*,*)
! Orthogonalize eigenvectors
S1 = + matmul(transpose(Z1),matmul(M,Z1))
S2 = - matmul(transpose(Z2),matmul(M,Z2))
if(ortho_eigvec) then
call orthogonalization_matrix(1,nVV,S1,O1)
call orthogonalization_matrix(1,nOO,S2,O2)
S1 = + matmul(transpose(Z1),matmul(M,Z1))
S2 = - matmul(transpose(Z2),matmul(M,Z2))
Z1 = matmul(Z1,O1)
Z2 = matmul(Z2,O2)
if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
Z1 = matmul(Z1,O1)
Z2 = matmul(Z2,O2)
end if
X1(1:nVV,1:nVV) = Z1( 1: nVV,1:nVV)
Y1(1:nOO,1:nVV) = Z1(nVV+1:nOO+nVV,1:nVV)

19
src/eDFT/Makefile
View File

@ -14,7 +14,7 @@ ifeq ($(PROFIL),1)
FC += -p -fno-inline
endif
LIBS = ../../lib/*.a
LIBS = -lstdc++ ../../lib/*.a
LIBS += -lblas -llapack
SRCF90 = $(wildcard *.f90)
@ -23,15 +23,18 @@ SRC = $(wildcard *.F)
OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.F,$(ODIR)/%.o,$(SRC))
$(ODIR)/%.o: %.f90
$(FC) -c -o $@ $< $(FFLAGS)
$(ODIR)/%.o: %.F
$(FC) -c -o $@ $< $(FFLAGS)
$(BDIR)/eDFT: $(OBJ)
$(FC) -o $@ $^ $(FFLAGS) $(LIBS)
numgrid.mod: numgrid.f90
$(FC) -c -o $@ $< $(FFLAGS)
$(ODIR)/%.o: %.f90 numgrid.mod
$(FC) -c -o $@ $< $(FFLAGS)
$(ODIR)/%.o: %.F numgrid.mod
$(FC) -c -o $@ $< $(FFLAGS)
debug:
DEBUG=1 make $(BDIR)/eDFT
@ -39,4 +42,4 @@ profil:
PROFIL=1 make $(BDIR)/eDFT
clean:
rm -f $(ODIR)/*.o $(BDIR)/eDFT $(BDIR)/debug
rm -f $(ODIR)/*.o $(BDIR)/eDFT $(BDIR)/debug numgrid.mod

113
src/eDFT/numgrid.f90 Normal file
View File

@ -0,0 +1,113 @@
module numgrid
use, intrinsic :: iso_c_binding, only: c_ptr, c_double, c_int
implicit none
public numgrid_new_atom_grid
public numgrid_free_atom_grid
public numgrid_get_num_grid_points
public numgrid_get_num_radial_grid_points
public numgrid_get_grid
public numgrid_get_radial_grid
public numgrid_get_angular_grid
private
interface numgrid_new_atom_grid
function numgrid_new_atom_grid(radial_precision, &
min_num_angular_points, &
max_num_angular_points, &
proton_charge, &
alpha_max, &
max_l_quantum_number, &
alpha_min) result(context) bind (c)
import :: c_ptr, c_double, c_int
type(c_ptr) :: context
real(c_double), intent(in), value :: radial_precision
integer(c_int), intent(in), value :: min_num_angular_points
integer(c_int), intent(in), value :: max_num_angular_points
integer(c_int), intent(in), value :: proton_charge
real(c_double), intent(in), value :: alpha_max
integer(c_int), intent(in), value :: max_l_quantum_number
real(c_double), intent(in) :: alpha_min(*)
end function
end interface
interface numgrid_free_atom_grid
subroutine numgrid_free_atom_grid(context) bind (c)
import :: c_ptr
type(c_ptr), value :: context
end subroutine
end interface
interface numgrid_get_num_grid_points
function numgrid_get_num_grid_points(context) result(num_grid_points) bind (c)
import :: c_ptr, c_int
type(c_ptr), value :: context
integer(c_int) :: num_grid_points
end function
end interface
interface numgrid_get_num_radial_grid_points
function numgrid_get_num_radial_grid_points(context) result(num_radial_grid_points) bind (c)
import :: c_ptr, c_int
type(c_ptr), value :: context
integer(c_int) :: num_radial_grid_points
end function
end interface
interface numgrid_get_grid
subroutine numgrid_get_grid(context, &
num_centers, &
center_index, &
x_coordinates_bohr, &
y_coordinates_bohr, &
z_coordinates_bohr, &
proton_charges, &
grid_x_bohr, &
grid_y_bohr, &
grid_z_bohr, &
grid_w) bind (c)
import :: c_ptr, c_double, c_int
type(c_ptr), value :: context
integer(c_int), intent(in), value :: num_centers
integer(c_int), intent(in), value :: center_index
real(c_double), intent(in) :: x_coordinates_bohr(*)
real(c_double), intent(in) :: y_coordinates_bohr(*)
real(c_double), intent(in) :: z_coordinates_bohr(*)
integer(c_int), intent(in) :: proton_charges(*)
real(c_double), intent(inout) :: grid_x_bohr(*)
real(c_double), intent(inout) :: grid_y_bohr(*)
real(c_double), intent(inout) :: grid_z_bohr(*)
real(c_double), intent(inout) :: grid_w(*)
end subroutine
end interface
interface numgrid_get_radial_grid
subroutine numgrid_get_radial_grid(context, &
radial_grid_r_bohr,&
radial_grid_w) bind (c)
import :: c_ptr, c_double, c_int
type(c_ptr), value :: context
real(c_double), intent(inout) :: radial_grid_r_bohr(*)
real(c_double), intent(inout) :: radial_grid_w(*)
end subroutine
end interface
interface numgrid_get_angular_grid
subroutine numgrid_get_angular_grid(num_angular_grid_points, &
angular_grid_x_bohr, &
angular_grid_y_bohr, &
angular_grid_z_bohr, &
angular_grid_w) bind (c)
import :: c_double, c_int
integer(c_int), intent(in), value :: num_angular_grid_points
real(c_double), intent(inout) :: angular_grid_x_bohr(*)
real(c_double), intent(inout) :: angular_grid_y_bohr(*)
real(c_double), intent(inout) :: angular_grid_z_bohr(*)
real(c_double), intent(inout) :: angular_grid_w(*)
end subroutine
end interface
end module

156
src/eDFT/test_explicit.f90 Normal file
View File

@ -0,0 +1,156 @@
subroutine test_numgrid
! in this test we compute the grid for O
! in the presence of two H in the H2O molecule
use numgrid
use, intrinsic :: iso_c_binding, only: c_ptr
implicit none
real(8) :: radial_precision
integer :: min_num_angular_points
integer :: max_num_angular_points
integer :: num_centers
integer :: proton_charges(3)
real(8) :: x_coordinates_bohr(3)
real(8) :: y_coordinates_bohr(3)
real(8) :: z_coordinates_bohr(3)
real(8) :: alpha_max
integer :: max_l_quantum_number
real(8) :: alpha_min(3)
integer :: num_points
integer :: num_radial_points
integer :: num_angular_points
real(8), allocatable :: grid_x_bohr(:)
real(8), allocatable :: grid_y_bohr(:)
real(8), allocatable :: grid_z_bohr(:)
real(8), allocatable :: grid_w(:)
real(8), allocatable :: angular_grid_x_bohr(:)
real(8), allocatable :: angular_grid_y_bohr(:)
real(8), allocatable :: angular_grid_z_bohr(:)
real(8), allocatable :: angular_grid_w(:)
real(8), allocatable :: radial_grid_r_bohr(:)
real(8), allocatable :: radial_grid_w(:)
integer :: center_index
integer, parameter :: io_unit = 13
real(8) :: ref(4), own(4)
integer :: ipoint
integer :: j
real(8) :: error
type(c_ptr) :: context
radial_precision = 1.0d-12
min_num_angular_points = 86
max_num_angular_points = 302
num_centers = 3
proton_charges(1) = 8
proton_charges(2) = 1
proton_charges(3) = 1
x_coordinates_bohr(1) = 0.0d0
x_coordinates_bohr(2) = 1.43d0
x_coordinates_bohr(3) = -1.43d0
y_coordinates_bohr(1) = 0.0d0
y_coordinates_bohr(2) = 0.0d0
y_coordinates_bohr(3) = 0.0d0
z_coordinates_bohr(1) = 0.0d0
z_coordinates_bohr(2) = 1.1d0
z_coordinates_bohr(3) = 1.1d0
! oxygen
alpha_max = 11720.0d0
max_l_quantum_number = 2
alpha_min(1) = 0.3023d0
alpha_min(2) = 0.2753d0
alpha_min(3) = 1.185d0
open(unit=io_unit, file='../test/reference_grid/cc-pVDZ.txt', access='sequential', action='read')
context = numgrid_new_atom_grid(radial_precision, &
min_num_angular_points, &
max_num_angular_points, &
proton_charges(1), &
alpha_max, &
max_l_quantum_number, &
alpha_min)
num_points = numgrid_get_num_grid_points(context)
if (num_points /= 16364) stop 1
allocate(grid_x_bohr(num_points))
allocate(grid_y_bohr(num_points))
allocate(grid_z_bohr(num_points))
allocate(grid_w(num_points))
center_index = 0
call numgrid_get_grid(context, &
num_centers, &
center_index, &
x_coordinates_bohr, &
y_coordinates_bohr, &
z_coordinates_bohr, &
proton_charges, &
grid_x_bohr, &
grid_y_bohr, &
grid_z_bohr, &
grid_w)
do ipoint = 1, num_points
read(io_unit, *) ref(1), ref(2), ref(3), ref(4)
own(1) = grid_x_bohr(ipoint)
own(2) = grid_y_bohr(ipoint)
own(3) = grid_z_bohr(ipoint)
own(4) = grid_w(ipoint)
do j = 1, 4
error = own(j) - ref(j)
if (dabs(ref(j)) > 1.0e-15) error = error/ref(j)
if (dabs(error) > 1.0e-5) stop 1
end do
end do
close(io_unit)
deallocate(grid_x_bohr)
deallocate(grid_y_bohr)
deallocate(grid_z_bohr)
deallocate(grid_w)
num_radial_points = numgrid_get_num_radial_grid_points(context)
if (num_radial_points /= 106) stop 1
allocate(radial_grid_r_bohr(num_radial_points))
allocate(radial_grid_w(num_radial_points))
call numgrid_get_radial_grid(context, &
radial_grid_r_bohr, &
radial_grid_w)
deallocate(radial_grid_r_bohr)
deallocate(radial_grid_w)
num_angular_points = 14
allocate(angular_grid_x_bohr(num_angular_points))
allocate(angular_grid_y_bohr(num_angular_points))
allocate(angular_grid_z_bohr(num_angular_points))
allocate(angular_grid_w(num_angular_points))
call numgrid_get_angular_grid(num_angular_points, &
angular_grid_x_bohr, &
angular_grid_y_bohr, &
angular_grid_z_bohr, &
angular_grid_w)
deallocate(angular_grid_x_bohr)
deallocate(angular_grid_y_bohr)
deallocate(angular_grid_z_bohr)
deallocate(angular_grid_w)
call numgrid_free_atom_grid(context)
end subroutine