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

introduce nBas_MOs in ROHF

Browse Source
This commit is contained in:
Abdallah Ammar 2024-08-28 18:39:51 +02:00
parent de867f2b54
commit d823cdcd20
18 changed files with 171 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
src/AOtoMO/AOtoMO.f90
View File

@ -1,14 +1,15 @@
subroutine AOtoMO(nBas,C,A,B)
subroutine AOtoMO(nBas_AOs, nBas_MOs, C, M_AOs, M_MOs)
! Perform AO to MO transformation of a matrix A for given coefficients c
! Perform AO to MO transformation of a matrix M_AOs for given coefficients c
! M_MOs = C.T M_AOs C
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: C(nBas,nBas)
double precision,intent(in) :: A(nBas,nBas)
integer,intent(in) :: nBas_AOs, nBas_MOs
double precision,intent(in) :: C(nBas_AOs,nBas_MOs)
double precision,intent(in) :: M_AOs(nBas_AOs,nBas_AOs)
! Local variables
@ -16,11 +17,21 @@ subroutine AOtoMO(nBas,C,A,B)
! Output variables
double precision,intent(out) :: B(nBas,nBas)
double precision,intent(out) :: M_MOs(nBas_MOs,nBas_MOs)
allocate(AC(nBas,nBas))
allocate(AC(nBas_AOs,nBas_MOs))
AC = matmul(A,C)
B = matmul(transpose(C),AC)
!AC = matmul(M_AOs, C)
!M_MOs = matmul(transpose(C), AC)
call dgemm("N", "N", nBas_AOs, nBas_MOs, nBas_AOs, 1.d0, &
M_AOs(1,1), nBas_AOs, C(1,1), nBas_AOs, &
0.d0, AC(1,1), nBas_AOs)
call dgemm("T", "N", nBas_MOs, nBas_MOs, nBas_AOs, 1.d0, &
C(1,1), nBas_AOs, AC(1,1), nBas_AOs, &
0.d0, M_MOs(1,1), nBas_MOs)
deallocate(AC)
end subroutine

39
src/AOtoMO/MOtoAO.f90
View File

@ -1,16 +1,19 @@
subroutine MOtoAO(nBas,S,C,B,A)
subroutine MOtoAO(nBas_AOs, nBas_MOs, S, C, M_MOs, M_AOs)
! Perform MO to AO transformation of a matrix A for a given metric S
! Perform MO to AO transformation of a matrix M_AOs for a given metric S
! and coefficients c
!
! M_AOs = S C M_MOs (S C).T
!
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: C(nBas,nBas)
double precision,intent(in) :: B(nBas,nBas)
integer,intent(in) :: nBas_AOs, nBas_MOs
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: C(nBas_AOs,nBas_MOs)
double precision,intent(in) :: M_MOs(nBas_MOs,nBas_MOs)
! Local variables
@ -18,14 +21,28 @@ subroutine MOtoAO(nBas,S,C,B,A)
! Output variables
double precision,intent(out) :: A(nBas,nBas)
double precision,intent(out) :: M_AOs(nBas_AOs,nBas_AOs)
! Memory allocation
allocate(SC(nBas,nBas),BSC(nBas,nBas))
allocate(SC(nBas_AOs,nBas_MOs), BSC(nBas_MOs,nBas_AOs))
SC = matmul(S,C)
BSC = matmul(B,transpose(SC))
A = matmul(SC,BSC)
!SC = matmul(S, C)
!BSC = matmul(M_MOs, transpose(SC))
!M_AOs = matmul(SC, BSC)
call dgemm("N", "N", nBas_AOs, nBas_MOs, nBas_AOs, 0.d0, &
S(1,1), nBas_AOs, C(1,1), nBas_AOs, &
1.d0, SC(1,1), nBas_AOs)
call dgemm("N", "T", nBas_MOs, nBas_AOs, nBas_MOs, 0.d0, &
M_MOs(1,1), nBas_MOs, SC(1,1), nBas_AOs, &
1.d0, BSC(1,1), nBas_MOs)
call dgemm("N", "N", nBas_AOs, nBas_AOs, nBas_MOs, 0.d0, &
SC(1,1), nBas_AOs, BSC(1,1), nBas_MOs, &
1.d0, M_AOs(1,1), nBas_AOs)
deallocate(SC, BSC)
end subroutine

2
src/GF/qsRGF2.f90
View File

@ -161,7 +161,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si
SigC = 0.5d0*(SigC + transpose(SigC))
call MOtoAO(nBas,S,c,SigC,SigCp)
call MOtoAO(nBas,nBas,S,c,SigC,SigCp)
! Solve the quasi-particle equation

2
src/GF/qsUGF2.f90
View File

@ -197,7 +197,7 @@ subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved
end do
do is=1,nspin
call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
end do
! Solve the quasi-particle equation

2
src/GT/qsRGTeh.f90
View File

@ -192,7 +192,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO(nBas,S,c,Sig,Sigp)
call MOtoAO(nBas,nBas,S,c,Sig,Sigp)
! Solve the quasi-particle equation

2
src/GT/qsRGTpp.f90
View File

@ -235,7 +235,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d
Sig = 0.5d0*(Sig + transpose(Sig))
call MOtoAO(nBas,S,c,Sig,Sigp)
call MOtoAO(nBas,nBas,S,c,Sig,Sigp)
! Solve the quasi-particle equation

2
src/GT/qsUGTpp.f90
View File

@ -280,7 +280,7 @@ subroutine qsUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B
end do
do ispin=1,nspin
call MOtoAO(nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin))
call MOtoAO(nBas,nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin))
end do
! Solve the quasi-particle equation

6
src/GW/SRG_qsGW.f90
View File

@ -175,7 +175,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
call wall_time(tao1)
do ixyz=1,ncart
call AOtoMO(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
call AOtoMO(nBas,nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
@ -218,7 +218,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
! Make correlation self-energy Hermitian and transform it back to AO basis
call wall_time(tmo1)
call MOtoAO(nBas,S,c,SigC,SigCp)
call MOtoAO(nBas,nBas,S,c,SigC,SigCp)
call wall_time(tmo2)
tmo = tmo + tmo2 - tmo1
! Solve the quasi-particle equation
@ -245,7 +245,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
call AOtoMO(nBas,c,SigCp,SigC)
call AOtoMO(nBas,nBas,c,SigCp,SigC)
! Compute new density matrix in the AO basis

8
src/GW/SRG_qsUGW.f90
View File

@ -184,8 +184,8 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS
!--------------------------------------------------
do ixyz=1,ncart
call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block
@ -228,7 +228,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS
! Make correlation self-energy Hermitian and transform it back to AO basis
do is=1,nspin
call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
end do
! Solve the quasi-particle equation
@ -279,7 +279,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS
! Back-transform self-energy
do is=1,nspin
call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
end do
! Compute density matrix

6
src/GW/qsRGW.f90
View File

@ -166,7 +166,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
! AO to MO transformation of two-electron integrals
do ixyz=1,ncart
call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
call AOtoMO(nBas,nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
@ -191,7 +191,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
SigC = 0.5d0*(SigC + transpose(SigC))
call MOtoAO(nBas,S,c,SigC,SigCp)
call MOtoAO(nBas,nBas,S,c,SigC,SigCp)
! Solve the quasi-particle equation
@ -238,7 +238,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eGW)
c = matmul(X,cp)
call AOtoMO(nBas,c,SigCp,SigC)
call AOtoMO(nBas,nBas,c,SigCp,SigC)
! Density matrix

8
src/GW/qsUGW.f90
View File

@ -184,8 +184,8 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
!--------------------------------------------------
do ixyz=1,ncart
call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block
@ -232,7 +232,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
end do
do is=1,nspin
call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is))
end do
! Solve the quasi-particle equation
@ -283,7 +283,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE
! Back-transform self-energy
do is=1,nspin
call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is))
end do
! Compute density matrix

2
src/HF/RHF_search.f90
View File

@ -109,7 +109,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
do ixyz=1,ncart
call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
call AOtoMO(nBas,nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO)
call wall_time(end_AOtoMO)

111
src/HF/ROHF.f90
View File

@ -1,5 +1,8 @@
subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EROHF,eHF,c,Ptot)
! ---
subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, &
nBas_AOs, nBas_MOs, nO, S, T, V, Hc, ERI, dipole_int, X, EROHF, eHF, c, Ptot)
! Perform restricted open-shell Hartree-Fock calculation
@ -16,7 +19,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
double precision,intent(in) :: mix
double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
@ -24,18 +27,18 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
double precision,intent(in) :: ENuc
integer,intent(in) :: nO(nspin)
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 :: Conv
double precision :: rcond
@ -62,9 +65,9 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Output variables
double precision,intent(out) :: EROHF
double precision,intent(out) :: eHF(nBas)
double precision,intent(inout):: c(nBas,nBas)
double precision,intent(out) :: Ptot(nBas,nBas)
double precision,intent(out) :: eHF(nBas_MOs)
double precision,intent(inout):: c(nBas_AOs,nBas_MOs)
double precision,intent(out) :: Ptot(nBas_AOs,nBas_AOs)
! Hello world
@ -76,19 +79,30 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Useful stuff
nBasSq = nBas*nBas
nBas_AOs_Sq = nBas_AOs*nBas_AOs
! Memory allocation
allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), &
P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), &
err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(J(nBas_AOs,nBas_AOs,nspin))
allocate(K(nBas_AOs,nBas_AOs,nspin))
allocate(F(nBas_AOs,nBas_AOs,nspin))
allocate(Ftot(nBas_AOs,nBas_AOs))
allocate(P(nBas_AOs,nBas_AOs,nspin))
allocate(err(nBas_AOs,nBas_AOs))
allocate(Fp(nBas_MOs,nBas_MOs))
allocate(cp(nBas_MOs,nBas_MOs))
allocate(err_diis(nBas_AOs_Sq,max_diis))
allocate(F_diis(nBas_AOs_Sq,max_diis))
! Guess coefficients and demsity matrices
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))))
call mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c)
do ispin = 1, nspin
!P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin))))
call dgemm('N', 'T', nBas_AOs, nBas_AOs, nO(ispin), 1.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P(1,1,ispin), nBas_AOs)
end do
Ptot(:,:) = P(:,:,1) + P(:,:,2)
@ -120,51 +134,51 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Build Hartree repulsion
do ispin=1,nspin
call Hartree_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin))
do ispin = 1, nspin
call Hartree_matrix_AO_basis(nBas_AOs, P(:,:,ispin), ERI(:,:,:,:), J(:,:,ispin))
end do
! Compute exchange potential
do ispin=1,nspin
call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),K(:,:,ispin))
do ispin = 1, nspin
call exchange_matrix_AO_basis(nBas_AOs, P(:,:,ispin), ERI(:,:,:,:), K(:,:,ispin))
end do
! Build Fock operator
do ispin=1,nspin
do ispin = 1, nspin
F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin)
end do
call ROHF_fock_matrix(nBas,nO(1),nO(2),S,c,F(:,:,1),F(:,:,2),Ftot)
call ROHF_fock_matrix(nBas_AOs, nBas_MOs, nO(1), nO(2), S, c, F(:,:,1), F(:,:,2), Ftot)
! 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(:,:)))
! Kinetic energy
do ispin=1,nspin
ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:)))
do ispin = 1, nspin
ET(ispin) = trace_matrix(nBas_AOs, matmul(P(:,:,ispin), T(:,:)))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
do ispin = 1, nspin
EV(ispin) = trace_matrix(nBas_AOs, matmul(P(:,:,ispin), V(:,:)))
end do
! Hartree energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
EJ(1) = 0.5d0*trace_matrix(nBas_AOs, matmul(P(:,:,1), J(:,:,1)))
EJ(2) = trace_matrix(nBas_AOs, matmul(P(:,:,1), J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas_AOs, matmul(P(:,:,2), J(:,:,2)))
! Exchange energy
do ispin=1,nspin
EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
do ispin = 1, nspin
EK(ispin) = 0.5d0*trace_matrix(nBas_AOs, matmul(P(:,:,ispin), K(:,:,ispin)))
end do
! Total energy
@ -176,7 +190,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
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,Ftot)
call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,err_diis,F_diis,err,Ftot)
end if
@ -185,28 +199,29 @@ 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,nBas,maxval(nO),S,c,Ftot)
call level_shifting(level_shift, nBas_AOs, nBas_MOs, maxval(nO), S, c, Ftot)
end do
end if
! 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
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,eHF)
call diagonalize_matrix(nBas_MOs, cp, eHF)
! Back-transform eigenvectors in non-orthogonal basis
c(:,:) = matmul(X(:,:),cp(:,:))
c(:,:) = matmul(X(:,:), cp(:,:))
! Compute density matrix
do ispin=1,nspin
P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin))))
do ispin = 1, nspin
!P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin))))
call dgemm('N', 'T', nBas_AOs, nBas_AOs, nO(ispin), 1.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P(1,1,ispin), nBas_AOs)
end do
Ptot(:,:) = P(:,:,1) + P(:,:,2)
@ -231,14 +246,16 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis)
stop
end if
! Compute final UHF energy
call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,EK,EROHF,dipole)
call dipole_moment(nBas_AOs,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_ROHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, EK, EROHF, dipole)
! Print test values
@ -248,4 +265,6 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
end if
deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis)
end subroutine

35
src/HF/ROHF_fock_matrix.f90
View File

@ -1,4 +1,7 @@
subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
! ---
subroutine ROHF_fock_matrix(nBas_AOs, nBas_MOs, nOa, nOb, S, c, FaAO, FbAO, FAO)
! Construct the ROHF Fock matrix in the AO basis
! For open shells, the ROHF Fock matrix in the MO basis reads
@ -17,14 +20,14 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nOa
integer,intent(in) :: nOb
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(inout):: FaAO(nBas,nBas)
double precision,intent(inout):: FbAO(nBas,nBas)
double precision,intent(in) :: S(nBas_AOs,nBas_AOs)
double precision,intent(in) :: c(nBas_AOs,nBas_MOs)
double precision,intent(inout):: FaAO(nBas_AOs,nBas_AOs)
double precision,intent(inout):: FbAO(nBas_AOs,nBas_AOs)
! Local variables
@ -42,11 +45,11 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
! Output variables
double precision,intent(out) :: FAO(nBas,nBas)
double precision,intent(out) :: FAO(nBas_AOs,nBas_AOs)
! Memory allocation
allocate(F(nBas,nBas),Fa(nBas,nBas),Fb(nBas,nBas))
allocate(F(nBas_MOs,nBas_MOs), Fa(nBas_MOs,nBas_MOs), Fb(nBas_MOs,nBas_MOs))
! Roothan canonicalization parameters
@ -61,14 +64,14 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
! Number of closed, open, and virtual orbitals
nC = min(nOa,nOb)
nC = min(nOa, nOb)
nO = abs(nOa - nOb)
nV = nBas - nC - nO
nV = nBas_AOs - nC - nO
! Block-by-block Fock matrix
call AOtoMO(nBas,c,FaAO,Fa)
call AOtoMO(nBas,c,FbAO,Fb)
call AOtoMO(nBas_AOs, nBas_MOs, c, FaAO, Fa)
call AOtoMO(nBas_AOs, nBas_MOs, c, FbAO, Fb)
F(1:nC, 1:nC ) = aC*Fa(1:nC, 1:nC ) + bC*Fb(1:nC, 1:nC )
F(1:nC, nC+1:nC+nO ) = Fb(1:nC, nC+1:nC+nO )
@ -82,8 +85,10 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO)
F(nO+nC+1:nC+nO+nV, nC+1:nC+nO ) = Fa(nO+nC+1:nC+nO+nV, nC+1:nC+nO )
F(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) = aV*Fa(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) + bV*Fb(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV)
call MOtoAO(nBas,S,c,F,FAO)
call MOtoAO(nBas,S,c,Fa,FaAO)
call MOtoAO(nBas,S,c,Fb,FbAO)
call MOtoAO(nBas_AOs, nBas_MOs, S, c, F, FAO)
call MOtoAO(nBas_AOs, nBas_MOs, S, c, Fa, FaAO)
call MOtoAO(nBas_AOs, nBas_MOs, S, c, Fb, FbAO)
deallocate(F, Fa, Fb)
end subroutine

4
src/HF/UHF_search.f90
View File

@ -124,8 +124,8 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block

17
src/HF/print_ROHF.f90
View File

@ -1,14 +1,17 @@
subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
! ---
subroutine print_ROHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, Ex, EROHF, dipole)
! Print one- and two-electron energies and other stuff for RoHF calculation
implicit none
include 'parameters.h'
integer,intent(in) :: nBas
integer,intent(in) :: nBas_AOs, nBas_MOs
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: eHF(nBas_MOs)
double precision,intent(in) :: c(nBas_AOs,nBas_MOs)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET(nspin)
double precision,intent(in) :: EV(nspin)
@ -31,7 +34,7 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
do ispin=1,nspin
if(nO(ispin) > 0) then
HOMO(ispin) = eHF(nO(ispin))
if(nO(ispin) < nBas) then
if(nO(ispin) < nBas_MOs) then
LUMO(ispin) = eHF(nO(ispin)+1)
else
LUMO(ispin) = 0d0
@ -102,13 +105,13 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') 'ROHF orbital coefficients '
write(*,'(A50)') '-----------------------------------------'
call matout(nBas,nBas,c)
call matout(nBas_AOs, nBas_MOs, c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' ROHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------'
call vecout(nBas,eHF)
call vecout(nBas_MOs, eHF)
write(*,*)
end subroutine

9
src/QuAcK/RQuAcK.f90
View File

@ -112,7 +112,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
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(dipole_int_MO(nBas_MOs,nBas_MOs,ncart))
allocate(ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs))
!---------------------!
@ -122,7 +122,6 @@ 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_AOs, nBas_MOs, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF)
call wall_time(end_HF)
@ -137,8 +136,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
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_AOs,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF)
call ROHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, &
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
@ -161,7 +160,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d
do ixyz = 1, ncart
! TODO
call AOtoMO(nBas_AOs,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
call AOtoMO(nBas_AOs,nBas_MOs,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz))
end do
! 4-index transform

4
src/QuAcK/UQuAcK.f90
View File

@ -163,8 +163,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do
! Read and transform dipole-related integrals
do ixyz=1,ncart
call AOtoMO(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
call AOtoMO(nBas,nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz))
call AOtoMO(nBas,nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block