diff --git a/input/options b/input/options index f61a0ff..c310e86 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.000001 T 5 1 1 F + 128 0.000001 T 5 1 1 T # MP: # CC: maxSCF thresh DIIS n_diis @@ -13,6 +13,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T T T F + F T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/h2.xyz b/mol/h2.xyz index 7ab70eb..dfa658c 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.741 +H 0.0 0.0 2.645875 diff --git a/src/LR/linear_response_B_pp.f90 b/src/LR/linear_response_B_pp.f90 index fe2515d..59e7a3d 100644 --- a/src/LR/linear_response_B_pp.f90 +++ b/src/LR/linear_response_B_pp.f90 @@ -86,4 +86,7 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) end if + print*,'B pp-matrix' + call matout(nVV,nOO,B_pp) + end subroutine linear_response_B_pp diff --git a/src/LR/linear_response_C_pp.f90 b/src/LR/linear_response_C_pp.f90 index 82f00f7..d22b11d 100644 --- a/src/LR/linear_response_C_pp.f90 +++ b/src/LR/linear_response_C_pp.f90 @@ -96,4 +96,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) end if + print*,'C pp-matrix' + call matout(nVV,nVV,C_pp) + end subroutine linear_response_C_pp diff --git a/src/LR/linear_response_D_pp.f90 b/src/LR/linear_response_D_pp.f90 index 4257795..d06a119 100644 --- a/src/LR/linear_response_D_pp.f90 +++ b/src/LR/linear_response_D_pp.f90 @@ -96,4 +96,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) end if + print*,'D pp-matrix' + call matout(nOO,nOO,D_pp) + end subroutine linear_response_D_pp diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index e0a75d6..710c9e4 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -316,7 +316,7 @@ program QuAcK if(doKS) then call cpu_time(start_KS) - call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, & + call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, & nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO) diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 93f7816..ab6b5f0 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -1,4 +1,4 @@ -subroutine eDFT(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, & +subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, & nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int) @@ -14,6 +14,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,n integer,intent(in) :: maxSCF integer,intent(in) :: max_diis integer,intent(in) :: guess_type + logical,intent(in) :: mix double precision,intent(in) :: thresh integer,intent(in) :: nNuc @@ -201,13 +202,13 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,n end if !------------------------------------------------------------------------ -! Compute N-centered UKS energy (UNBROKEN) +! Compute N-centered UKS energy !------------------------------------------------------------------------ if(method == 'eDFT-UKS') then call cpu_time(start_KS) - call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type, & + call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),maxSCF,thresh,max_diis,guess_type,mix, & nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered) call cpu_time(end_KS) diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index fa8b158..f045c3b 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -1,4 +1,4 @@ -subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh,max_diis,guess_type, & +subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, & nBas,AO,dAO,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice,doNcentered) ! Perform unrestricted Kohn-Sham calculation for ensembles @@ -17,6 +17,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: maxSCF,max_diis,guess_type + logical,intent(in) :: mix double precision,intent(in) :: thresh integer,intent(in) :: nBas double precision,intent(in) :: AO(nBas,nGrid) @@ -39,6 +40,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig logical :: LDA_centered = .false. integer :: nSCF,nBasSq integer :: n_diis + integer :: nO(nspin) double precision :: conv double precision :: rcond(nspin) double precision :: ET(nspin) @@ -75,7 +77,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig double precision :: E(nEns) double precision :: Om(nEns) - integer :: ispin,iEns + integer :: ispin,iEns,iBas ! Hello world @@ -136,6 +138,15 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin)) end do + ! Mix guess to enforce symmetry breaking + + nO(:) = 0 + do ispin=1,nspin + nO(ispin) = int(sum(occnum(:,ispin,1))) + end do + + if(mix) call mix_guess(nBas,nO,c) + else if(guess_type == 2) then do ispin=1,nspin @@ -260,7 +271,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig err(:,:,ispin) = matmul(F(:,:,ispin),matmul(Pw(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:,ispin)),F(:,:,ispin)) end do - conv = maxval(abs(err(:,:,:))) + if(nSCF > 1) conv = maxval(abs(err(:,:,:))) ! DIIS extrapolation