From 16f7366e363d3c60ced7624fd2ee7a51c2d5a92e Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 26 Oct 2023 19:27:30 +0200 Subject: [PATCH] cleaning GHF code --- src/HF/GHF.f90 | 69 +++++++++++++++++++------------------------- src/QuAcK/GQuAcK.f90 | 33 +++++++++++++-------- src/QuAcK/QuAcK.f90 | 21 +++++--------- src/QuAcK/RQuAcK.f90 | 14 +++++---- src/QuAcK/UQuAcK.f90 | 6 ++-- 5 files changed, 71 insertions(+), 72 deletions(-) diff --git a/src/HF/GHF.f90 b/src/HF/GHF.f90 index 2e73046..8d8f38f 100644 --- a/src/HF/GHF.f90 +++ b/src/HF/GHF.f90 @@ -22,7 +22,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nO(nspin) + integer,intent(in) :: nO double precision,intent(in) :: Ov(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) @@ -36,9 +36,6 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, integer :: nSCF integer :: nBasSq integer :: nBas2Sq - integer :: nOa - integer :: nOb - integer :: nOcc integer :: n_diis double precision :: Conv double precision :: rcond @@ -48,7 +45,6 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, double precision :: Ex,Exaaaa,Exabba,Exbaab,Exbbbb double precision :: dipole(ncart) - double precision,allocatable :: Ca(:,:),Cb(:,:) double precision,allocatable :: Jaa(:,:),Jbb(:,:) double precision,allocatable :: Kaa(:,:),Kab(:,:),Kba(:,:),Kbb(:,:) double precision,allocatable :: Faa(:,:),Fab(:,:),Fba(:,:),Fbb(:,:) @@ -83,21 +79,17 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Useful stuff - nBasSq = nBas*nBas + nBasSq = nBas *nBas nBas2Sq = nBas2*nBas2 - nOa = nO(1) - nOb = nO(2) - nOcc = nOa + nOb - ! Memory allocation - allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),Jaa(nBas,nBas),Jbb(nBas,nBas), & - Kaa(nBas,nBas),Kab(nBas,nBas),Kba(nBas,nBas),Kbb(nBas,nBas), & - Faa(nBas,nBas),Fab(nBas,nBas),Fba(nBas,nBas),Fbb(nBas,nBas), & - Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas), & - F(nBas2,nBas2),Fp(nBas2,nBas2),Cp(nBas2,nBas2), & - H(nBas2,nBas2),S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2), & + allocate(Jaa(nBas,nBas),Jbb(nBas,nBas), & + Kaa(nBas,nBas),Kab(nBas,nBas),Kba(nBas,nBas),Kbb(nBas,nBas), & + Faa(nBas,nBas),Fab(nBas,nBas),Fba(nBas,nBas),Fbb(nBas,nBas), & + Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas), & + F(nBas2,nBas2),Fp(nBas2,nBas2),Cp(nBas2,nBas2),H(nBas2,nBas2), & + S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2), & err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis)) ! Initialization @@ -131,12 +123,14 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, call mo_guess(nBas2,guess_type,S,H,X,C) - P(:,:) = matmul(C(:,1:nOcc),transpose(C(:,1:nOcc))) +! Construct super density matrix - Paa(:,:) = P(1:nBas,1:nBas) - Pab(:,:) = P(1:nBas,nBas+1:nBas2) - Pba(:,:) = P(nBas+1:nBas2,1:nBas) - Pbb(:,:) = P(nBas+1:nBAs2,nBAs+1:nBas2) + P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO))) + + Paa(:,:) = P( 1:nBas , 1:nBas ) + Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) + Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) + Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) !------------------------------------------------------------------------ ! Main SCF loop @@ -189,18 +183,14 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! DIIS extrapolation if(max_diis > 1) then - n_diis = min(n_diis+1,max_diis) call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis(:,1:n_diis),F_diis(:,1:n_diis),err,F) - end if ! Level-shifting if(level_shift > 0d0 .and. Conv > thresh) then - - call level_shifting(level_shift,nBas,nOa+nOb,S,C,F) - + call level_shifting(level_shift,nBas,nO,S,C,F) end if ! Transform Fock matrix in orthogonal basis @@ -216,23 +206,21 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, C(:,:) = matmul(X,Cp) -! Form individual coefficient matrices - - Ca(1:nBas,1:nBas2) = C( 1:nBas ,1:nBas2) - Cb(1:nBas,1:nBas2) = C(nBas+1:nBas2,1:nBas2) - ! Mix guess for UHF solution in singlet states ! if(nSCF == 1) call mix_guess(nBas,nO,mix,c) -! Compute individual density matrices - P(:,:) = matmul(C(:,1:nOcc),transpose(C(:,1:nOcc))) +! Form super density matrix - Paa(:,:) = P(1:nBas,1:nBas) - Pab(:,:) = P(1:nBas,nBas+1:nBas2) - Pba(:,:) = P(nBas+1:nBas2,1:nBas) + P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO))) + +! Compute individual density matrices + + Paa(:,:) = P( 1:nBas , 1:nBas ) + Pab(:,:) = P( 1:nBas ,nBas+1:nBas2) + Pba(:,:) = P(nBas+1:nBas2, 1:nBas ) Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2) - + !------------------------------------------------------------------------ ! Compute UHF energy !------------------------------------------------------------------------ @@ -298,9 +286,12 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, end if -! Compute final UHF energy +! Compute dipole moments + + call dipole_moment(nBas2,P,nNuc,ZNuc,rNuc,dipole_int,dipole) + +! Compute final GHF energy -! call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) ! call print_GHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) end subroutine diff --git a/src/QuAcK/GQuAcK.f90 b/src/QuAcK/GQuAcK.f90 index 3c207f4..82fbe24 100644 --- a/src/QuAcK/GQuAcK.f90 +++ b/src/QuAcK/GQuAcK.f90 @@ -1,5 +1,5 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & - nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & + nNuc,nBas,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, & spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & @@ -16,11 +16,10 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs logical,intent(in) :: doG0W0,doevGW,doqsGW integer,intent(in) :: nNuc,nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: ENuc double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart) @@ -77,6 +76,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs integer :: ixyz integer :: nBas2 + integer :: nS write(*,*) write(*,*) '*******************************' @@ -131,17 +131,24 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs ! 4-index transform allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),ERI_tmp(nBas2,nBas2,nBas2,nBas2)) + Ca(:,:) = cHF(1:nBas,1:nBas2) Cb(:,:) = cHF(nBas+1:nBas2,1:nBas2) + call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Ca,Ca,Ca,ERI_AO,ERI_tmp) ERI_MO(:,:,:,:) = ERI_tmp(:,:,:,:) + call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Cb,Ca,Cb,ERI_AO,ERI_tmp) ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:) + call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Ca,Cb,Ca,ERI_AO,ERI_tmp) ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:) + call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Cb,Cb,Cb,ERI_AO,ERI_tmp) ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:) + deallocate(Ca,Cb,ERI_tmp) + call wall_time(end_AOtoMO) t_AOtoMO = end_AOtoMO - start_AOtoMO @@ -152,10 +159,12 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs ! Stability analysis of HF solution ! !-----------------------------------! + nS = (nO - nC)*(nV - nR) + if(dostab) then call wall_time(start_stab) - call GHF_stability(nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV),epsHF,ERI_MO) + call GHF_stability(nBas2,nC,nO,nV,nR,nS,epsHF,ERI_MO) call wall_time(end_stab) t_stab = end_stab - start_stab @@ -173,7 +182,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs if(doMP) then call wall_time(start_MP) - call GMP(doMP2,reg_MP,nBas2,sum(nC),sum(nO),sum(nV),sum(nR),ERI_MO,ENuc,EHF,epsHF) + call GMP(doMP2,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) call wall_time(end_MP) t_MP = end_MP - start_MP @@ -192,7 +201,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs call wall_time(start_RPA) call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel, & - nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV),ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S) + nBas2,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -231,9 +240,9 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs if(doGW) then call wall_time(start_GW) - call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, & - exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, & - lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV), & + call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, & + exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, & + lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS, & EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) call wall_time(end_GW) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 210a429..18b7005 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -20,7 +20,6 @@ program QuAcK integer :: nO(nspin) integer :: nV(nspin) integer :: nR(nspin) - integer :: nS(nspin) double precision :: ENuc double precision,allocatable :: ZNuc(:),rNuc(:,:) @@ -126,9 +125,6 @@ program QuAcK ! nV = number of virtual orbitals (see below) ! ! nR = number of Rydberg orbitals ! ! nBas = number of basis functions (see below) ! -! = nO + nV ! -! nS = number of single excitation ! -! = nO*nV ! !------------------------------------------------! call read_molecule(nNuc,nO,nC,nR) @@ -142,8 +138,7 @@ program QuAcK ! Read basis set information from PySCF ! !---------------------------------------! - call read_basis_pyscf (nBas,nO,nV) - nS(:) = (nO(:) - nC(:))*(nV(:) - nR(:)) + call read_basis_pyscf(nBas,nO,nV) !--------------------------------------! ! Read one- and two-electron integrals ! @@ -193,7 +188,7 @@ program QuAcK call RQuAcK(doRHF,doROHF,dostab,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,doufG0W0,doufGW,doSRGqsGW, & - doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc, & + doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,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, & @@ -208,7 +203,7 @@ program QuAcK call UQuAcK(doUHF,dostab,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,doufG0W0,doufGW,doSRGqsGW, & - doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc, & + doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,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, & @@ -220,11 +215,11 @@ program QuAcK !--------------------------! if(doGQuAcK) & - call GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & - nNuc,nBas,nC,nO,nV,nR,nS,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, & - spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & - maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & + call GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & + nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(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, & + spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & + maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS) !--------------! diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 5c65a06..a50dab6 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -1,7 +1,7 @@ subroutine RQuAcK(doRHF,doROHF,dostab,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,doufG0W0,doufGW,doSRGqsGW, & - doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc, & + doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,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, & @@ -25,11 +25,10 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh integer,intent(in) :: nNuc,nBas - integer,intent(in) :: nC(nspin) - integer,intent(in) :: nO(nspin) - integer,intent(in) :: nV(nspin) - integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR double precision,intent(in) :: ENuc double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart) @@ -93,6 +92,7 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: ERI_MO(:,:,:,:) integer :: ixyz + integer :: nS write(*,*) write(*,*) '******************************' @@ -168,6 +168,8 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC ! Stability analysis of HF solution ! !-----------------------------------! + nS = (nO - nC)*(nV - nR) + if(dostab) then call wall_time(start_stab) diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index e837367..5eea1e8 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -1,7 +1,7 @@ subroutine UQuAcK(doUHF,dostab,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,doufG0W0,doufGW,doSRGqsGW, & - doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc, & + doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,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, & @@ -27,7 +27,6 @@ subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, integer,intent(in) :: nO(nspin) integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) - integer,intent(in) :: nS(nspin) double precision,intent(in) :: ENuc double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart) @@ -91,6 +90,7 @@ subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, double precision,allocatable :: dipole_int_aa(:,:,:),dipole_int_bb(:,:,:) double precision,allocatable :: ERI_aaaa(:,:,:,:),ERI_aabb(:,:,:,:),ERI_bbbb(:,:,:,:) integer :: ixyz + integer :: nS(nspin) write(*,*) write(*,*) '********************************' @@ -186,6 +186,8 @@ subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, ! Stability analysis of HF solution ! !-----------------------------------! + nS(:) = (nO(:) - nC(:))*(nV(:) - nR(:)) + if(dostab) then call wall_time(start_stab)