diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index 60fd9d2..5659a04 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -123,7 +123,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Guess coefficients and density matrices - call mo_guess(nBas2,guess_type,S,H,X,C) + call mo_guess(nBas2,nBas2,guess_type,S,H,X,C) ! Construct super density matrix @@ -227,7 +227,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Level-shifting - if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,C,F) + if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nBas,nO,S,C,F) ! Transform Fock matrix in orthogonal basis diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 8b1e441..5756744 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,5 +1,8 @@ -subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI,dipole_int,X,ERHF,eHF,c,P) + +! --- + +subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas_AOs, nBas_MOs, nO, S, T, V, Hc, ERI, dipole_int, X, ERHF, eHF, c, P) ! Perform restricted Hartree-Fock calculation @@ -16,24 +19,24 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision,intent(in) :: thresh double precision,intent(in) :: level_shift - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nO integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - 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 :: ET double precision :: EV @@ -56,9 +59,9 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Output variables double precision,intent(out) :: ERHF - double precision,intent(out) :: eHF(nBas) - double precision,intent(inout):: c(nBas,nBas) - double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: eHF(nBas_MOs) + double precision,intent(inout):: c(nBas_AOs,nBas_MOs) + double precision,intent(out) :: P(nBas_AOs,nBas_AOs) ! Hello world @@ -70,21 +73,32 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Useful quantities - nBasSq = nBas*nBas + nBas_AOs_Sq = nBas_AOs*nBas_AOs ! Memory allocation - allocate(J(nBas,nBas),K(nBas,nBas),err(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), & - Fp(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(J(nBas_AOs,nBas_AOs)) + allocate(K(nBas_AOs,nBas_AOs)) + + allocate(err(nBas_AOs,nBas_AOs)) + allocate(F(nBas_AOs,nBas_AOs)) + + allocate(cp(nBas_MOs,nBas_MOs)) + allocate(Fp(nBas_MOs,nBas_MOs)) + + allocate(err_diis(nBas_AOs_Sq,max_diis)) + allocate(F_diis(nBas_AOs_Sq,max_diis)) ! Guess coefficients and density matrix - call mo_guess(nBas,guess_type,S,Hc,X,c) - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + call mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c) + + !P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) + call dgemm('N', 'T', nBas_AOs, nBas_AOs, nO, 2.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P, nBas_AOs) ! Initialization - n_diis = 0 + n_diis = 0 F_diis(:,:) = 0d0 err_diis(:,:) = 0d0 rcond = 0d0 @@ -110,31 +124,31 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Build Fock matrix - call Hartree_matrix_AO_basis(nBas,P,ERI,J) - call exchange_matrix_AO_basis(nBas,P,ERI,K) + call Hartree_matrix_AO_basis(nBas_AOs, P, ERI, J) + call exchange_matrix_AO_basis(nBas_AOs, P, ERI, K) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) ! 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)) ! Kinetic energy - ET = trace_matrix(nBas,matmul(P,T)) + ET = trace_matrix(nBas_AOs, matmul(P, T)) ! Potential energy - EV = trace_matrix(nBas,matmul(P,V)) + EV = trace_matrix(nBas_AOs, matmul(P, V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + EJ = 0.5d0*trace_matrix(nBas_AOs, matmul(P, J)) ! Exchange energy - EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) + EK = 0.25d0*trace_matrix(nBas_AOs, matmul(P, K)) ! Total energy @@ -144,25 +158,28 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN 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,F) + n_diis = min(n_diis+1, max_diis) + call DIIS_extrapolation(rcond, nBas_AOs_Sq, nBas_AOs_Sq, n_diis, err_diis, F_diis, err, F) end if ! Level shift - if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F) + if(level_shift > 0d0 .and. Conv > thresh) then + call level_shifting(level_shift, nBas_AOs, nBas_MOs, nO, S, c, F) + endif ! Diagonalize Fock matrix - Fp = matmul(transpose(X),matmul(F,X)) + Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eHF) - c = matmul(X,cp) + call diagonalize_matrix(nBas_MOs, cp, eHF) + c = matmul(X, cp) ! Density matrix - 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_AOs, nBas_AOs, nO, 2.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P, nBas_AOs) ! Dump results @@ -185,14 +202,16 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(J, K, err, F, cp, Fp, err_diis, F_diis) + stop end if ! Compute dipole moments - call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_RHF(nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole) + call dipole_moment(nBas_AOs, P, nNuc, ZNuc, rNuc, dipole_int, dipole) + call print_RHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, EK, ERHF, dipole) ! Testing zone @@ -205,4 +224,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN end if + deallocate(J, K, err, F, cp, Fp, err_diis, F_diis) + end subroutine diff --git a/src/HF/RMOM.f90 b/src/HF/RMOM.f90 index a49746d..fc43b76 100644 --- a/src/HF/RMOM.f90 +++ b/src/HF/RMOM.f90 @@ -190,6 +190,6 @@ subroutine RMOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) EK = 0.5d0*trace_matrix(nBas,matmul(P,K)) ERHF = ET + EV + EJ + EK - call print_RHF(nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF) + call print_RHF(nBas,nBas,nO,e,c,ENuc,ET,EV,EJ,EK,ERHF) end subroutine diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index 7d86359..0cd1f5f 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -86,7 +86,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Guess coefficients and demsity matrices - call mo_guess(nBas,guess_type,S,Hc,X,c) + 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)))) end do @@ -185,7 +185,7 @@ 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,maxval(nO),S,c,Ftot) + call level_shifting(level_shift,nBas,nBas,maxval(nO),S,c,Ftot) end do end if diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 919cb0f..faa7de1 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -85,7 +85,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Guess coefficients and demsity matrices do ispin=1,nspin - call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin)) + call mo_guess(nBas,nBas,guess_type,S,Hc,X,c(:,:,ispin)) P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) end do @@ -186,7 +186,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu if(level_shift > 0d0 .and. Conv > thresh) then do ispin=1,nspin - call level_shifting(level_shift,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin)) + call level_shifting(level_shift,nBas,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin)) end do end if diff --git a/src/HF/cRHF.f90 b/src/HF/cRHF.f90 index 9640a48..36cfa41 100644 --- a/src/HF/cRHF.f90 +++ b/src/HF/cRHF.f90 @@ -94,7 +94,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! Guess coefficients and density matrix - call mo_guess(nBas,guess_type,S,Hc,X,c) + call mo_guess(nBas,nBas,guess_type,S,Hc,X,c) P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) ! Initialization @@ -166,7 +166,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! Level shift - if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F) + if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nBas,nO,S,c,F) ! Diagonalize Fock matrix @@ -207,7 +207,7 @@ subroutine cRHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,r ! Compute dipole moments call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_RHF(nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole) + call print_RHF(nBas,nBas,nO,eHF,C,ENuc,ET,EV,EJ,EK,ERHF,dipole) ! Testing zone diff --git a/src/HF/core_guess.f90 b/src/HF/core_guess.f90 index 7e45e8e..c48c3ae 100644 --- a/src/HF/core_guess.f90 +++ b/src/HF/core_guess.f90 @@ -1,4 +1,4 @@ -subroutine core_guess(nBas,Hc,X,c) +subroutine core_guess(nBas_AOs, nBas_MOs, Hc, X, c) ! Core guess of the molecular orbitals for HF calculation @@ -6,9 +6,9 @@ subroutine core_guess(nBas,Hc,X,c) ! Input variables - integer,intent(in) :: nBas - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + integer,intent(in) :: nBas_AOs, nBas_MOs + double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(in) :: X(nBas_AOs,nBas_MOs) ! Local variables @@ -18,16 +18,19 @@ subroutine core_guess(nBas,Hc,X,c) ! Output variables - double precision,intent(out) :: c(nBas,nBas) + double precision,intent(out) :: c(nBas_AOs,nBas_MOs) ! Memory allocation - allocate(cp(nBas,nBas),e(nBas)) + allocate(cp(nBas_MOs,nBas_MOs), e(nBas_MOs)) ! Core guess - cp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) - call diagonalize_matrix(nBas,cp,e) - c(:,:) = matmul(X(:,:),cp(:,:)) + cp(:,:) = matmul(transpose(X(:,:)), matmul(Hc(:,:), X(:,:))) + + call diagonalize_matrix(nBas_MOs, cp, e) + c(:,:) = matmul(X(:,:), cp(:,:)) + + deallocate(cp, e) end subroutine diff --git a/src/HF/huckel_guess.f90 b/src/HF/huckel_guess.f90 index 658c853..a8e7e52 100644 --- a/src/HF/huckel_guess.f90 +++ b/src/HF/huckel_guess.f90 @@ -1,4 +1,4 @@ -subroutine huckel_guess(nBas,S,Hc,X,c) +subroutine huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c) ! Hickel guess @@ -6,10 +6,10 @@ subroutine huckel_guess(nBas,S,Hc,X,c) ! Input variables - integer,intent(in) :: nBas - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + integer,intent(in) :: nBas_AOs, nBas_MOs + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(in) :: X(nBas_AOs,nBas_MOs) ! Local variables @@ -20,11 +20,11 @@ subroutine huckel_guess(nBas,S,Hc,X,c) ! Output variables - double precision,intent(out) :: c(nBas,nBas) + double precision,intent(out) :: c(nBas_AOs,nBas_MOs) ! Memory allocation - allocate(F(nBas,nBas)) + allocate(F(nBas_AOs,nBas_AOs)) ! Extended Huckel parameter @@ -32,9 +32,9 @@ subroutine huckel_guess(nBas,S,Hc,X,c) ! GWH approximation - do mu=1,nBas + do mu = 1, nBas_AOs F(mu,mu) = Hc(mu,mu) - do nu=mu+1,nBas + do nu = mu+1, nBas_AOs F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu)) F(nu,mu) = F(mu,nu) @@ -42,6 +42,8 @@ subroutine huckel_guess(nBas,S,Hc,X,c) end do end do - call core_guess(nBas,F,X,c) + call core_guess(nBas_AOs, nBas_MOs, F, X, c) + + deallocate(F) end subroutine diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90 index fef8d87..fd1f20d 100644 --- a/src/HF/mo_guess.f90 +++ b/src/HF/mo_guess.f90 @@ -1,4 +1,7 @@ -subroutine mo_guess(nBas,guess_type,S,Hc,X,c) + +! --- + +subroutine mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c) ! Guess of the molecular orbitals for HF calculation @@ -6,15 +9,15 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: guess_type - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(in) :: X(nBas_AOs,nBas_MOs) ! Output variables - double precision,intent(inout) :: c(nBas,nBas) + double precision,intent(inout) :: c(nBas_AOs,nBas_MOs) if(guess_type == 0) then @@ -24,12 +27,12 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c) elseif(guess_type == 1) then write(*,*) 'Core guess...' - call core_guess(nBas,Hc,X,c) + call core_guess(nBas_AOs, nBas_MOs, Hc, X, c) elseif(guess_type == 2) then write(*,*) 'Huckel guess...' - call huckel_guess(nBas,S,Hc,X,c) + call huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c) elseif(guess_type == 3) then diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 3a06d31..6790fcf 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -1,4 +1,7 @@ -subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) + +! --- + +subroutine print_RHF(nBas_AOs, nBas_MOs, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipole) ! Print one-electron energies and other stuff for G0W0 @@ -7,10 +10,10 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nO - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: eHF(nBas_MOs) + double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) double precision,intent(in) :: ENuc double precision,intent(in) :: ET double precision,intent(in) :: EV @@ -75,13 +78,13 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' RHF orbital coefficients ' write(*,'(A50)') '---------------------------------------' - call matout(nBas,nBas,cHF) + call matout(nBas_AOs, nBas_MOs, cHF) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' RHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eHF) + call vecout(nBas_MOs, eHF) write(*,*) end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index d1aa007..75514f4 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -15,7 +15,7 @@ program QuAcK logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh - integer :: nNuc,nBas_AOs + integer :: nNuc, nBas_AOs, nBas_MOs integer :: nC(nspin) integer :: nO(nspin) integer :: nV(nspin) @@ -31,6 +31,7 @@ program QuAcK double precision,allocatable :: X(:,:) double precision,allocatable :: dipole_int_AO(:,:,:) double precision,allocatable :: ERI_AO(:,:,:,:) + double precision,allocatable :: Uvec(:,:), Uval(:) double precision :: start_QuAcK,end_QuAcK,t_QuAcK double precision :: start_int ,end_int ,t_int @@ -68,6 +69,8 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest + integer :: i, j + !-------------! ! Hello World ! !-------------! @@ -128,6 +131,7 @@ program QuAcK ! nV = number of virtual orbitals (see below) ! ! nR = number of Rydberg orbitals ! ! nBas_AOs = number of basis functions in AOs ! +! nBas_MOs = number of basis functions in MOs ! !---------------------------------------------------! call read_molecule(nNuc,nO,nC,nR) @@ -153,7 +157,6 @@ program QuAcK allocate(T(nBas_AOs,nBas_AOs)) allocate(V(nBas_AOs,nBas_AOs)) allocate(Hc(nBas_AOs,nBas_AOs)) - allocate(X(nBas_AOs,nBas_AOs)) allocate(ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)) allocate(dipole_int_AO(nBas_AOs,nBas_AOs,ncart)) @@ -173,7 +176,39 @@ program QuAcK ! Compute orthogonalization matrix - call orthogonalization_matrix(nBas_AOs, S, X) + !call orthogonalization_matrix(nBas_AOs, S, X) + + allocate(Uvec(nBas_AOs,nBas_AOs), Uval(nBas_AOs)) + + Uvec(1:nBas_AOs,1:nBas_AOs) = S(1:nBas_AOs,1:nBas_AOs) + call diagonalize_matrix(nBas_AOs, Uvec, Uval) + + nBas_MOs = 0 + do i = 1, nBas_AOs + if(Uval(i) > 1d-6) then + Uval(i) = 1d0 / dsqrt(Uval(i)) + nBas_MOs = nBas_MOs + 1 + else + write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' + end if + end do + + write(*,'(A38)') '--------------------------------------' + write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs + write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nBas_MOs + write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nBas_MOs)/dble(nBas_AOs)) + write(*,'(A38)') '--------------------------------------' + write(*,*) + + allocate(X(nBas_AOs,nBas_MOs)) + do j = 1, nBas_MOs + do i = 1, nBas_AOs + X(i,j) = Uvec(i,j) * Uval(j) + enddo + enddo + + deallocate(Uvec, Uval) + !---------------------! ! Choose QuAcK branch ! @@ -205,7 +240,7 @@ program QuAcK dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - nNuc,nBas_AOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + nNuc,nBas_AOs,nBas_MOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 3897f37..1c9545b 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -2,7 +2,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + nNuc,nBas_AOs,nBas_MOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & @@ -29,7 +29,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp logical,intent(in) :: doG0T0eh,doevGTeh,doqsGTeh - integer,intent(in) :: nNuc,nBas + integer,intent(in) :: nNuc,nBas_AOs,nBas_MOs integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -38,13 +38,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart) - 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) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + 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) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart) + double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) integer,intent(in) :: maxSCF_HF,max_diis_HF double precision,intent(in) :: thresh_HF,level_shift,mix @@ -109,8 +109,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! Memory allocation ! !-------------------! - allocate(cHF(nBas,nBas),eHF(nBas),PHF(nBas,nBas), & - dipole_int_MO(nBas,nBas,ncart),ERI_MO(nBas,nBas,nBas,nBas)) + 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(ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs)) !---------------------! ! Hartree-Fock module ! @@ -119,12 +122,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doRHF) then call wall_time(start_HF) - call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF) + ! 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) t_HF = end_HF - start_HF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds' write(*,*) end if @@ -132,12 +136,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doROHF) then 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,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF) + nBas_AOs,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 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ROHF = ',t_HF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for ROHF = ',t_HF,' seconds' write(*,*) end if @@ -154,18 +159,20 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! Read and transform dipole-related integrals - do ixyz=1,ncart - call AOtoMO(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + do ixyz = 1, ncart + ! TODO + call AOtoMO(nBas_AOs,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do ! 4-index transform - call AOtoMO_ERI_RHF(nBas,cHF,ERI_AO,ERI_MO) + ! TODO + call AOtoMO_ERI_RHF(nBas_AOs,cHF,ERI_AO,ERI_MO) call wall_time(end_AOtoMO) t_AOtoMO = end_AOtoMO - start_AOtoMO - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds' write(*,*) !-----------------------------------! @@ -177,11 +184,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(dostab) then call wall_time(start_stab) - call RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO) + call RHF_stability(nBas_AOs,nC,nO,nV,nR,nS,eHF,ERI_MO) call wall_time(end_stab) t_stab = end_stab - start_stab - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds' write(*,*) end if @@ -189,12 +196,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(dosearch) then call wall_time(start_stab) + ! TODO call RHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X,ERHF,eHF,cHF,PHF) + nBas_AOs,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X,ERHF,eHF,cHF,PHF) call wall_time(end_stab) t_stab = end_stab - start_stab - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for stability analysis = ',t_stab,' seconds' write(*,*) end if @@ -208,11 +216,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doMP) then call wall_time(start_MP) - call RMP(dotest,doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + ! TODO + call RMP(dotest,doMP2,doMP3,reg_MP,nBas_AOs,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_MP) t_MP = end_MP - start_MP - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP = ',t_MP,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP = ',t_MP,' seconds' write(*,*) end if @@ -227,12 +236,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doCC) then call wall_time(start_CC) + ! TODO call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,Hc,ERI_MO,ENuc,ERHF,eHF,cHF) + maxSCF_CC,thresh_CC,max_diis_CC,nBas_AOs,nC,nO,nV,nR,Hc,ERI_MO,ENuc,ERHF,eHF,cHF) call wall_time(end_CC) t_CC = end_CC - start_CC - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CC = ',t_CC,' seconds' write(*,*) end if @@ -246,12 +256,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doCI) then call wall_time(start_CI) - call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & + ! TODO + call RCI(dotest,doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas_AOs,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & eHF,ERHF,cHF,S) call wall_time(end_CI) t_CI = end_CI - start_CI - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CI = ',t_CI,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for CI = ',t_CI,' seconds' write(*,*) end if @@ -265,12 +276,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doRPA) then call wall_time(start_RPA) + ! TODO call RRPA(dotest,dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S) + nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF,cHF,S) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RPA = ',t_RPA,' seconds' write(*,*) end if @@ -284,14 +296,15 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGF) then call wall_time(start_GF) + ! TODO call RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF, & - nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & + nNuc,ZNuc,rNuc,ENuc,nBas_AOs,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GF) t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GF2 = ',t_GF,' seconds' write(*,*) end if @@ -305,14 +318,15 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGW) then call wall_time(start_GW) + ! TODO call RGW(dotest,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW, & doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & - lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & + lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas_AOs,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GW = ',t_GW,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GW = ',t_GW,' seconds' write(*,*) end if @@ -326,14 +340,15 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doGT) then call wall_time(start_GT) + ! TODO call RGT(dotest,doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & maxSCF_GT,thresh_GT,max_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE, & TDA_T,TDA,dBSE,dTDA,singlet,triplet,lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + nBas_AOs,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GT = ',t_GT,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for GT = ',t_GT,' seconds' write(*,*) end if diff --git a/src/utils/level_shifting.f90 b/src/utils/level_shifting.f90 index 8435b32..a006622 100644 --- a/src/utils/level_shifting.f90 +++ b/src/utils/level_shifting.f90 @@ -1,4 +1,4 @@ -subroutine level_shifting(level_shift,nBas,nO,S,c,F) +subroutine level_shifting(level_shift, nBas_AOs, nBas_MOs, nO, S, c, F) ! Perform level-shifting on the Fock matrix @@ -7,10 +7,10 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F) ! Input variables double precision,intent(in) :: level_shift - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nO - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: c(nBas_AOs,nBas_MOs) ! Local variables @@ -21,17 +21,19 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F) ! Output variables - double precision,intent(inout):: F(nBas,nBas) + double precision,intent(inout):: F(nBas_AOs,nBas_AOs) - allocate(F_MO(nBas,nBas),Sc(nBas,nBas)) + allocate(F_MO(nBas_MOs,nBas_MOs), Sc(nBas_AOs,nBas_MOs)) - F_MO(:,:) = matmul(transpose(c),matmul(F,c)) + F_MO(:,:) = matmul(transpose(c), matmul(F, c)) - do a=nO+1,nBas + do a = nO+1, nBas_MOs F_MO(a,a) = F_MO(a,a) + level_shift end do - Sc(:,:) = matmul(S,c) - F(:,:) = matmul(Sc,matmul(F_MO,transpose(Sc))) + Sc(:,:) = matmul(S, c) + F(:,:) = matmul(Sc, matmul(F_MO, transpose(Sc))) + + deallocate(F_MO, Sc) end subroutine diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90 index d7e0089..af781e1 100644 --- a/src/utils/orthogonalization_matrix.f90 +++ b/src/utils/orthogonalization_matrix.f90 @@ -1,4 +1,7 @@ -subroutine orthogonalization_matrix(nBas,S,X) + +! --- + +subroutine orthogonalization_matrix(nBas, S, X) ! Compute the orthogonalization matrix X @@ -35,14 +38,32 @@ subroutine orthogonalization_matrix(nBas,S,X) if(ortho_type == 1) then + ! + ! S V = V s where + ! + ! V.T V = 1 and s > 0 (S is positive def) + ! + ! S = V s V.T + ! = V s^0.5 s^0.5 V.T + ! = V s^0.5 V.T V s^0.5 V.T + ! = S^0.5 S^0.5 + ! + ! where + ! + ! S^0.5 = V s^0.5 V.T + ! + ! X = S^(-0.5) + ! = V s^(-0.5) V.T + ! + ! write(*,*) ! write(*,*) ' Lowdin orthogonalization' ! write(*,*) Uvec = S - call diagonalize_matrix(nBas,Uvec,Uval) + call diagonalize_matrix(nBas, Uvec, Uval) - do i=1,nBas + do i = 1, nBas if(Uval(i) < thresh) then @@ -50,7 +71,7 @@ subroutine orthogonalization_matrix(nBas,S,X) end if - Uval(i) = 1d0/sqrt(Uval(i)) + Uval(i) = 1d0 / dsqrt(Uval(i)) end do @@ -63,13 +84,13 @@ subroutine orthogonalization_matrix(nBas,S,X) ! write(*,*) Uvec = S - call diagonalize_matrix(nBas,Uvec,Uval) + call diagonalize_matrix(nBas, Uvec, Uval) - do i=1,nBas + do i = 1, nBas if(Uval(i) > thresh) then - Uval(i) = 1d0/sqrt(Uval(i)) + Uval(i) = 1d0 / dsqrt(Uval(i)) else @@ -79,7 +100,7 @@ subroutine orthogonalization_matrix(nBas,S,X) end do - call AD(nBas,Uvec,Uval) + call AD(nBas, Uvec, Uval) X = Uvec elseif(ortho_type == 3) then @@ -117,3 +138,6 @@ subroutine orthogonalization_matrix(nBas,S,X) end if end subroutine + +! --- + diff --git a/src/utils/read_basis_pyscf.f90 b/src/utils/read_basis_pyscf.f90 index 027fcac..42dfde2 100644 --- a/src/utils/read_basis_pyscf.f90 +++ b/src/utils/read_basis_pyscf.f90 @@ -24,10 +24,10 @@ subroutine read_basis_pyscf(nBas_AOs, nO, nV) read(3, *) nBas_AOs close(unit=3) - write(*,'(A38)') '--------------------------------------' - write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs - write(*,'(A38)') '--------------------------------------' - write(*,*) +! write(*,'(A38)') '--------------------------------------' +! write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs +! write(*,'(A38)') '--------------------------------------' +! write(*,*) ! Number of virtual orbitals