From fde7d422023ea29c49cc6314521844440d790fef Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 27 Jan 2025 16:52:08 +0100 Subject: [PATCH] Starting chem pot --- src/HF/HFB.f90 | 17 +++++++++++------ src/HF/print_HFB.f90 | 5 ++++- src/QuAcK/BQuAcK.f90 | 12 +++++------- 3 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 44ee47d..a764e88 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -1,5 +1,5 @@ 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,Panom,F) + nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EHFB,eHF,c,P,Panom,F,Delta) ! Perform Hartree-Fock Bogoliubov calculation @@ -40,6 +40,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision :: EJ double precision :: EK double precision :: EL + double precision :: chem_pot double precision :: dipole(ncart) double precision :: Conv @@ -62,6 +63,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision,intent(out) :: P(nBas,nBas) double precision,intent(out) :: Panom(nBas,nBas) double precision,intent(out) :: F(nBas,nBas) + double precision,intent(out) :: Delta(nBas,nBas) ! Hello world @@ -89,10 +91,11 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & allocate(err_diis(nBas_Sq,max_diis)) allocate(F_diis(nBas_Sq,max_diis)) -! Guess coefficients and density matrix +! Guess coefficients, density matrix, and chem. pot + chem_pot = 0d0 P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) - Panom(:,:) = -P(:,:) + Panom(:,:) = -P(:,:) ! Do sth TODO ! Initialization @@ -127,6 +130,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & call anomalous_matrix_AO_basis(nBas,Panom,ERI,L) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + Delta(:,:) = L(:,:) ! Check convergence @@ -151,11 +155,11 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Anomalous energy - EL = 0.25d0*trace_matrix(nBas,matmul(-Panom,L)) + EL = -0.25d0*trace_matrix(nBas,matmul(Panom,L)) ! Total energy - EHFB = ET + EV + EJ + EK + EHFB = ET + EV + EJ + EK + EL ! DIIS extrapolation @@ -182,6 +186,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Density matrix P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + Panom(:,:) = -P(:,:) ! Do sth TODO ! Dump results @@ -213,7 +218,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! 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,EL,EHFB,dipole) + call print_HFB(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole) ! Testing zone diff --git a/src/HF/print_HFB.f90 b/src/HF/print_HFB.f90 index e3aa5e7..ef25a8e 100644 --- a/src/HF/print_HFB.f90 +++ b/src/HF/print_HFB.f90 @@ -1,7 +1,7 @@ ! --- -subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, dipole) +subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_pot, dipole) ! Print one-electron energies and other stuff for G0W0 @@ -21,6 +21,7 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, d double precision,intent(in) :: EK double precision,intent(in) :: EL double precision,intent(in) :: ERHF + double precision,intent(in) :: chem_pot double precision,intent(in) :: dipole(ncart) ! Local variables @@ -61,6 +62,8 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, d 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.10,A3)') ' Chemical potential = ',chem_pot,' au' + 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 diff --git a/src/QuAcK/BQuAcK.f90 b/src/QuAcK/BQuAcK.f90 index edd38f0..01ff9dd 100644 --- a/src/QuAcK/BQuAcK.f90 +++ b/src/QuAcK/BQuAcK.f90 @@ -43,10 +43,9 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, double precision,allocatable :: PHF(:,:) double precision,allocatable :: PanomHF(:,:) double precision,allocatable :: FHF(:,:) + double precision,allocatable :: Delta(:,:) double precision :: ERHF,EHFB -! double precision,allocatable :: dipole_int_MO(:,:,:) double precision,allocatable :: ERI_AO(:,:,:,:) -! double precision,allocatable :: ERI_MO(:,:,:,:) write(*,*) write(*,*) '******************************' @@ -63,8 +62,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, allocate(PHF(nBas,nBas)) allocate(PanomHF(nBas,nBas)) allocate(FHF(nBas,nBas)) -! allocate(dipole_int_MO(nOrb,nOrb,ncart)) -! allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) + allocate(Delta(nBas,nBas)) allocate(ERI_AO(nBas,nBas,nBas,nBas)) call wall_time(start_int) @@ -94,7 +92,8 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, ! 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,PanomHF,FHF) + nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHFB,eHF,cHF,PHF,PanomHF, & + FHF,Delta) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -110,8 +109,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, deallocate(PHF) deallocate(PanomHF) deallocate(FHF) -! deallocate(dipole_int_MO) -! deallocate(ERI_MO) + deallocate(Delta) deallocate(ERI_AO) end subroutine