PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

passing Fock matrix and others in R branch

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-09-02 22:18:53 +02:00
parent e37becbba4
commit fda6a45276
5 changed files with 129 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/CC/RCC.f90
View File

@ -1,8 +1,6 @@
subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
! --- maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc, &
ERHF,eHF,cHF,PHF,FHF)
subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, &
maxSCF, thresh, max_diis, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_AO, ERI_MO, ENuc, ERHF, eHF, cHF)
! Coupled-cluster module ! Coupled-cluster module
@ -27,7 +25,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: nBas, nOrb integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nC integer,intent(in) :: nC
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nV integer,intent(in) :: nV
@ -36,6 +35,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
double precision,intent(in) :: ERHF double precision,intent(in) :: ERHF
double precision,intent(in) :: eHF(nOrb) double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: FHF(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb)
@ -166,9 +167,8 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d
if(dopCCD) then if(dopCCD) then
call wall_time(start_CC) call wall_time(start_CC)
call pCCD(dotest, maxSCF, thresh, max_diis, nBas, nOrb, & call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, &
nC, nO, nV, nR, Hc, ERI_AO, ENuc, ERHF, eHF, cHF) Hc,ERI_AO,ENuc,ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_CC) call wall_time(end_CC)
t_CC = end_CC - start_CC t_CC = end_CC - start_CC

69
src/CC/pCCD.f90
View File

@ -1,8 +1,5 @@
subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nOrb,nC,nO,nV,nR, &
! --- Hc,ERI_AO,ENuc,ERHF,eHF,cHF,PHF,FHF)
subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
nC, nO, nV, nR, Hc, ERI_AO, ENuc, ERHF, eHF, cHF)
! pair CCD module ! pair CCD module
@ -16,15 +13,23 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: nBas, nOrb, nC, nO, nV, nR 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,ERHF double precision,intent(in) :: ENuc,ERHF
double precision,intent(in) :: eHF(nOrb) double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: FHF(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
! Local variables ! Local variables
integer :: mu,nu
integer :: p,q,r,s,t,u,w integer :: p,q,r,s,t,u,w
integer :: pq,rs integer :: pq,rs
integer :: i,j,a,b integer :: i,j,a,b
@ -35,6 +40,7 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
double precision :: CvgOrb double precision :: CvgOrb
double precision :: ECC double precision :: ECC
double precision :: EcCC double precision :: EcCC
double precision :: dECC
double precision,allocatable :: eO(:) double precision,allocatable :: eO(:)
double precision,allocatable :: eV(:) double precision,allocatable :: eV(:)
@ -93,18 +99,18 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
O = nO - nC O = nO - nC
V = nV - nR V = nV - nR
N = O + V ! nOrb - nC - nR N = O + V
!------------------------------------! !------------------------------------!
! Star Loop for orbital optimization ! ! Star Loop for orbital optimization !
!------------------------------------! !------------------------------------!
allocate(ERI_MO(N,N,N,N)) allocate(ERI_MO(N,N,N,N))
allocate(c(nBas,N), h(N,N)) allocate(c(nBas,N),h(N,N))
allocate(eO(O), eV(V), delta_OV(O,V)) allocate(eO(O),eV(V),delta_OV(O,V))
allocate(OOOO(O,O), OOVV(O,V), OVOV(O,V), OVVO(O,V), VVVV(V,V)) allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V))
do i = 1, N do i=1,N
c(:,i) = cHF(:,nC+i) c(:,i) = cHF(:,nC+i)
enddo enddo
@ -116,20 +122,34 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
write(*,*)'| Orbital Optimization for pCCD |' write(*,*)'| Orbital Optimization for pCCD |'
write(*,*)'----------------------------------------------------' write(*,*)'----------------------------------------------------'
do while(CvgOrb > thresh .and. nItOrb < 1) do while(CvgOrb > thresh .and. nItOrb < maxIt)
nItOrb = nItOrb + 1 nItOrb = nItOrb + 1
! Transform integrals ! Transform integrals
h = matmul(transpose(c), matmul(Hc, c)) h = matmul(transpose(c),matmul(Hc,c))
call AOtoMO_ERI_RHF(nBas, N, c(1,1), ERI_AO(1,1,1,1), ERI_MO(1,1,1,1)) call AOtoMO_ERI_RHF(nBas,N,c(1,1),ERI_AO(1,1,1,1),ERI_MO(1,1,1,1))
! Form energy denominator ! Form energy denominator
eO(:) = eHF(nC+1:nO) eO(:) = 0d0
eV(:) = eHF(nO+1:nOrb-nR) eV(:) = 0d0
do mu=1,nBas
do nu=1,nBas
do i=1,O
eO(i) = eO(i) + c(mu,i)*FHF(mu,nu)*c(nu,i)
end do
do a=1,V
eV(a) = eV(a) + c(mu,O+a)*FHF(mu,nu)*c(nu,O+a)
end do
end do
end do
do i=1,O do i=1,O
do a=1,V do a=1,V
@ -170,6 +190,7 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
nItAmp = 0 nItAmp = 0
ECC = ERHF ECC = ERHF
EcCC = 0d0 EcCC = 0d0
dECC = ECC
n_diis = 0 n_diis = 0
t2(:,:) = 0d0 t2(:,:) = 0d0
@ -573,9 +594,13 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
! Check convergence of orbital optimization ! Check convergence of orbital optimization
CvgOrb = maxval(abs(grad)) CvgOrb = maxval(abs(grad))
write(*,*) ' Iteration',nItOrb,'for pCCD orbital optimization' write(*,'(A10,I4,A30)') ' Iteration',nItOrb,'for pCCD orbital optimization'
write(*,*) ' Convergence of orbital gradient = ',CvgOrb write(*,*) '----------------------------------------------------------'
write(*,'(A40,F16.10,A3)') ' Convergence of orbital gradient = ',CvgOrb,' au'
write(*,'(A40,F16.10,A3)') ' Energy difference = ',ECC-dECC,' au'
write(*,*) write(*,*)
dECC = ECC
!-------------------------! !-------------------------!
! Compute orbital Hessian ! ! Compute orbital Hessian !
@ -703,18 +728,18 @@ subroutine pCCD(dotest, maxIt, thresh, max_diis, nBas, nOrb, &
deallocate(Kap) deallocate(Kap)
write(*,*) 'e^kappa' write(*,*) 'e^kappa'
call matout(N, N, ExpKap) call matout(N,N,ExpKap)
write(*,*) write(*,*)
write(*,*) 'Old orbitals' write(*,*) 'Old orbitals'
call matout(nBas, N, c) call matout(nBas,N,c)
write(*,*) write(*,*)
c = matmul(c, ExpKap) c = matmul(c,ExpKap)
deallocate(ExpKap) deallocate(ExpKap)
write(*,*) 'Rotated orbitals' write(*,*) 'Rotated orbitals'
call matout(nBas, N, c) call matout(nBas,N,c)
write(*,*) write(*,*)
end do end do

59
src/HF/RHF.f90
View File

@ -1,8 +1,5 @@
subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
! --- nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,eHF,c,P,F)
subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI, dipole_int, X, ERHF, eHF, c, P)
! Perform restricted Hartree-Fock calculation ! Perform restricted Hartree-Fock calculation
@ -19,7 +16,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift double precision,intent(in) :: level_shift
integer,intent(in) :: nBas, nOrb integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO integer,intent(in) :: nO
integer,intent(in) :: nNuc integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: ZNuc(nNuc)
@ -53,7 +51,6 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
double precision,allocatable :: J(:,:) double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:) double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:) double precision,allocatable :: cp(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: Fp(:,:) double precision,allocatable :: Fp(:,:)
! Output variables ! Output variables
@ -62,6 +59,7 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
double precision,intent(out) :: eHF(nOrb) double precision,intent(out) :: eHF(nOrb)
double precision,intent(inout):: c(nBas,nOrb) double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: P(nBas,nBas) double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: F(nBas,nBas)
! Hello world ! Hello world
@ -81,7 +79,6 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
allocate(K(nBas,nBas)) allocate(K(nBas,nBas))
allocate(err(nBas,nBas)) allocate(err(nBas,nBas))
allocate(F(nBas,nBas))
allocate(cp(nOrb,nOrb)) allocate(cp(nOrb,nOrb))
allocate(Fp(nOrb,nOrb)) allocate(Fp(nOrb,nOrb))
@ -91,7 +88,7 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
! Guess coefficients and density matrix ! Guess coefficients and density matrix
call mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c)
!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, nBas, nO, 2.d0, & call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
@ -105,8 +102,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
err_diis(:,:) = 0d0 err_diis(:,:) = 0d0
rcond = 0d0 rcond = 0d0
Conv = 1d0 Conv = 1d0
nSCF = 0 nSCF = 0
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Main SCF loop ! Main SCF loop
@ -126,8 +123,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
! Build Fock matrix ! Build Fock matrix
call Hartree_matrix_AO_basis(nBas, P, ERI, J) call Hartree_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas, P, ERI, K) call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
if(nBas .ne. nOrb) then if(nBas .ne. nOrb) then
@ -137,24 +134,24 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
! Check convergence ! Check convergence
err = matmul(F, matmul(P, S)) - matmul(matmul(S, P), F) err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
if(nSCF > 1) Conv = maxval(abs(err)) if(nSCF > 1) Conv = maxval(abs(err))
! Kinetic energy ! Kinetic energy
ET = trace_matrix(nBas, matmul(P, T)) ET = trace_matrix(nBas,matmul(P,T))
! Potential energy ! Potential energy
EV = trace_matrix(nBas, matmul(P, V)) EV = trace_matrix(nBas,matmul(P,V))
! Hartree energy ! Hartree energy
EJ = 0.5d0*trace_matrix(nBas, matmul(P, J)) EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
! Exchange energy ! Exchange energy
EK = 0.25d0*trace_matrix(nBas, matmul(P, K)) EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
! Total energy ! Total energy
@ -164,29 +161,29 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
if(max_diis > 1) then if(max_diis > 1) then
n_diis = min(n_diis+1, max_diis) n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond, nBas_Sq, nBas_Sq, n_diis, err_diis, F_diis, err, F) call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
end if end if
! Level shift ! Level shift
if(level_shift > 0d0 .and. Conv > thresh) then if(level_shift > 0d0 .and. Conv > thresh) then
call level_shifting(level_shift, nBas, nOrb, nO, S, c, F) call level_shifting(level_shift,nBas,nOrb,nO,S,c,F)
endif endif
! Diagonalize Fock matrix ! Diagonalize Fock matrix
if(nBas .eq. nOrb) then if(nBas .eq. nOrb) then
Fp = matmul(transpose(X), matmul(F, X)) Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eHF) call diagonalize_matrix(nOrb,cp,eHF)
c = matmul(X, cp) c = matmul(X,cp)
else else
Fp = matmul(transpose(c), matmul(F, c)) Fp = matmul(transpose(c),matmul(F,c))
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eHF) call diagonalize_matrix(nOrb,cp,eHF)
c = matmul(c, cp) c = matmul(c,cp)
endif endif
! Density matrix ! Density matrix
@ -217,7 +214,7 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*) write(*,*)
deallocate(J, K, err, F, cp, Fp, err_diis, F_diis) deallocate(J,K,err,cp,Fp,err_diis,F_diis)
stop stop
@ -225,8 +222,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
! Compute dipole moments ! Compute dipole moments
call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int, dipole) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_RHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, EK, ERHF, dipole) call print_RHF(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,ERHF,dipole)
! Testing zone ! Testing zone
@ -239,6 +236,6 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc,
end if end if
deallocate(J, K, err, F, cp, Fp, err_diis, F_diis) deallocate(J,K,err,cp,Fp,err_diis,F_diis)
end subroutine end subroutine

65
src/HF/ROHF.f90
View File

@ -1,8 +1,5 @@
subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
! --- nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EROHF,eHF,c,Ptot,Ftot)
subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI, dipole_int, X, EROHF, eHF, c, Ptot)
! Perform restricted open-shell Hartree-Fock calculation ! Perform restricted open-shell Hartree-Fock calculation
@ -19,7 +16,8 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
double precision,intent(in) :: mix double precision,intent(in) :: mix
double precision,intent(in) :: level_shift double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: nBas, nOrb integer,intent(in) :: nOrb
integer,intent(in) :: nBas
integer,intent(in) :: nNuc integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: ZNuc(nNuc)
@ -52,7 +50,6 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
double precision,allocatable :: J(:,:,:) double precision,allocatable :: J(:,:,:)
double precision,allocatable :: F(:,:,:) double precision,allocatable :: F(:,:,:)
double precision,allocatable :: Fp(:,:) double precision,allocatable :: Fp(:,:)
double precision,allocatable :: Ftot(:,:)
double precision,allocatable :: P(:,:,:) double precision,allocatable :: P(:,:,:)
double precision,allocatable :: K(:,:,:) double precision,allocatable :: K(:,:,:)
double precision,allocatable :: err(:,:) double precision,allocatable :: err(:,:)
@ -68,6 +65,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
double precision,intent(out) :: eHF(nOrb) double precision,intent(out) :: eHF(nOrb)
double precision,intent(inout):: c(nBas,nOrb) double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: Ptot(nBas,nBas) double precision,intent(out) :: Ptot(nBas,nBas)
double precision,intent(out) :: Ftot(nBas,nBas)
! Hello world ! Hello world
@ -86,7 +84,6 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
allocate(J(nBas,nBas,nspin)) allocate(J(nBas,nBas,nspin))
allocate(K(nBas,nBas,nspin)) allocate(K(nBas,nBas,nspin))
allocate(F(nBas,nBas,nspin)) allocate(F(nBas,nBas,nspin))
allocate(Ftot(nBas,nBas))
allocate(P(nBas,nBas,nspin)) allocate(P(nBas,nBas,nspin))
allocate(err(nBas,nBas)) allocate(err(nBas,nBas))
@ -98,9 +95,9 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
! Guess coefficients and demsity matrices ! Guess coefficients and demsity matrices
call mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) call mo_guess(nBas,nOrb,guess_type,S,Hc,X,c)
do ispin = 1, nspin do ispin = 1,nspin
!P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin)))) !P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin))))
call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas) call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas)
end do end do
@ -134,51 +131,51 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
! Build Hartree repulsion ! Build Hartree repulsion
do ispin = 1, nspin do ispin=1,nspin
call Hartree_matrix_AO_basis(nBas, P(:,:,ispin), ERI(:,:,:,:), J(:,:,ispin)) call Hartree_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin))
end do end do
! Compute exchange potential ! Compute exchange potential
do ispin = 1, nspin do ispin=1,nspin
call exchange_matrix_AO_basis(nBas, P(:,:,ispin), ERI(:,:,:,:), K(:,:,ispin)) call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),K(:,:,ispin))
end do end do
! Build Fock operator ! Build Fock operator
do ispin = 1, nspin do ispin=1,nspin
F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin) F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin)
end do end do
call ROHF_fock_matrix(nBas, nOrb, nO(1), nO(2), S, c, F(:,:,1), F(:,:,2), Ftot) call ROHF_fock_matrix(nBas,nOrb,nO(1),nO(2),S,c,F(:,:,1),F(:,:,2),Ftot)
! Check convergence ! Check convergence
err(:,:) = matmul(Ftot, matmul(Ptot, S)) - matmul(matmul(S, Ptot), Ftot) err(:,:) = matmul(Ftot,matmul(Ptot,S)) - matmul(matmul(S,Ptot),Ftot)
if(nSCF > 1) Conv = maxval(abs(err(:,:))) if(nSCF > 1) Conv = maxval(abs(err(:,:)))
! Kinetic energy ! Kinetic energy
do ispin = 1, nspin do ispin=1,nspin
ET(ispin) = trace_matrix(nBas, matmul(P(:,:,ispin), T(:,:))) ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:)))
end do end do
! Potential energy ! Potential energy
do ispin = 1, nspin do ispin=1,nspin
EV(ispin) = trace_matrix(nBas, matmul(P(:,:,ispin), V(:,:))) EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
end do end do
! Hartree energy ! Hartree energy
EJ(1) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,1), J(:,:,1))) EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas, matmul(P(:,:,1), J(:,:,2))) EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,2), J(:,:,2))) EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
! Exchange energy ! Exchange energy
do ispin = 1, nspin do ispin = 1,nspin
EK(ispin) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,ispin), K(:,:,ispin))) EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
end do end do
! Total energy ! Total energy
@ -199,27 +196,27 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
if(level_shift > 0d0 .and. Conv > thresh) then if(level_shift > 0d0 .and. Conv > thresh) then
do ispin=1,nspin do ispin=1,nspin
call level_shifting(level_shift, nBas, nOrb, maxval(nO), S, c, Ftot) call level_shifting(level_shift,nBas,nOrb,maxval(nO),S,c,Ftot)
end do end do
end if end if
! Transform Fock matrix in orthogonal basis ! Transform Fock matrix in orthogonal basis
Fp(:,:) = matmul(transpose(X(:,:)), matmul(Ftot(:,:), X(:,:))) Fp(:,:) = matmul(transpose(X(:,:)),matmul(Ftot(:,:),X(:,:)))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues ! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb, cp, eHF) call diagonalize_matrix(nOrb,cp,eHF)
! Back-transform eigenvectors in non-orthogonal basis ! Back-transform eigenvectors in non-orthogonal basis
c(:,:) = matmul(X(:,:), cp(:,:)) c(:,:) = matmul(X(:,:),cp(:,:))
! Compute density matrix ! Compute density matrix
do ispin = 1, nspin do ispin=1,nspin
!P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin)))) !P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin))))
call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas) call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas)
end do end do
@ -246,7 +243,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*) write(*,*)
deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis) deallocate(J,K,F,P,err,Fp,cp,err_diis,F_diis)
stop stop
@ -255,7 +252,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
! Compute final UHF energy ! Compute final UHF energy
call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_ROHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, EK, EROHF, dipole) call print_ROHF(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,EROHF,dipole)
! Print test values ! Print test values
@ -265,6 +262,6 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift,
end if end if
deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis) deallocate(J,K,F,P,err,Fp,cp,err_diis,F_diis)
end subroutine end subroutine

23
src/QuAcK/RQuAcK.f90
View File

@ -92,7 +92,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
double precision :: start_GW ,end_GW ,t_GW double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT double precision :: start_GT ,end_GT ,t_GT
double precision,allocatable :: cHF(:,:),eHF(:),PHF(:,:) double precision,allocatable :: eHF(:)
double precision,allocatable :: cHF(:,:)
double precision,allocatable :: PHF(:,:)
double precision,allocatable :: FHF(:,:)
double precision :: ERHF double precision :: ERHF
double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:)
@ -109,9 +112,10 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
! Memory allocation ! ! Memory allocation !
!-------------------! !-------------------!
allocate(cHF(nBas,nOrb))
allocate(eHF(nOrb)) allocate(eHF(nOrb))
allocate(cHF(nBas,nOrb))
allocate(PHF(nBas,nBas)) allocate(PHF(nBas,nBas))
allocate(FHF(nBas,nBas))
allocate(dipole_int_MO(nOrb,nOrb,ncart)) allocate(dipole_int_MO(nOrb,nOrb,ncart))
allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
@ -122,8 +126,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doRHF) then if(doRHF) then
call wall_time(start_HF) call wall_time(start_HF)
call RHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF) nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_HF) call wall_time(end_HF)
t_HF = end_HF - start_HF t_HF = end_HF - start_HF
@ -135,8 +139,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doROHF) then if(doROHF) then
call wall_time(start_HF) call wall_time(start_HF)
call ROHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, & call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas, nOrb, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF) nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_HF) call wall_time(end_HF)
t_HF = end_HF - start_HF t_HF = end_HF - start_HF
@ -163,7 +167,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
! 4-index transform ! 4-index transform
call AOtoMO_ERI_RHF(nBas, nOrb, cHF, ERI_AO, ERI_MO) call AOtoMO_ERI_RHF(nBas,nOrb,cHF,ERI_AO,ERI_MO)
call wall_time(end_AOtoMO) call wall_time(end_AOtoMO)
@ -231,8 +235,9 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
if(doCC) then if(doCC) then
call wall_time(start_CC) call wall_time(start_CC)
call RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, & call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
maxSCF_CC, thresh_CC, max_diis_CC, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_AO, ERI_MO, ENuc, ERHF, eHF, cHF) maxSCF_CC,thresh_CC,max_diis_CC,nBas,nOrb,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO, &
ENuc,ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_CC) call wall_time(end_CC)
t_CC = end_CC - start_CC t_CC = end_CC - start_CC