mirror of
https://github.com/pfloos/quack
synced 2025-05-06 15:24:43 +02:00
small modifs in GParquet
This commit is contained in:
parent
ea26528afc
commit
1c279e528e
@ -1,4 +1,4 @@
|
||||
2
|
||||
|
||||
Mg 0.0000 0.0000 0.0000
|
||||
O 0.0000 0.0000 1.728
|
||||
O 0.0000 0.0000 1.749
|
||||
|
||||
@ -7,8 +7,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
! Hard-coded parameters
|
||||
|
||||
logical :: linearize = .true.
|
||||
logical :: TDA = .true.
|
||||
logical :: TDA = .true.
|
||||
logical :: linearize = .true.
|
||||
logical :: print_phLR = .false.
|
||||
logical :: print_ppLR = .false.
|
||||
|
||||
@ -25,10 +25,10 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
integer :: n_it_1b,n_it_2b
|
||||
double precision :: err_1b,err_2b
|
||||
double precision :: err_eh, err_pp
|
||||
double precision :: err_eig_eh, err_eig_hh, err_eig_ee
|
||||
double precision :: start_t, end_t, t
|
||||
double precision :: start_1b, end_1b, t_1b
|
||||
double precision :: start_2b, end_2b, t_2b
|
||||
double precision :: err_eig_eh,err_eig_pp,err_eig_hh,err_eig_ee
|
||||
double precision :: start_t,end_t,t
|
||||
double precision :: start_1b,end_1b,t_1b
|
||||
double precision :: start_2b,end_2b,t_2b
|
||||
|
||||
integer :: nOO,nVV
|
||||
|
||||
@ -100,11 +100,11 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
write(*,*)'---------------------------------------------------------------'
|
||||
write(*,*)' Parquet parameters for one-body and two-body self-consistency '
|
||||
write(*,*)'---------------------------------------------------------------'
|
||||
write(*,'(1X,A50,1X,I5)') 'Maximum number for one-body self-consistency:', max_it_1b
|
||||
write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for one-body energies:', conv_1b
|
||||
write(*,'(1X,A50,1X,I5)') 'Maximum number of one-body iteration:',max_it_1b
|
||||
write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for one-body energies:',conv_1b
|
||||
write(*,*)'---------------------------------------------------------------'
|
||||
write(*,'(1X,A50,1X,I5)') 'Maximum number for two-body self-consistency:', max_it_2b
|
||||
write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for two-body energies:', conv_2b
|
||||
write(*,'(1X,A50,1X,I5)') 'Maximum number of two-body iteration:',max_it_2b
|
||||
write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for two-body energies:',conv_2b
|
||||
write(*,*)'---------------------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
@ -177,14 +177,10 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
! eh channel !
|
||||
!-----------------!
|
||||
|
||||
write(*,*)' ------------------------------'
|
||||
write(*,*)' | Diagonalizing ehBSE |'
|
||||
write(*,*)' ------------------------------'
|
||||
write(*,*)
|
||||
write(*,*) 'Diagonalizing eh BSE problem...'
|
||||
|
||||
allocate(Aph(nS,nS),Bph(nS,nS),eh_Om(nS),XpY(nS,nS),XmY(nS,nS),eh_Gam_A(nS,nS),eh_Gam_B(nS,nS))
|
||||
|
||||
|
||||
Aph(:,:) = 0d0
|
||||
Bph(:,:) = 0d0
|
||||
|
||||
@ -208,14 +204,12 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
Aph(:,:) = Aph(:,:) + eh_Gam_A(:,:)
|
||||
Bph(:,:) = Bph(:,:) + eh_Gam_B(:,:)
|
||||
|
||||
|
||||
|
||||
call phGLR(TDA,nS,Aph,Bph,EcRPA,eh_Om,XpY,XmY)
|
||||
|
||||
call wall_time(end_t)
|
||||
|
||||
t = end_t - start_t
|
||||
write(*,'(A50,1X,F9.3,A8)') 'Wall time for phBSE =',t,' seconds'
|
||||
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for phBSE problem =',t,' seconds'
|
||||
write(*,*)
|
||||
|
||||
if(print_phLR) call print_excitation_energies('phBSE@Parquet','eh generalized',nS,eh_Om)
|
||||
@ -228,10 +222,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
! pp channel !
|
||||
!-----------------!
|
||||
|
||||
write(*,*)' ------------------------------'
|
||||
write(*,*)' | Diagonalizing ppBSE |'
|
||||
write(*,*)' ------------------------------'
|
||||
write(*,*)
|
||||
write(*,*) 'Diagonalizing pp BSE problem...'
|
||||
|
||||
|
||||
allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), &
|
||||
ee_Om(nVV),X1(nVV,nVV),Y1(nOO,nVV), &
|
||||
@ -256,8 +248,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
else
|
||||
|
||||
if(.not.TDA) call G_pp_Gamma_B(nOrb,nC,nO,nR,nOO,nVV,old_eh_Phi,pp_Gam_B)
|
||||
call G_pp_Gamma_C(nOrb,nO,nR,nVV,old_eh_Phi,pp_Gam_C)
|
||||
call G_pp_Gamma_D(nOrb,nC,nO,nOO,old_eh_Phi,pp_Gam_D)
|
||||
call G_pp_Gamma_C(nOrb,nO,nR,nVV,old_eh_Phi,pp_Gam_C)
|
||||
call G_pp_Gamma_D(nOrb,nC,nO,nOO,old_eh_Phi,pp_Gam_D)
|
||||
|
||||
end if
|
||||
|
||||
@ -269,7 +261,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
call wall_time(end_t)
|
||||
t = end_t - start_t
|
||||
|
||||
write(*,'(A50,1X,F9.3,A8)') 'Wall time for ppBSE =',t,' seconds'
|
||||
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for ppBSE problem =',t,' seconds'
|
||||
write(*,*)
|
||||
|
||||
if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2p generalized',nVV,ee_Om)
|
||||
@ -277,18 +269,10 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
err_eig_ee = maxval(abs(old_ee_Om - ee_Om))
|
||||
err_eig_hh = maxval(abs(old_hh_Om - hh_Om))
|
||||
err_eig_pp = max(err_eig_ee,err_eig_hh)
|
||||
|
||||
deallocate(Bpp,Cpp,Dpp,pp_Gam_B,pp_Gam_C,pp_Gam_D)
|
||||
|
||||
|
||||
write(*,*) '----------------------------------------'
|
||||
write(*,*) ' Two-body (eigenvalue) convergence '
|
||||
write(*,*) '----------------------------------------'
|
||||
write(*,'(1X,A30,F10.6)')'Error for eh channel = ',err_eig_eh
|
||||
write(*,'(1X,A30,F10.6)')'Error for pp channel = ',max(err_eig_ee,err_eig_hh)
|
||||
write(*,*) '----------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
!----------!
|
||||
! Updating !
|
||||
!----------!
|
||||
@ -344,7 +328,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
allocate(pp_Phi(nOrb,nOrb,nOrb,nOrb))
|
||||
|
||||
! Build eh reducible kernels
|
||||
write(*,*) 'Computing eh reducible kernel ...'
|
||||
write(*,*) 'Computing eh reducible kernel...'
|
||||
|
||||
call wall_time(start_t)
|
||||
call G_eh_Phi(nOrb,nC,nR,nS,old_eh_Om,eh_rho,eh_Phi)
|
||||
@ -354,7 +338,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
write(*,*)
|
||||
|
||||
! Build pp reducible kernels
|
||||
write(*,*) 'Computing pp reducible kernel ...'
|
||||
write(*,*) 'Computing pp reducible kernel...'
|
||||
|
||||
call wall_time(start_t)
|
||||
call G_pp_Phi(nOrb,nC,nR,nOO,nVV,old_ee_Om,ee_rho,old_hh_Om,hh_rho,pp_Phi)
|
||||
@ -363,15 +347,15 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for pp reducible kernel =',t,' seconds'
|
||||
write(*,*)
|
||||
|
||||
! alpha = 0.05d0
|
||||
! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:)
|
||||
! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:)
|
||||
|
||||
err_eh = maxval(abs(eh_Phi - old_eh_Phi))
|
||||
err_pp = maxval(abs(pp_Phi - old_pp_Phi))
|
||||
|
||||
call matout(nOrb**2,nOrb**2,eh_Phi - old_eh_Phi)
|
||||
call matout(nOrb**2,nOrb**2,pp_Phi - old_pp_Phi)
|
||||
! alpha = 0.05d0
|
||||
! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:)
|
||||
! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:)
|
||||
|
||||
! call matout(nOrb**2,nOrb**2,eh_Phi - old_eh_Phi)
|
||||
! call matout(nOrb**2,nOrb**2,pp_Phi - old_pp_Phi)
|
||||
|
||||
!--------------------!
|
||||
! DIIS extrapolation !
|
||||
@ -425,12 +409,12 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
deallocate(eh_Phi,pp_Phi)
|
||||
|
||||
write(*,*) '----------------------------------------'
|
||||
write(*,*) ' Two-body (kernel) convergence '
|
||||
write(*,*) '----------------------------------------'
|
||||
write(*,'(1X,A30,F10.6)')'Error for eh channel = ',err_eh
|
||||
write(*,'(1X,A30,F10.6)')'Error for pp channel = ',err_pp
|
||||
write(*,*) '----------------------------------------'
|
||||
write(*,*) '----------------------------------------------'
|
||||
write(*,*) ' Two-body (frequency/kernel) convergence '
|
||||
write(*,*) '----------------------------------------------'
|
||||
write(*,'(1X,A24,F10.6,1X,F10.6)')'Error for eh channel = ',err_eig_eh,err_eh
|
||||
write(*,'(1X,A24,F10.6,1X,F10.6)')'Error for pp channel = ',err_eig_pp,err_pp
|
||||
write(*,*) '----------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
! Convergence criteria
|
||||
@ -438,7 +422,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
|
||||
|
||||
call wall_time(end_2b)
|
||||
t_2b = end_2b - start_2b
|
||||
write(*,'(A50,1X,F9.3,A8)') 'Wall time for two-body iteration =',t_2b,' seconds'
|
||||
write(*,'(A50,1X,I4,A2,F9.3,A8)') 'Wall time for two-body iteration #',n_it_2b,' =',t_2b,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end do
|
||||
|
||||
@ -37,9 +37,9 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho
|
||||
Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb))
|
||||
|
||||
rho(p,q,ia) = rho(p,q,ia) &
|
||||
+ (ERI(q,j,p,b) - ERI(q,j,b,p) &
|
||||
+ (1d0*ERI(q,j,p,b) - 1d0*ERI(q,j,b,p) &
|
||||
- 1d0*eh_Phi(q,j,b,p) + 1d0*pp_Phi(q,j,p,b)) * X &
|
||||
+ (ERI(q,b,p,j) - ERI(q,b,j,p) &
|
||||
+ (1d0*ERI(q,b,p,j) - 1d0*ERI(q,b,j,p) &
|
||||
- 1d0*eh_Phi(q,b,j,p) + 1d0*pp_Phi(q,b,p,j)) * Y
|
||||
|
||||
end do
|
||||
@ -107,7 +107,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
|
||||
do c=nO+1,nOrb-nR
|
||||
do d=c+1,nOrb-nR
|
||||
cd = cd + 1
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,c,d) - ERI(p,q,d,c) &
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ( 1d0*ERI(p,q,c,d) - 1d0*ERI(p,q,d,c) &
|
||||
+ 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X1(cd,ab)
|
||||
end do
|
||||
end do
|
||||
@ -116,7 +116,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,k,l) - ERI(p,q,l,k) &
|
||||
rho1(p,q,ab) = rho1(p,q,ab) + ( 1d0*ERI(p,q,k,l) - 1d0*ERI(p,q,l,k) &
|
||||
+ 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab)
|
||||
end do
|
||||
end do
|
||||
@ -132,7 +132,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
|
||||
do c=nO+1,nOrb-nR
|
||||
do d=c+1,nOrb-nR
|
||||
cd = cd + 1
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,c,d) - ERI(p,q,d,c) &
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ( 1d0*ERI(p,q,c,d) - 1d0*ERI(p,q,d,c) &
|
||||
+ 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X2(cd,ij)
|
||||
end do
|
||||
end do
|
||||
@ -141,7 +141,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
|
||||
do k=nC+1,nO
|
||||
do l=k+1,nO
|
||||
kl = kl + 1
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,k,l) - ERI(p,q,l,k) &
|
||||
rho2(p,q,ij) = rho2(p,q,ij) + ( 1d0*ERI(p,q,k,l) - 1d0*ERI(p,q,l,k) &
|
||||
+ 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij)
|
||||
end do
|
||||
end do
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user