diff --git a/input/methods b/input/methods index 2452cdc..738a281 100644 --- a/input/methods +++ b/input/methods @@ -13,7 +13,7 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW* - T F T + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index eef2cae..5f8d15d 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ -# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess - 128 0.0000001 T 5 1 1 F +# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability + 128 0.0000001 T 5 1 1 T T # MP: # CC: maxSCF thresh DIIS n_diis @@ -9,7 +9,7 @@ # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 - 256 0.0000001 T 5 T 0.00367493 F F F F F + 256 0.00001 T 5 T 0.0 F F F F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn diff --git a/mol/h2.xyz b/mol/h2.xyz index ce5f117..a302682 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0.0 0.0 0.0 -H 0.0 0.0 0.7408481486 +H 0.0 0.0 1.4 diff --git a/src/HF/RHF_stability.f90 b/src/HF/RHF_stability.f90 new file mode 100644 index 0000000..bade757 --- /dev/null +++ b/src/HF/RHF_stability.f90 @@ -0,0 +1,147 @@ +subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) + +! Perform a stability analysis of the RHF solution + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer,parameter :: maxS = 20 + integer :: ia + integer :: ispin + + double precision,allocatable :: A(:,:) + double precision,allocatable :: B(:,:) + double precision,allocatable :: AB(:,:) + double precision,allocatable :: Om(:) + +! Memory allocation + + allocate(A(nS,nS),B(nS,nS),AB(nS,nS),Om(nS)) + +!-------------------------------------------------------------! +! Stability analysis: Real RHF -> Real RHF +!-------------------------------------------------------------! + + ispin = 1 + + call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A) + call linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B) + + AB(:,:) = A(:,:) + B(:,:) + + call diagonalize_matrix(nS,AB,Om) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real RHF -> Real RHF |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + do ia=1,min(nS,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' + enddo + write(*,*)'-------------------------------------------------------------' + + if(minval(Om(:)) < 0d0) then + + write(*,'(1X,A40,1X)') 'Too bad, RHF solution is unstable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au' + + else + + write(*,'(1X,A40,1X)') 'Well done, RHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au' + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + +!-------------------------------------------------------------! +! Stability analysis: Real RHF -> Complex RHF +!-------------------------------------------------------------! + + AB(:,:) = A(:,:) - B(:,:) + + call diagonalize_matrix(nS,AB,Om) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real RHF -> Complex RHF |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + do ia=1,min(nS,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' + enddo + write(*,*)'-------------------------------------------------------------' + + if(minval(Om(:)) < 0d0) then + + write(*,'(1X,A40,1X)') 'Too bad, RHF solution is unstable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au' + + else + + write(*,'(1X,A40,1X)') 'Well done, RHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au' + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + +!-------------------------------------------------------------! +! Stability analysis: Real RHF -> Real UHF +!-------------------------------------------------------------! + + ispin = 2 + + call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A) + call linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B) + + AB(:,:) = A(:,:) + B(:,:) + + call diagonalize_matrix(nS,AB,Om) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real RHF -> Real UHF |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + do ia=1,min(nS,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Om(ia),'|',Om(ia)*HaToeV,'|' + enddo + write(*,*)'-------------------------------------------------------------' + + if(minval(Om(:)) < 0d0) then + + write(*,'(1X,A40,1X)') 'Too bad, RHF solution is unstable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om(1),' au' + + else + + write(*,'(1X,A40,1X)') 'Well done, RHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om(1),' au' + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + +end subroutine RHF_stability diff --git a/src/HF/UHF_stability.f90 b/src/HF/UHF_stability.f90 new file mode 100644 index 0000000..db77387 --- /dev/null +++ b/src/HF/UHF_stability.f90 @@ -0,0 +1,174 @@ +subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) + +! Perform a stability analysis of the UHF solution + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: 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) + + double precision,intent(in) :: eHF(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + +! Local variables + + integer,parameter :: maxS = 20 + integer :: ia + integer :: ispin + + integer :: nS_aa,nS_bb,nS_sc + double precision,allocatable :: Om_sc(:) + double precision,allocatable :: A_sc(:,:) + double precision,allocatable :: B_sc(:,:) + double precision,allocatable :: AB_sc(:,:) + + integer :: nS_ab,nS_ba,nS_sf + double precision,allocatable :: Om_sf(:) + double precision,allocatable :: A_sf(:,:) + double precision,allocatable :: B_sf(:,:) + double precision,allocatable :: AB_sf(:,:) + +! Menory allocation + + nS_aa = nS(1) + nS_bb = nS(2) + nS_sc = nS_aa + nS_bb + + allocate(Om_sc(nS_sc),A_sc(nS_sc,nS_sc),B_sc(nS_sc,nS_sc),AB_sc(nS_sc,nS_sc)) + +!-------------------------------------------------------------! +! Stability analysis: Real UHF -> Real UHF +!-------------------------------------------------------------! + + ispin = 1 + + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0,eHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,A_sc) + call unrestricted_linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, & + ERI_aaaa,ERI_aabb,ERI_bbbb,B_sc) + + AB_sc(:,:) = A_sc(:,:) + B_sc(:,:) + + call diagonalize_matrix(nS_sc,AB_sc,Om_sc) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real UHF -> Real UHF |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + do ia=1,min(nS_sc,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Om_sc(ia),'|',Om_sc(ia)*HaToeV,'|' + enddo + write(*,*)'-------------------------------------------------------------' + + if(minval(Om_sc(:)) < 0d0) then + + write(*,'(1X,A40,1X)') 'Too bad, UHF solution is unstable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om_sc(1),' au' + + else + + write(*,'(1X,A40,1X)') 'Well done, UHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om_sc(1),' au' + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + +!-------------------------------------------------------------! +! Stability analysis: Real UHF -> Complex UHF +!-------------------------------------------------------------! + + AB_sc(:,:) = A_sc(:,:) - B_sc(:,:) + + call diagonalize_matrix(nS_sc,AB_sc,Om_sc) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real UHF -> Complex UHF |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + do ia=1,min(nS_sc,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Om_sc(ia),'|',Om_sc(ia)*HaToeV,'|' + enddo + write(*,*)'-------------------------------------------------------------' + + if(minval(Om_sc(:)) < 0d0) then + + write(*,'(1X,A40,1X)') 'Too bad, UHF solution is unstable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om_sc(1),' au' + + else + + write(*,'(1X,A40,1X)') 'Well done, UHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om_sc(1),' au' + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + +! Menory (de)allocation + + + nS_ab = (nO(1) - nC(1))*(nV(2) - nR(2)) + nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) + nS_sf = nS_ab + nS_ba + + deallocate(Om_sc,A_sc,B_sc,AB_sc) + allocate(Om_sf(nS_sf),A_sf(nS_sf,nS_sf),B_sf(nS_sf,nS_sf),AB_sf(nS_sf,nS_sf)) + +!-------------------------------------------------------------! +! Stability analysis: Real UHF -> Real GHF +!-------------------------------------------------------------! + + ispin = 2 + + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,1d0,eHF, & + ERI_aaaa,ERI_aabb,ERI_bbbb,A_sf) + call unrestricted_linear_response_B_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,1d0, & + ERI_aaaa,ERI_aabb,ERI_bbbb,B_sf) + + AB_sf(:,:) = A_sf(:,:) + B_sf(:,:) + + call diagonalize_matrix(nS_sf,AB_sf,Om_sf) + + write(*,*)'-------------------------------------------------------------' + write(*,*)'| Stability analysis: Real UHF -> Real GHF |' + write(*,*)'-------------------------------------------------------------' + write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & + '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' + write(*,*)'-------------------------------------------------------------' + do ia=1,min(nS_sf,maxS) + write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & + '|',ia,'|',Om_sf(ia),'|',Om_sf(ia)*HaToeV,'|' + enddo + write(*,*)'-------------------------------------------------------------' + + if(minval(Om_sf(:)) < 0d0) then + + write(*,'(1X,A40,1X)') 'Too bad, UHF solution is unstable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Largest negative eigenvalue: ',Om_sf(1),' au' + + else + + write(*,'(1X,A40,1X)') 'Well done, UHF solution is stable!' + write(*,'(1X,A40,1X,F15.10,A3)') 'Smallest eigenvalue: ',Om_sf(1),' au' + + end if + write(*,*)'-------------------------------------------------------------' + write(*,*) + +end subroutine UHF_stability diff --git a/src/MBPT/print_qsUGW.f90 b/src/MBPT/print_qsUGW.f90 index 144b2d9..32dc633 100644 --- a/src/MBPT/print_qsUGW.f90 +++ b/src/MBPT/print_qsUGW.f90 @@ -113,6 +113,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, & ! Dump results for final iteration if(Conv < thresh) then +! if(.true.) then write(*,*) write(*,'(A60)') '-------------------------------------------------' diff --git a/src/MP/MP2.f90 b/src/MP/MP2.f90 index 2c6d75c..8b028df 100644 --- a/src/MP/MP2.f90 +++ b/src/MP/MP2.f90 @@ -6,9 +6,15 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) ! Input variables - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: ENuc,EHF - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas),e(nBas) + integer,intent(in) :: nBas + 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) :: EHF + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index ad67e96..177f7e4 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -6,6 +6,7 @@ program QuAcK logical :: doSph logical :: unrestricted = .false. logical :: doRHF,doUHF,doMOM + logical :: dostab logical :: doKS logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doDCD,doCCSD,doCCSDT @@ -78,6 +79,7 @@ program QuAcK double precision :: start_QuAcK ,end_QuAcK ,t_QuAcK double precision :: start_int ,end_int ,t_int double precision :: start_HF ,end_HF ,t_HF + double precision :: start_stab ,end_stab ,t_stab double precision :: start_KS ,end_KS ,t_KS double precision :: start_MOM ,end_MOM ,t_MOM double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO @@ -169,14 +171,14 @@ program QuAcK ! Read options for methods - call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix, & - maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - TDA,singlet,triplet,spin_conserved,spin_flip, & - maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & - COHSEX,SOSEX,TDA_W,G0W,GW0, & - doACFDT,exchange_kernel,doXBS, & - BSE,dBSE,dTDA,evDyn, & + call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, & + maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & + TDA,singlet,triplet,spin_conserved,spin_flip, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & + COHSEX,SOSEX,TDA_W,G0W,GW0, & + doACFDT,exchange_kernel,doXBS, & + BSE,dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Weird stuff @@ -439,6 +441,32 @@ program QuAcK write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' write(*,*) +!------------------------------------------------------------------------ +! Stability analysis of HF solution +!------------------------------------------------------------------------ + + if(dostab) then + + call cpu_time(start_stab) + + if(unrestricted) then + + call UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb) + + else + + call RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_MO) + + end if + + call cpu_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(*,*) + + end if + !------------------------------------------------------------------------ ! Compute MP2 energy !------------------------------------------------------------------------ diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 1492197..2d7aa4f 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,11 +1,11 @@ -subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix, & - maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - TDA,singlet,triplet,spin_conserved,spin_flip, & - maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & - COHSEX,SOSEX,TDA_W,G0W,GW0, & - doACFDT,exchange_kernel,doXBS, & - BSE,dBSE,dTDA,evDyn, & +subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, & + maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & + TDA,singlet,triplet,spin_conserved,spin_flip, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & + COHSEX,SOSEX,TDA_W,G0W,GW0, & + doACFDT,exchange_kernel,doXBS, & + BSE,dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Read desired methods @@ -21,6 +21,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t integer,intent(out) :: guess_type integer,intent(out) :: ortho_type logical,intent(out) :: mix + logical,intent(out) :: dostab integer,intent(out) :: maxSCF_CC double precision,intent(out) :: thresh_CC @@ -87,12 +88,14 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t guess_type = 1 ortho_type = 1 mix = .false. + dostab = .false. read(1,*) - read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2 + read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,answer3 if(answer1 == 'T') DIIS_HF = .true. if(answer2 == 'T') mix = .true. + if(answer3 == 'T') dostab = .true. if(.not.DIIS_HF) n_diis_HF = 1