10
1
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00

introduce nBas_MOs in RHF

This commit is contained in:
Abdallah Ammar 2024-08-28 17:27:46 +02:00
parent 567e77dd4a
commit de867f2b54
15 changed files with 251 additions and 143 deletions

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

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)
nBas_AOs, nBas_MOs, nO, S, T, V, Hc, ERI, dipole_int, X, ERHF, eHF, c, P)
! Perform restricted Hartree-Fock calculation
@ -16,24 +19,24 @@ 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_AOs, nBas_MOs
integer,intent(in) :: nO
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
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) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: T(nBas_AOs,nBas_AOs)
double precision,intent(in) :: V(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
double precision,intent(in) :: ERI(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
double precision,intent(in) :: dipole_int(nBas_AOs,nBas_AOs,ncart)
! Local variables
integer :: nSCF
integer :: nBasSq
integer :: nBas_AOs_Sq
integer :: n_diis
double precision :: ET
double precision :: EV
@ -56,9 +59,9 @@ 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) :: P(nBas,nBas)
double precision,intent(out) :: eHF(nBas_MOs)
double precision,intent(inout):: c(nBas_AOs,nBas_MOs)
double precision,intent(out) :: P(nBas_AOs,nBas_AOs)
! Hello world
@ -70,17 +73,28 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Useful quantities
nBasSq = nBas*nBas
nBas_AOs_Sq = nBas_AOs*nBas_AOs
! 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_AOs,nBas_AOs))
allocate(K(nBas_AOs,nBas_AOs))
allocate(err(nBas_AOs,nBas_AOs))
allocate(F(nBas_AOs,nBas_AOs))
allocate(cp(nBas_MOs,nBas_MOs))
allocate(Fp(nBas_MOs,nBas_MOs))
allocate(err_diis(nBas_AOs_Sq,max_diis))
allocate(F_diis(nBas_AOs_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_AOs, nBas_MOs, guess_type, S, Hc, X, c)
!P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO)))
call dgemm('N', 'T', nBas_AOs, nBas_AOs, nO, 2.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P, nBas_AOs)
! Initialization
@ -110,8 +124,8 @@ 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_AOs, P, ERI, J)
call exchange_matrix_AO_basis(nBas_AOs, P, ERI, K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
@ -122,19 +136,19 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
ET = trace_matrix(nBas_AOs, matmul(P, T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
EV = trace_matrix(nBas_AOs, matmul(P, V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EJ = 0.5d0*trace_matrix(nBas_AOs, matmul(P, J))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
EK = 0.25d0*trace_matrix(nBas_AOs, matmul(P, K))
! Total energy
@ -145,24 +159,27 @@ 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)
call DIIS_extrapolation(rcond, nBas_AOs_Sq, nBas_AOs_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_AOs, nBas_MOs, nO, S, c, F)
endif
! Diagonalize Fock matrix
Fp = matmul(transpose(X), matmul(F, X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eHF)
call diagonalize_matrix(nBas_MOs, cp, eHF)
c = matmul(X, cp)
! 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_AOs, nBas_AOs, nO, 2.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P, nBas_AOs)
! Dump results
@ -185,14 +202,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_AOs, P, nNuc, ZNuc, rNuc, dipole_int, dipole)
call print_RHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, EK, ERHF, dipole)
! Testing zone
@ -205,4 +224,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

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

View File

@ -86,7 +86,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Guess coefficients and demsity matrices
call mo_guess(nBas,guess_type,S,Hc,X,c)
call mo_guess(nBas,nBas,guess_type,S,Hc,X,c)
do ispin=1,nspin
P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin))))
end do
@ -185,7 +185,7 @@ 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,nBas,maxval(nO),S,c,Ftot)
end do
end if

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

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

View File

@ -1,4 +1,4 @@
subroutine core_guess(nBas,Hc,X,c)
subroutine core_guess(nBas_AOs, nBas_MOs, 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
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
integer,intent(in) :: nBas_AOs, nBas_MOs
double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
! 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_AOs,nBas_MOs)
! Memory allocation
allocate(cp(nBas,nBas),e(nBas))
allocate(cp(nBas_MOs,nBas_MOs), e(nBas_MOs))
! Core guess
cp(:,:) = matmul(transpose(X(:,:)), matmul(Hc(:,:), X(:,:)))
call diagonalize_matrix(nBas,cp,e)
call diagonalize_matrix(nBas_MOs, cp, e)
c(:,:) = matmul(X(:,:), cp(:,:))
deallocate(cp, e)
end subroutine

View File

@ -1,4 +1,4 @@
subroutine huckel_guess(nBas,S,Hc,X,c)
subroutine huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c)
! Hickel guess
@ -6,10 +6,10 @@ subroutine huckel_guess(nBas,S,Hc,X,c)
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
integer,intent(in) :: nBas_AOs, nBas_MOs
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
! Local variables
@ -20,11 +20,11 @@ subroutine huckel_guess(nBas,S,Hc,X,c)
! Output variables
double precision,intent(out) :: c(nBas,nBas)
double precision,intent(out) :: c(nBas_AOs,nBas_MOs)
! Memory allocation
allocate(F(nBas,nBas))
allocate(F(nBas_AOs,nBas_AOs))
! Extended Huckel parameter
@ -32,9 +32,9 @@ subroutine huckel_guess(nBas,S,Hc,X,c)
! GWH approximation
do mu=1,nBas
do mu = 1, nBas_AOs
F(mu,mu) = Hc(mu,mu)
do nu=mu+1,nBas
do nu = mu+1, nBas_AOs
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_AOs, nBas_MOs, F, X, c)
deallocate(F)
end subroutine

View File

@ -1,4 +1,7 @@
subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
! ---
subroutine mo_guess(nBas_AOs, nBas_MOs, 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_AOs, nBas_MOs
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) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
! Output variables
double precision,intent(inout) :: c(nBas,nBas)
double precision,intent(inout) :: c(nBas_AOs,nBas_MOs)
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_AOs, nBas_MOs, Hc, X, c)
elseif(guess_type == 2) then
write(*,*) 'Huckel guess...'
call huckel_guess(nBas,S,Hc,X,c)
call huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c)
elseif(guess_type == 3) then

View File

@ -1,4 +1,7 @@
subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! ---
subroutine print_RHF(nBas_AOs, nBas_MOs, 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_AOs, nBas_MOs
integer,intent(in) :: nO
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: eHF(nBas_MOs)
double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs)
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_AOs, nBas_MOs, cHF)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' RHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eHF)
call vecout(nBas_MOs, eHF)
write(*,*)
end subroutine

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_AOs
integer :: nNuc, nBas_AOs, nBas_MOs
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,8 @@ program QuAcK
logical :: dotest,doRtest,doUtest,doGtest
integer :: i, j
!-------------!
! Hello World !
!-------------!
@ -128,6 +131,7 @@ program QuAcK
! nV = number of virtual orbitals (see below) !
! nR = number of Rydberg orbitals !
! nBas_AOs = number of basis functions in AOs !
! nBas_MOs = number of basis functions in MOs !
!---------------------------------------------------!
call read_molecule(nNuc,nO,nC,nR)
@ -153,7 +157,6 @@ program QuAcK
allocate(T(nBas_AOs,nBas_AOs))
allocate(V(nBas_AOs,nBas_AOs))
allocate(Hc(nBas_AOs,nBas_AOs))
allocate(X(nBas_AOs,nBas_AOs))
allocate(ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs))
allocate(dipole_int_AO(nBas_AOs,nBas_AOs,ncart))
@ -173,7 +176,39 @@ program QuAcK
! Compute orthogonalization matrix
call orthogonalization_matrix(nBas_AOs, S, X)
!call orthogonalization_matrix(nBas_AOs, S, X)
allocate(Uvec(nBas_AOs,nBas_AOs), Uval(nBas_AOs))
Uvec(1:nBas_AOs,1:nBas_AOs) = S(1:nBas_AOs,1:nBas_AOs)
call diagonalize_matrix(nBas_AOs, Uvec, Uval)
nBas_MOs = 0
do i = 1, nBas_AOs
if(Uval(i) > 1d-6) then
Uval(i) = 1d0 / dsqrt(Uval(i))
nBas_MOs = nBas_MOs + 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_AOs
write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nBas_MOs
write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nBas_MOs)/dble(nBas_AOs))
write(*,'(A38)') '--------------------------------------'
write(*,*)
allocate(X(nBas_AOs,nBas_MOs))
do j = 1, nBas_MOs
do i = 1, nBas_AOs
X(i,j) = Uvec(i,j) * Uval(j)
enddo
enddo
deallocate(Uvec, Uval)
!---------------------!
! Choose QuAcK branch !
@ -205,7 +240,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_AOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
nNuc,nBas_AOs,nBas_MOs,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, &

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_AOs,nBas_MOs,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_AOs,nBas_MOs
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -38,13 +38,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart)
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) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: T(nBas_AOs,nBas_AOs)
double precision,intent(in) :: V(nBas_AOs,nBas_AOs)
double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs)
double precision,intent(in) :: X(nBas_AOs,nBas_MOs)
double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart)
double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)
integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix
@ -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_AOs,nBas_MOs))
allocate(eHF(nBas_MOs))
allocate(PHF(nBas_AOs,nBas_AOs))
allocate(dipole_int_MO(nBas_AOs,nBas_AOs,ncart))
allocate(ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs))
!---------------------!
! Hartree-Fock module !
@ -119,12 +122,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doRHF) then
call wall_time(start_HF)
! TODO
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)
nBas_AOs, nBas_MOs, 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 +136,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doROHF) then
call wall_time(start_HF)
! TODO
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)
nBas_AOs,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
@ -155,17 +160,19 @@ 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))
! TODO
call AOtoMO(nBas_AOs,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
! 4-index transform
call AOtoMO_ERI_RHF(nBas,cHF,ERI_AO,ERI_MO)
! TODO
call AOtoMO_ERI_RHF(nBas_AOs,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 +184,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(nBas_AOs,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 +196,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(dosearch) then
call wall_time(start_stab)
! TODO
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)
nBas_AOs,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 +216,12 @@ 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)
! TODO
call RMP(dotest,doMP2,doMP3,reg_MP,nBas_AOs,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 +236,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doCC) then
call wall_time(start_CC)
! TODO
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_MO,ENuc,ERHF,eHF,cHF)
maxSCF_CC,thresh_CC,max_diis_CC,nBas_AOs,nC,nO,nV,nR,Hc,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 +256,13 @@ 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, &
! TODO
call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas_AOs,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, &
eHF,ERHF,cHF,S)
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 +276,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doRPA) then
call wall_time(start_RPA)
! TODO
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)
nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S)
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 +296,15 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGF) then
call wall_time(start_GF)
! TODO
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, &
nNuc,ZNuc,rNuc,ENuc,nBas_AOs,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 +318,15 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGW) then
call wall_time(start_GW)
! TODO
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, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas_AOs,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 +340,15 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doGT) then
call wall_time(start_GT)
! TODO
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)
nBas_AOs,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

View File

@ -1,4 +1,4 @@
subroutine level_shifting(level_shift,nBas,nO,S,c,F)
subroutine level_shifting(level_shift, nBas_AOs, nBas_MOs, 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_AOs, nBas_MOs
integer,intent(in) :: nO
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: c(nBas_AOs,nBas_MOs)
! Local variables
@ -21,17 +21,19 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F)
! Output variables
double precision,intent(inout):: F(nBas,nBas)
double precision,intent(inout):: F(nBas_AOs,nBas_AOs)
allocate(F_MO(nBas,nBas),Sc(nBas,nBas))
allocate(F_MO(nBas_MOs,nBas_MOs), Sc(nBas_AOs,nBas_MOs))
F_MO(:,:) = matmul(transpose(c), matmul(F, c))
do a=nO+1,nBas
do a = nO+1, nBas_MOs
F_MO(a,a) = F_MO(a,a) + level_shift
end do
Sc(:,:) = matmul(S, c)
F(:,:) = matmul(Sc, matmul(F_MO, transpose(Sc)))
deallocate(F_MO, Sc)
end subroutine

View File

@ -1,3 +1,6 @@
! ---
subroutine orthogonalization_matrix(nBas, S, X)
! Compute the orthogonalization matrix X
@ -35,6 +38,24 @@ 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(*,*)
@ -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
@ -69,7 +90,7 @@ subroutine orthogonalization_matrix(nBas,S,X)
if(Uval(i) > thresh) then
Uval(i) = 1d0/sqrt(Uval(i))
Uval(i) = 1d0 / dsqrt(Uval(i))
else
@ -117,3 +138,6 @@ subroutine orthogonalization_matrix(nBas,S,X)
end if
end subroutine
! ---

View File

@ -24,10 +24,10 @@ subroutine read_basis_pyscf(nBas_AOs, nO, nV)
read(3, *) nBas_AOs
close(unit=3)
write(*,'(A38)') '--------------------------------------'
write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs
write(*,'(A38)') '--------------------------------------'
write(*,*)
! write(*,'(A38)') '--------------------------------------'
! write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs
! write(*,'(A38)') '--------------------------------------'
! write(*,*)
! Number of virtual orbitals