LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-10-19 22:41:50 +02:00
Code Releases Activity

Merge pull request #10 from pfloos/mo_num

Browse Source 
Mo num
This commit is contained in:
AbdAmmar 2024-09-01 14:33:53 +02:00 committed by GitHub
commit fe1abef827
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
51 changed files with 1137 additions and 726 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
PyDuck.py
View File

@ -18,6 +18,7 @@ parser.add_argument('-b', '--basis', type=str, required=True, help='Name of the
parser.add_argument('--bohr', default='Angstrom', action='store_const', const='Bohr', help='By default QuAcK assumes that the xyz files are in Angstrom. Add this argument if your xyz file is in Bohr.')
parser.add_argument('-c', '--charge', type=int, default=0, help='Total charge of the molecule. Specify negative charges with "m" instead of the minus sign, for example m1 instead of -1. Default is 0')
parser.add_argument('--cartesian', default=False, action='store_true', help='Add this option if you want to use cartesian basis functions.')
parser.add_argument('--print_2e', default=False, action='store_true', help='Add this option if you want to print 2e-integrals.')
parser.add_argument('-fc', '--frozen_core', type=bool, default=False, help='Freeze core MOs. Default is false')
parser.add_argument('-m', '--multiplicity', type=int, default=1, help='Spin multiplicity. Default is 1 therefore singlet')
parser.add_argument('--working_dir', type=str, default=QuAcK_dir, help='Set a working directory to run the calculation.')
@ -32,6 +33,7 @@ frozen_core=args.frozen_core
multiplicity=args.multiplicity
xyz=args.xyz + '.xyz'
cartesian=args.cartesian
print_2e=args.print_2e
working_dir=args.working_dir
#Read molecule
@ -90,11 +92,10 @@ t1e = mol.intor('int1e_kin') #Kinetic energy matrix elements
dipole = mol.intor('int1e_r') #Matrix elements of the x, y, z operators
x,y,z = dipole[0],dipole[1],dipole[2]
norb = len(ovlp)
norb = len(ovlp) # nBAS_AOs
subprocess.call(['rm', working_dir + '/int/nBas.dat'])
f = open(working_dir+'/int/nBas.dat','w')
f.write(str(norb))
f.write(' ')
f.write(" {} ".format(str(norb)))
f.close()
@ -122,7 +123,6 @@ write_matrix_to_file(y,norb,working_dir+'/int/y.dat')
subprocess.call(['rm', working_dir + '/int/z.dat'])
write_matrix_to_file(z,norb,working_dir+'/int/z.dat')
#Write two-electron integrals
eri_ao = mol.intor('int2e')
def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
@ -132,12 +132,23 @@ def write_tensor_to_file(tensor,size,file,cutoff=1e-15):
for k in range(i,size):
for l in range(j,size):
if abs(tensor[i][k][j][l]) > cutoff:
#f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l]))
f.write(str(i+1)+' '+str(j+1)+' '+str(k+1)+' '+str(l+1)+' '+"{:.16E}".format(tensor[i][k][j][l]))
f.write('\n')
f.close()
subprocess.call(['rm', working_dir + '/int/ERI.dat'])
write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat')
# Write two-electron integrals
if print_2e:
# (formatted)
subprocess.call(['rm', working_dir + '/int/ERI.dat'])
write_tensor_to_file(eri_ao,norb,working_dir+'/int/ERI.dat')
else:
# (binary)
subprocess.call(['rm', working_dir + '/int/ERI.bin'])
# chem -> phys notation
eri_ao = eri_ao.transpose(0, 2, 1, 3)
eri_ao.tofile('int/ERI.bin')
#Execute the QuAcK fortran program
subprocess.call(QuAcK_dir+'/bin/QuAcK')

31
src/AOtoMO/AOtoMO.f90
View File

@ -1,26 +1,31 @@
subroutine AOtoMO(nBas,C,A,B)
subroutine AOtoMO(nBas, nOrb, C, M_AOs, M_MOs)
! Perform AO to MO transformation of a matrix A for given coefficients c
! Perform AO to MO transformation of a matrix M_AOs for given coefficients c
! M_MOs = C.T M_AOs C
implicit none
! Input variables
integer, intent(in) :: nBas, nOrb
double precision, intent(in) :: C(nBas,nOrb)
double precision, intent(in) :: M_AOs(nBas,nBas)
integer,intent(in) :: nBas
double precision,intent(in) :: C(nBas,nBas)
double precision,intent(in) :: A(nBas,nBas)
double precision, intent(out) :: M_MOs(nOrb,nOrb)
! Local variables
double precision, allocatable :: AC(:,:)
double precision,allocatable :: AC(:,:)
allocate(AC(nBas,nOrb))
! Output variables
!AC = matmul(M_AOs, C)
!M_MOs = matmul(transpose(C), AC)
double precision,intent(out) :: B(nBas,nBas)
call dgemm("N", "N", nBas, nOrb, nBas, 1.d0, &
M_AOs(1,1), nBas, C(1,1), nBas, &
0.d0, AC(1,1), nBas)
allocate(AC(nBas,nBas))
call dgemm("T", "N", nOrb, nOrb, nBas, 1.d0, &
C(1,1), nBas, AC(1,1), nBas, &
0.d0, M_MOs(1,1), nOrb)
AC = matmul(A,C)
B = matmul(transpose(C),AC)
deallocate(AC)
end subroutine

48
src/AOtoMO/AOtoMO_ERI_RHF.f90
View File

@ -1,4 +1,7 @@
subroutine AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
! ---
subroutine AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
@ -7,32 +10,51 @@ subroutine AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: c(nBas,nOrb)
! Local variables
double precision,allocatable :: scr(:,:,:,:)
integer :: mu,nu,la,si
integer :: i,j,k,l
double precision,allocatable :: a1(:,:,:,:)
double precision,allocatable :: a2(:,:,:,:)
! Output variables
double precision,intent(out) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(out) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
! Memory allocation
allocate(scr(nBas,nBas,nBas,nBas))
allocate(a2(nBas,nBas,nBas,nOrb))
allocate(a1(nBas,nBas,nOrb,nOrb))
! Four-index transform via semi-direct O(N^5) algorithm
call dgemm('T','N',nBas**3,nBas,nBas,1d0,ERI_AO,nBas,c(1,1),size(c,1),0d0,scr,nBas**3)
call dgemm('T','N',nBas**3,nBas,nBas,1d0,scr,nBas,c(1,1),size(c,1),0d0,ERI_MO,nBas**3)
call dgemm( 'T', 'N', nBas*nBas*nBas, nOrb, nBas, 1.d0 &
, ERI_AO(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, a2(1,1,1,1), nBas*nBas*nBas)
call dgemm('T','N',nBas**3,nBas,nBas,1d0,ERI_MO,nBas,c(1,1),size(c,1),0d0,scr,nBas**3)
call dgemm( 'T', 'N', nBas*nBas*nOrb, nOrb, nBas, 1.d0 &
, a2(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, a1(1,1,1,1), nBas*nBas*nOrb)
call dgemm('T','N',nBas**3,nBas,nBas,1d0,scr,nBas,c(1,1),size(c,1),0d0,ERI_MO,nBas**3)
deallocate(a2)
allocate(a2(nBas,nOrb,nOrb,nOrb))
call dgemm( 'T', 'N', nBas*nOrb*nOrb, nOrb, nBas, 1.d0 &
, a1(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, a2(1,1,1,1), nBas*nOrb*nOrb)
deallocate(a1)
call dgemm( 'T', 'N', nOrb*nOrb*nOrb, nOrb, nBas, 1.d0 &
, a2(1,1,1,1), nBas, c(1,1), nBas &
, 0.d0, ERI_MO(1,1,1,1), nOrb*nOrb*nOrb)
deallocate(a2)
end subroutine

42
src/AOtoMO/MOtoAO.f90
View File

@ -1,31 +1,39 @@
subroutine MOtoAO(nBas,S,C,B,A)
subroutine MOtoAO(nBas, nOrb, S, C, M_MOs, M_AOs)
! Perform MO to AO transformation of a matrix A for a given metric S
! and coefficients c
! Perform MO to AO transformation of a matrix M_AOs for a given metric S
! and coefficients c
!
! M_AOs = S C M_MOs (S C).T
implicit none
! Input variables
integer, intent(in) :: nBas, nOrb
double precision, intent(in) :: S(nBas,nBas)
double precision, intent(in) :: C(nBas,nOrb)
double precision, intent(in) :: M_MOs(nOrb,nOrb)
double precision, intent(out) :: M_AOs(nBas,nBas)
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: C(nBas,nBas)
double precision,intent(in) :: B(nBas,nBas)
double precision, allocatable :: SC(:,:),BSC(:,:)
! Local variables
double precision,allocatable :: SC(:,:),BSC(:,:)
allocate(SC(nBas,nOrb), BSC(nOrb,nBas))
! Output variables
!SC = matmul(S, C)
!BSC = matmul(M_MOs, transpose(SC))
!M_AOs = matmul(SC, BSC)
double precision,intent(out) :: A(nBas,nBas)
call dgemm("N", "N", nBas, nOrb, nBas, 1.d0, &
S(1,1), nBas, C(1,1), nBas, &
0.d0, SC(1,1), nBas)
! Memory allocation
call dgemm("N", "T", nOrb, nBas, nOrb, 1.d0, &
M_MOs(1,1), nOrb, SC(1,1), nBas, &
0.d0, BSC(1,1), nOrb)
allocate(SC(nBas,nBas),BSC(nBas,nBas))
call dgemm("N", "N", nBas, nBas, nOrb, 1.d0, &
SC(1,1), nBas, BSC(1,1), nOrb, &
0.d0, M_AOs(1,1), nBas)
SC = matmul(S,C)
BSC = matmul(B,transpose(SC))
A = matmul(SC,BSC)
deallocate(SC, BSC)
end subroutine

48
src/CC/RCC.f90
View File

@ -1,5 +1,8 @@
subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc,ERHF,eHF,cHF)
! ---
subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, &
maxSCF, thresh, max_diis, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_AO, ERI_MO, ENuc, ERHF, eHF, cHF)
! Coupled-cluster module
@ -24,18 +27,18 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
! Local variables
@ -48,11 +51,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(doCCD) then
call wall_time(start_CC)
call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call CCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CCD = ',t_CC,' seconds'
write(*,*)
end if
@ -64,12 +67,12 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(doDCD) then
call wall_time(start_CC)
call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, &
call DCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR, &
ERI_MO,ENuc,ERHF,eHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for DCD = ',t_CC,' seconds'
write(*,*)
end if
@ -83,11 +86,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(doCCSD) then
call wall_time(start_CC)
call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CCSD or CCSD(T)= ',t_CC,' seconds'
write(*,*)
end if
@ -99,11 +102,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(dodrCCD) then
call wall_time(start_CC)
call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call drCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for direct ring CCD = ',t_CC,' seconds'
write(*,*)
end if
@ -115,11 +118,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(dorCCD) then
call wall_time(start_CC)
call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call rCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for rCCD = ',t_CC,' seconds'
write(*,*)
end if
@ -131,11 +134,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(docrCCD) then
call wall_time(start_CC)
call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call crCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for crossed-ring CCD = ',t_CC,' seconds'
write(*,*)
end if
@ -147,11 +150,11 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(dolCCD) then
call wall_time(start_CC)
call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call lCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ladder CCD = ',t_CC,' seconds'
write(*,*)
end if
@ -163,12 +166,13 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d
if(dopCCD) then
call wall_time(start_CC)
call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERHF,eHF,cHF)
call pCCD(dotest, maxSCF, thresh, max_diis, nBas, nOrb, &
nC, nO, nV, nR, Hc, ERI_AO, ENuc, ERHF, eHF, cHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pair CCD = ',t_CC,' seconds'
write(*,*)
end if

46
src/CC/pCCD.f90
View File

@ -1,4 +1,8 @@
subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERHF,eHF,cHF)
! ---
subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
nC, nO, nV, nR, Hc, ERI_AO, ENuc, ERHF, eHF, cHF)
! pair CCD module
@ -12,15 +16,10 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERH
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
integer,intent(in) :: nBas, nOrb, nC, nO, nV, nR
double precision,intent(in) :: ENuc,ERHF
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
@ -94,18 +93,20 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERH
O = nO - nC
V = nV - nR
N = O + V
N = O + V ! nOrb - nC - nR
!------------------------------------!
! Star Loop for orbital optimization !
!------------------------------------!
allocate(ERI_MO(N,N,N,N))
allocate(c(N,N),h(N,N))
allocate(eO(O),eV(V),delta_OV(O,V))
allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V))
allocate(c(nBas,N), h(N,N))
allocate(eO(O), eV(V), delta_OV(O,V))
allocate(OOOO(O,O), OOVV(O,V), OVOV(O,V), OVVO(O,V), VVVV(V,V))
c(:,:) = cHF(nC+1:nBas-nR,nC+1:nBas-nR)
do i = 1, N
c(:,i) = cHF(:,nC+i)
enddo
CvgOrb = 1d0
nItOrb = 0
@ -121,13 +122,14 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERH
! Transform integrals
h = matmul(transpose(c),matmul(Hc(nC+1:nBas-nR,nC+1:nBas-nR),c))
call AOtoMO_ERI_RHF(N,c,ERI_AO(nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR),ERI_MO)
h = matmul(transpose(c), matmul(Hc, c))
call AOtoMO_ERI_RHF(nBas, N, c(1,1), ERI_AO(1,1,1,1), ERI_MO(1,1,1,1))
! Form energy denominator
eO(:) = eHF(nC+1:nO)
eV(:) = eHF(nO+1:nBas-nR)
eV(:) = eHF(nO+1:nOrb-nR)
do i=1,O
do a=1,V
@ -701,18 +703,18 @@ subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERH
deallocate(Kap)
write(*,*) 'e^kappa'
call matout(N,N,ExpKap)
call matout(N, N, ExpKap)
write(*,*)
write(*,*) 'Old orbitals'
call matout(N,N,c)
call matout(nBas, N, c)
write(*,*)
c = matmul(c,ExpKap)
c = matmul(c, ExpKap)
deallocate(ExpKap)
write(*,*) 'Rotated orbitals'
call matout(N,N,c)
call matout(nBas, N, c)
write(*,*)
end do

33
src/CI/RCI.f90
View File

@ -1,5 +1,8 @@
subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, &
epsHF,EHF,cHF,S)
! ---
subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nOrb, &
nC, nO, nV, nR, nS, ERI, dipole_int, epsHF, EHF)
! Configuration interaction module
@ -18,18 +21,16 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n
logical,intent(in) :: singlet
logical,intent(in) :: triplet
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: epsHF(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart)
! Local variables
@ -42,11 +43,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n
if(doCIS) then
call wall_time(start_CI)
call RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF)
call RCIS(dotest,singlet,triplet,doCIS_D,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CI,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CIS = ',t_CI,' seconds'
write(*,*)
end if
@ -58,11 +59,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n
if(doCID) then
call wall_time(start_CI)
call CID(dotest,singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF)
call CID(dotest,singlet,triplet,nOrb,nC,nO,nV,nR,ERI,epsHF,EHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CI,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CID = ',t_CI,' seconds'
write(*,*)
end if
@ -74,11 +75,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n
if(doCISD) then
call wall_time(start_CI)
call CISD(dotest,singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF)
call CISD(dotest,singlet,triplet,nOrb,nC,nO,nV,nR,ERI,epsHF,EHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CI,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CISD = ',t_CI,' seconds'
write(*,*)
end if
@ -91,11 +92,11 @@ subroutine RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,n
call wall_time(start_CI)
write(*,*) ' FCI is not yet implemented! Sorry.'
! call FCI(nBas,nC,nO,nV,nR,ERI,epsHF)
! call FCI(nOrb,nC,nO,nV,nR,ERI,epsHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_CI,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for FCI = ',t_CI,' seconds'
write(*,*)
end if

4
src/CI/RCIS.f90
View File

@ -41,7 +41,7 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in
! Memory allocation
allocate(A(nS,nS),Om(nS))
allocate(A(nS,nS), Om(nS))
! Compute CIS matrix
@ -117,4 +117,6 @@ subroutine RCIS(dotest,singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_in
end if
deallocate(A, Om)
end subroutine

4
src/GF/RG0F2.f90
View File

@ -51,7 +51,7 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
! Memory allocation
allocate(SigC(nBas),Z(nBas),eGFlin(nBas),eGF(nBas))
allocate(SigC(nBas), Z(nBas), eGFlin(nBas), eGF(nBas))
! Frequency-dependent second-order contribution
@ -133,4 +133,6 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
end if
deallocate(SigC, Z, eGFlin, eGF)
end subroutine

56
src/GF/RGF.f90
View File

@ -1,7 +1,10 @@
subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis, &
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
! ---
subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm, maxSCF, &
thresh, max_diis, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, linearize, &
eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, EHF, &
S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF)
! Green's function module
@ -39,7 +42,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -47,18 +50,18 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
integer,intent(in) :: nS
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: epsHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
@ -71,12 +74,13 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doG0F2) then
call wall_time(start_GF)
call RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call RG0F2(dotest, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, &
linearize, eta, regularize, nOrb, nC, nO, nV, nR, nS, &
ENuc, EHF, ERI_MO, dipole_int_MO, epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
@ -89,12 +93,12 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
call wall_time(start_GF)
call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI,dipole_int,epsHF)
singlet,triplet,linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI_MO,dipole_int_MO,epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
@ -106,12 +110,14 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doqsGF2) then
call wall_time(start_GF)
call qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
call qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, &
dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, &
rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, EHF, S, &
X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGF2 = ',t_GF,' seconds'
write(*,*)
end if
@ -123,11 +129,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doufG0F02) then
call wall_time(start_GF)
call ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF)
call ufRG0F02(dotest, nOrb, nC, nO, nV, nR, nS, ENuc, EHF, ERI_MO, epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0F02 = ',t_GF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0F02 = ',t_GF,' seconds'
write(*,*)
end if
@ -139,11 +145,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doG0F3) then
call wall_time(start_GF)
call RG0F3(dotest,renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
call RG0F3(dotest, renorm, nOrb, nC, nO, nV, nR, ERI_MO, epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF3 = ',t_GF,' seconds'
write(*,*)
end if
@ -155,11 +161,11 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm,max
if(doevGF3) then
call wall_time(start_GF)
call evRGF3(dotest,maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
call evRGF3(dotest, maxSCF, thresh, max_diis, renorm, nOrb, nC, nO, nV, nR, ERI_MO, epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF3 = ',t_GF,' seconds'
write(*,*)
end if

4
src/GF/evRGF2.f90
View File

@ -62,7 +62,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
! Memory allocation
allocate(SigC(nBas),Z(nBas),eGF(nBas),eOld(nBas),error_diis(nBas,max_diis),e_diis(nBas,max_diis))
allocate(SigC(nBas), Z(nBas), eGF(nBas), eOld(nBas), error_diis(nBas,max_diis), e_diis(nBas,max_diis))
! Initialization
@ -189,4 +189,6 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
end if
deallocate(SigC, Z, eGF, eOld, error_diis, e_diis)
end subroutine

24
src/GF/print_qsRGF2.f90
View File

@ -1,4 +1,8 @@
subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF,dipole)
! ---
subroutine print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, c, &
SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF, dipole)
! Print one-electron energies and other stuff for qsGF2
@ -7,17 +11,17 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
double precision,intent(in) :: Conv
double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGF(nBas)
double precision,intent(in) :: c(nBas)
double precision,intent(in) :: SigC(nBas,nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: eGF(nOrb)
double precision,intent(in) :: c(nBas,nOrb)
double precision,intent(in) :: SigC(nOrb,nOrb)
double precision,intent(in) :: Z(nOrb)
double precision,intent(in) :: ET
double precision,intent(in) :: EV
double precision,intent(in) :: EJ
@ -53,7 +57,7 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,
'|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do q=1,nBas
do q = 1, nOrb
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF(q)*HaToeV,'|'
end do
@ -102,12 +106,12 @@ subroutine print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigC,Z,ENuc,ET,EV,EJ,
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGF2 MO coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
call matout(nBas, nOrb, c)
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGF2 MO energies'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,eGF)
call matout(nOrb, 1, eGF)
write(*,*)
end if

103
src/GF/qsRGF2.f90
View File

@ -1,6 +1,10 @@
subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! ---
subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, &
dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, &
rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, &
S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Perform a quasiparticle self-consistent GF2 calculation
@ -29,25 +33,25 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
integer,intent(in) :: nBas,nOrb,nC,nO,nV,nR,nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_Sq
integer :: ispin
integer :: n_diis
double precision :: EqsGF2
@ -94,7 +98,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Stuff
nBasSq = nBas*nBas
nBas_Sq = nBas*nBas
! TDA
@ -105,9 +109,27 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Memory allocation
allocate(eGF(nBas),eOld(nbas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(eGF(nOrb))
allocate(eOld(nOrb))
allocate(c(nBas,nOrb))
allocate(cp(nOrb,nOrb))
allocate(Fp(nOrb,nOrb))
allocate(P(nBas,nBas))
allocate(F(nBas,nBas))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(error(nBas,nBas))
allocate(Z(nOrb))
allocate(SigC(nOrb,nOrb))
allocate(SigCp(nBas,nBas))
allocate(error_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Initialization
@ -117,7 +139,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
Conv = 1d0
P(:,:) = PHF(:,:)
eOld(:) = eHF(:)
eGF(:) = eHF(:)
eGF(:) = eHF(:)
c(:,:) = cHF(:,:)
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
@ -135,25 +157,25 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Buid Hartree matrix
call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J)
call Hartree_matrix_AO_basis(nBas, P, ERI_AO, J)
! Compute exchange part of the self-energy
call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
call exchange_matrix_AO_basis(nBas, P, ERI_AO, K)
! AO to MO transformation of two-electron integrals
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
! Compute self-energy and renormalization factor
if(regularize) then
call GF2_reg_self_energy(eta,nBas,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z)
call GF2_reg_self_energy(eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
else
call GF2_self_energy(eta,nBas,nC,nO,nV,nR,eGF,ERI_MO,SigC,Z)
call GF2_self_energy(eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
end if
@ -161,7 +183,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
SigC = 0.5d0*(SigC + transpose(SigC))
call MOtoAO(nBas,S,c,SigC,SigCp)
call MOtoAO(nBas, nOrb, S, c, SigC, SigCp)
! Solve the quasi-particle equation
@ -169,28 +191,27 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Compute commutator and convergence criteria
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
error = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
n_diis = min(n_diis+1, max_diis)
if(abs(rcond) > 1d-7) then
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F)
else
n_diis = 0
end if
! Diagonalize Hamiltonian in AO basis
Fp = matmul(transpose(X),matmul(F,X))
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGF)
c = matmul(X,cp)
SigCp = matmul(transpose(c),matmul(SigCp,c))
call diagonalize_matrix(nOrb, cp, eGF)
c = matmul(X, cp)
! Compute new density matrix in the AO basis
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
! Save quasiparticles energy for next cycle
@ -203,23 +224,23 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
ET = trace_matrix(nBas, matmul(P, T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
EV = trace_matrix(nBas, matmul(P, V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EJ = 0.5d0*trace_matrix(nBas, matmul(P, J))
! Exchange energy
Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
Ex = 0.25d0*trace_matrix(nBas, matmul(P, K))
! Correlation energy
call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EqsGF2,eGF,Ec)
call RMP2(.false., regularize, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec)
! Total energy
@ -230,8 +251,9 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
! Print results
!------------------------------------------------------------------------
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole)
call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole)
call print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, &
c, SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF2, dipole)
end do
!------------------------------------------------------------------------
@ -248,19 +270,21 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
stop
end if
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,error,error_diis,F_diis)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
! Perform BSE calculation
if(dophBSE) then
call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGF,EcBSE)
call GF2_phBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, &
nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -278,7 +302,8 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
if(doppBSE) then
call GF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI_MO,dipole_int_MO,eGF,EcBSE)
call GF2_ppBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

2
src/GF/qsUGF2.f90
View File

@ -197,7 +197,7 @@ subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved
end do
do is=1,nspin
call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
end do
! Solve the quasi-particle equation

61
src/GT/RGT.f90
View File

@ -1,7 +1,11 @@
subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis, &
doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int,PHF,cHF,eHF)
! ---
subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, &
maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, &
doppBSE, TDA_T, TDA, dBSE, dTDA, singlet, triplet, linearize, eta, regularize, &
nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, &
V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! T-matrix module
@ -44,7 +48,7 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -52,18 +56,18 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
@ -78,11 +82,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
call wall_time(start_GT)
call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0pp = ',t_GT,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0T0pp = ',t_GT,' seconds'
write(*,*)
end if
@ -95,11 +99,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
call wall_time(start_GT)
call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGTpp = ',t_GT,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGTpp = ',t_GT,' seconds'
write(*,*)
end if
@ -111,13 +115,13 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
if(doqsGTpp) then
call wall_time(start_GT)
call qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int, &
PHF,cHF,eHF)
call qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, dBSE, &
dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, &
nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGTpp = ',t_GT,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGTpp = ',t_GT,' seconds'
write(*,*)
end if
@ -129,11 +133,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
if(doufG0T0pp) then
call wall_time(start_GT)
call ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call ufG0T0pp(dotest,TDA_T,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0T0pp = ',t_GT,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0T0pp = ',t_GT,' seconds'
write(*,*)
end if
@ -146,11 +150,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
call wall_time(start_GT)
call RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0eh = ',t_GT,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0T0eh = ',t_GT,' seconds'
write(*,*)
end if
@ -163,11 +167,11 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
call wall_time(start_GT)
call evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
singlet,triplet,linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGTeh = ',t_GT,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGTeh = ',t_GT,' seconds'
write(*,*)
end if
@ -179,13 +183,14 @@ subroutine RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,do
if(doqsGTeh) then
call wall_time(start_GT)
call qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int, &
PHF,cHF,eHF)
call qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, &
dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, &
ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, &
Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGTeh = ',t_GT,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGTeh = ',t_GT,' seconds'
write(*,*)
end if

24
src/GT/print_qsRGTeh.f90
View File

@ -1,4 +1,8 @@
subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
! ---
subroutine print_qsRGTeh(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, &
Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole)
! Print one-electron energies and other stuff for qsGTeh
@ -7,7 +11,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
@ -19,11 +23,11 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
double precision,intent(in) :: EcRPA(nspin)
double precision,intent(in) :: Conv
double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: c(nBas)
double precision,intent(in) :: SigC(nBas,nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: eGT(nOrb)
double precision,intent(in) :: c(nBas,nOrb)
double precision,intent(in) :: SigC(nOrb,nOrb)
double precision,intent(in) :: Z(nOrb)
double precision,intent(in) :: EqsGT
double precision,intent(in) :: dipole(ncart)
@ -58,7 +62,7 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
'|','#','|','e_HF (eV)','|','Sig_GTeh (eV)','|','Z','|','e_GTeh (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do p=1,nBas
do p=1,nOrb
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|'
end do
@ -109,13 +113,13 @@ subroutine print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTeh MO coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
call matout(nBas, nOrb, c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTeh MO energies'
write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eGT)
call vecout(nOrb, eGT)
write(*,*)
end if

24
src/GT/print_qsRGTpp.f90
View File

@ -1,4 +1,8 @@
subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
! ---
subroutine print_qsRGTpp(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, Z, &
ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole)
! Print one-electron energies and other stuff for qsGT
@ -7,7 +11,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
@ -19,11 +23,11 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
double precision,intent(in) :: EcRPA(nspin)
double precision,intent(in) :: Conv
double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: c(nBas)
double precision,intent(in) :: SigC(nBas,nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: eGT(nOrb)
double precision,intent(in) :: c(nBas,nOrb)
double precision,intent(in) :: SigC(nOrb,nOrb)
double precision,intent(in) :: Z(nOrb)
double precision,intent(in) :: EqsGT
double precision,intent(in) :: dipole(ncart)
@ -58,7 +62,7 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
'|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do p=1,nBas
do p=1,nOrb
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|'
end do
@ -109,13 +113,13 @@ subroutine print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTpp MO coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
call matout(nBas, nOrb, c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' qsGTpp MO energies'
write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eGT)
call vecout(nOrb, eGT)
write(*,*)
end if

78
src/GT/qsRGTeh.f90
View File

@ -1,6 +1,10 @@
subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA, &
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! ---
subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, &
dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, &
ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, &
Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Perform a quasiparticle self-consistent GTeh calculation
@ -33,31 +37,31 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
logical :: dRPA = .false.
integer :: nSCF
integer :: nBasSq
integer :: nBas_Sq
integer :: ispin
integer :: n_diis
double precision :: ET
@ -113,7 +117,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Stuff
nBasSq = nBas*nBas
nBas_Sq = nBas*nBas
! TDA for T
@ -131,9 +135,29 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Memory allocation
allocate(Aph(nS,nS),Bph(nS,nS),eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS), &
rhoL(nBas,nBas,nS),rhoR(nBas,nBas,nS),err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(Aph(nS,nS), Bph(nS,nS), Om(nS), XpY(nS,nS), XmY(nS,nS))
allocate(eGT(nOrb))
allocate(eOld(nOrb))
allocate(Z(nOrb))
allocate(c(nBas,nOrb))
allocate(cp(nOrb,nOrb))
allocate(Fp(nOrb,nOrb))
allocate(Sig(nOrb,nOrb))
allocate(P(nBas,nBas))
allocate(F(nBas,nBas))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(Sigp(nBas,nBas))
allocate(err(nBas,nBas))
allocate(err_diis(nBas_Sq,max_diis), F_diis(nBas_Sq,max_diis))
allocate(rhoL(nOrb,nOrb,nS), rhoR(nOrb,nOrb,nS))
! Initialization
@ -169,12 +193,12 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! AO to MO transformation of two-electron integrals
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
! Compute linear response
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph)
if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -182,17 +206,17 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Compute correlation part of the self-energy
call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR)
call GTeh_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR)
if(regularize) call GTeh_regularization(nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR)
if(regularize) call GTeh_regularization(nOrb,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR)
call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z)
call GTeh_self_energy(eta,nOrb,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z)
! Make correlation self-energy Hermitian and transform it back to AO basis
Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO(nBas,S,c,Sig,Sigp)
call MOtoAO(nBas, nOrb, S, c, Sig, Sigp)
! Solve the quasi-particle equation
@ -207,7 +231,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
end if
@ -215,9 +239,8 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGT)
call diagonalize_matrix(nOrb, cp, eGT)
c = matmul(X,cp)
Sigp = matmul(transpose(c),matmul(Sigp,c))
! Compute new density matrix in the AO basis
@ -255,7 +278,8 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
call print_qsRGTeh(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, Sig, &
Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole)
end do
!------------------------------------------------------------------------
@ -272,13 +296,15 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, Om, XpY, XmY, rhoL, rhoR, err, err_diis, F_diis)
stop
end if
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,Om,XpY,XmY,rhoL,rhoR,err,err_diis,F_diis)
deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, Om, XpY, XmY, rhoL, rhoR, err, err_diis, F_diis)
! Perform BSE calculation

119
src/GT/qsRGTpp.f90
View File

@ -1,6 +1,9 @@
subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA, &
dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! ---
subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, &
dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, &
nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Perform a quasiparticle self-consistent GT calculation
@ -31,25 +34,26 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC,nO,nV,nR,nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_Sq
integer :: ispin
integer :: iblock
integer :: n_diis
@ -119,7 +123,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Stuff
nBasSq = nBas*nBas
nBas_Sq = nBas*nBas
! TDA for T
@ -137,16 +141,30 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Memory allocation
allocate(eGT(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),Sig(nBas,nBas),Sigp(nBas,nBas),Z(nBas), &
error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(eGT(nOrb))
allocate(eOld(nOrb))
allocate(Z(nOrb))
allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), &
Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
allocate(c(nBas,nOrb))
allocate(Fp(nOrb,nOrb))
allocate(cp(nOrb,nOrb))
allocate(Sig(nOrb,nOrb))
allocate(P(nBas,nBas))
allocate(F(nBas,nBas))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(error(nBas,nBas))
allocate(Sigp(nBas,nBas))
allocate(error_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
allocate(Om1s(nVVs), X1s(nVVs,nVVs), Y1s(nOOs,nVVs), rho1s(nOrb,nOrb,nVVs))
allocate(Om2s(nOOs), X2s(nVVs,nOOs), Y2s(nOOs,nOOs), rho2s(nOrb,nOrb,nOOs))
allocate(Om1t(nVVt), X1t(nVVt,nVVt), Y1t(nOOt,nVVt), rho1t(nOrb,nOrb,nVVt))
allocate(Om2t(nOOt), X2t(nVVt,nOOt), Y2t(nOOt,nOOt), rho2t(nOrb,nOrb,nOOt))
! Initialization
@ -182,7 +200,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! AO to MO transformation of two-electron integrals
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
! Compute linear response
@ -191,9 +209,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp)
call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin))
@ -204,9 +222,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp)
call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp)
call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp)
if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp)
call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin))
@ -217,25 +235,25 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Compute correlation part of the self-energy
iblock = 3
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
iblock = 3
call GTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
iblock = 4
call GTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
iblock = 4
call GTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
if(regularize) then
call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s)
call GTpp_regularization(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t)
call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s)
call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t)
end if
call GTpp_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, &
call GTpp_self_energy(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, &
Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z)
! Make correlation self-energy Hermitian and transform it back to AO basis
Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO(nBas,S,c,Sig,Sigp)
call MOtoAO(nBas, nOrb, S, c, Sig, Sigp)
! Solve the quasi-particle equation
@ -249,22 +267,21 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
n_diis = min(n_diis+1,max_diis)
if(abs(rcond) > 1d-7) then
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F)
else
n_diis = 0
end if
! Diagonalize Hamiltonian in AO basis
Fp = matmul(transpose(X),matmul(F,X))
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGT)
c = matmul(X,cp)
Sigp = matmul(transpose(c),matmul(Sigp,c))
call diagonalize_matrix(nOrb, cp, eGT)
c = matmul(X, cp)
! Compute new density matrix in the AO basis
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
! Save quasiparticles energy for next cycle
@ -298,7 +315,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGTpp(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,Sigp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole)
call print_qsRGTpp(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, &
eGT, c, Sig, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, &
EqsGT, dipole)
end do
!------------------------------------------------------------------------
@ -315,19 +334,24 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis)
deallocate(Om1s, X1s, Y1s, rho1s)
deallocate(Om2s, X2s, Y2s, rho2s)
deallocate(Om1t, X1t, Y1t, rho1t)
deallocate(Om2t, X2t, Y2t, rho2t)
stop
end if
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,Sig,Sigp,Z,error,error_diis,F_diis)
deallocate(c, cp, P, F, Fp, J, K, Sig, Sigp, Z, error, error_diis, F_diis)
! Perform BSE calculation
if(dophBSE) then
call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, &
Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, &
ERI_MO,dipole_int_MO,eGT,eGT,EcBSE)
@ -363,7 +387,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
end if
call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, &
nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, &
Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE)
@ -391,4 +415,9 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
end if
deallocate(Om1s, X1s, Y1s, rho1s)
deallocate(Om2s, X2s, Y2s, rho2s)
deallocate(Om1t, X1t, Y1t, rho1t)
deallocate(Om2t, X2t, Y2t, rho2t)
end subroutine

2
src/GT/qsUGTpp.f90
View File

@ -280,7 +280,7 @@ subroutine qsUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
end do
do ispin=1,nspin
call MOtoAO(nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin))
call MOtoAO(nBas,nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin))
end do
! Solve the quasi-particle equation

56
src/GW/RGW.f90
View File

@ -1,7 +1,10 @@
subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int,PHF,cHF,eHF)
! ---
subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF, thresh, max_diis, doACFDT, &
exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, &
linearize, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, &
S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Restricted GW module
@ -43,7 +46,7 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -51,18 +54,18 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
@ -76,11 +79,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for G0W0 = ',t_GW,' seconds'
write(*,*)
end if
@ -93,11 +96,11 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
call wall_time(start_GW)
call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF)
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGW = ',t_GW,' seconds'
write(*,*)
end if
@ -109,13 +112,14 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doqsGW) then
call wall_time(start_GW)
call qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, &
TDA_W, TDA, dBSE, dTDA, doppBSE, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, &
ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, &
dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds'
write(*,*)
end if
@ -127,13 +131,15 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doSRGqsGW) then
call wall_time(start_GW)
call SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, &
singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
dipole_int_AO,dipole_int,PHF,cHF,eHF)
call SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, &
dophBSE, dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, &
nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, &
ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, &
PHF, cHF, eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds'
write(*,*)
end if
@ -145,11 +151,12 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doufG0W0) then
call wall_time(start_GW)
! TODO
call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufG0W0 = ',t_GW,' seconds'
write(*,*)
end if
@ -161,11 +168,12 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thre
if(doufGW) then
call wall_time(start_GW)
! TODO
call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ufGW = ',t_GW,' seconds'
write(*,*)
end if

98
src/GW/SRG_qsGW.f90
View File

@ -1,6 +1,10 @@
subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA, &
dBSE,dTDA,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, &
S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! ---
subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, &
BSE, BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nNuc, &
ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, &
X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Perform a quasiparticle self-consistent GW calculation
@ -32,30 +36,30 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_Sq
integer :: ispin
integer :: ixyz
integer :: n_diis
@ -114,7 +118,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Stuff
nBasSq = nBas*nBas
nBas_Sq = nBas*nBas
! TDA for W
@ -132,9 +136,32 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Memory allocation
allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas),Aph(nS,nS),Bph(nS,nS), &
Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(eGW(nOrb))
allocate(eOld(nOrb))
allocate(Z(nOrb))
allocate(c(nBas,nOrb))
allocate(cp(nOrb,nOrb))
allocate(Fp(nOrb,nOrb))
allocate(SigC(nOrb,nOrb))
allocate(P(nBas,nBas))
allocate(F(nBas,nBas))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(error(nBas,nBas))
allocate(SigCp(nBas,nBas))
allocate(Aph(nS,nS))
allocate(Bph(nS,nS))
allocate(Om(nS))
allocate(XpY(nS,nS))
allocate(XmY(nS,nS))
allocate(rho(nOrb,nOrb,nS))
allocate(error_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Initialization
@ -174,22 +201,22 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
call wall_time(tao1)
do ixyz=1,ncart
call AOtoMO(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
do ixyz = 1, ncart
call AOtoMO(nBas, nOrb, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
call wall_time(tao2)
tao = tao + tao2 -tao1
tao = tao + tao2 - tao1
! Compute linear response
call wall_time(tlr1)
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
@ -203,13 +230,13 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
call wall_time(tex1)
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
call GW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho)
call wall_time(tex2)
tex=tex+tex2-tex1
call wall_time(tsrg1)
call SRG_self_energy(flow,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call SRG_self_energy(flow,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call wall_time(tsrg2)
@ -218,7 +245,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Make correlation self-energy Hermitian and transform it back to AO basis
call wall_time(tmo1)
call MOtoAO(nBas,S,c,SigC,SigCp)
call MOtoAO(nBas, nOrb, S, c, SigC, SigCp)
call wall_time(tmo2)
tmo = tmo + tmo2 - tmo1
! Solve the quasi-particle equation
@ -234,18 +261,18 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F)
end if
! Diagonalize Hamiltonian in AO basis
Fp = matmul(transpose(X),matmul(F,X))
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
call diagonalize_matrix(nOrb, cp, eGW)
c = matmul(X, cp)
call AOtoMO(nBas,c,SigCp,SigC)
call AOtoMO(nBas, nOrb, c, SigCp, SigC)
! Compute new density matrix in the AO basis
@ -283,7 +310,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole)
call print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, &
SigC, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGW, dipole)
end do
!------------------------------------------------------------------------
@ -300,6 +328,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis)
stop
end if
@ -313,17 +343,18 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Cumulant expansion
call RGWC(dotest,eta,nBas,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z)
call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z)
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,error,error_diis,F_diis)
deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis)
! Perform BSE calculation
if(BSE) then
call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE)
call GW_phBSE(BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE)
if(exchange_kernel) then
@ -357,7 +388,8 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
end if
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE)
call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, BSE, singlet, triplet, &
eta, nOrb, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

8
src/GW/SRG_qsUGW.f90
View File

@ -184,8 +184,8 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS
!--------------------------------------------------
do ixyz=1,ncart
call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block
@ -228,7 +228,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS
! Make correlation self-energy Hermitian and transform it back to AO basis
do is=1,nspin
call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
end do
! Solve the quasi-particle equation
@ -279,7 +279,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS
! Back-transform self-energy
do is=1,nspin
call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
end do
! Compute density matrix

24
src/GW/print_qsRGW.f90
View File

@ -1,4 +1,8 @@
subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
! ---
subroutine print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, &
Z, ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole)
! Print useful information about qsRGW calculation
@ -7,7 +11,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
integer,intent(in) :: nSCF
double precision,intent(in) :: ENuc
@ -19,11 +23,11 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
double precision,intent(in) :: EcRPA
double precision,intent(in) :: Conv
double precision,intent(in) :: thresh
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGW(nBas)
double precision,intent(in) :: c(nBas)
double precision,intent(in) :: SigC(nBas,nBas)
double precision,intent(in) :: Z(nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: eGW(nOrb)
double precision,intent(in) :: c(nBas,nOrb)
double precision,intent(in) :: SigC(nOrb,nOrb)
double precision,intent(in) :: Z(nOrb)
double precision,intent(in) :: EqsGW
double precision,intent(in) :: dipole(ncart)
@ -59,7 +63,7 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
'|','#','|','e_HF (eV)','|','Sig_GW (eV)','|','Z','|','e_GW (eV)','|'
write(*,*)'-------------------------------------------------------------------------------'
do p=1,nBas
do p=1,nOrb
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|'
end do
@ -110,13 +114,13 @@ subroutine print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigC,Z,ENuc,ET,EV,EJ,E
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Restricted qsGW orbital coefficients'
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,c)
call matout(nBas, nOrb, c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Restricted qsGW orbital energies (au) '
write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eGW)
call vecout(nOrb, eGW)
write(*,*)
end if

114
src/GW/qsRGW.f90
View File

@ -1,6 +1,10 @@
subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO, &
ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
! ---
subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, &
TDA_W, TDA, dBSE, dTDA, doppBSE, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, &
ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, &
ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
! Perform a quasiparticle self-consistent GW calculation
@ -34,15 +38,15 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
integer,intent(in) :: nS
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
@ -50,14 +54,14 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_Sq
integer :: ispin
integer :: ixyz
integer :: n_diis
@ -112,7 +116,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Stuff
nBasSq = nBas*nBas
nBas_Sq = nBas*nBas
! TDA for W
@ -130,10 +134,31 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Memory allocation
allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), &
J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),Z(nBas), &
Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), &
err(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(eGW(nOrb))
allocate(Z(nOrb))
allocate(c(nBas,nOrb))
allocate(cp(nOrb,nOrb))
allocate(Fp(nOrb,nOrb))
allocate(SigC(nOrb,nOrb))
allocate(P(nBas,nBas))
allocate(F(nBas,nBas))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(err(nBas,nBas))
allocate(SigCp(nBas,nBas))
allocate(Aph(nS,nS))
allocate(Bph(nS,nS))
allocate(Om(nS))
allocate(XpY(nS,nS))
allocate(XmY(nS,nS))
allocate(rho(nOrb,nOrb,nS))
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Initialization
@ -160,38 +185,38 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Build Hartree-exchange matrix
call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J)
call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
call Hartree_matrix_AO_basis(nBas, P, ERI_AO, J)
call exchange_matrix_AO_basis(nBas, P, ERI_AO, K)
! AO to MO transformation of two-electron integrals
do ixyz=1,ncart
call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
do ixyz = 1, ncart
call AOtoMO(nBas, nOrb, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
! Compute linear response
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph)
if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR_A(ispin, dRPA, nOrb, nC, nO, nV, nR, nS, 1d0, eGW, ERI_MO, Aph)
if(.not.TDA_W) call phLR_B(ispin, dRPA, nOrb, nC, nO, nV, nR, nS, 1d0, ERI_MO, Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
call phLR(TDA_W, nS, Aph, Bph, EcRPA, Om, XpY, XmY)
if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om)
! Compute correlation part of the self-energy
call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho)
call GW_excitation_density(nOrb, nC, nO, nR, nS, ERI_MO, XpY, rho)
if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eGW,Om,rho)
if(regularize) call GW_regularization(nOrb, nC, nO, nV, nR, nS, eGW, Om, rho)
call GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
call GW_self_energy(eta, nOrb, nC, nO, nV, nR, nS, eGW, Om, rho, EcGM, SigC, Z)
! Make correlation self-energy Hermitian and transform it back to AO basis
SigC = 0.5d0*(SigC + transpose(SigC))
call MOtoAO(nBas,S,c,SigC,SigCp)
call MOtoAO(nBas, nOrb, S, c, SigC, SigCp)
! Solve the quasi-particle equation
@ -205,19 +230,19 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
ET = trace_matrix(nBas, matmul(P, T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
EV = trace_matrix(nBas, matmul(P, V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EJ = 0.5d0*trace_matrix(nBas, matmul(P, J))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
EK = 0.25d0*trace_matrix(nBas, matmul(P, K))
! Total energy
@ -228,26 +253,27 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
end if
! Diagonalize Hamiltonian in AO basis
Fp = matmul(transpose(X),matmul(F,X))
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
call AOtoMO(nBas,c,SigCp,SigC)
call diagonalize_matrix(nOrb, cp, eGW)
c = matmul(X, cp)
call AOtoMO(nBas, nOrb, c, SigCp, SigC)
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
! Print results
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
call print_qsRGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,SigCp,Z,ENuc,ET,EV,EJ,EK,EcGM,EcRPA,EqsGW,dipole)
call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole)
call print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, Z, &
ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole)
end do
!------------------------------------------------------------------------
@ -264,19 +290,21 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, Om, XpY, XmY, rho, err, err_diis, F_diis)
stop
end if
! Deallocate memory
deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,Om,XpY,XmY,rho,err,err_diis,F_diis)
deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, Om, XpY, XmY, rho, err, err_diis, F_diis)
! Perform BSE calculation
if(dophBSE) then
call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE)
call GW_phBSE(dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, &
nOrb, nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE)
if(exchange_kernel) then
@ -310,7 +338,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
end if
call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcBSE)
call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, dophBSE, singlet, triplet, &
eta, nOrb, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
@ -327,7 +356,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
if(doppBSE) then
call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE)
call GW_ppBSE(TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, eGW, EcBSE)
EcBSE(2) = 3d0*EcBSE(2)

8
src/GW/qsUGW.f90
View File

@ -184,8 +184,8 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
!--------------------------------------------------
do ixyz=1,ncart
call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block
@ -232,7 +232,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
end do
do is=1,nspin
call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
end do
! Solve the quasi-particle equation
@ -283,7 +283,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
! Back-transform self-energy
do is=1,nspin
call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
end do
! Compute density matrix

4
src/HF/GHF.f90
View File

@ -123,7 +123,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Guess coefficients and density matrices
call mo_guess(nBas2,guess_type,S,H,X,C)
call mo_guess(nBas2,nBas2,guess_type,S,H,X,C)
! Construct super density matrix
@ -227,7 +227,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Level-shifting
if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,C,F)
if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nBas,nO,S,C,F)
! Transform Fock matrix in orthogonal basis

98
src/HF/RHF.f90
View File

@ -1,5 +1,8 @@
subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,eHF,c,P)
! ---
subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI, dipole_int, X, ERHF, eHF, c, P)
! Perform restricted Hartree-Fock calculation
@ -16,7 +19,7 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -26,14 +29,14 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_Sq
integer :: n_diis
double precision :: ET
double precision :: EV
@ -56,8 +59,8 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Output variables
double precision,intent(out) :: ERHF
double precision,intent(out) :: eHF(nBas)
double precision,intent(inout):: c(nBas,nBas)
double precision,intent(out) :: eHF(nOrb)
double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: P(nBas,nBas)
! Hello world
@ -70,24 +73,37 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Useful quantities
nBasSq = nBas*nBas
nBas_Sq = nBas*nBas
! Memory allocation
allocate(J(nBas,nBas),K(nBas,nBas),err(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), &
Fp(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(err(nBas,nBas))
allocate(F(nBas,nBas))
allocate(cp(nOrb,nOrb))
allocate(Fp(nOrb,nOrb))
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Guess coefficients and density matrix
call mo_guess(nBas,guess_type,S,Hc,X,c)
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
call mo_guess(nBas, nOrb, guess_type, S, Hc, X, c)
!P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO)))
call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
c(1,1), nBas, c(1,1), nBas, &
0.d0, P(1,1), nBas)
! Initialization
n_diis = 0
n_diis = 0
F_diis(:,:) = 0d0
err_diis(:,:) = 0d0
rcond = 0d0
rcond = 0d0
Conv = 1d0
nSCF = 0
@ -110,31 +126,35 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Build Fock matrix
call Hartree_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
call Hartree_matrix_AO_basis(nBas, P, ERI, J)
call exchange_matrix_AO_basis(nBas, P, ERI, K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
if(nBas .ne. nOrb) then
call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1))
call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1))
endif
! Check convergence
err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F)
if(nSCF > 1) Conv = maxval(abs(err))
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
ET = trace_matrix(nBas, matmul(P, T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
EV = trace_matrix(nBas, matmul(P, V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EJ = 0.5d0*trace_matrix(nBas, matmul(P, J))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
EK = 0.25d0*trace_matrix(nBas, matmul(P, K))
! Total energy
@ -144,25 +164,37 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F)
n_diis = min(n_diis+1, max_diis)
call DIIS_extrapolation(rcond, nBas_Sq, nBas_Sq, n_diis, err_diis, F_diis, err, F)
end if
! Level shift
if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F)
if(level_shift > 0d0 .and. Conv > thresh) then
call level_shifting(level_shift, nBas, nOrb, nO, S, c, F)
endif
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eHF)
c = matmul(X,cp)
if(nBas .eq. nOrb) then
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eHF)
c = matmul(X, cp)
else
Fp = matmul(transpose(c), matmul(F, c))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eHF)
c = matmul(c, cp)
endif
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
!P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO)))
call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
c(1,1), nBas, c(1,1), nBas, &
0.d0, P(1,1), nBas)
! Dump results
@ -185,14 +217,16 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(J, K, err, F, cp, Fp, err_diis, F_diis)
stop
end if
! Compute dipole moments
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_RHF(nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole)
call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int, dipole)
call print_RHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, EK, ERHF, dipole)
! Testing zone
@ -205,4 +239,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
end if
deallocate(J, K, err, F, cp, Fp, err_diis, F_diis)
end subroutine

53
src/HF/RHF_search.f90
View File

@ -1,6 +1,9 @@
subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO, &
X,ERHF,e,c,P)
! ---
subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nC, nO, nV, nR, S, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, &
X, ERHF, e, c, P)
! Search for RHF solutions
@ -13,7 +16,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -26,11 +29,11 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart)
double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart)
! Local variables
@ -59,8 +62,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Output variables
double precision,intent(out) :: ERHF
double precision,intent(out) :: e(nBas)
double precision,intent(inout):: c(nBas,nBas)
double precision,intent(out) :: e(nOrb)
double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: P(nBas,nBas)
! Memory allocation
@ -76,7 +79,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
!-------------------!
nS = (nO - nC)*(nV - nR)
allocate(Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS),R(nBas,nBas),ExpR(nBas,nBas))
allocate(Aph(nS,nS), Bph(nS,nS), AB(nS,nS), Om(nS))
allocate(R(nOrb,nOrb), ExpR(nOrb,nOrb))
!------------------!
! Search algorithm !
@ -92,12 +96,12 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
!---------------------!
call wall_time(start_HF)
call RHF(.false.,maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,e,c,P)
call RHF(.false., maxSCF, thresh, max_diis, guess, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, e, c, P)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
write(*,*)
!----------------------------------!
@ -108,14 +112,14 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
write(*,*)
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
do ixyz=1,ncart
call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
do ixyz = 1, ncart
call AOtoMO(nBas, nOrb, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
call wall_time(end_AOtoMO)
t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*)
!-------------------------------------------------------------!
@ -124,12 +128,12 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
ispin = 1
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phLR_A(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
call phLR_B(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
AB(:,:) = Aph(:,:) + Bph(:,:)
call diagonalize_matrix(nS,AB,Om)
call diagonalize_matrix(nS, AB, Om)
Om(:) = 2d0*Om(:)
write(*,*)'-------------------------------------------------------------'
@ -156,6 +160,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
if(eig < 0 .or. eig > nS) then
write(*,'(1X,A40,1X,A10)') 'Invalid option...','Stop...'
write(*,*)
deallocate(Aph, Bph, AB, Om, R, ExpR)
stop
end if
@ -164,15 +169,15 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
R(:,:) = 0d0
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
do a=nO+1,nOrb-nR
ia = ia + 1
R(a,i) = +AB(ia,eig)
R(i,a) = -AB(ia,eig)
end do
end do
call matrix_exponential(nBas,R,ExpR)
c = matmul(c,ExpR)
call matrix_exponential(nOrb, R, ExpR)
c = matmul(c, ExpR)
else
@ -191,4 +196,6 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
!---------------!
end do
deallocate(Aph, Bph, AB, Om, R, ExpR)
end subroutine

4
src/HF/RHF_stability.f90
View File

@ -30,7 +30,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
! Memory allocation
allocate(A(nS,nS),B(nS,nS),AB(nS,nS),Om(nS))
allocate(A(nS,nS), B(nS,nS), AB(nS,nS), Om(nS))
!-------------------------------------------------------------!
! Stability analysis: Real RHF -> Real RHF
@ -148,5 +148,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
end if
write(*,*)'-------------------------------------------------------------'
write(*,*)
deallocate(A, B, AB, Om)
end subroutine

2
src/HF/RMOM.f90
View File

@ -190,6 +190,6 @@ subroutine RMOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P)
EK = 0.5d0*trace_matrix(nBas,matmul(P,K))
ERHF = ET + EV + EJ + EK
call print_RHF(nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF)
call print_RHF(nBas,nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF)
end subroutine

95
src/HF/ROHF.f90
View File

@ -1,5 +1,8 @@
subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EROHF,eHF,c,Ptot)
! ---
subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI, dipole_int, X, EROHF, eHF, c, Ptot)
! Perform restricted open-shell Hartree-Fock calculation
@ -16,7 +19,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
double precision,intent(in) :: mix
double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -28,14 +31,14 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_Sq
integer :: n_diis
double precision :: Conv
double precision :: rcond
@ -62,8 +65,8 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Output variables
double precision,intent(out) :: EROHF
double precision,intent(out) :: eHF(nBas)
double precision,intent(inout):: c(nBas,nBas)
double precision,intent(out) :: eHF(nOrb)
double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: Ptot(nBas,nBas)
! Hello world
@ -76,19 +79,30 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Useful stuff
nBasSq = nBas*nBas
nBas_Sq = nBas*nBas
! Memory allocation
allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), &
P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), &
err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(J(nBas,nBas,nspin))
allocate(K(nBas,nBas,nspin))
allocate(F(nBas,nBas,nspin))
allocate(Ftot(nBas,nBas))
allocate(P(nBas,nBas,nspin))
allocate(err(nBas,nBas))
allocate(Fp(nOrb,nOrb))
allocate(cp(nOrb,nOrb))
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Guess coefficients and demsity matrices
call mo_guess(nBas,guess_type,S,Hc,X,c)
do ispin=1,nspin
P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin))))
call mo_guess(nBas, nOrb, guess_type, S, Hc, X, c)
do ispin = 1, nspin
!P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin))))
call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas)
end do
Ptot(:,:) = P(:,:,1) + P(:,:,2)
@ -120,51 +134,51 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Build Hartree repulsion
do ispin=1,nspin
call Hartree_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin))
do ispin = 1, nspin
call Hartree_matrix_AO_basis(nBas, P(:,:,ispin), ERI(:,:,:,:), J(:,:,ispin))
end do
! Compute exchange potential
do ispin=1,nspin
call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),K(:,:,ispin))
do ispin = 1, nspin
call exchange_matrix_AO_basis(nBas, P(:,:,ispin), ERI(:,:,:,:), K(:,:,ispin))
end do
! Build Fock operator
do ispin=1,nspin
do ispin = 1, nspin
F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin)
end do
call ROHF_fock_matrix(nBas,nO(1),nO(2),S,c,F(:,:,1),F(:,:,2),Ftot)
call ROHF_fock_matrix(nBas, nOrb, nO(1), nO(2), S, c, F(:,:,1), F(:,:,2), Ftot)
! Check convergence
err(:,:) = matmul(Ftot,matmul(Ptot,S)) - matmul(matmul(S,Ptot),Ftot)
err(:,:) = matmul(Ftot, matmul(Ptot, S)) - matmul(matmul(S, Ptot), Ftot)
if(nSCF > 1) Conv = maxval(abs(err(:,:)))
! Kinetic energy
do ispin=1,nspin
ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:)))
do ispin = 1, nspin
ET(ispin) = trace_matrix(nBas, matmul(P(:,:,ispin), T(:,:)))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
do ispin = 1, nspin
EV(ispin) = trace_matrix(nBas, matmul(P(:,:,ispin), V(:,:)))
end do
! Hartree energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
EJ(1) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,1), J(:,:,1)))
EJ(2) = trace_matrix(nBas, matmul(P(:,:,1), J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,2), J(:,:,2)))
! Exchange energy
do ispin=1,nspin
EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
do ispin = 1, nspin
EK(ispin) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,ispin), K(:,:,ispin)))
end do
! Total energy
@ -176,7 +190,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,Ftot)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,Ftot)
end if
@ -185,28 +199,29 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
if(level_shift > 0d0 .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(level_shift,nBas,maxval(nO),S,c,Ftot)
call level_shifting(level_shift, nBas, nOrb, maxval(nO), S, c, Ftot)
end do
end if
! Transform Fock matrix in orthogonal basis
Fp(:,:) = matmul(transpose(X(:,:)),matmul(Ftot(:,:),X(:,:)))
Fp(:,:) = matmul(transpose(X(:,:)), matmul(Ftot(:,:), X(:,:)))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eHF)
call diagonalize_matrix(nOrb, cp, eHF)
! Back-transform eigenvectors in non-orthogonal basis
c(:,:) = matmul(X(:,:),cp(:,:))
c(:,:) = matmul(X(:,:), cp(:,:))
! Compute density matrix
do ispin=1,nspin
P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin))))
do ispin = 1, nspin
!P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin))))
call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas)
end do
Ptot(:,:) = P(:,:,1) + P(:,:,2)
@ -231,6 +246,8 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis)
stop
end if
@ -238,7 +255,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Compute final UHF energy
call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,EK,EROHF,dipole)
call print_ROHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, EK, EROHF, dipole)
! Print test values
@ -248,4 +265,6 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
end if
deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis)
end subroutine

25
src/HF/ROHF_fock_matrix.f90
View File

@ -1,4 +1,7 @@
subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
! ---
subroutine ROHF_fock_matrix(nBas, nOrb, nOa, nOb, S, c, FaAO, FbAO, FAO)
! Construct the ROHF Fock matrix in the AO basis
! For open shells, the ROHF Fock matrix in the MO basis reads
@ -17,12 +20,12 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nOa
integer,intent(in) :: nOb
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: c(nBas,nOrb)
double precision,intent(inout):: FaAO(nBas,nBas)
double precision,intent(inout):: FbAO(nBas,nBas)
@ -46,7 +49,7 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
! Memory allocation
allocate(F(nBas,nBas),Fa(nBas,nBas),Fb(nBas,nBas))
allocate(F(nOrb,nOrb), Fa(nOrb,nOrb), Fb(nOrb,nOrb))
! Roothan canonicalization parameters
@ -61,14 +64,14 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
! Number of closed, open, and virtual orbitals
nC = min(nOa,nOb)
nC = min(nOa, nOb)
nO = abs(nOa - nOb)
nV = nBas - nC - nO
! Block-by-block Fock matrix
call AOtoMO(nBas,c,FaAO,Fa)
call AOtoMO(nBas,c,FbAO,Fb)
call AOtoMO(nBas, nOrb, c, FaAO, Fa)
call AOtoMO(nBas, nOrb, c, FbAO, Fb)
F(1:nC, 1:nC ) = aC*Fa(1:nC, 1:nC ) + bC*Fb(1:nC, 1:nC )
F(1:nC, nC+1:nC+nO ) = Fb(1:nC, nC+1:nC+nO )
@ -82,8 +85,10 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
F(nO+nC+1:nC+nO+nV, nC+1:nC+nO ) = Fa(nO+nC+1:nC+nO+nV, nC+1:nC+nO )
F(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) = aV*Fa(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) + bV*Fb(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV)
call MOtoAO(nBas,S,c,F,FAO)
call MOtoAO(nBas,S,c,Fa,FaAO)
call MOtoAO(nBas,S,c,Fb,FbAO)
call MOtoAO(nBas, nOrb, S, c, F, FAO)
call MOtoAO(nBas, nOrb, S, c, Fa, FaAO)
call MOtoAO(nBas, nOrb, S, c, Fb, FbAO)
deallocate(F, Fa, Fb)
end subroutine

4
src/HF/UHF.f90
View File

@ -85,7 +85,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Guess coefficients and demsity matrices
do ispin=1,nspin
call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin))
call mo_guess(nBas,nBas,guess_type,S,Hc,X,c(:,:,ispin))
P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
end do
@ -186,7 +186,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
if(level_shift > 0d0 .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(level_shift,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin))
call level_shifting(level_shift,nBas,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin))
end do
end if

4
src/HF/UHF_search.f90
View File

@ -124,8 +124,8 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block

6
src/HF/cRHF.f90
View File

@ -94,7 +94,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
! Guess coefficients and density matrix
call mo_guess(nBas,guess_type,S,Hc,X,c)
call mo_guess(nBas,nBas,guess_type,S,Hc,X,c)
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization
@ -166,7 +166,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
! Level shift
if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F)
if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nBas,nO,S,c,F)
! Diagonalize Fock matrix
@ -207,7 +207,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r
! Compute dipole moments
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_RHF(nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole)
call print_RHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! Testing zone

19
src/HF/core_guess.f90
View File

@ -1,4 +1,4 @@
subroutine core_guess(nBas,Hc,X,c)
subroutine core_guess(nBas, nOrb, Hc, X, c)
! Core guess of the molecular orbitals for HF calculation
@ -6,9 +6,9 @@ subroutine core_guess(nBas,Hc,X,c)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
! Local variables
@ -18,16 +18,19 @@ subroutine core_guess(nBas,Hc,X,c)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
double precision,intent(out) :: c(nBas,nOrb)
! Memory allocation
allocate(cp(nBas,nBas),e(nBas))
allocate(cp(nOrb,nOrb), e(nOrb))
! Core guess
cp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:)))
call diagonalize_matrix(nBas,cp,e)
c(:,:) = matmul(X(:,:),cp(:,:))
cp(:,:) = matmul(transpose(X(:,:)), matmul(Hc(:,:), X(:,:)))
call diagonalize_matrix(nOrb, cp, e)
c(:,:) = matmul(X(:,:), cp(:,:))
deallocate(cp, e)
end subroutine

16
src/HF/huckel_guess.f90
View File

@ -1,4 +1,4 @@
subroutine huckel_guess(nBas,S,Hc,X,c)
subroutine huckel_guess(nBas, nOrb, S, Hc, X, c)
! Hickel guess
@ -6,10 +6,10 @@ subroutine huckel_guess(nBas,S,Hc,X,c)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
! Local variables
@ -20,7 +20,7 @@ subroutine huckel_guess(nBas,S,Hc,X,c)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
double precision,intent(out) :: c(nBas,nOrb)
! Memory allocation
@ -32,9 +32,9 @@ subroutine huckel_guess(nBas,S,Hc,X,c)
! GWH approximation
do mu=1,nBas
do mu = 1, nBas
F(mu,mu) = Hc(mu,mu)
do nu=mu+1,nBas
do nu = mu+1, nBas
F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu))
F(nu,mu) = F(mu,nu)
@ -42,6 +42,8 @@ subroutine huckel_guess(nBas,S,Hc,X,c)
end do
end do
call core_guess(nBas,F,X,c)
call core_guess(nBas, nOrb, F, X, c)
deallocate(F)
end subroutine

15
src/HF/mo_guess.f90
View File

@ -1,4 +1,7 @@
subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
! ---
subroutine mo_guess(nBas, nOrb, guess_type, S, Hc, X, c)
! Guess of the molecular orbitals for HF calculation
@ -6,15 +9,15 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: guess_type
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
! Output variables
double precision,intent(inout) :: c(nBas,nBas)
double precision,intent(inout) :: c(nBas,nOrb)
if(guess_type == 0) then
@ -24,12 +27,12 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
elseif(guess_type == 1) then
write(*,*) 'Core guess...'
call core_guess(nBas,Hc,X,c)
call core_guess(nBas, nOrb, Hc, X, c)
elseif(guess_type == 2) then
write(*,*) 'Huckel guess...'
call huckel_guess(nBas,S,Hc,X,c)
call huckel_guess(nBas, nOrb, S, Hc, X, c)
elseif(guess_type == 3) then

15
src/HF/print_RHF.f90
View File

@ -1,4 +1,7 @@
subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! ---
subroutine print_RHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipole)
! Print one-electron energies and other stuff for G0W0
@ -7,10 +10,10 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET
double precision,intent(in) :: EV
@ -75,13 +78,13 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' RHF orbital coefficients '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,cHF)
call matout(nBas, nOrb, cHF)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' RHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eHF)
call vecout(nOrb, eHF)
write(*,*)
end subroutine

17
src/HF/print_ROHF.f90
View File

@ -1,14 +1,17 @@
subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
! ---
subroutine print_ROHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, Ex, EROHF, dipole)
! Print one- and two-electron energies and other stuff for RoHF calculation
implicit none
include 'parameters.h'
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: c(nBas,nOrb)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET(nspin)
double precision,intent(in) :: EV(nspin)
@ -31,7 +34,7 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
do ispin=1,nspin
if(nO(ispin) > 0) then
HOMO(ispin) = eHF(nO(ispin))
if(nO(ispin) < nBas) then
if(nO(ispin) < nOrb) then
LUMO(ispin) = eHF(nO(ispin)+1)
else
LUMO(ispin) = 0d0
@ -102,13 +105,13 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'ROHF orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,c)
call matout(nBas, nOrb, c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' ROHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eHF)
call vecout(nOrb, eHF)
write(*,*)
end subroutine

4
src/MP/RMP.f90
View File

@ -41,7 +41,7 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP2 = ',t_MP,' seconds'
write(*,*)
end if
@ -57,7 +57,7 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP2 = ',t_MP,' seconds'
write(*,*)
end if

116
src/QuAcK/QuAcK.f90
View File

@ -15,7 +15,7 @@ program QuAcK
logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW
logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh
integer :: nNuc,nBas
integer :: nNuc, nBas, nOrb
integer :: nC(nspin)
integer :: nO(nspin)
integer :: nV(nspin)
@ -31,6 +31,7 @@ program QuAcK
double precision,allocatable :: X(:,:)
double precision,allocatable :: dipole_int_AO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: Uvec(:,:), Uval(:)
double precision :: start_QuAcK,end_QuAcK,t_QuAcK
double precision :: start_int ,end_int ,t_int
@ -68,6 +69,10 @@ program QuAcK
logical :: dotest,doRtest,doUtest,doGtest
integer :: i, j, j0
double precision :: acc_d, acc_nd
double precision, allocatable :: tmp1(:,:), tmp2(:,:)
!-------------!
! Hello World !
!-------------!
@ -120,15 +125,16 @@ program QuAcK
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
!------------------------------------------------!
! Read input information !
!------------------------------------------------!
! nC = number of core orbitals !
! nO = number of occupied orbitals !
! nV = number of virtual orbitals (see below) !
! nR = number of Rydberg orbitals !
! nBas = number of basis functions (see below) !
!------------------------------------------------!
!---------------------------------------------------!
! Read input information !
!---------------------------------------------------!
! nC = number of core orbitals !
! nO = number of occupied orbitals !
! nV = number of virtual orbitals (see below) !
! nR = number of Rydberg orbitals !
! nBas = number of basis functions in AOs !
! nOrb = number of basis functions in MOs !
!---------------------------------------------------!
call read_molecule(nNuc,nO,nC,nR)
allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
@ -141,7 +147,7 @@ program QuAcK
! Read basis set information from PySCF !
!---------------------------------------!
call read_basis_pyscf(nBas,nO,nV)
call read_basis_pyscf(nBas, nO, nV)
!--------------------------------------!
! Read one- and two-electron integrals !
@ -149,26 +155,84 @@ program QuAcK
! Memory allocation for one- and two-electron integrals
allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), &
ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart))
allocate(S(nBas,nBas))
allocate(T(nBas,nBas))
allocate(V(nBas,nBas))
allocate(Hc(nBas,nBas))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
allocate(dipole_int_AO(nBas,nBas,ncart))
! Read integrals
call wall_time(start_int)
call read_integrals(nBas,S,T,V,Hc,ERI_AO)
call read_dipole_integrals(nBas,dipole_int_AO)
call read_integrals(nBas, S(1,1), T(1,1), V(1,1), Hc(1,1), ERI_AO(1,1,1,1))
call read_dipole_integrals(nBas, dipole_int_AO)
call wall_time(end_int)
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for reading integrals = ',t_int,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading integrals = ',t_int,' seconds'
write(*,*)
! Compute orthogonalization matrix
call orthogonalization_matrix(nBas,S,X)
!call orthogonalization_matrix(nBas, S, X)
allocate(Uvec(nBas,nBas), Uval(nBas))
Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas)
call diagonalize_matrix(nBas, Uvec, Uval)
nOrb = 0
do i = 1, nBas
if(Uval(i) > 1d-6) then
Uval(i) = 1d0 / dsqrt(Uval(i))
nOrb = nOrb + 1
else
write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization'
end if
end do
write(*,'(A38)') '--------------------------------------'
write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas
write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nOrb
write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nOrb)/dble(nBas))
write(*,'(A38)') '--------------------------------------'
write(*,*)
j0 = nBas - nOrb
allocate(X(nBas,nOrb))
do j = j0+1, nBas
do i = 1, nBas
X(i,j-j0) = Uvec(i,j) * Uval(j)
enddo
enddo
deallocate(Uvec, Uval)
!! check if X.T S X = 1_(nOrb,nOrb)
!allocate(tmp1(nOrb,nBas), tmp2(nOrb,nOrb))
!call dgemm("T", "N", nOrb, nBas, nBas, 1.d0, &
! X(1,1), nBas, S(1,1), nBas, &
! 0.d0, tmp1(1,1), nOrb)
!call dgemm("N", "N", nOrb, nOrb, nBas, 1.d0, &
! tmp1(1,1), nOrb, X(1,1), nBas, &
! 0.d0, tmp2(1,1), nOrb)
!acc_d = 0.d0
!acc_nd = 0.d0
!do i = 1, nOrb
! !write(*,'(1000(F15.7,2X))') (tmp2(i,j), j = 1, nOrb)
! acc_d = acc_d + tmp2(i,i)
! do j = 1, nOrb
! if(j == i) cycle
! acc_nd = acc_nd + dabs(tmp2(j,i))
! enddo
!enddo
!print*, ' diag part: ', dabs(acc_d - dble(nOrb)) / dble(nOrb)
!print*, ' non-diag part: ', acc_nd
!deallocate(tmp1, tmp2)
!---------------------!
! Choose QuAcK branch !
@ -200,7 +264,7 @@ program QuAcK
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
@ -216,7 +280,7 @@ program QuAcK
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
@ -228,13 +292,13 @@ program QuAcK
!--------------------------!
if(doGQuAcK) &
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, &
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, &
maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, &
maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, &
maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
!-----------!
@ -256,7 +320,7 @@ program QuAcK
call wall_time(end_QuAcK)
t_QuAcK = end_QuAcK - start_QuAcK
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for QuAcK = ',t_QuAcK,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for QuAcK = ',t_QuAcK,' seconds'
write(*,*)
end program

97
src/QuAcK/RQuAcK.f90
View File

@ -2,7 +2,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
@ -29,7 +29,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp
logical,intent(in) :: doG0T0eh,doevGTeh,doqsGTeh
integer,intent(in) :: nNuc,nBas
integer,intent(in) :: nNuc,nBas,nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -42,7 +42,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
@ -109,8 +109,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
! Memory allocation !
!-------------------!
allocate(cHF(nBas,nBas),eHF(nBas),PHF(nBas,nBas), &
dipole_int_MO(nBas,nBas,ncart),ERI_MO(nBas,nBas,nBas,nBas))
allocate(cHF(nBas,nOrb))
allocate(eHF(nOrb))
allocate(PHF(nBas,nBas))
allocate(dipole_int_MO(nOrb,nOrb,ncart))
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
!---------------------!
! Hartree-Fock module !
@ -119,12 +122,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doRHF) then
call wall_time(start_HF)
call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
call RHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
write(*,*)
end if
@ -132,12 +135,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doROHF) then
call wall_time(start_HF)
call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
call ROHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ROHF = ',t_HF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ROHF = ',t_HF,' seconds'
write(*,*)
end if
@ -154,18 +157,18 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
! Read and transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
do ixyz = 1, ncart
call AOtoMO(nBas, nOrb, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz))
end do
! 4-index transform
call AOtoMO_ERI_RHF(nBas,cHF,ERI_AO,ERI_MO)
call AOtoMO_ERI_RHF(nBas, nOrb, cHF, ERI_AO, ERI_MO)
call wall_time(end_AOtoMO)
t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*)
!-----------------------------------!
@ -177,11 +180,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(dostab) then
call wall_time(start_stab)
call RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO)
call RHF_stability(nOrb, nC, nO, nV, nR, nS, eHF, ERI_MO)
call wall_time(end_stab)
t_stab = end_stab - start_stab
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds'
write(*,*)
end if
@ -189,12 +192,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(dosearch) then
call wall_time(start_stab)
call RHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X,ERHF,eHF,cHF,PHF)
call RHF_search(maxSCF_HF, thresh_HF, max_diis_HF, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nC, nO, nV, nR, S, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, &
dipole_int_MO, X, ERHF, eHF, cHF, PHF)
call wall_time(end_stab)
t_stab = end_stab - start_stab
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds'
write(*,*)
end if
@ -208,11 +212,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doMP) then
call wall_time(start_MP)
call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
call RMP(dotest, doMP2, doMP3, reg_MP, nOrb, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, ERHF, eHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP = ',t_MP,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP = ',t_MP,' seconds'
write(*,*)
end if
@ -227,12 +231,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doCC) then
call wall_time(start_CC)
call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc,ERHF,eHF,cHF)
call RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, &
maxSCF_CC, thresh_CC, max_diis_CC, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_AO, ERI_MO, ENuc, ERHF, eHF, cHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CC = ',t_CC,' seconds'
write(*,*)
end if
@ -246,12 +250,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doCI) then
call wall_time(start_CI)
call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, &
eHF,ERHF,cHF,S)
call RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nOrb, &
nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, ERHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CI = ',t_CI,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CI = ',t_CI,' seconds'
write(*,*)
end if
@ -265,12 +269,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doRPA) then
call wall_time(start_RPA)
call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S)
call RRPA(dotest, dophRPA, dophRPAx, docrRPA, doppRPA, TDA, doACFDT, exchange_kernel, singlet, triplet, &
nOrb, nC, nO, nV, nR, nS, ENuc, ERHF, ERI_MO, dipole_int_MO, eHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPA = ',t_RPA,' seconds'
write(*,*)
end if
@ -284,14 +288,14 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGF) then
call wall_time(start_GF)
call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, &
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm_GF, maxSCF_GF, &
thresh_GF, max_diis_GF, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, lin_GF, &
eta_GF, reg_GF, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, &
S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
@ -305,14 +309,14 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGW) then
call wall_time(start_GW)
call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW, &
doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF_GW, thresh_GW, max_diis_GW, &
doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, &
lin_GW, eta_GW, reg_GW, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, &
V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GW = ',t_GW,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GW = ',t_GW,' seconds'
write(*,*)
end if
@ -326,14 +330,15 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGT) then
call wall_time(start_GT)
call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, &
maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, &
TDA_T,TDA,dBSE,dTDA,singlet,triplet,lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
call RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, &
maxSCF_GT, thresh_GT, max_diis_GT, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, &
TDA_T, TDA, dBSE, dTDA, singlet, triplet, lin_GT, eta_GT, reg_GT, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, &
dipole_int_MO, PHF, cHF, eHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GT = ',t_GT,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GT = ',t_GT,' seconds'
write(*,*)
end if

4
src/QuAcK/UQuAcK.f90
View File

@ -163,8 +163,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
! Read and transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block

12
src/RPA/RRPA.f90
View File

@ -1,5 +1,5 @@
subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,cHF,S)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
! Random-phase approximation module
@ -29,8 +29,6 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
@ -49,7 +47,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPA = ',t_RPA,' seconds'
write(*,*)
end if
@ -65,7 +63,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPAx = ',t_RPA,' seconds'
write(*,*)
end if
@ -81,7 +79,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pp-RPA = ',t_RPA,' seconds'
write(*,*)
end if
@ -97,7 +95,7 @@ subroutine RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_ker
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for pp-RPA = ',t_RPA,' seconds'
write(*,*)
end if

18
src/utils/level_shifting.f90
View File

@ -1,4 +1,4 @@
subroutine level_shifting(level_shift,nBas,nO,S,c,F)
subroutine level_shifting(level_shift, nBas, nOrb, nO, S, c, F)
! Perform level-shifting on the Fock matrix
@ -7,10 +7,10 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F)
! Input variables
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: c(nBas,nOrb)
! Local variables
@ -23,15 +23,17 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F)
double precision,intent(inout):: F(nBas,nBas)
allocate(F_MO(nBas,nBas),Sc(nBas,nBas))
allocate(F_MO(nOrb,nOrb), Sc(nBas,nOrb))
F_MO(:,:) = matmul(transpose(c),matmul(F,c))
F_MO(:,:) = matmul(transpose(c), matmul(F, c))
do a=nO+1,nBas
do a = nO+1, nOrb
F_MO(a,a) = F_MO(a,a) + level_shift
end do
Sc(:,:) = matmul(S,c)
F(:,:) = matmul(Sc,matmul(F_MO,transpose(Sc)))
Sc(:,:) = matmul(S, c)
F(:,:) = matmul(Sc, matmul(F_MO, transpose(Sc)))
deallocate(F_MO, Sc)
end subroutine

40
src/utils/orthogonalization_matrix.f90
View File

@ -1,4 +1,7 @@
subroutine orthogonalization_matrix(nBas,S,X)
! ---
subroutine orthogonalization_matrix(nBas, S, X)
! Compute the orthogonalization matrix X
@ -35,14 +38,32 @@ subroutine orthogonalization_matrix(nBas,S,X)
if(ortho_type == 1) then
!
! S V = V s where
!
! V.T V = 1 and s > 0 (S is positive def)
!
! S = V s V.T
! = V s^0.5 s^0.5 V.T
! = V s^0.5 V.T V s^0.5 V.T
! = S^0.5 S^0.5
!
! where
!
! S^0.5 = V s^0.5 V.T
!
! X = S^(-0.5)
! = V s^(-0.5) V.T
!
! write(*,*)
! write(*,*) ' Lowdin orthogonalization'
! write(*,*)
Uvec = S
call diagonalize_matrix(nBas,Uvec,Uval)
call diagonalize_matrix(nBas, Uvec, Uval)
do i=1,nBas
do i = 1, nBas
if(Uval(i) < thresh) then
@ -50,7 +71,7 @@ subroutine orthogonalization_matrix(nBas,S,X)
end if
Uval(i) = 1d0/sqrt(Uval(i))
Uval(i) = 1d0 / dsqrt(Uval(i))
end do
@ -63,13 +84,13 @@ subroutine orthogonalization_matrix(nBas,S,X)
! write(*,*)
Uvec = S
call diagonalize_matrix(nBas,Uvec,Uval)
call diagonalize_matrix(nBas, Uvec, Uval)
do i=1,nBas
do i = 1, nBas
if(Uval(i) > thresh) then
Uval(i) = 1d0/sqrt(Uval(i))
Uval(i) = 1d0 / dsqrt(Uval(i))
else
@ -79,7 +100,7 @@ subroutine orthogonalization_matrix(nBas,S,X)
end do
call AD(nBas,Uvec,Uval)
call AD(nBas, Uvec, Uval)
X = Uvec
elseif(ortho_type == 3) then
@ -117,3 +138,6 @@ subroutine orthogonalization_matrix(nBas,S,X)
end if
end subroutine
! ---

16
src/utils/read_basis_pyscf.f90
View File

@ -1,4 +1,4 @@
subroutine read_basis_pyscf(nBas,nO,nV)
subroutine read_basis_pyscf(nBas_AOs, nO, nV)
! Read basis set information from PySCF
@ -14,23 +14,23 @@ subroutine read_basis_pyscf(nBas,nO,nV)
! Output variables
integer,intent(out) :: nV(nspin)
integer,intent(out) :: nBas
integer,intent(out) :: nBas_AOs
!------------------------------------------------------------------------
! Primary basis set information
!------------------------------------------------------------------------
open(unit=3,file='int/nBas.dat')
read(3,*) nBas
read(3, *) nBas_AOs
close(unit=3)
write(*,'(A28)') '------------------'
write(*,'(A28,1X,I16)') 'Number of basis functions',nBas
write(*,'(A28)') '------------------'
write(*,*)
! write(*,'(A38)') '--------------------------------------'
! write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs
! write(*,'(A38)') '--------------------------------------'
! write(*,*)
! Number of virtual orbitals
nV(:) = nBas - nO(:)
nV(:) = nBas_AOs - nO(:)
end subroutine

71
src/utils/read_integrals.f90
View File

@ -1,4 +1,4 @@
subroutine read_integrals(nBas,S,T,V,Hc,G)
subroutine read_integrals(nBas_AOs, S, T, V, Hc, G)
! Read one- and two-electron integrals from files
@ -7,7 +7,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas_AOs
! Local variables
@ -18,11 +18,11 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
! Output variables
double precision,intent(out) :: S(nBas,nBas)
double precision,intent(out) :: T(nBas,nBas)
double precision,intent(out) :: V(nBas,nBas)
double precision,intent(out) :: Hc(nBas,nBas)
double precision,intent(out) :: G(nBas,nBas,nBas,nBas)
double precision,intent(out) :: S(nBas_AOs,nBas_AOs)
double precision,intent(out) :: T(nBas_AOs,nBas_AOs)
double precision,intent(out) :: V(nBas_AOs,nBas_AOs)
double precision,intent(out) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(out) :: G(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
! Open file with integrals
@ -35,7 +35,6 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
open(unit=8 ,file='int/Ov.dat')
open(unit=9 ,file='int/Kin.dat')
open(unit=10,file='int/Nuc.dat')
open(unit=11,file='int/ERI.dat')
open(unit=21,file='int/x.dat')
open(unit=22,file='int/y.dat')
@ -75,31 +74,29 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
Hc(:,:) = T(:,:) + V(:,:)
! Read nuclear integrals
! Read 2e-integrals
G(:,:,:,:) = 0d0
do
read(11,*,end=11) mu,nu,la,si,ERI
! ! formatted file
! open(unit=11, file='int/ERI.dat')
! G(:,:,:,:) = 0d0
! do
! read(11,*,end=11) mu, nu, la, si, ERI
! ERI = lambda*ERI
! G(mu,nu,la,si) = ERI ! <12|34>
! G(la,nu,mu,si) = ERI ! <32|14>
! G(mu,si,la,nu) = ERI ! <14|32>
! G(la,si,mu,nu) = ERI ! <34|12>
! G(si,mu,nu,la) = ERI ! <41|23>
! G(nu,la,si,mu) = ERI ! <23|41>
! G(nu,mu,si,la) = ERI ! <21|43>
! G(si,la,nu,mu) = ERI ! <43|21>
! end do
! 11 close(unit=11)
ERI = lambda*ERI
! <12|34>
G(mu,nu,la,si) = ERI
! <32|14>
G(la,nu,mu,si) = ERI
! <14|32>
G(mu,si,la,nu) = ERI
! <34|12>
G(la,si,mu,nu) = ERI
! <41|23>
G(si,mu,nu,la) = ERI
! <23|41>
G(nu,la,si,mu) = ERI
! <21|43>
G(nu,mu,si,la) = ERI
! <43|21>
G(si,la,nu,mu) = ERI
end do
11 close(unit=11)
! binary file
open(unit=11, file='int/ERI.bin', form='unformatted', access='stream')
read(11) G
close(11)
! Print results
@ -107,24 +104,24 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Overlap integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,S)
call matout(nBas_AOs,nBas_AOs,S)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Kinetic integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,T)
call matout(nBas_AOs,nBas_AOs,T)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Nuclear integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,V)
call matout(nBas_AOs,nBas_AOs,V)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Electron repulsion integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,G(1,1,la,si))
do la=1,nBas_AOs
do si=1,nBas_AOs
call matout(nBas_AOs, nBas_AOs, G(1,1,la,si))
end do
end do
write(*,*)

24
src/utils/utils.f90
View File

@ -65,7 +65,7 @@ subroutine diagonal_matrix(N,D,A)
end subroutine
!------------------------------------------------------------------------
subroutine matrix_exponential(N,A,ExpA)
subroutine matrix_exponential(N, A, ExpA)
! Compute Exp(A)
@ -81,7 +81,7 @@ subroutine matrix_exponential(N,A,ExpA)
! Memory allocation
allocate(W(N,N),tau(N),t(N,N))
allocate(W(N,N), tau(N), t(N,N))
! Initialize
@ -89,8 +89,8 @@ subroutine matrix_exponential(N,A,ExpA)
! Diagonalize
W(:,:) = - matmul(A,A)
call diagonalize_matrix(N,W,tau)
W(:,:) = - matmul(A, A)
call diagonalize_matrix(N, W, tau)
! do i=1,N
! tau(i) = max(abs(tau(i)),1d-14)
@ -99,16 +99,18 @@ subroutine matrix_exponential(N,A,ExpA)
! Construct cos part
call diagonal_matrix(N,cos(tau),t)
t(:,:) = matmul(t,transpose(W))
ExpA(:,:) = ExpA(:,:) + matmul(W,t)
call diagonal_matrix(N, cos(tau), t)
t(:,:) = matmul(t, transpose(W))
ExpA(:,:) = ExpA(:,:) + matmul(W, t)
! Construct sin part
call diagonal_matrix(N,sin(tau)/tau,t)
t(:,:) = matmul(t,transpose(W))
t(:,:) = matmul(t,A)
ExpA(:,:) = ExpA(:,:) + matmul(W,t)
call diagonal_matrix(N, sin(tau)/tau, t)
t(:,:) = matmul(t, transpose(W))
t(:,:) = matmul(t, A)
ExpA(:,:) = ExpA(:,:) + matmul(W, t)
deallocate(W, tau, t)
end subroutine