mirror of
https://github.com/pfloos/quack
synced 2024-12-22 20:34:46 +01:00
ROHF not yet debugged
This commit is contained in:
parent
ea1fd8eebe
commit
99cbd0e95d
@ -1,5 +1,5 @@
|
||||
# RHF UHF RMOM UMOM KS
|
||||
T F F F F
|
||||
# RHF UHF ROHF RMOM UMOM KS
|
||||
T F T F F F
|
||||
# MP2* MP3
|
||||
F F
|
||||
# CCD pCCD DCD CCSD CCSD(T)
|
||||
@ -15,5 +15,5 @@
|
||||
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
|
||||
F F F F F F
|
||||
# G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
|
||||
F F F T F F
|
||||
F F F F F F
|
||||
# * unrestricted version available
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability
|
||||
512 0.0000001 T 5 1 1 F 0.0 F
|
||||
64 0.0000001 F 5 1 1 F 0.0 F
|
||||
# MP: reg
|
||||
F
|
||||
# CC: maxSCF thresh DIIS n_diis
|
||||
@ -9,9 +9,9 @@
|
||||
# GF: maxSCF thresh DIIS n_diis lin eta renorm reg
|
||||
256 0.00001 T 5 F 0.0 0 F
|
||||
# GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg
|
||||
256 0.00001 T 5 F 0.0 F F
|
||||
256 0.00001 T 5 F 0.001 F F
|
||||
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
|
||||
256 0.00001 T 5 F 0.0 F F
|
||||
256 0.00001 T 5 F 0.001 F F
|
||||
# ACFDT: AC Kx XBS
|
||||
F F T
|
||||
# BSE: phBSE phBSE2 ppBSE dBSE dTDA
|
||||
|
||||
@ -159,6 +159,7 @@ subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,do_drCCD,do_rCCD,do_crCCD,do_lCC
|
||||
|
||||
call wall_time(start_CC)
|
||||
call pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
call pCCDso(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
|
||||
call wall_time(end_CC)
|
||||
|
||||
t_CC = end_CC - start_CC
|
||||
|
||||
@ -149,7 +149,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,
|
||||
|
||||
end if
|
||||
|
||||
call GTeh_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rhoL,rhoR)
|
||||
! call GTeh_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rhoL,rhoR)
|
||||
|
||||
! Compute the RPA correlation energy based on the G0T0eh quasiparticle energies
|
||||
|
||||
|
||||
@ -197,8 +197,8 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,dopp
|
||||
|
||||
end if
|
||||
|
||||
call GTpp_plot_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, &
|
||||
Om1t,rho1t,Om2t,rho2t)
|
||||
! call GTpp_plot_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, &
|
||||
! Om1t,rho1t,Om2t,rho2t)
|
||||
|
||||
!----------------------------------------------
|
||||
! Dump results
|
||||
|
||||
@ -33,7 +33,7 @@ subroutine GTeh_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGT,Om,rhoL,rhoR)
|
||||
|
||||
! Construct grid
|
||||
|
||||
nGrid = 1000
|
||||
nGrid = 5000
|
||||
allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid))
|
||||
|
||||
! Initialize
|
||||
|
||||
@ -37,7 +37,7 @@ subroutine GTpp_plot_self_energy(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eG
|
||||
|
||||
! Construct grid
|
||||
|
||||
nGrid = 1000
|
||||
nGrid = 5000
|
||||
allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid))
|
||||
|
||||
! Initialize
|
||||
|
||||
@ -140,7 +140,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT
|
||||
|
||||
end if
|
||||
|
||||
call GW_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
|
||||
! call GW_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eHF,Om,rho)
|
||||
|
||||
! Compute the RPA correlation energy
|
||||
|
||||
|
||||
@ -32,7 +32,7 @@ subroutine GW_plot_self_energy(eta,nBas,nC,nO,nV,nR,nS,eHF,eGW,Om,rho)
|
||||
|
||||
! Construct grid
|
||||
|
||||
nGrid = 1000
|
||||
nGrid = 5000
|
||||
allocate(w(nGrid),SigC(nBas,nGrid),Z(nBas,nGrid),S(nBas,nGrid))
|
||||
|
||||
! Initialize
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,guess_type,mix,level_shift, &
|
||||
subroutine HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,guess_type,mix,level_shift, &
|
||||
nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
|
||||
|
||||
! Hartree-Fock module
|
||||
@ -10,6 +10,7 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues
|
||||
|
||||
logical,intent(in) :: doRHF
|
||||
logical,intent(in) :: doUHF
|
||||
logical,intent(in) :: doROHF
|
||||
logical,intent(in) :: doRMOM
|
||||
logical,intent(in) :: doUMOM
|
||||
|
||||
@ -87,7 +88,7 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues
|
||||
call wall_time(end_HF)
|
||||
|
||||
t_HF = end_HF - start_HF
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for HF = ',t_HF,' seconds'
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
@ -107,7 +108,27 @@ subroutine HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,gues
|
||||
call wall_time(end_HF)
|
||||
|
||||
t_HF = end_HF - start_HF
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for HF = ',t_HF,' seconds'
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for UHF = ',t_HF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute ROHF energy
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doROHF) then
|
||||
|
||||
! Switch on the unrestricted flag
|
||||
unrestricted = .true.
|
||||
|
||||
call wall_time(start_HF)
|
||||
call ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
|
||||
call wall_time(end_HF)
|
||||
|
||||
t_HF = end_HF - start_HF
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ROHF = ',t_HF,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
244
src/HF/ROHF.f90
Normal file
244
src/HF/ROHF.f90
Normal file
@ -0,0 +1,244 @@
|
||||
subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P)
|
||||
|
||||
! Perform unrestricted Hartree-Fock calculation
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: maxSCF
|
||||
integer,intent(in) :: max_diis
|
||||
integer,intent(in) :: guess_type
|
||||
logical,intent(in) :: mix
|
||||
double precision,intent(in) :: level_shift
|
||||
double precision,intent(in) :: thresh
|
||||
integer,intent(in) :: nBas
|
||||
|
||||
integer,intent(in) :: nNuc
|
||||
double precision,intent(in) :: ZNuc(nNuc)
|
||||
double precision,intent(in) :: rNuc(nNuc,ncart)
|
||||
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)
|
||||
|
||||
! Local variables
|
||||
|
||||
integer :: nSCF
|
||||
integer :: nBasSq
|
||||
integer :: n_diis
|
||||
double precision :: conv
|
||||
double precision :: rcond
|
||||
double precision :: ET(nspin)
|
||||
double precision :: EV(nspin)
|
||||
double precision :: EJ(nsp)
|
||||
double precision :: Ex(nspin)
|
||||
double precision :: dipole(ncart)
|
||||
|
||||
double precision,allocatable :: cp(:,:)
|
||||
double precision,allocatable :: J(:,:,:)
|
||||
double precision,allocatable :: F(:,:,:)
|
||||
double precision,allocatable :: Fp(:,:)
|
||||
double precision,allocatable :: Ftot(:,:)
|
||||
double precision,allocatable :: Ptot(:,:)
|
||||
double precision,allocatable :: K(:,:,:)
|
||||
double precision,allocatable :: err(:,:)
|
||||
double precision,allocatable :: err_diis(:,:)
|
||||
double precision,allocatable :: F_diis(:,:)
|
||||
double precision,external :: trace_matrix
|
||||
|
||||
integer :: ispin
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EHF
|
||||
double precision,intent(out) :: e(nBas)
|
||||
double precision,intent(out) :: c(nBas,nBas)
|
||||
double precision,intent(out) :: P(nBas,nBas,nspin)
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)'* Restricted Open-Shell Hartree-Fock *'
|
||||
write(*,*)'************************************************'
|
||||
write(*,*)
|
||||
|
||||
! Useful stuff
|
||||
|
||||
nBasSq = nBas*nBas
|
||||
|
||||
! Memory allocation
|
||||
|
||||
allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), &
|
||||
Ptot(nBas,nBas),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), &
|
||||
err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
|
||||
|
||||
! Guess coefficients and demsity matrices
|
||||
|
||||
call mo_guess(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
|
||||
|
||||
! Initialization
|
||||
|
||||
n_diis = 0
|
||||
F_diis(:,:) = 0d0
|
||||
err_diis(:,:) = 0d0
|
||||
rcond = 0d0
|
||||
|
||||
nSCF = 0
|
||||
conv = 1d0
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Main SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'----------------------------------------------------------'
|
||||
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
|
||||
'|','#','|','E(ROHF)','|','Ex(ROHF)','|','Conv','|'
|
||||
write(*,*)'----------------------------------------------------------'
|
||||
|
||||
do while(conv > thresh .and. nSCF < maxSCF)
|
||||
|
||||
! Increment
|
||||
|
||||
nSCF = nSCF + 1
|
||||
|
||||
! Build Coulomb repulsion
|
||||
|
||||
do ispin=1,nspin
|
||||
call Coulomb_matrix_AO_basis(nBas,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))
|
||||
end do
|
||||
|
||||
! Build Fock operator
|
||||
|
||||
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)
|
||||
|
||||
! Check convergence
|
||||
|
||||
err(:,:) = matmul(Ftot(:,:),matmul(Ptot(:,:),S(:,:))) - matmul(matmul(S(:,:),Ptot(:,:)),Ftot(:,:))
|
||||
if(nSCF > 1) conv = maxval(abs(err(:,:)))
|
||||
|
||||
! DIIS extrapolation
|
||||
|
||||
if(max_diis > 1) then
|
||||
|
||||
n_diis = min(n_diis+1,max_diis)
|
||||
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis),F_diis(:,1:n_diis),err,Ftot)
|
||||
|
||||
end if
|
||||
|
||||
! Level-shifting
|
||||
|
||||
if(level_shift > 0d0 .and. Conv > thresh) then
|
||||
|
||||
do ispin=1,nspin
|
||||
call level_shifting(level_shift,nBas,maxval(nO),S,c,Ftot)
|
||||
end do
|
||||
|
||||
end if
|
||||
|
||||
! Transform Fock matrix in orthogonal basis
|
||||
|
||||
Fp(:,:) = matmul(transpose(X(:,:)),matmul(Ftot(:,:),X(:,:)))
|
||||
|
||||
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
|
||||
|
||||
cp(:,:) = Ftot(:,:)
|
||||
call diagonalize_matrix(nBas,cp,e)
|
||||
|
||||
! Back-transform eigenvectors in non-orthogonal basis
|
||||
|
||||
c(:,:) = matmul(X(:,:),cp(:,:))
|
||||
|
||||
! Compute density matrix
|
||||
|
||||
do ispin=1,nspin
|
||||
P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin))))
|
||||
end do
|
||||
Ptot(:,:) = P(:,:,1) + P(:,:,2)
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute ROHF energy
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Kinetic energy
|
||||
|
||||
do ispin=1,nspin
|
||||
ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:)))
|
||||
end do
|
||||
|
||||
! Potential energy
|
||||
|
||||
do ispin=1,nspin
|
||||
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
|
||||
end do
|
||||
|
||||
! Coulomb 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)))
|
||||
|
||||
! Exchange energy
|
||||
|
||||
do ispin=1,nspin
|
||||
Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
|
||||
end do
|
||||
|
||||
! Total energy
|
||||
|
||||
EHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
|
||||
|
||||
! Dump results
|
||||
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') &
|
||||
'|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',conv,'|'
|
||||
|
||||
end do
|
||||
write(*,*)'----------------------------------------------------------'
|
||||
!------------------------------------------------------------------------
|
||||
! End of SCF loop
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
! Did it actually converge?
|
||||
|
||||
if(nSCF == maxSCF) then
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)' Convergence failed '
|
||||
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
write(*,*)
|
||||
|
||||
stop
|
||||
|
||||
end if
|
||||
|
||||
! Compute final UHF energy
|
||||
|
||||
! call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
|
||||
! call print_ROHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
|
||||
|
||||
end subroutine
|
||||
71
src/HF/ROHF_fock_matrix.f90
Normal file
71
src/HF/ROHF_fock_matrix.f90
Normal file
@ -0,0 +1,71 @@
|
||||
subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,Fa,Fb,F)
|
||||
|
||||
! Construct the ROHF Fock matrix
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
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):: Fa(nBas,nBas)
|
||||
double precision,intent(inout):: Fb(nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: aC,bC
|
||||
double precision :: aO,bO
|
||||
double precision :: aV,bV
|
||||
|
||||
integer :: nC
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: F(nBas,nBas)
|
||||
|
||||
! Roothan canonicalization parameters
|
||||
|
||||
aC = -0.5d0
|
||||
bC = +1.5d0
|
||||
|
||||
aO = +0.5d0
|
||||
bO = +0.5d0
|
||||
|
||||
aV = +1.5d0
|
||||
bV = -0.5d0
|
||||
|
||||
! Number of closed, open, and virtual orbitals
|
||||
|
||||
nC = min(nOa,nOb)
|
||||
nO = abs(nOa - nOb)
|
||||
nV = nBas - nC - nO
|
||||
|
||||
! Block-by-block Fock matrix
|
||||
|
||||
call AOtoMO_transform(nBas,c,Fa)
|
||||
call AOtoMO_transform(nBas,c,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 )
|
||||
F(1:nC,nO+nC+1:nC+nO+nV) = 0.5d0*Fa(1:nC,nO+nC+1:nC+nO+nV) + 0.5d0*Fb(1:nC,nO+nC+1:nC+nO+nV)
|
||||
|
||||
F(nC+1:nO+nC, 1:nC ) = Fb(nC+1:nO+nC, 1:nC )
|
||||
F(nC+1:nO+nC, nC+1:nC+nO ) = aO*Fa(nC+1:nC+nO, nC+1:nC+nO ) + bO*Fb(nC+1:nO+nC,nC+1:nC+nO)
|
||||
F(nC+1:nO+nC,nO+nC+1:nC+nO+nV) = Fa(nC+1:nC+nO,nO+nC+1:nC+nO+nV)
|
||||
|
||||
F(nO+nC+1:nC+nO+nV, 1:nC ) = 0.5d0*Fa(nO+nC+1:nC+nO+nV, 1:nC ) + 0.5d0*Fb(nO+nC+1:nC+nO+nV, 1:nC )
|
||||
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_transform(nBas,S,c,F)
|
||||
call MOtoAO_transform(nBas,S,c,Fa)
|
||||
call MOtoAO_transform(nBas,S,c,Fb)
|
||||
|
||||
end subroutine
|
||||
@ -1,5 +1,5 @@
|
||||
subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
|
||||
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P)
|
||||
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P)
|
||||
|
||||
! Perform unrestricted Hartree-Fock calculation
|
||||
|
||||
@ -57,7 +57,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: EUHF
|
||||
double precision,intent(out) :: EHF
|
||||
double precision,intent(out) :: e(nBas,nspin)
|
||||
double precision,intent(out) :: c(nBas,nBas,nspin)
|
||||
double precision,intent(out) :: P(nBas,nBas,nspin)
|
||||
@ -219,12 +219,12 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
|
||||
! Total energy
|
||||
|
||||
EUHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
|
||||
EHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
|
||||
|
||||
! Dump results
|
||||
|
||||
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') &
|
||||
'|',nSCF,'|',EUHF + ENuc,'|',sum(Ex(:)),'|',conv,'|'
|
||||
'|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',conv,'|'
|
||||
|
||||
end do
|
||||
write(*,*)'----------------------------------------------------------'
|
||||
@ -249,6 +249,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
|
||||
! Compute final UHF energy
|
||||
|
||||
call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
|
||||
call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
|
||||
call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
|
||||
|
||||
end subroutine
|
||||
|
||||
@ -4,7 +4,7 @@ program QuAcK
|
||||
include 'parameters.h'
|
||||
|
||||
logical :: unrestricted = .false.
|
||||
logical :: doHF,doRHF,doUHF,doRMOM,doUMOM
|
||||
logical :: doHF,doRHF,doUHF,doROHF,doRMOM,doUMOM
|
||||
logical :: dostab
|
||||
logical :: doKS
|
||||
logical :: doMP,doMP2,doMP3
|
||||
@ -123,7 +123,7 @@ program QuAcK
|
||||
! Method selection !
|
||||
!------------------!
|
||||
|
||||
call read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, &
|
||||
call read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
|
||||
doMP2,doMP3, &
|
||||
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
|
||||
dodrCCD,dorCCD,docrCCD,dolCCD, &
|
||||
@ -209,12 +209,12 @@ program QuAcK
|
||||
! Hartree-Fock module !
|
||||
!---------------------!
|
||||
|
||||
doHF = doRHF .or. doUHF .or. doRMOM .or. doUMOM
|
||||
doHF = doRHF .or. doUHF .or. doROHF .or. doRMOM .or. doUMOM
|
||||
|
||||
if(doHF) then
|
||||
|
||||
call wall_time(start_HF)
|
||||
call HF(doRHF,doUHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,max_diis_HF, &
|
||||
call HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,max_diis_HF, &
|
||||
guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO, &
|
||||
ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
|
||||
call wall_time(end_HF)
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
subroutine read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, &
|
||||
subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
|
||||
doMP2,doMP3, &
|
||||
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
|
||||
do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
|
||||
@ -17,7 +17,7 @@ subroutine read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, &
|
||||
|
||||
! Input variables
|
||||
|
||||
logical,intent(out) :: doRHF,doUHF,doRMOM,doUMOM,doKS
|
||||
logical,intent(out) :: doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS
|
||||
logical,intent(out) :: doMP2,doMP3
|
||||
logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
|
||||
logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
|
||||
@ -40,17 +40,19 @@ subroutine read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, &
|
||||
|
||||
doRHF = .false.
|
||||
doUHF = .false.
|
||||
doROHF = .false.
|
||||
doRMOM = .false.
|
||||
doUMOM = .false.
|
||||
doKS = .false.
|
||||
|
||||
read(1,*)
|
||||
read(1,*) answer1,answer2,answer3,answer4,answer5
|
||||
read(1,*) answer1,answer2,answer3,answer4,answer5,answer6
|
||||
if(answer1 == 'T') doRHF = .true.
|
||||
if(answer2 == 'T') doUHF = .true.
|
||||
if(answer3 == 'T') doRMOM = .true.
|
||||
if(answer4 == 'T') doUMOM = .true.
|
||||
if(answer5 == 'T') doKS = .true.
|
||||
if(answer3 == 'T') doROHF = .true.
|
||||
if(answer4 == 'T') doRMOM = .true.
|
||||
if(answer5 == 'T') doUMOM = .true.
|
||||
if(answer6 == 'T') doKS = .true.
|
||||
|
||||
! Read MPn methods
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user