mirror of
https://github.com/pfloos/quack
synced 2024-10-20 06:48:15 +02:00
clean up in ccG0W0
This commit is contained in:
parent
42f5439924
commit
0e43294211
@ -21,10 +21,10 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H)
|
||||
|
||||
H(:,:) = 0d0
|
||||
|
||||
do mu=1,nBas
|
||||
do si=1,nBas
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do si=1,nBas
|
||||
do mu=1,nBas
|
||||
H(mu,nu) = H(mu,nu) + P(la,si)*G(mu,la,nu,si)
|
||||
end do
|
||||
end do
|
||||
|
||||
@ -167,7 +167,6 @@ subroutine RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,maxSCF,thresh,max_dii
|
||||
|
||||
call wall_time(start_GW)
|
||||
call ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
! call ccRG0W0_mat(maxSCF,thresh,nBas,nOrb,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF)
|
||||
call wall_time(end_GW)
|
||||
|
||||
t_GW = end_GW - start_GW
|
||||
|
||||
@ -31,9 +31,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
integer :: nSCF
|
||||
double precision :: Conv
|
||||
|
||||
double precision :: x
|
||||
|
||||
double precision,allocatable :: eGW(:)
|
||||
double precision,allocatable :: Sig(:)
|
||||
double precision,allocatable :: Z(:)
|
||||
|
||||
double precision,allocatable :: del(:,:,:)
|
||||
@ -55,50 +53,29 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
write(*,*)'*****************************'
|
||||
write(*,*)
|
||||
|
||||
! Form energy denominator and guess amplitudes
|
||||
! Memory allocation
|
||||
|
||||
allocate(del(nO,nV,nOrb))
|
||||
allocate(vec(nO,nV,nOrb))
|
||||
allocate(res(nO,nV,nOrb))
|
||||
allocate(amp(nO,nV,nOrb))
|
||||
|
||||
allocate(eGW(nOrb),Z(nOrb))
|
||||
allocate(Sig(nOrb))
|
||||
allocate(Z(nOrb))
|
||||
|
||||
allocate(r_diis(nO*nV*nOrb,max_diis))
|
||||
allocate(t_diis(nO*nV*nOrb,max_diis))
|
||||
allocate(r_diis(nO*nV*nOrb,max_diis))
|
||||
allocate(t_diis(nO*nV*nOrb,max_diis))
|
||||
|
||||
! Initialization
|
||||
|
||||
eGW(:) = eHF(:)
|
||||
Sig(:) = 0d0
|
||||
Z(:) = 1d0
|
||||
|
||||
! Compute energy differences
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
|
||||
del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
!-------------------------!
|
||||
! Main loop over orbitals !
|
||||
!-------------------------!
|
||||
|
||||
do p=nC+1,nO+1
|
||||
do p=nO,nO+1
|
||||
|
||||
! Initialization
|
||||
|
||||
@ -115,6 +92,26 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
|
||||
! Compute energy differences
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
|
||||
del(i,a,j) = eHF(i) + eHF(j) - eHF(nO+a) - eHF(p)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
del(i,a,nO+b) = eHF(nO+a) + eHF(nO+b) - eHF(i) - eHF(p)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do j=nC+1,nO
|
||||
@ -144,7 +141,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
write(*,*)'| CC-based G0W0 calculation |'
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
|
||||
'|','#','|','HF','|','G0W0','|','Conv','|'
|
||||
'|','#','|','Sig_c (eV)','|','e_GW (eV)','|','Conv','|'
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
|
||||
do while(Conv > thresh .and. nSCF < maxSCF)
|
||||
@ -155,7 +152,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
|
||||
! Compute residual for 2h1p sector
|
||||
|
||||
res(:,:,:) = vec(:,:,:) + (del(:,:,:) - eGW(p))*amp(:,:,:)
|
||||
res(:,:,:) = vec(:,:,:) + (del(:,:,:) - Sig(p))*amp(:,:,:)
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
@ -197,7 +194,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
|
||||
! Update amplitudes
|
||||
|
||||
amp(:,:,:) = amp(:,:,:) - res(:,:,:)/(del(:,:,:) - eHF(p))
|
||||
amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:)
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
@ -210,13 +207,13 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
|
||||
! Compute quasiparticle energy
|
||||
|
||||
eGW(p) = eHF(p)
|
||||
Sig(p) = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do q=nC+1,nOrb-nR
|
||||
|
||||
eGW(p) = eGW(p) + vec(i,a,q)*amp(i,a,q)
|
||||
Sig(p) = Sig(p) + vec(i,a,q)*amp(i,a,q)
|
||||
|
||||
end do
|
||||
end do
|
||||
@ -225,11 +222,12 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
! Dump results
|
||||
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') &
|
||||
'|',nSCF,'|',eHF(p)*HaToeV,'|',eGW(p)*HaToeV,'|',Conv,'|'
|
||||
'|',nSCF,'|',Sig(p)*HaToeV,'|',(eHF(p)+Sig(p))*HaToeV,'|',Conv,'|'
|
||||
|
||||
end do
|
||||
|
||||
write(*,*)'-------------------------------------------------------------'
|
||||
write(*,*)
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
@ -244,19 +242,17 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
end if
|
||||
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,*)' CC-G0W0 calculation '
|
||||
write(*,*)'| CC-based G0W0 calculation |'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
|
||||
'|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|'
|
||||
'|','Orb','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X,F15.10,1X,A1,1X)') &
|
||||
'|',p,'|',eHF(p)*HaToeV,'|',(eGW(p)-eHF(p))*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|'
|
||||
'|',p,'|',eHF(p)*HaToeV,'|',Sig(p)*HaToeV,'|',Z(p),'|',(eHF(p)+Sig(p))*HaToeV,'|'
|
||||
write(*,*)'-------------------------------------------------------------------------------'
|
||||
|
||||
end do
|
||||
|
||||
@ -90,10 +90,10 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
|
||||
|
||||
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)
|
||||
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
|
||||
|
||||
@ -127,10 +127,10 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
|
||||
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
|
||||
! 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
|
||||
|
||||
@ -174,24 +174,24 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
|
||||
|
||||
! Diagonalize Fock matrix
|
||||
|
||||
if(nBas .eq. nOrb) then
|
||||
! 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
|
||||
! 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)))
|
||||
call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
|
||||
c(1,1), nBas, c(1,1), nBas, &
|
||||
0.d0, P(1,1), nBas)
|
||||
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
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user