From 60269ea3cec2f8c4a49143bf08b83ff0e195b59d Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 19 Mar 2019 09:30:14 +0100 Subject: [PATCH] UHF code --- input/basis | 29 +++-- input/methods | 6 +- input/molecule | 4 +- input/weight | 29 +++-- src/MCQC/CCD.f90 | 14 ++- src/MCQC/EOM_CCSD.f90 | 78 -------------- src/MCQC/UHF.f90 | 237 +++++++++++++++++++++++++++++++++++++++++ src/MCQC/form_u.f90 | 71 ++++++++++++ src/MCQC/print_UHF.f90 | 102 ++++++++++++++++++ 9 files changed, 467 insertions(+), 103 deletions(-) delete mode 100644 src/MCQC/EOM_CCSD.f90 create mode 100644 src/MCQC/UHF.f90 create mode 100644 src/MCQC/form_u.f90 create mode 100644 src/MCQC/print_UHF.f90 diff --git a/input/basis b/input/basis index b9ca7b5..c12e4d6 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,20 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 -S 1 1.00 - 0.2976000 1.0000000 -P 1 1.00 - 1.2750000 1.0000000 +1 5 +S 6 1.00 + 1264.5857000 0.0019448 + 189.9368100 0.0148351 + 43.1590890 0.0720906 + 12.0986630 0.2371542 + 3.8063232 0.4691987 + 1.2728903 0.3565202 +S 3 1.00 + 3.1964631 -0.1126487 + 0.7478133 -0.2295064 + 0.2199663 1.1869167 +S 1 1.00 + 0.0823099 1.0000000 +P 3 1.00 + 3.1964631 0.0559802 + 0.7478133 0.2615506 + 0.2199663 0.7939723 +P 1 1.00 + 0.0823099 1.0000000 diff --git a/input/methods b/input/methods index 722184a..10d4455 100644 --- a/input/methods +++ b/input/methods @@ -1,14 +1,14 @@ # HF MOM T F # MP2 MP3 MP2-F12 - F F F + T T F # CCD CCSD CCSD(T) - F F F + T F F # CIS TDHF ADC F F F # GF2 GF3 F F # G0W0 evGW qsGW - T F F + F F F # MCMP2 F diff --git a/input/molecule b/input/molecule index 7d017f4..67ee9ac 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEl nCore nRyd - 1 2 0 0 + 1 4 0 0 # Znuc x y z - He 0.0 0.0 0.0 + Be 0.0 0.0 0.0 diff --git a/input/weight b/input/weight index b9ca7b5..c12e4d6 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,20 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 -S 1 1.00 - 0.2976000 1.0000000 -P 1 1.00 - 1.2750000 1.0000000 +1 5 +S 6 1.00 + 1264.5857000 0.0019448 + 189.9368100 0.0148351 + 43.1590890 0.0720906 + 12.0986630 0.2371542 + 3.8063232 0.4691987 + 1.2728903 0.3565202 +S 3 1.00 + 3.1964631 -0.1126487 + 0.7478133 -0.2295064 + 0.2199663 1.1869167 +S 1 1.00 + 0.0823099 1.0000000 +P 3 1.00 + 3.1964631 0.0559802 + 0.7478133 0.2615506 + 0.2199663 0.7939723 +P 1 1.00 + 0.0823099 1.0000000 diff --git a/src/MCQC/CCD.f90 b/src/MCQC/CCD.f90 index de3ff2c..9caf576 100644 --- a/src/MCQC/CCD.f90 +++ b/src/MCQC/CCD.f90 @@ -22,7 +22,7 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) integer :: nV integer :: nSCF double precision :: Conv - double precision :: EcMP2 + double precision :: EcMP2,EcMP3,EcMP4 double precision :: ECCD,EcCCD double precision,allocatable :: seHF(:) double precision,allocatable :: sERI(:,:,:,:) @@ -108,7 +108,7 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) - write(*,'(1X,A10,1X,F10.6)') 'Ec(MP2) = ',EcMP2 + EcMP4 = 0d0 ! Initialization @@ -163,6 +163,8 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.)) + ! Dump results ECCD = ERHF + EcCCD @@ -190,4 +192,12 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) endif +! Moller-Plesset energies + + write(*,*) + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2 + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3 + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4 + write(*,*) + end subroutine CCD diff --git a/src/MCQC/EOM_CCSD.f90 b/src/MCQC/EOM_CCSD.f90 deleted file mode 100644 index 2be73aa..0000000 --- a/src/MCQC/EOM_CCSD.f90 +++ /dev/null @@ -1,78 +0,0 @@ -subroutine EOM_CCSD(ispin,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,e,ERI) - -! Compute EOM-CCSD excitation energies: see Stanton & Bartlett JCP 98 7029 (1993) - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: ispin - integer,intent(in) :: maxSCF - double precision,intent(in) :: thresh - integer,intent(in) :: max_diis - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) - -! Local variables - - integer :: nH,nP,nHH,nPP,nSCF,n_diis - double precision :: Conv - double precision,external :: Kronecker_delta - double precision,allocatable :: B_ADC(:,:),X_ADC(:,:),e_ADC(:),SigInf(:,:),G_ADC(:,:) - double precision,allocatable :: db_ERI(:,:,:,:),eOld(:),error_diis(:,:),e_diis(:,:) - - integer :: i,j,k,l - integer :: a,b,c,d - integer :: p,q,r,s - integer :: nADC,iADC,jADC - - -! Hello world - - write(*,*) - write(*,*)'***********************************' - write(*,*)'| EOM-CCSD calculation |' - write(*,*)'***********************************' - write(*,*) - -! Number of holes - - nH = nO - nHH = nH*nH - -! Number of particles - - nP = nV - nPP = nP*nP - - write(*,*) 'Total states: ',nH + nP - write(*,*) 'Hole states: ',nH - write(*,*) 'Particle states: ',nP - -! Size of EOM-CCSD matrices - - nEOM = nH + nP + nH*nPP + nHH*nP + nHH*nPP - write(*,'(1X,A25,I3,A6,I6)') 'Size of EOM-CCSD matrix: ',nEOM,' x ',nEOM - -! Memory allocation - - allocate() - -! Construct EOM-CCSD matrix - - H(:,:) = 0d0 - - iEOM = 1 - jEOM = 1 - - H(iEOM,jEOM) = ECCSD - - do p=1,nO - - jADC = jADC + 1 - B_ADC(jADC,jADC) = e(p) - - enddo - -end subroutine EOM_CCSD diff --git a/src/MCQC/UHF.f90 b/src/MCQC/UHF.f90 new file mode 100644 index 0000000..79dac5d --- /dev/null +++ b/src/MCQC/UHF.f90 @@ -0,0 +1,237 @@ +subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,nV,S,T,V,Hc,ERI,X,ENuc,EUHF) + +! 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 + double precision,intent(in) :: thresh + integer,intent(in) :: nBas + + integer,intent(in) :: nO(nspin),nV(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) :: ENuc + +! Local variables + + integer :: nSCF + integer :: nBasSq + integer :: n_diis + double precision :: conv + double precision :: rcond(nspin) + double precision :: ET(nspin) + double precision :: EV(nspin) + double precision :: EJ(nsp) + double precision :: Ex(nspin) + double precision :: Ec(nsp) + double precision :: EUHF + + double precision,allocatable :: eps(:,:) + double precision,allocatable :: c(:,:,:) + double precision,allocatable :: cp(:,:,:) + double precision,allocatable :: J(:,:,:) + double precision,allocatable :: F(:,:,:) + double precision,allocatable :: Fp(:,:,:) + double precision,allocatable :: Fx(:,:,:) + double precision,allocatable :: err(:,:,:) + double precision,allocatable :: err_diis(:,:,:) + double precision,allocatable :: F_diis(:,:,:) + double precision,external :: trace_matrix + + double precision,allocatable :: P(:,:,:) + + integer :: ispin + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'* Unrestricted Hartree-Fock calculation *' + write(*,*)'************************************************' + write(*,*) + +! Useful stuff + + nBasSq = nBas*nBas + +! Memory allocation + + allocate(eps(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin), & + J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), & + Fx(nBas,nBas,nspin),err(nBas,nBas,nspin),P(nBas,nBas,nspin), & + err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) + +! Guess coefficients and eigenvalues + + if(guess_type == 1) then + + do ispin=1,nspin + F(:,:,ispin) = Hc(:,:) + end do + + else if(guess_type == 2) then + + do ispin=1,nspin + call random_number(F(:,:,ispin)) + end do + + end if + +! Initialization + + nSCF = 0 + conv = 1d0 + + n_diis = 0 + F_diis(:,:,:) = 0d0 + err_diis(:,:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + + write(*,*) + write(*,*)'------------------------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|','nEl','|' + write(*,*)'------------------------------------------------------------------------------------------' + + do while(conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Transform Fock matrix in orthogonal basis + + do ispin=1,nspin + Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:))) + end do + +! Diagonalize Fock matrix to get eigenvectors and eigenvalues + + cp(:,:,:) = Fp(:,:,:) + do ispin=1,nspin + call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin)) + end do + +! Back-transform eigenvectors in non-orthogonal basis + + do ispin=1,nspin + c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin)) + end do + +! Compute density matrix + + do ispin=1,nspin + P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) + end do + +! 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(:,:,:,:),Fx(:,:,ispin)) + end do + +! Build Fock operator + do ispin=1,nspin + F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + Fx(:,:,ispin) + end do + +! Check convergence + + do ispin=1,nspin + err(:,:,ispin) = matmul(F(:,:,ispin),matmul(P(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),P(:,:,ispin)),F(:,:,ispin)) + end do + + conv = maxval(abs(err(:,:,:))) + +! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + do ispin=1,nspin + call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin)) + end do + +! Reset DIIS if required + + if(minval(rcond(:)) < 1d-15) n_diis = 0 + +!------------------------------------------------------------------------ +! Compute UHF 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) = trace_matrix(nBas,matmul(P(:,:,ispin),Fx(:,:,ispin))) + end do + +! Total energy + + EUHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)) + +! Dump results + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',EUHF + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',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 print_UHF(nBas,nO(:),eps(:,:),c(:,:,:),ENuc,ET(:),EV(:),EJ(:),Ex(:),Ec(:),EUHF) + +end subroutine UHF diff --git a/src/MCQC/form_u.f90 b/src/MCQC/form_u.f90 new file mode 100644 index 0000000..c8fc76d --- /dev/null +++ b/src/MCQC/form_u.f90 @@ -0,0 +1,71 @@ +subroutine form_u(nO,nV,OOOO,VVVV,OVOV,t2,u) + +! Form linear array in CCD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: t2(nO,nO,nV,nV) + double precision,intent(in) :: OOOO(nO,nO,nO,nO) + double precision,intent(in) :: VVVV(nV,nV,nV,nV) + double precision,intent(in) :: OVOV(nO,nV,nO,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: u(nO,nO,nV,nV) + + u(:,:,:,:) = 0d0 + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + do c=1,nV + do d=1,nV + u(i,j,a,b) = u(i,j,a,b) + 0.5d0*VVVV(a,b,c,d)*t2(i,j,c,d) + enddo + enddo + enddo + enddo + enddo + enddo + + do i=1,nO + do j=1,nO + do k=1,nO + do l=1,nO + do a=1,nV + do b=1,nV + u(i,j,a,b) = u(i,j,a,b) + 0.5d0*OOOO(k,l,i,j)*t2(k,l,a,b) + enddo + enddo + enddo + enddo + enddo + enddo + + do i=1,nO + do j=1,nO + do k=1,nO + do a=1,nV + do b=1,nV + do c=1,nV + u(i,j,a,b) = u(i,j,a,b) - OVOV(k,b,j,c)*t2(i,k,a,c) & + + OVOV(k,a,j,c)*t2(i,k,b,c) & + - OVOV(k,a,i,c)*t2(j,k,b,c) & + + OVOV(k,b,i,c)*t2(j,k,a,c) + enddo + enddo + enddo + enddo + enddo + enddo + +end subroutine form_u diff --git a/src/MCQC/print_UHF.f90 b/src/MCQC/print_UHF.f90 new file mode 100644 index 0000000..d0a16c1 --- /dev/null +++ b/src/MCQC/print_UHF.f90 @@ -0,0 +1,102 @@ +subroutine print_UHF(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) + +! Print one- and two-electron energies and other stuff for UHF calculation + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas + integer,intent(in) :: nO(nspin) + double precision,intent(in) :: eps(nBas,nspin) + double precision,intent(in) :: c(nBas,nBas,nspin) + double precision,intent(in) :: ENuc + double precision,intent(in) :: ET(nspin) + double precision,intent(in) :: EV(nspin) + double precision,intent(in) :: EJ(nsp) + double precision,intent(in) :: Ex(nspin) + double precision,intent(in) :: Ec(nsp) + double precision,intent(in) :: Ew + + integer :: HOMO(nspin) + integer :: LUMO(nspin) + double precision :: Gap(nspin) + +! HOMO and LUMO + + HOMO(:) = nO(:) + + LUMO(:) = HOMO(:) + 1 + + Gap(1) = eps(LUMO(1),1) - eps(HOMO(1),1) + Gap(2) = eps(LUMO(2),2) - eps(HOMO(2),2) + +! Dump results + + + write(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40)') ' Summary ' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' One-electron energy: ',sum(ET(:)) + sum(EV(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' One-electron a energy: ',ET(1) + EV(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' One-electron b energy: ',ET(2) + EV(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic energy: ',sum(ET(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic a energy: ',ET(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Kinetic b energy: ',ET(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential energy: ',sum(EV(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential a energy: ',EV(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Potential b energy: ',EV(2),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron a energy: ',sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron aa energy: ',EJ(1) + Ex(1) + Ec(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron ab energy: ',EJ(2) + Ec(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Two-electron bb energy: ',EJ(3) + Ex(2) + Ec(3),' au' + write(*,'(A40,1X,F16.10,A3)') ' Coulomb energy: ',sum(EJ(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Coulomb aa energy: ',EJ(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Coulomb ab energy: ',EJ(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Coulomb bb energy: ',EJ(3),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange energy: ',sum(Ex(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange a energy: ',Ex(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Exchange b energy: ',Ex(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation energy: ',sum(Ec(:)),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation aa energy: ',Ec(1),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation ab energy: ',Ec(2),' au' + write(*,'(A40,1X,F16.10,A3)') ' Correlation bb energy: ',Ec(3),' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,1X,F16.10,A3)') ' Electronic energy: ',Ew,' au' + write(*,'(A40,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' + write(*,'(A40,1X,F16.10,A3)') ' UHF energy: ',Ew + ENuc,' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,F13.6,A3)') ' UHF HOMO a energy:',eps(HOMO(1),1)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF LUMO a energy:',eps(LUMO(1),1)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF HOMOa-LUMOa gap:',Gap(1)*HatoeV,' eV' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,F13.6,A3)') ' UHF HOMO b energy:',eps(HOMO(2),2)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',eps(LUMO(2),2)*HatoeV,' eV' + write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV' + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +! Print results + + write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') 'UHF spin-up orbital coefficients ' + write(*,'(A50)') '-----------------------------------------' + call matout(nBas,nBas,c(:,:,1)) + write(*,'(A50)') '-----------------------------------------' + write(*,'(A50)') 'UHF spin-down orbital coefficients ' + write(*,'(A50)') '-----------------------------------------' + call matout(nBas,nBas,c(:,:,2)) + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A50)') ' UHF spin-up orbital energies ' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,1,eps(:,1)) + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A50)') ' UHF spin-down orbital energies ' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,1,eps(:,2)) + write(*,*) + +end subroutine print_UHF