PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-01 06:21:37 +02:00
Code Releases Activity

diis on Phis

Browse Source
This commit is contained in:
Pierre-Francois Loos 2025-03-27 11:14:13 +01:00
parent b3178e8a08
commit ea26528afc
4 changed files with 70 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
mol/H2O.xyz
View File

@ -1,5 +1,5 @@
3
O 0.0000 0.0000 0.0000
H 0.9591 0.0000 0.0000
H -0.2373 0.9293 0.0000
H 0.7571 0.0000 0.5861
H -0.7571 0.0000 0.5861

97
src/Parquet/GParquet.f90
View File

@ -59,11 +59,13 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
integer :: max_diis,n_diis
double precision :: rcond
double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: Om_diis(:,:)
double precision,allocatable :: Phi_diis(:,:)
double precision,allocatable :: err(:)
double precision,allocatable :: Om(:)
double precision,allocatable :: Phi(:)
double precision :: alpha
integer ::p,q,r,s,pqrs
! Output variables
! None
@ -75,15 +77,15 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
! DIIS parameters
! max_diis = 10
! n_diis = 0
! rcond = 0d0
max_diis = 1
n_diis = 0
rcond = 1d0
! allocate(err_diis(nS+nOO+nVV,max_diis),Om_diis(nS+nOO+nVV,max_diis))
! allocate(err(nS+nOO+nVV),Om(nS+nOO+nVV))
allocate(err_diis(2*nOrb**4,max_diis),Phi_diis(2*nOrb**4,max_diis))
allocate(err(2*nOrb**4),Phi(2*nOrb**4))
! err_diis(:,:) = 0d0
! Om_diis(:,:) = 0d0
err_diis(:,:) = 0d0
Phi_diis(:,:) = 0d0
! Start
@ -365,38 +367,64 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:)
! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:)
err_eh = maxval(abs(old_eh_Phi - eh_Phi))
err_pp = maxval(abs(old_pp_Phi - 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)
!--------------------!
! DIIS extrapolation !
!--------------------!
pqrs = 0
do p=1,nOrb
do q=1,nOrb
do r=1,nOrb
do s=1,nOrb
pqrs = pqrs + 1
err( pqrs) = eh_Phi(p,q,r,s) - old_eh_Phi(p,q,r,s)
err(nOrb**4+pqrs) = pp_Phi(p,q,r,s) - old_pp_Phi(p,q,r,s)
Phi( pqrs) = eh_Phi(p,q,r,s)
Phi(nOrb**4+pqrs) = pp_Phi(p,q,r,s)
end do
end do
end do
end do
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,2*nOrb**4,2*nOrb**4,n_diis,err_diis,Phi_diis,err,Phi)
print*,rcond
end if
pqrs = 0
do p=1,nOrb
do q=1,nOrb
do r=1,nOrb
do s=1,nOrb
pqrs = pqrs + 1
eh_Phi(p,q,r,s) = Phi( pqrs)
pp_Phi(p,q,r,s) = Phi(nOrb**4+pqrs)
end do
end do
end do
end do
old_eh_Phi(:,:,:,:) = eh_Phi(:,:,:,:)
old_pp_Phi(:,:,:,:) = pp_Phi(:,:,:,:)
! Free memory
deallocate(eh_Phi,pp_Phi)
!--------------------!
! DIIS extrapolation !
!--------------------!
! err( 1:nS ) = eh_Om(:) - old_eh_Om(:)
! err( nS+1:nS+nVV ) = ee_Om(:) - old_ee_Om(:)
! err(nVV+nS+1:nS+nVV+nOO) = hh_Om(:) - old_hh_Om(:)
! Om( 1:nS ) = eh_Om(:)
! Om( nS+1:nS+nVV ) = ee_Om(:)
! Om(nVV+nS+1:nS+nVV+nOO) = hh_Om(:)
! if(max_diis > 1) then
! n_diis = min(n_diis+1,max_diis)
! call DIIS_extrapolation(rcond,nS+nOO+nVV,nS+nOO+nVV,n_diis,err_diis,Om_diis,err,Om)
! end if
! eh_Om(:) = Om( 1:nS )
! ee_Om(:) = Om( nS+1:nS+nVV )
! hh_Om(:) = Om(nVV+nS+1:nS+nVV+nOO)
write(*,*) '----------------------------------------'
write(*,*) ' Two-body (kernel) convergence '
write(*,*) '----------------------------------------'
@ -404,7 +432,6 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF,
write(*,'(1X,A30,F10.6)')'Error for pp channel = ',err_pp
write(*,*) '----------------------------------------'
write(*,*)
! Convergence criteria
err_2b = max(err_eh,err_pp)

12
src/Parquet/G_screened_integrals.f90
View File

@ -38,9 +38,9 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho
rho(p,q,ia) = rho(p,q,ia) &
+ (ERI(q,j,p,b) - ERI(q,j,b,p) &
- eh_Phi(q,j,b,p) + pp_Phi(q,j,p,b)) * X &
- 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) &
- eh_Phi(q,b,j,p) + pp_Phi(q,b,p,j)) * Y
- 1d0*eh_Phi(q,b,j,p) + 1d0*pp_Phi(q,b,p,j)) * Y
end do
end do
@ -108,7 +108,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
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) &
+ eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X1(cd,ab)
+ 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X1(cd,ab)
end do
end do
@ -117,7 +117,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
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) &
+ eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y1(kl,ab)
+ 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab)
end do
end do
@ -133,7 +133,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
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) &
+ eh_Phi(p,q,c,d) - eh_Phi(p,q,d,c) ) * X2(cd,ij)
+ 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X2(cd,ij)
end do
end do
@ -142,7 +142,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2
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) &
+ eh_Phi(p,q,k,l) - eh_Phi(p,q,l,k) ) * Y2(kl,ij)
+ 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij)
end do
end do
end do

2
src/QuAcK/QuAcK.f90
View File

@ -198,8 +198,6 @@ program QuAcK
call read_dipole_integrals(working_dir,nBas,dipole_int_AO)
call wall_time(end_int)
call matout(nBas,nBas,dipole_int_AO(:,:,1))
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 1e-integrals = ',t_int,' seconds'