From cf1ad337e594a5ba31ba16523204e1348fb1acfe Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 27 Jan 2025 15:24:43 +0100 Subject: [PATCH] Starting HFB scf --- src/HF/HFB.f90 | 229 +++++++++++++++++++++++++++++++++++++++++++ src/HF/print_HFB.f90 | 90 +++++++++++++++++ src/QuAcK/BQuAcK.f90 | 80 ++++++++++++++- src/QuAcK/QuAcK.f90 | 3 +- 4 files changed, 396 insertions(+), 6 deletions(-) create mode 100644 src/HF/HFB.f90 create mode 100644 src/HF/print_HFB.f90 diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 new file mode 100644 index 0000000..5fab200 --- /dev/null +++ b/src/HF/HFB.f90 @@ -0,0 +1,229 @@ +subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EHFB,eHF,c,P,F) + +! Perform Hartree-Fock Bogoliubov calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + double precision,intent(in) :: level_shift + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + 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,nOrb) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: nSCF + integer :: nBas_Sq + integer :: n_diis + double precision :: ET + double precision :: EV + double precision :: EJ + double precision :: EK + double precision :: dipole(ncart) + + double precision :: Conv + double precision :: rcond + double precision,external :: trace_matrix + double precision,allocatable :: err(:,:) + double precision,allocatable :: err_diis(:,:) + double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: J(:,:) + double precision,allocatable :: K(:,:) + double precision,allocatable :: cp(:,:) + double precision,allocatable :: Fp(:,:) + +! Output variables + + double precision,intent(out) :: EHFB + double precision,intent(out) :: eHF(nOrb) + double precision,intent(inout):: c(nBas,nOrb) + double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: F(nBas,nBas) + +! Hello world + + write(*,*) + write(*,*)'*****************************' + write(*,*)'* HF Bogoliubov Calculation *' + write(*,*)'*****************************' + write(*,*) + +! Useful quantities + + nBas_Sq = nBas*nBas + +! Memory allocation + + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + + allocate(err(nBas,nBas)) + + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + + allocate(err_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) + +! Guess coefficients and density matrix + + P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) +! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & +! c(1,1), nBas, c(1,1), nBas, & +! 0.d0, P(1,1), nBas) + +! Initialization + + n_diis = 0 + F_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 + rcond = 0d0 + + Conv = 1d0 + nSCF = 0 + +!------------------------------------------------------------------------ +! 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)') & + '|','#','|','E(HFB)','|','EJ(HFB)','|','EK(HFB)','|','Conv','|' + write(*,*)'-----------------------------------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Build Fock matrix + + call Hartree_matrix_AO_basis(nBas,P,ERI,J) + call exchange_matrix_AO_basis(nBas,P,ERI,K) + + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + + ! Check convergence + + 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)) + + ! Potential energy + + EV = trace_matrix(nBas,matmul(P,V)) + + ! Hartree energy + + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + + ! Exchange energy + + EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) + + ! Total energy + + EHFB = ET + EV + EJ + EK + + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) + + end if + + ! Level shift + + if(level_shift > 0d0 .and. Conv > thresh) then + call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) + endif + + ! Diagonalize Fock matrix + + Fp = matmul(transpose(X),matmul(F,X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nOrb,cp,eHF) + c = matmul(X,cp) + + ! Density matrix + + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) +! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & +! c(1,1), nBas, c(1,1), nBas, & +! 0.d0, P(1,1), nBas) + + ! 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,E10.2,1X,A1,1X)') & + '|',nSCF,'|',EHFB + ENuc,'|',EJ,'|',EK,'|',Conv,'|' + + end do + write(*,*)'-----------------------------------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + deallocate(J,K,err,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_HFB(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,EHFB,dipole) + +! Testing zone + + if(dotest) then + + call dump_test_value('R','HFB energy',EHFB) + call dump_test_value('R','HFB HOMO energy',eHF(nO)) + call dump_test_value('R','HFB LUMO energy',eHF(nO+1)) + call dump_test_value('R','HFB dipole moment',norm2(dipole)) + + end if + +! Memory deallocation + + deallocate(J,K,err,cp,Fp,err_diis,F_diis) + +end subroutine diff --git a/src/HF/print_HFB.f90 b/src/HF/print_HFB.f90 new file mode 100644 index 0000000..75efb4c --- /dev/null +++ b/src/HF/print_HFB.f90 @@ -0,0 +1,90 @@ + +! --- + +subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipole) + +! Print one-electron energies and other stuff for G0W0 + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas, nOrb + integer,intent(in) :: nO + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: ENuc + double precision,intent(in) :: ET + double precision,intent(in) :: EV + double precision,intent(in) :: EJ + double precision,intent(in) :: EK + double precision,intent(in) :: ERHF + double precision,intent(in) :: dipole(ncart) + +! Local variables + + integer :: ixyz + integer :: HOMO + integer :: LUMO + double precision :: Gap + double precision :: S,S2 + + logical :: dump_orb = .false. + +! HOMO and LUMO + + HOMO = nO + LUMO = HOMO + 1 + Gap = eHF(LUMO)-eHF(HOMO) + + S2 = 0d0 + S = 0d0 + +! Dump results + + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A33)') ' Summary ' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',ET + EV,' au' + write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',ET,' au' + write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',EV,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',EJ + EK,' au' + write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',EJ,' au' + write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',EK,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',ERHF,' au' + write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au' + write(*,'(A33,1X,F16.10,A3)') ' HFB energy = ',ERHF + ENuc,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.6,A3)') ' HFB HOMO energy = ',eHF(HOMO)*HaToeV,' eV' + write(*,'(A33,1X,F16.6,A3)') ' HFB LUMO energy = ',eHF(LUMO)*HaToeV,' eV' + write(*,'(A33,1X,F16.6,A3)') ' HFB HOMO-LUMO gap = ',Gap*HaToeV,' eV' + write(*,'(A50)') '---------------------------------------' + write(*,'(A33,1X,F16.6)') ' = ',S + write(*,'(A33,1X,F16.6)') ' = ',S2 + write(*,'(A50)') '---------------------------------------' + write(*,'(A36)') ' Dipole moment (Debye) ' + write(*,'(10X,4A10)') 'X','Y','Z','Tot.' + write(*,'(10X,4F10.4)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + write(*,'(A50)') '---------------------------------------' + write(*,*) + +! Print results + + if(dump_orb) then + write(*,'(A50)') '---------------------------------------' + write(*,'(A50)') ' HFB orbital coefficients ' + write(*,'(A50)') '---------------------------------------' + call matout(nBas, nOrb, cHF) + write(*,*) + end if + write(*,'(A50)') '---------------------------------------' + write(*,'(A50)') ' HFB orbital energies (au) ' + write(*,'(A50)') '---------------------------------------' + call vecout(nOrb, eHF) + write(*,*) + +end subroutine diff --git a/src/QuAcK/BQuAcK.f90 b/src/QuAcK/BQuAcK.f90 index 889b518..a6a2749 100644 --- a/src/QuAcK/BQuAcK.f90 +++ b/src/QuAcK/BQuAcK.f90 @@ -1,4 +1,6 @@ -subroutine BQuAcK(working_dir,dotest,doHFB) +subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & + guess_type,mix) ! Restricted branch of QuAcK @@ -11,10 +13,40 @@ subroutine BQuAcK(working_dir,dotest,doHFB) logical,intent(in) :: doHFB + integer,intent(in) :: nNuc,nBas,nOrb + 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) + + 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,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + + integer,intent(in) :: maxSCF_HF,max_diis_HF + double precision,intent(in) :: thresh_HF,level_shift,mix + integer,intent(in) :: guess_type + ! Local variables double precision :: start_HF ,end_HF ,t_HF + double precision :: start_int, end_int, t_int + double precision,allocatable :: eHF(:) + double precision,allocatable :: cHF(:,:) + double precision,allocatable :: PHF(:,:) + double precision,allocatable :: FHF(:,:) + double precision :: ERHF,EHFB +! double precision,allocatable :: dipole_int_MO(:,:,:) + double precision,allocatable :: ERI_AO(:,:,:,:) +! double precision,allocatable :: ERI_MO(:,:,:,:) + write(*,*) write(*,*) '******************************' write(*,*) '* Bogoliubov Branch of QuAcK *' @@ -25,14 +57,42 @@ subroutine BQuAcK(working_dir,dotest,doHFB) ! Memory allocation ! !-------------------! -!---------------------! -! Hartree-Fock module ! -!---------------------! + allocate(eHF(nOrb)) + allocate(cHF(nBas,nOrb)) + allocate(PHF(nBas,nBas)) + allocate(FHF(nBas,nBas)) +! allocate(dipole_int_MO(nOrb,nOrb,ncart)) +! allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) + + allocate(ERI_AO(nBas,nBas,nBas,nBas)) + call wall_time(start_int) + call read_2e_integrals(working_dir,nBas,ERI_AO) + call wall_time(end_int) + t_int = end_int - start_int + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds' + write(*,*) + +!--------------------------------! +! Hartree-Fock Bogoliubov module ! +!--------------------------------! if(doHFB) then + ! Run first a RHF calculation call wall_time(start_HF) -! call HFB(dotest) + call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF) + call wall_time(end_HF) + + t_HF = end_HF - start_HF + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds' + write(*,*) + + ! Continue with a HFB calculation + call wall_time(start_HF) + call HFB(dotest,maxSCF_HF,thresh_HF,max_diis_HF,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHFB,eHF,cHF,PHF,FHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -41,4 +101,14 @@ subroutine BQuAcK(working_dir,dotest,doHFB) end if +! Memory deallocation + + deallocate(eHF) + deallocate(cHF) + deallocate(PHF) + deallocate(FHF) +! deallocate(dipole_int_MO) +! deallocate(ERI_MO) + deallocate(ERI_AO) + end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 516ca40..0342309 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -289,7 +289,8 @@ program QuAcK ! Bogoliubov QuAcK branch ! !--------------------------! if(doBQuAcK) & - call BQuAcK(working_dir,doGtest,doHFB) + call BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix) !-----------! ! Stop Test !