From 99cbd0e95d5200a5bd4ccb2bf1284b7e1961e4b3 Mon Sep 17 00:00:00 2001 From: pfloos Date: Tue, 5 Sep 2023 17:23:40 +0200 Subject: [PATCH] ROHF not yet debugged --- input/methods | 6 +- input/options | 6 +- src/CC/CC.f90 | 1 + src/GT/G0T0eh.f90 | 2 +- src/GT/G0T0pp.f90 | 4 +- src/GT/GTeh_plot_self_energy.f90 | 2 +- src/GT/GTpp_plot_self_energy.f90 | 2 +- src/GW/G0W0.f90 | 2 +- src/GW/GW_plot_self_energy.f90 | 2 +- src/HF/HF.f90 | 27 +++- src/HF/ROHF.f90 | 244 +++++++++++++++++++++++++++++++ src/HF/ROHF_fock_matrix.f90 | 71 +++++++++ src/HF/UHF.f90 | 10 +- src/QuAcK/QuAcK.f90 | 30 ++-- src/QuAcK/read_methods.f90 | 34 +++-- 15 files changed, 391 insertions(+), 52 deletions(-) create mode 100644 src/HF/ROHF.f90 create mode 100644 src/HF/ROHF_fock_matrix.f90 diff --git a/input/methods b/input/methods index 16a717e..9843b39 100644 --- a/input/methods +++ b/input/methods @@ -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 diff --git a/input/options b/input/options index acdfcbb..0a181b9 100644 --- a/input/options +++ b/input/options @@ -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 diff --git a/src/CC/CC.f90 b/src/CC/CC.f90 index 2a90da6..411b7c5 100644 --- a/src/CC/CC.f90 +++ b/src/CC/CC.f90 @@ -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 diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index 5de9200..483b39c 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -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 diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index 63c0ddb..56a7a40 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -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 diff --git a/src/GT/GTeh_plot_self_energy.f90 b/src/GT/GTeh_plot_self_energy.f90 index f1a07f3..54d6bca 100644 --- a/src/GT/GTeh_plot_self_energy.f90 +++ b/src/GT/GTeh_plot_self_energy.f90 @@ -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 diff --git a/src/GT/GTpp_plot_self_energy.f90 b/src/GT/GTpp_plot_self_energy.f90 index 62f418d..668a7a3 100644 --- a/src/GT/GTpp_plot_self_energy.f90 +++ b/src/GT/GTpp_plot_self_energy.f90 @@ -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 diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 6827e01..3b2ad07 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -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 diff --git a/src/GW/GW_plot_self_energy.f90 b/src/GW/GW_plot_self_energy.f90 index 0f079f1..9574a73 100644 --- a/src/GW/GW_plot_self_energy.f90 +++ b/src/GW/GW_plot_self_energy.f90 @@ -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 diff --git a/src/HF/HF.f90 b/src/HF/HF.f90 index 7e074af..1b62bda 100644 --- a/src/HF/HF.f90 +++ b/src/HF/HF.f90 @@ -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 diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 new file mode 100644 index 0000000..2981141 --- /dev/null +++ b/src/HF/ROHF.f90 @@ -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 diff --git a/src/HF/ROHF_fock_matrix.f90 b/src/HF/ROHF_fock_matrix.f90 new file mode 100644 index 0000000..8c6de5e --- /dev/null +++ b/src/HF/ROHF_fock_matrix.f90 @@ -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 diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index dad1b69..4f2a900 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -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 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 097010a..40d2e79 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -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,17 +123,17 @@ program QuAcK ! Method selection ! !------------------! - call read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, & - doMP2,doMP3, & - doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - dodrCCD,dorCCD,docrCCD,dolCCD, & - doCIS,doCIS_D,doCID,doCISD,doFCI, & - dophRPA,dophRPAx,docrRPA,doppRPA, & - doG0F2,doevGF2,doqsGF2, & - doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW,doSRGqsGW, & - doufG0W0,doufGW, & - doG0T0pp,doevGTpp,doqsGTpp, & + call read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & + doMP2,doMP3, & + doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + dodrCCD,dorCCD,docrCCD,dolCCD, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & + dophRPA,dophRPAx,docrRPA,doppRPA, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW,doSRGqsGW, & + doufG0W0,doufGW, & + doG0T0pp,doevGTpp,doqsGTpp, & doG0T0eh,doevGTeh,doqsGTeh) !--------------------------! @@ -209,13 +209,13 @@ 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, & - guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO, & + 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) diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 1626798..0ccd26d 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,14 +1,14 @@ -subroutine read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, & - doMP2,doMP3, & - doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - do_drCCD,do_rCCD,do_crCCD,do_lCCD, & - doCIS,doCIS_D,doCID,doCISD,doFCI, & - dophRPA,dophRPAx,docrRPA,doppRPA, & - doG0F2,doevGF2,doqsGF2, & - doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW,doSRGqsGW, & - doufG0W0,doufGW, & - doG0T0pp,doevGTpp,doqsGTpp, & +subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & + doMP2,doMP3, & + doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & + dophRPA,dophRPAx,docrRPA,doppRPA, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW,doSRGqsGW, & + doufG0W0,doufGW, & + doG0T0pp,doevGTpp,doqsGTpp, & doG0T0eh,doevGTeh,doqsGTeh) ! Read desired methods @@ -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