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

saving wok in ccGW

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-17 17:17:57 +02:00
parent 8b7c715811
commit 2fd282fcce
2 changed files with 256 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
src/GW/ccRG0W0.f90
View File

@ -31,8 +31,6 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
integer :: nSCF integer :: nSCF
double precision :: Conv double precision :: Conv
double precision :: x
double precision,allocatable :: eGW(:) double precision,allocatable :: eGW(:)
double precision,allocatable :: Z(:) double precision,allocatable :: Z(:)
@ -96,7 +94,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
do j=nC+1,nO do j=nC+1,nO
do a=nO+1,nOrb-nR do a=nO+1,nOrb-nR
del(i,j,a) = eHF(i) + eHF(j) - eHF(a) - eHF(p) del(i,j,a) = eHF(i) + eHF(j) - eHF(a)
end do end do
end do end do
@ -106,7 +104,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
do a=nO+1,nOrb-nR do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR do b=nO+1,nOrb-nR
del(b,a,i) = eHF(a) + eHF(b) - eHF(i) - eHF(p) del(b,a,i) = eHF(a) + eHF(b) - eHF(i)
end do end do
end do end do
@ -130,27 +128,13 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
nSCF = nSCF + 1 nSCF = nSCF + 1
! Compute intermediates x_2h1p and x_2p1h
x = 0d0
do q=nC+1,nOrb-nR
do r=nC+1,nOrb-nR
do s=nC+1,nOrb-nR
x = x + sqrt(2d0)*ERI(p,s,q,r)*amp(q,r,s)
end do
end do
end do
! Compute residual for 2h1p sector ! Compute residual for 2h1p sector
do i=nC+1,nO do i=nC+1,nO
do j=nC+1,nO do j=nC+1,nO
do a=nO+1,nOrb-nR do a=nO+1,nOrb-nR
res(i,j,a) = sqrt(2d0)*ERI(p,a,i,j) + (del(i,j,a) - x)*amp(i,j,a) res(i,j,a) = sqrt(2d0)*ERI(p,a,i,j) + (del(i,j,a) - eGW(p))*amp(i,j,a)
do k=nC+1,nO do k=nC+1,nO
do c=nO+1,nOrb-nR do c=nO+1,nOrb-nR
@ -170,7 +154,7 @@ subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eH
do a=nO+1,nOrb-nR do a=nO+1,nOrb-nR
do b=nO+1,nOrb-nR do b=nO+1,nOrb-nR
res(b,a,i) = sqrt(2d0)*ERI(p,i,b,a) + (del(b,a,i) - x)*amp(b,a,i) res(b,a,i) = sqrt(2d0)*ERI(p,i,b,a) + (del(b,a,i) - eGW(p))*amp(b,a,i)
do k=nC+1,nO do k=nC+1,nO
do c=nO+1,nOrb-nR do c=nO+1,nOrb-nR

252
src/GW/ccRG0W0.f90.bak Normal file
View File

@ -0,0 +1,252 @@
subroutine ccRG0W0(maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
! CC-based GW module
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF
double precision,intent(in) :: thresh
integer,intent(in) :: max_diis
integer,intent(in) :: nBas
integer,intent(in) :: 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(nOrb)
double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb)
! Local variables
integer :: p,q,r,s
integer :: i,j,k,l
integer :: a,b,c,d
integer :: nSCF
double precision :: Conv
double precision :: x
double precision,allocatable :: eGW(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: del(:,:,:)
double precision,allocatable :: res(:,:,:)
double precision,allocatable :: amp(:,:,:)
integer :: n_diis
double precision :: rcond
double precision,allocatable :: r_diis(:,:)
double precision,allocatable :: t_diis(:,:)
! Hello world
write(*,*)
write(*,*)'*****************************'
write(*,*)'* CC-based G0W0 Calculation *'
write(*,*)'*****************************'
write(*,*)
! Form energy denominator and guess amplitudes
allocate(del(nO,nV,nOrb))
allocate(res(nO,nV,nOrb))
allocate(amp(nO,nV,nOrb))
allocate(eGW(nOrb),Z(nOrb))
allocate(r_diis(nO*nV*nOrb,max_diis))
allocate(t_diis(nO*nV*nOrb,max_diis))
! Initialization
eGW(:) = eHF(:)
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=nO,nO
! Initialization
Conv = 1d0
nSCF = 0
n_diis = 0
t_diis(:,:) = 0d0
r_diis(:,:) = 0d0
rcond = 0d0
amp(:,:,:) = 0d0
res(:,:,:) = 0d0
!----------------------!
! Loop over amplitudes !
!----------------------!
write(*,*)
write(*,*)'-------------------------------------------------------------'
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','|'
write(*,*)'-------------------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Compute residual for 2h1p sector
do i=nC+1,nO
do a=1,nV-nR
do j=nC+1,nO
res(i,a,j) = sqrt(2d0)*ERI(p,nO+a,i,j) + (del(i,a,j) - eGW(p))*amp(i,a,j)
do k=nC+1,nO
do c=1,nV-nR
res(i,a,j) = res(i,a,j) - 2d0*ERI(j,nO+c,nO+a,k)*amp(i,c,k)
end do
end do
end do
end do
end do
! Compute residual for 2p1h sector
do i=nC+1,nO
do a=1,nV-nR
do b=1,nV-nR
res(i,a,nO+b) = sqrt(2d0)*ERI(p,i,nO+b,nO+a) + (del(i,a,nO+b) - eGW(p))*amp(i,a,nO+b)
do k=nC+1,nO
do c=1,nV-nR
res(i,a,nO+b) = res(i,a,nO+b) + 2d0*ERI(nO+a,k,i,nO+c)*amp(k,c,nO+b)
end do
end do
end do
end do
end do
! Check convergence
Conv = maxval(abs(res))
! Update amplitudes
amp(:,:,:) = amp(:,:,:) - res(:,:,:)/del(:,:,:)
! DIIS extrapolation
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nO*nV*nOrb,nO*nV*nOrb,n_diis,r_diis,t_diis,res,amp)
end if
! Compute quasiparticle energy
eGW(p) = eHF(p)
do i=nC+1,nOrb-nR
do a=1,nV-nR
do j=nC+1,nO
eGW(p) = eGW(p) + sqrt(2d0)*ERI(p,nO+a,i,j)*amp(i,a,j)
end do
end do
end do
do i=nC+1,nOrb-nR
do a=1,nV-nR
do b=1,nV-nR
eGW(p) = eGW(p) + sqrt(2d0)*ERI(p,i,nO+b,nO+a)*amp(i,a,nO+b)
end do
end do
end do
! 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,'|'
end do
write(*,*)'-------------------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
end if
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)' CC-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)','|'
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,'|'
write(*,*)'-------------------------------------------------------------------------------'
end do
end subroutine