From 12044b800dd54b56155d0aa09407cbadf832f891 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Sat, 25 Jan 2025 23:15:59 +0100 Subject: [PATCH 01/32] Starting bogoliubov --- PyDuck.py | 2 +- src/QuAcK/QuAcK.f90 | 2 +- src/QuAcK/read_methods.f90 | 8 +++++--- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/PyDuck.py b/PyDuck.py index 4779566b..b8706d68 100644 --- a/PyDuck.py +++ b/PyDuck.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 import os, sys import argparse diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 0954b000..516ca40d 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -108,7 +108,7 @@ program QuAcK !------------------! call read_methods(working_dir, & - doRHF,doUHF,doGHF,doROHF, & + doRHF,doUHF,doGHF,doROHF,doHFB, & doMP2,doMP3, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & dodrCCD,dorCCD,docrCCD,dolCCD, & diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index bef33057..4869262e 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,5 +1,5 @@ subroutine read_methods(working_dir, & - doRHF,doUHF,doGHF,doROHF, & + doRHF,doUHF,doGHF,doROHF,doHFB, & doMP2,doMP3, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_crCCD,do_lCCD, & @@ -22,7 +22,7 @@ subroutine read_methods(working_dir, & ! Output variables - logical,intent(out) :: doRHF,doUHF,doGHF,doROHF + logical,intent(out) :: doRHF,doUHF,doGHF,doROHF,doHFB logical,intent(out) :: doMP2,doMP3 logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD @@ -58,13 +58,15 @@ subroutine read_methods(working_dir, & doUHF = .false. doGHF = .false. doROHF = .false. + doHFB = .false. read(1,*) - read(1,*) ans1,ans2,ans3,ans4 + read(1,*) ans1,ans2,ans3,ans4,ans5 if(ans1 == 'T') doRHF = .true. if(ans2 == 'T') doUHF = .true. if(ans3 == 'T') doGHF = .true. if(ans4 == 'T') doROHF = .true. + if(ans5 == 'T') doHFB = .true. ! Read MPn methods From cf1ad337e594a5ba31ba16523204e1348fb1acfe Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 27 Jan 2025 15:24:43 +0100 Subject: [PATCH 02/32] 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 00000000..5fab200b --- /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 00000000..75efb4c2 --- /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 889b518d..a6a27496 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 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) + 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 516ca40d..0342309a 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 ! From 79fada93da8addc58825088e1b048ecacba912db Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 27 Jan 2025 16:21:35 +0100 Subject: [PATCH 03/32] Starting Panomalous --- src/AOtoMO/Hartree_matrix_AO_basis.f90 | 6 ++-- src/AOtoMO/anomalous_matrix_AO_basis.f90 | 37 +++++++++++++++++++++++ src/AOtoMO/exchange_matrix_AO_basis.f90 | 3 +- src/HF/HFB.f90 | 38 +++++++++++++----------- src/HF/print_HFB.f90 | 11 ++----- src/QuAcK/BQuAcK.f90 | 5 +++- 6 files changed, 70 insertions(+), 30 deletions(-) create mode 100644 src/AOtoMO/anomalous_matrix_AO_basis.f90 diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index 63141cc1..8b75a2bb 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_AO_basis.f90 @@ -1,4 +1,4 @@ -subroutine Hartree_matrix_AO_basis(nBas,P,G,H) +subroutine Hartree_matrix_AO_basis(nBas,P,ERI,H) ! Compute Hartree matrix in the AO basis @@ -9,7 +9,7 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H) integer,intent(in) :: nBas double precision,intent(in) :: P(nBas,nBas) - double precision,intent(in) :: G(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables @@ -25,7 +25,7 @@ subroutine Hartree_matrix_AO_basis(nBas,P,G,H) do nu=1,nBas do la=1,nBas do mu=1,nBas - H(mu,nu) = H(mu,nu) + P(la,si)*G(mu,la,nu,si) + H(mu,nu) = H(mu,nu) + P(la,si)*ERI(mu,la,nu,si) end do end do end do diff --git a/src/AOtoMO/anomalous_matrix_AO_basis.f90 b/src/AOtoMO/anomalous_matrix_AO_basis.f90 new file mode 100644 index 00000000..4f2583ef --- /dev/null +++ b/src/AOtoMO/anomalous_matrix_AO_basis.f90 @@ -0,0 +1,37 @@ +subroutine anomalous_matrix_AO_basis(nBas,Pa,ERI,L) + +! Compute anomalous L matrix in the AO basis + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: Pa(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: mu,nu,la,si + +! Output variables + + double precision,intent(out) :: L(nBas,nBas) + + L(:,:) = 0d0 + + do nu=1,nBas + do si=1,nBas + do la=1,nBas + do mu=1,nBas + L(mu,nu) = L(mu,nu) + Pa(la,si)*ERI(la,si,mu,nu) + end do + end do + end do + end do + +end subroutine + +! --- + diff --git a/src/AOtoMO/exchange_matrix_AO_basis.f90 b/src/AOtoMO/exchange_matrix_AO_basis.f90 index 8042420f..7e04096d 100644 --- a/src/AOtoMO/exchange_matrix_AO_basis.f90 +++ b/src/AOtoMO/exchange_matrix_AO_basis.f90 @@ -19,7 +19,8 @@ subroutine exchange_matrix_AO_basis(nBas,P,ERI,K) double precision,intent(out) :: K(nBas,nBas) - K = 0d0 + K(:,:) = 0d0 + do nu=1,nBas do si=1,nBas do la=1,nBas diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 5fab200b..44ee47db 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,F) + nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EHFB,eHF,c,P,Panom,F) ! Perform Hartree-Fock Bogoliubov calculation @@ -39,6 +39,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision :: EV double precision :: EJ double precision :: EK + double precision :: EL double precision :: dipole(ncart) double precision :: Conv @@ -49,6 +50,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision,allocatable :: F_diis(:,:) double precision,allocatable :: J(:,:) double precision,allocatable :: K(:,:) + double precision,allocatable :: L(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: Fp(:,:) @@ -58,6 +60,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision,intent(out) :: eHF(nOrb) double precision,intent(inout):: c(nBas,nOrb) double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: Panom(nBas,nBas) double precision,intent(out) :: F(nBas,nBas) ! Hello world @@ -76,6 +79,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & allocate(J(nBas,nBas)) allocate(K(nBas,nBas)) + allocate(L(nBas,nBas)) allocate(err(nBas,nBas)) @@ -88,9 +92,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! 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) + Panom(:,:) = -P(:,:) ! Initialization @@ -107,10 +109,10 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & !------------------------------------------------------------------------ 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(*,*)'-----------------------------------------------------------------------------' + write(*,*)'-----------------------------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1A16,1X,A1,1X,A10,2X,A1,1X)') & + '|','#','|','E(HFB)','|','EJ(HFB)','|','EK(HFB)','|','EL(HFB)','|','Conv','|' + write(*,*)'-----------------------------------------------------------------------------------------------' do while(Conv > thresh .and. nSCF < maxSCF) @@ -122,6 +124,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & call Hartree_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) + call anomalous_matrix_AO_basis(nBas,Panom,ERI,L) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) @@ -146,6 +149,10 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) + ! Anomalous energy + + EL = 0.25d0*trace_matrix(nBas,matmul(-Panom,L)) + ! Total energy EHFB = ET + EV + EJ + EK @@ -175,17 +182,14 @@ 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))) -! 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,'|' + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1XF16.10,1X,A1,1X,E10.2,1X,A1,1X)') & + '|',nSCF,'|',EHFB + ENuc,'|',EJ,'|',EK,'|',EL,'|',Conv,'|' end do - write(*,*)'-----------------------------------------------------------------------------' + write(*,*)'-----------------------------------------------------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ @@ -200,7 +204,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,err,cp,Fp,err_diis,F_diis) + deallocate(J,K,L,err,cp,Fp,err_diis,F_diis) stop @@ -209,7 +213,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,EHFB,dipole) + call print_HFB(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,EL,EHFB,dipole) ! Testing zone @@ -224,6 +228,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Memory deallocation - deallocate(J,K,err,cp,Fp,err_diis,F_diis) + deallocate(J,K,L,err,cp,Fp,err_diis,F_diis) end subroutine diff --git a/src/HF/print_HFB.f90 b/src/HF/print_HFB.f90 index 75efb4c2..e3aa5e7b 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, ERHF, dipole) +subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, dipole) ! Print one-electron energies and other stuff for G0W0 @@ -19,6 +19,7 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipol double precision,intent(in) :: EV double precision,intent(in) :: EJ double precision,intent(in) :: EK + double precision,intent(in) :: EL double precision,intent(in) :: ERHF double precision,intent(in) :: dipole(ncart) @@ -28,7 +29,6 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipol integer :: HOMO integer :: LUMO double precision :: Gap - double precision :: S,S2 logical :: dump_orb = .false. @@ -38,9 +38,6 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipol LUMO = HOMO + 1 Gap = eHF(LUMO)-eHF(HOMO) - S2 = 0d0 - S = 0d0 - ! Dump results write(*,*) @@ -54,6 +51,7 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipol 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(*,'(A33,1X,F16.10,A3)') ' Anomalous energy = ',EL,' au' write(*,'(A50)') '---------------------------------------' write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',ERHF,' au' write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au' @@ -63,9 +61,6 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipol 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 diff --git a/src/QuAcK/BQuAcK.f90 b/src/QuAcK/BQuAcK.f90 index a6a27496..edd38f0b 100644 --- a/src/QuAcK/BQuAcK.f90 +++ b/src/QuAcK/BQuAcK.f90 @@ -41,6 +41,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, double precision,allocatable :: eHF(:) double precision,allocatable :: cHF(:,:) double precision,allocatable :: PHF(:,:) + double precision,allocatable :: PanomHF(:,:) double precision,allocatable :: FHF(:,:) double precision :: ERHF,EHFB ! double precision,allocatable :: dipole_int_MO(:,:,:) @@ -60,6 +61,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, allocate(eHF(nOrb)) allocate(cHF(nBas,nOrb)) 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)) @@ -92,7 +94,7 @@ 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,FHF) + nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHFB,eHF,cHF,PHF,PanomHF,FHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -106,6 +108,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, deallocate(eHF) deallocate(cHF) deallocate(PHF) + deallocate(PanomHF) deallocate(FHF) ! deallocate(dipole_int_MO) ! deallocate(ERI_MO) From fde7d422023ea29c49cc6314521844440d790fef Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 27 Jan 2025 16:52:08 +0100 Subject: [PATCH 04/32] 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 44ee47db..a764e88c 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 e3aa5e7b..ef25a8ed 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 edd38f0b..01ff9dda 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 From f688dd6cec477c470536f69a5068a34cfad37cdd Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Tue, 28 Jan 2025 16:52:17 +0100 Subject: [PATCH 05/32] Incl up to reshufle after MOM --- src/HF/HFB.f90 | 109 ++++++++++++++++++++++++++++------------ src/HF/reshufle_hfb.f90 | 90 +++++++++++++++++++++++++++++++++ 2 files changed, 166 insertions(+), 33 deletions(-) create mode 100644 src/HF/reshufle_hfb.f90 diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index a764e88c..e5419bc6 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -32,7 +32,10 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Local variables + integer :: iorb integer :: nSCF + integer :: nOrb2 + integer :: nBas2 integer :: nBas_Sq integer :: n_diis double precision :: ET @@ -46,19 +49,25 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision :: Conv double precision :: rcond double precision,external :: trace_matrix + double precision,allocatable :: eHFB_(:) + double precision,allocatable :: project(:) 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 :: L(:,:) double precision,allocatable :: cp(:,:) - double precision,allocatable :: Fp(:,:) + double precision,allocatable :: cOld_T(:,:) + double precision,allocatable :: S_hfb(:,:) + double precision,allocatable :: X_hfb(:,:) + double precision,allocatable :: H_hfb(:,:) + double precision,allocatable :: c_hfb(:,:) + double precision,allocatable :: mom(:,:) ! Output variables double precision,intent(out) :: EHFB - double precision,intent(out) :: eHF(nOrb) + double precision,intent(inout):: eHF(nOrb) double precision,intent(inout):: c(nBas,nOrb) double precision,intent(out) :: P(nBas,nBas) double precision,intent(out) :: Panom(nBas,nBas) @@ -76,26 +85,32 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Useful quantities nBas_Sq = nBas*nBas + nBas2 = nBas+nBas + nOrb2 = nOrb+nOrb ! Memory allocation allocate(J(nBas,nBas)) allocate(K(nBas,nBas)) - allocate(L(nBas,nBas)) allocate(err(nBas,nBas)) - allocate(cp(nOrb,nOrb)) - allocate(Fp(nOrb,nOrb)) + allocate(cp(nOrb2,nOrb2)) + allocate(H_hfb(nOrb2,nOrb2)) + allocate(cOld_T(nOrb2,nBas2)) + allocate(S_hfb(nBas2,nBas2)) + allocate(X_hfb(nBas2,nOrb2)) + allocate(c_hfb(nBas2,nOrb2)) + allocate(mom(nOrb2,nOrb2)) + allocate(eHFB_(nOrb2)) + allocate(project(nOrb2)) allocate(err_diis(nBas_Sq,max_diis)) allocate(F_diis(nBas_Sq,max_diis)) -! Guess coefficients, density matrix, and chem. pot +! Guess coefficients chem. pot - chem_pot = 0d0 - P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) - Panom(:,:) = -P(:,:) ! Do sth TODO + chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1)) ! Initialization @@ -104,6 +119,18 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & err_diis(:,:) = 0d0 rcond = 0d0 + cOld_T(:,:) = 0d0 + cOld_T(1:nOrb,1:nbas) = transpose(c) + S_hfb(:,:) = 0d0 + S_hfb(1:nBas ,1:nBas ) = S(1:nBas,1:nBas) + S_hfb(nBas+1:nBas2,nBas+1:nBas2) = S(1:nBas,1:nBas) + X_hfb(:,:) = 0d0 + X_hfb(1:nBas ,1:nOrb) = X(1:nBas,1:nOrb) + X_hfb(nBas+1:nBas2,nOrb+1:nOrb2) = X(1:nBas,1:nOrb) + + P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) + Panom(:,:) = 0d0 !-P(:,:) ! Do sth TODO + Conv = 1d0 nSCF = 0 @@ -127,10 +154,9 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & call Hartree_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) - call anomalous_matrix_AO_basis(nBas,Panom,ERI,L) + call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - Delta(:,:) = L(:,:) + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) ! Check convergence @@ -155,38 +181,53 @@ 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,Delta)) ! Total energy - EHFB = ET + EV + EJ + EK + EL + EHFB = ET + EV + EJ + EK !+ EL - ! DIIS extrapolation + ! ! DIIS extrapolation TODO check and adapt + ! + ! 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 - if(max_diis > 1) then + ! ! Level shift TODO check and adapt + ! + ! if(level_shift > 0d0 .and. Conv > thresh) then + ! call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) + ! endif - 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 + ! Diagonalize Fock matrix - ! Level shift + H_hfb(:,:) = 0d0 + H_hfb(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X)) + H_hfb(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_hfb(1:nOrb,1:nOrb) + H_hfb(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) + H_hfb(nOrb+1:nOrb2,1:nOrb ) = transpose(H_hfb(1:nOrb,nOrb+1:nOrb2)) + + cp(:,:) = H_hfb(:,:) + call diagonalize_matrix(nOrb2,cp,eHFB_) + c_hfb = matmul(X_hfb,cp) - if(level_shift > 0d0 .and. Conv > thresh) then - call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) - endif + ! Compute MOM and reshufle the states + mom=matmul(cOld_T,matmul(S_hfb,c_hfb)) + do iorb=1,nOrb2 + project(iorb)=sqrt(dot_product(mom(:,iorb),mom(:,iorb))) + enddo + call reshufle_hfb(nBas2,nOrb2,c_hfb,eHFB_,project) - ! Diagonalize Fock matrix - Fp = matmul(transpose(X),matmul(F,X)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nOrb,cp,eHF) - c = matmul(X,cp) + !TODO extract c ! Density matrix P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - Panom(:,:) = -P(:,:) ! Do sth TODO + Panom(:,:) = 0d0 ! Do sth TODO ! Dump results @@ -209,7 +250,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,L,err,cp,Fp,err_diis,F_diis) + deallocate(J,K,err,cp,H_hfb,X_hfb,c_hfb,S_hfb,mom,cOld_T,eHFB_,project,err_diis,F_diis) stop @@ -217,6 +258,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Compute dipole moments +eHF(:)=eHF(:)-chem_pot!TODO remove + 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,chem_pot,dipole) @@ -233,6 +276,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Memory deallocation - deallocate(J,K,L,err,cp,Fp,err_diis,F_diis) + deallocate(J,K,err,cp,H_hfb,X_hfb,c_hfb,S_hfb,mom,cOld_T,eHFB_,project,err_diis,F_diis) end subroutine diff --git a/src/HF/reshufle_hfb.f90 b/src/HF/reshufle_hfb.f90 new file mode 100644 index 00000000..7ff97e73 --- /dev/null +++ b/src/HF/reshufle_hfb.f90 @@ -0,0 +1,90 @@ + +! --- + +subroutine reshufle_hfb(nBas2,nOrb2,c_hfb,eHFB,project) + +! Print one-electron energies and other stuff for G0W0 + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas2, nOrb2 + double precision,intent(inout) :: eHFB(nOrb2) + double precision,intent(inout) :: project(nOrb2) + double precision,intent(inout) :: c_hfb(nBas2,nOrb2) + +! Local variables + + integer :: iorb,iorb2,ipos + double precision :: val + double precision,allocatable :: eHFB_tmp(:) + double precision,allocatable :: project_tmp(:) + double precision,allocatable :: c_tmp(:,:) + + allocate(c_tmp(nBas2,nOrb2),project_tmp(nOrb2),eHFB_tmp(nOrb2)) + +!write(*,*) 'e_HFB' +!do iorb=1,nOrb2 +!write(*,*) eHFB(iorb),project(iorb) +!enddo + + ! Reshufle using MOM + do iorb=1,nOrb2 + val=-1d0 + ipos=0 + do iorb2=1,nOrb2 + if(project(iorb2)>val) then + val=project(iorb2) + ipos=iorb2 + endif + enddo + project(ipos)=-1.0d1 + project_tmp(iorb)=val + eHFB_tmp(iorb)=eHFB(ipos) + c_tmp(:,iorb)=c_hfb(:,ipos) + enddo + + ! Reshufle using energies + do iorb=1,nOrb2/2 ! electronic 1-RDM + val=1d10 + ipos=0 + do iorb2=1,nOrb2/2 + if(eHFB_tmp(iorb2) Date: Tue, 28 Jan 2025 21:40:17 +0100 Subject: [PATCH 06/32] Added one space, it is doing MOM --- src/HF/HFB.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index e5419bc6..29959062 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -215,6 +215,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & c_hfb = matmul(X_hfb,cp) ! Compute MOM and reshufle the states + mom=matmul(cOld_T,matmul(S_hfb,c_hfb)) do iorb=1,nOrb2 project(iorb)=sqrt(dot_product(mom(:,iorb),mom(:,iorb))) From 2cf4aa901bea7efc8e85485f895bd0d3d60a04f4 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Tue, 28 Jan 2025 23:16:17 +0100 Subject: [PATCH 07/32] Building R instead of MOM --- src/HF/HFB.f90 | 187 ++++++++++++++++++++++++++----------------------- 1 file changed, 100 insertions(+), 87 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 29959062..16be4a6f 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -32,8 +32,10 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Local variables + integer :: maxSCF_chem_pot integer :: iorb integer :: nSCF + integer :: nSCF_chem_pot integer :: nOrb2 integer :: nBas2 integer :: nBas_Sq @@ -48,21 +50,19 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision :: Conv double precision :: rcond + double precision :: delta_chem_pot + double precision :: trace_1rdm + double precision :: thrs_N double precision,external :: trace_matrix double precision,allocatable :: eHFB_(:) - double precision,allocatable :: project(:) 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 :: cOld_T(:,:) - double precision,allocatable :: S_hfb(:,:) - double precision,allocatable :: X_hfb(:,:) double precision,allocatable :: H_hfb(:,:) - double precision,allocatable :: c_hfb(:,:) - double precision,allocatable :: mom(:,:) + double precision,allocatable :: R(:,:) ! Output variables @@ -97,13 +97,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & allocate(cp(nOrb2,nOrb2)) allocate(H_hfb(nOrb2,nOrb2)) - allocate(cOld_T(nOrb2,nBas2)) - allocate(S_hfb(nBas2,nBas2)) - allocate(X_hfb(nBas2,nOrb2)) - allocate(c_hfb(nBas2,nOrb2)) - allocate(mom(nOrb2,nOrb2)) + allocate(R(nOrb2,nOrb2)) allocate(eHFB_(nOrb2)) - allocate(project(nOrb2)) allocate(err_diis(nBas_Sq,max_diis)) allocate(F_diis(nBas_Sq,max_diis)) @@ -114,22 +109,15 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Initialization - n_diis = 0 - F_diis(:,:) = 0d0 - err_diis(:,:) = 0d0 - rcond = 0d0 - - cOld_T(:,:) = 0d0 - cOld_T(1:nOrb,1:nbas) = transpose(c) - S_hfb(:,:) = 0d0 - S_hfb(1:nBas ,1:nBas ) = S(1:nBas,1:nBas) - S_hfb(nBas+1:nBas2,nBas+1:nBas2) = S(1:nBas,1:nBas) - X_hfb(:,:) = 0d0 - X_hfb(1:nBas ,1:nOrb) = X(1:nBas,1:nOrb) - X_hfb(nBas+1:nBas2,nOrb+1:nOrb2) = X(1:nBas,1:nOrb) + thrs_N = 1d-6 + maxSCF_chem_pot = 1000 + n_diis = 0 + F_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 + rcond = 0d0 P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) - Panom(:,:) = 0d0 !-P(:,:) ! Do sth TODO + Panom(:,:) = 0d0 ! Do sth TODO Conv = 1d0 nSCF = 0 @@ -138,11 +126,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Main SCF loop !------------------------------------------------------------------------ - write(*,*) - write(*,*)'-----------------------------------------------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1A16,1X,A1,1X,A10,2X,A1,1X)') & - '|','#','|','E(HFB)','|','EJ(HFB)','|','EK(HFB)','|','EL(HFB)','|','Conv','|' - write(*,*)'-----------------------------------------------------------------------------------------------' do while(Conv > thresh .and. nSCF < maxSCF) @@ -150,18 +133,63 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & nSCF = nSCF + 1 - ! Build Fock matrix - - call Hartree_matrix_AO_basis(nBas,P,ERI,J) - call exchange_matrix_AO_basis(nBas,P,ERI,K) - call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) + write(*,*) + write(*,*)'-------------------------------------' + write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') & + '|','Tr[1D]','|','Chem. Pot.','|' + write(*,*)'-------------------------------------' + + ! Loop to adjust chem_pot + delta_chem_pot=1d0 + nSCF_chem_pot=0 + do + + ! Build Fock matrix + + call Hartree_matrix_AO_basis(nBas,P,ERI,J) + call exchange_matrix_AO_basis(nBas,P,ERI,K) + call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) + + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) + + ! Diagonalize H_HFB matrix + + H_hfb(:,:) = 0d0 + H_hfb(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X)) + H_hfb(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_hfb(1:nOrb,1:nOrb) + H_hfb(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) + H_hfb(nOrb+1:nOrb2,1:nOrb ) = transpose(H_hfb(1:nOrb,nOrb+1:nOrb2)) + + cp(:,:) = H_hfb(:,:) + call diagonalize_matrix(nOrb2,cp,eHFB_) + + ! Build R and extract P and Panom + + trace_1rdm = 0d0 + R(:,:) = 0d0 + do iorb=1,nOrb + R(:,:) = R(:,:) + matmul(cp(:,iorb:iorb),transpose(cp(:,iorb:iorb))) + enddo + do iorb=1,nOrb + trace_1rdm = trace_1rdm + R(iorb,iorb) + enddo + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|' + + nSCF_chem_pot = nSCF_chem_pot + 1 + if( abs(trace_1rdm-nO) < thrs_N .or. nSCF_chem_pot > maxSCF_chem_pot ) exit + + + + enddo + + ! Extract P and Panom from R - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) - - ! Check convergence + P(:,:) = 0d0 + Panom(:,:) = 0d0 + P(:,:) = 2d0*matmul(X,matmul(R(1:nOrb,1:nOrb),transpose(X))) + Panom(:,:) = matmul(X,matmul(R(1:nOrb,nOrb+1:nOrb2),transpose(X))) - err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) - if(nSCF > 1) Conv = maxval(abs(err)) ! Kinetic energy @@ -187,56 +215,23 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & EHFB = ET + EV + EJ + EK !+ EL - ! ! DIIS extrapolation TODO check and adapt - ! - ! 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 TODO check and adapt - ! - ! if(level_shift > 0d0 .and. Conv > thresh) then - ! call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) - ! endif - - ! Diagonalize Fock matrix - - H_hfb(:,:) = 0d0 - H_hfb(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X)) - H_hfb(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_hfb(1:nOrb,1:nOrb) - H_hfb(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) - H_hfb(nOrb+1:nOrb2,1:nOrb ) = transpose(H_hfb(1:nOrb,nOrb+1:nOrb2)) - - cp(:,:) = H_hfb(:,:) - call diagonalize_matrix(nOrb2,cp,eHFB_) - c_hfb = matmul(X_hfb,cp) - - ! Compute MOM and reshufle the states - - mom=matmul(cOld_T,matmul(S_hfb,c_hfb)) - do iorb=1,nOrb2 - project(iorb)=sqrt(dot_product(mom(:,iorb),mom(:,iorb))) - enddo - call reshufle_hfb(nBas2,nOrb2,c_hfb,eHFB_,project) - - - !TODO extract c - - ! Density matrix - - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - Panom(:,:) = 0d0 ! Do sth TODO - ! Dump results + write(*,*)'-----------------------------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1A16,1X,A1,1X,A10,2X,A1,1X)') & + '|','#','|','E(HFB)','|','EJ(HFB)','|','EK(HFB)','|','EL(HFB)','|','Conv','|' + write(*,*)'-----------------------------------------------------------------------------------------------' write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1XF16.10,1X,A1,1X,E10.2,1X,A1,1X)') & '|',nSCF,'|',EHFB + ENuc,'|',EJ,'|',EK,'|',EL,'|',Conv,'|' + write(*,*)'-----------------------------------------------------------------------------------------------' + write(*,*) + ! Check convergence TODO HFB does not fulfill [F,1D]=0 ... + + err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + if(nSCF > 1) Conv = maxval(abs(err)) + end do - write(*,*)'-----------------------------------------------------------------------------------------------' !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ @@ -251,7 +246,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,err,cp,H_hfb,X_hfb,c_hfb,S_hfb,mom,cOld_T,eHFB_,project,err_diis,F_diis) + deallocate(J,K,err,cp,H_hfb,R,eHFB_,err_diis,F_diis) stop @@ -277,6 +272,24 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Memory deallocation - deallocate(J,K,err,cp,H_hfb,X_hfb,c_hfb,S_hfb,mom,cOld_T,eHFB_,project,err_diis,F_diis) + deallocate(J,K,err,cp,H_hfb,R,eHFB_,err_diis,F_diis) end subroutine + + + + ! ! DIIS extrapolation TODO check and adapt + ! + ! 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 TODO check and adapt + ! + ! if(level_shift > 0d0 .and. Conv > thresh) then + ! call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) + ! endif + From c85fa22a98f0b8a589b9182c8f17d6c8c0b6f429 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Wed, 29 Jan 2025 14:12:00 +0100 Subject: [PATCH 08/32] Incl. check conv R. TODO fix chem_pot --- src/HF/HFB.f90 | 154 ++++++++++++++++++++++++++++--------------------- 1 file changed, 87 insertions(+), 67 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 16be4a6f..e0051190 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -32,12 +32,9 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Local variables - integer :: maxSCF_chem_pot integer :: iorb integer :: nSCF - integer :: nSCF_chem_pot integer :: nOrb2 - integer :: nBas2 integer :: nBas_Sq integer :: n_diis double precision :: ET @@ -50,7 +47,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision :: Conv double precision :: rcond - double precision :: delta_chem_pot double precision :: trace_1rdm double precision :: thrs_N double precision,external :: trace_matrix @@ -63,6 +59,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision,allocatable :: cp(:,:) double precision,allocatable :: H_hfb(:,:) double precision,allocatable :: R(:,:) + double precision,allocatable :: P_old(:,:) + double precision,allocatable :: Panom_old(:,:) ! Output variables @@ -85,7 +83,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Useful quantities nBas_Sq = nBas*nBas - nBas2 = nBas+nBas nOrb2 = nOrb+nOrb ! Memory allocation @@ -98,6 +95,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & allocate(cp(nOrb2,nOrb2)) allocate(H_hfb(nOrb2,nOrb2)) allocate(R(nOrb2,nOrb2)) + allocate(P_old(nBas,nBas)) + allocate(Panom_old(nBas,nBas)) allocate(eHFB_(nOrb2)) allocate(err_diis(nBas_Sq,max_diis)) @@ -110,14 +109,16 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Initialization thrs_N = 1d-6 - maxSCF_chem_pot = 1000 n_diis = 0 F_diis(:,:) = 0d0 err_diis(:,:) = 0d0 rcond = 0d0 - P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) - Panom(:,:) = 0d0 ! Do sth TODO + P(:,:) = matmul(c(:,1:nO), transpose(c(:,1:nO))) + Panom(:,:) = 0d0 ! Do sth TODO + P_old(:,:) = P(:,:) + Panom_old(:,:) = Panom(:,:) + P(:,:) = 2d0*P(:,:) Conv = 1d0 nSCF = 0 @@ -125,64 +126,53 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ - - - do while(Conv > thresh .and. nSCF < maxSCF) + + !!do while(Conv > thresh .and. nSCF < maxSCF) +do while(Conv > 1e-7 .and. nSCF < maxSCF) ! TODO ! Increment nSCF = nSCF + 1 - write(*,*) - write(*,*)'-------------------------------------' - write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') & - '|','Tr[1D]','|','Chem. Pot.','|' - write(*,*)'-------------------------------------' - - ! Loop to adjust chem_pot - delta_chem_pot=1d0 - nSCF_chem_pot=0 - do - - ! Build Fock matrix - - call Hartree_matrix_AO_basis(nBas,P,ERI,J) - call exchange_matrix_AO_basis(nBas,P,ERI,K) - call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) - - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) - - ! Diagonalize H_HFB matrix - - H_hfb(:,:) = 0d0 - H_hfb(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X)) - H_hfb(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_hfb(1:nOrb,1:nOrb) - H_hfb(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) - H_hfb(nOrb+1:nOrb2,1:nOrb ) = transpose(H_hfb(1:nOrb,nOrb+1:nOrb2)) - - cp(:,:) = H_hfb(:,:) - call diagonalize_matrix(nOrb2,cp,eHFB_) - - ! Build R and extract P and Panom - - trace_1rdm = 0d0 - R(:,:) = 0d0 - do iorb=1,nOrb - R(:,:) = R(:,:) + matmul(cp(:,iorb:iorb),transpose(cp(:,iorb:iorb))) - enddo - do iorb=1,nOrb - trace_1rdm = trace_1rdm + R(iorb,iorb) - enddo - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_1rdm,'|',chem_pot,'|' - - nSCF_chem_pot = nSCF_chem_pot + 1 - if( abs(trace_1rdm-nO) < thrs_N .or. nSCF_chem_pot > maxSCF_chem_pot ) exit - - - + + ! Build Fock and Delta matrices + + call Hartree_matrix_AO_basis(nBas,P,ERI,J) + call exchange_matrix_AO_basis(nBas,P,ERI,K) + call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) + + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) + + ! Diagonalize H_HFB matrix + + H_hfb(:,:) = 0d0 + H_hfb(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X)) + H_hfb(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_hfb(1:nOrb,1:nOrb) + H_hfb(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) + H_hfb(nOrb+1:nOrb2,1:nOrb ) = transpose(H_hfb(1:nOrb,nOrb+1:nOrb2)) + + cp(:,:) = H_hfb(:,:) + call diagonalize_matrix(nOrb2,cp,eHFB_) + + ! Build R and extract P and Panom + + trace_1rdm = 0d0 + R(:,:) = 0d0 + do iorb=1,nOrb + R(:,:) = R(:,:) + matmul(cp(:,iorb:iorb),transpose(cp(:,iorb:iorb))) + enddo + do iorb=1,nOrb + trace_1rdm = trace_1rdm + R(iorb,iorb) enddo + ! Adjust the chemical potential + + if( abs(trace_1rdm-nO) > thrs_N ) then + +! call fix_chem_pot(nO,nOrb,nOrb2,nBas,thrs_N,chem_pot,F,Hc,J,K,Delta,X,H_hfb,cp,R,eHFB_) + + endif + ! Extract P and Panom from R P(:,:) = 0d0 @@ -190,6 +180,18 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & P(:,:) = 2d0*matmul(X,matmul(R(1:nOrb,1:nOrb),transpose(X))) Panom(:,:) = matmul(X,matmul(R(1:nOrb,nOrb+1:nOrb2),transpose(X))) +block +integer::iorb1 +write(*,*) 'Tr[1D] and chem_pot',trace_1rdm,chem_pot +write(*,*) 'iter',nSCF +do iorb1=1,nOrb2 + write(*,*) eHFB_(iorb1) +enddo +do iorb1=1,nOrb + write(*,'(7f10.5)') P(iorb1,:) +enddo +end block + ! Kinetic energy @@ -213,7 +215,20 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Total energy - EHFB = ET + EV + EJ + EK !+ EL + EHFB = ET + EV + EJ + EK + EL + + ! Check convergence + + if(nSCF > 1) then + + err = P_old - P + Conv = maxval(abs(err)) + err = Panom_old - Panom + Conv = Conv + maxval(abs(err)) + P_old(:,:) = P(:,:) + Panom_old(:,:) = Panom(:,:) + + endif ! Dump results write(*,*)'-----------------------------------------------------------------------------------------------' @@ -226,10 +241,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & write(*,*)'-----------------------------------------------------------------------------------------------' write(*,*) - ! Check convergence TODO HFB does not fulfill [F,1D]=0 ... - - err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) - if(nSCF > 1) Conv = maxval(abs(err)) end do !------------------------------------------------------------------------ @@ -246,7 +257,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,err,cp,H_hfb,R,eHFB_,err_diis,F_diis) + deallocate(J,K,err,cp,H_hfb,R,P_old,Panom_old,eHFB_,err_diis,F_diis) stop @@ -254,7 +265,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Compute dipole moments -eHF(:)=eHF(:)-chem_pot!TODO remove + eHF(1:nOrb)=eHF(1:nOrb)-chem_pot ! TODO fix it 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,chem_pot,dipole) @@ -272,7 +283,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Memory deallocation - deallocate(J,K,err,cp,H_hfb,R,eHFB_,err_diis,F_diis) + deallocate(J,K,err,cp,H_hfb,R,P_old,Panom_old,eHFB_,err_diis,F_diis) end subroutine @@ -293,3 +304,12 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) ! endif +! write(*,*) +! write(*,*)'-------------------------------------' +! write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') & +! '|','Tr[1D]','|','Chem. Pot.','|' +! write(*,*)'-------------------------------------' + +! write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & +! '|',trace_1rdm,'|',chem_pot,'|' + From 94062e66cdb07644b71ed4959547ac34582ab192 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Wed, 29 Jan 2025 15:52:04 +0100 Subject: [PATCH 09/32] Todo golden-section search for mu --- src/HF/HFB.f90 | 36 ++++------ src/HF/fix_chem_pot.f90 | 150 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 164 insertions(+), 22 deletions(-) create mode 100644 src/HF/fix_chem_pot.f90 diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index e0051190..039a572d 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -134,6 +134,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & nSCF = nSCF + 1 +! TODO remove +chem_pot=-10 ! Build Fock and Delta matrices @@ -169,7 +171,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & if( abs(trace_1rdm-nO) > thrs_N ) then -! call fix_chem_pot(nO,nOrb,nOrb2,nBas,thrs_N,chem_pot,F,Hc,J,K,Delta,X,H_hfb,cp,R,eHFB_) + call fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB_) endif @@ -180,18 +182,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & P(:,:) = 2d0*matmul(X,matmul(R(1:nOrb,1:nOrb),transpose(X))) Panom(:,:) = matmul(X,matmul(R(1:nOrb,nOrb+1:nOrb2),transpose(X))) -block -integer::iorb1 -write(*,*) 'Tr[1D] and chem_pot',trace_1rdm,chem_pot -write(*,*) 'iter',nSCF -do iorb1=1,nOrb2 - write(*,*) eHFB_(iorb1) -enddo -do iorb1=1,nOrb - write(*,'(7f10.5)') P(iorb1,:) -enddo -end block - ! Kinetic energy @@ -304,12 +294,14 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) ! endif -! write(*,*) -! write(*,*)'-------------------------------------' -! write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') & -! '|','Tr[1D]','|','Chem. Pot.','|' -! write(*,*)'-------------------------------------' - -! write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & -! '|',trace_1rdm,'|',chem_pot,'|' - +!block +!integer::iorb1 +!write(*,*) 'Tr[1D] and chem_pot',trace_1rdm,chem_pot +!write(*,*) 'iter',nSCF +!do iorb1=1,nOrb2 +! write(*,*) eHFB_(iorb1) +!enddo +!do iorb1=1,nOrb +! write(*,'(7f10.5)') P(iorb1,:) +!enddo +!end block diff --git a/src/HF/fix_chem_pot.f90 b/src/HF/fix_chem_pot.f90 new file mode 100644 index 00000000..b5c86101 --- /dev/null +++ b/src/HF/fix_chem_pot.f90 @@ -0,0 +1,150 @@ +subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB_) + +! Fix the chemical potential to integrate to the N for 2N electron systems + + implicit none + +! Input variables + + integer,intent(in) :: nO,nOrb,nOrb2 + double precision,intent(in) :: thrs_N + +! Local variables + + integer :: iorb + double precision :: trace_curr,trace_down,trace_up + double precision :: chem_pot_curr,chem_pot_down,chem_pot_up + double precision :: delta_chem_pot + +! Output variables + + double precision :: trace_1rdm + double precision,intent(inout):: chem_pot + double precision,intent(inout):: cp(nOrb2,nOrb2) + double precision,intent(inout):: R(nOrb2,nOrb2) + double precision,intent(inout):: eHFB_(nOrb2) + double precision,intent(inout):: H_hfb(nOrb2,nOrb2) + + ! Initialize delta_chem_pot + + delta_chem_pot = 1d0 + trace_down = 0d0 + trace_up = 0d0 + chem_pot_down = 0d0 + chem_pot_up = 0d0 + + ! Set H_HFB to its non-chemical potential dependent contribution + + do iorb=1,nOrb + H_hfb(iorb,iorb)=H_hfb(iorb,iorb)+chem_pot + H_hfb(iorb+nOrb,iorb+nOrb)=H_hfb(iorb+nOrb,iorb+nOrb)-chem_pot + enddo + + write(*,*) + write(*,*)'-------------------------------------' + write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') & + '|','Tr[1D]','|','Chem. Pot.','|' + write(*,*)'-------------------------------------' + + ! Set interval to search + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_curr,H_hfb,cp,R,eHFB_) + chem_pot_curr=chem_pot + if(trace_curr=nO+1) exit + enddo + chem_pot_up=chem_pot + else + ! Decrease chem_pot to occupy less orbs. + do + chem_pot = chem_pot - delta_chem_pot + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_down,H_hfb,cp,R,eHFB_) + if(trace_down<=nO-1) exit + enddo + chem_pot_down=chem_pot + endif + + if(abs(chem_pot_up)>1e-4) then + + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_curr,'|',chem_pot_curr,'|' + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_up ,'|',chem_pot_up ,'|' + + + write(*,*)'-------------------------------------' + write(*,*) + + endif + + if(abs(chem_pot_down)>1e-4) then + + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_curr,'|',chem_pot_curr,'|' + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_down,'|',chem_pot_down,'|' + + + write(*,*)'-------------------------------------' + write(*,*) + + endif + + ! Reset H_HFB to its chemical potential version + + do iorb=1,nOrb + H_hfb(iorb,iorb)=H_hfb(iorb,iorb)-chem_pot + H_hfb(iorb+nOrb,iorb+nOrb)=H_hfb(iorb+nOrb,iorb+nOrb)+chem_pot + enddo + +end subroutine + +subroutine diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) + +! Fix the chemical potential to integrate to the N for 2N electron systems + + implicit none + +! Input variables + + integer,intent(in) :: nOrb,nOrb2 + double precision,intent(in) :: H_hfb(nOrb2,nOrb2) + double precision,intent(in) :: chem_pot + +! Local variables + + integer :: iorb + +! Output variables + + double precision,intent(out) :: trace_1rdm + double precision,intent(inout):: cp(nOrb2,nOrb2) + double precision,intent(inout):: R(nOrb2,nOrb2) + double precision,intent(inout):: eHFB_(nOrb2) + + cp(:,:) = H_hfb(:,:) + do iorb=1,nOrb + cp(iorb,iorb) = cp(iorb,iorb) - chem_pot + cp(iorb+nOrb,iorb+nOrb) = cp(iorb+nOrb,iorb+nOrb) + chem_pot + enddo + + ! Diagonalize H_HFB matrix + + call diagonalize_matrix(nOrb2,cp,eHFB_) + + ! Build R and extract P and Panom + + trace_1rdm = 0d0 + R(:,:) = 0d0 + do iorb=1,nOrb + R(:,:) = R(:,:) + matmul(cp(:,iorb:iorb),transpose(cp(:,iorb:iorb))) + enddo + do iorb=1,nOrb + trace_1rdm = trace_1rdm + R(iorb,iorb) + enddo + +end subroutine + From e105665f8e6a878dd6d07f8ab48eb75482441f47 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Wed, 29 Jan 2025 18:04:21 +0100 Subject: [PATCH 10/32] Apparently golden-search works HF --- src/HF/HFB.f90 | 6 +-- src/HF/fix_chem_pot.f90 | 83 +++++++++++++++++++++++++++++++---------- 2 files changed, 64 insertions(+), 25 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 039a572d..9e8bee34 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -127,16 +127,12 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Main SCF loop !------------------------------------------------------------------------ - !!do while(Conv > thresh .and. nSCF < maxSCF) -do while(Conv > 1e-7 .and. nSCF < maxSCF) ! TODO + do while(Conv > thresh .and. nSCF < maxSCF) ! Increment nSCF = nSCF + 1 -! TODO remove -chem_pot=-10 - ! Build Fock and Delta matrices call Hartree_matrix_AO_basis(nBas,P,ERI,J) diff --git a/src/HF/fix_chem_pot.f90 b/src/HF/fix_chem_pot.f90 index b5c86101..0f4f7559 100644 --- a/src/HF/fix_chem_pot.f90 +++ b/src/HF/fix_chem_pot.f90 @@ -11,10 +11,15 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB ! Local variables + logical :: is_up integer :: iorb + integer :: isteps double precision :: trace_curr,trace_down,trace_up double precision :: chem_pot_curr,chem_pot_down,chem_pot_up double precision :: delta_chem_pot + double precision :: golden_ratio + double precision :: a + double precision :: c ! Output variables @@ -27,7 +32,10 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB ! Initialize delta_chem_pot - delta_chem_pot = 1d0 + is_up=.false. + isteps=0 + golden_ratio = 1.618033988 + delta_chem_pot = 5d-1 trace_down = 0d0 trace_up = 0d0 chem_pot_down = 0d0 @@ -46,52 +54,87 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB '|','Tr[1D]','|','Chem. Pot.','|' write(*,*)'-------------------------------------' - ! Set interval to search + ! Set interval to search (the good chem_pot is in [chem_pot_down,chem_pot_up]) call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_curr,H_hfb,cp,R,eHFB_) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_curr,'|',chem_pot,'|' chem_pot_curr=chem_pot if(trace_curr=nO+1) exit + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_up,'|',chem_pot,'|' + if(trace_up>=nO+1 .or. isteps>100) exit ! max 50 au steps for mu (is a lot) enddo chem_pot_up=chem_pot else ! Decrease chem_pot to occupy less orbs. do + isteps = isteps+1 chem_pot = chem_pot - delta_chem_pot call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_down,H_hfb,cp,R,eHFB_) - if(trace_down<=nO-1) exit + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_down,'|',chem_pot,'|' + if(trace_down<=nO-1 .or. isteps>100) exit ! max 50 au steps for mu (is a lot) enddo chem_pot_down=chem_pot endif - - if(abs(chem_pot_up)>1e-4) then - + if(is_up) then write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & '|',trace_curr,'|',chem_pot_curr,'|' write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & '|',trace_up ,'|',chem_pot_up ,'|' - - - write(*,*)'-------------------------------------' - write(*,*) - - endif - - if(abs(chem_pot_down)>1e-4) then - + trace_down=trace_curr + chem_pot_down=chem_pot_curr + else write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & '|',trace_curr,'|',chem_pot_curr,'|' write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & '|',trace_down,'|',chem_pot_down,'|' + trace_up=trace_curr + chem_pot_up=chem_pot_curr + endif + ! Use Golden-section search algorithm to find chem_pot + isteps=0 + do + isteps = isteps+1 + a=(chem_pot_up-chem_pot_down)/(1d0+golden_ratio) + c=a/golden_ratio + chem_pot_curr=chem_pot_down+a + chem_pot=chem_pot_curr+c + call diag_H_hfb(nOrb,nOrb2,chem_pot_curr,trace_curr,H_hfb,cp,R,eHFB_) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_curr,'|',chem_pot_curr,'|' + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|' + if(abs(trace_1rdm-nO)1000) exit ! 1000 steps for finding mu + if(is_up) then + if(trace_1rdm>=trace_curr .or. abs(trace_1rdm-trace_curr)<1d-8 ) then + chem_pot_down=chem_pot_curr + trace_down=trace_curr + else + chem_pot_up=chem_pot + trace_up=trace_1rdm + endif + else + if(trace_1rdm>=trace_curr .or. abs(trace_1rdm-trace_curr)<1d-8 ) then + chem_pot_up=chem_pot + trace_up=trace_1rdm + else + chem_pot_down=chem_pot_curr + trace_down=trace_curr + endif + endif + enddo - write(*,*)'-------------------------------------' - write(*,*) - - endif + write(*,*)'-------------------------------------' + write(*,*) ! Reset H_HFB to its chemical potential version From d0bdf0661e88dacc4bd9122bea78e5cdbac39fc2 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Wed, 29 Jan 2025 23:03:54 +0100 Subject: [PATCH 11/32] TODO FD occupancies and chem pot --- src/HF/HFB.f90 | 44 ++++++++++++++++++++++---------------- src/QuAcK/BQuAcK.f90 | 5 +++-- src/QuAcK/QuAcK.f90 | 7 ++++-- src/QuAcK/read_options.f90 | 11 +++++++++- 4 files changed, 44 insertions(+), 23 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 9e8bee34..9c8c1622 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -1,5 +1,6 @@ -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,Delta) +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,Delta, & + temperature) ! Perform Hartree-Fock Bogoliubov calculation @@ -22,6 +23,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc + double precision,intent(in) :: temperature double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) @@ -51,6 +53,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & double precision :: thrs_N double precision,external :: trace_matrix double precision,allocatable :: eHFB_(:) + double precision,allocatable :: Occ(:) double precision,allocatable :: err(:,:) double precision,allocatable :: err_diis(:,:) double precision,allocatable :: F_diis(:,:) @@ -102,7 +105,7 @@ 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 chem. pot +! Guess chem. pot. chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1)) @@ -115,7 +118,26 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & rcond = 0d0 P(:,:) = matmul(c(:,1:nO), transpose(c(:,1:nO))) - Panom(:,:) = 0d0 ! Do sth TODO + Panom(:,:) = 0d0 + + ! Use Fermi-Dirac occupancies to compute P, Panom, and chem_pot + + if(abs(temperature)>1e-4) then ! TODO + allocate(Occ(nOrb)) + Occ(:) = 0d0 + Occ(1:nO) = 1d0 +! call fermi_dirac_occ(nOrb,chem_pot,Occ,eHF) + P(:,:) = 0d0 + Panom(:,:) = 0d0 + do iorb=1,nOrb + P(:,:) = P(:,:) + Occ(iorb) * & + matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) + Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & + matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) + enddo + deallocate(Occ) + endif + P_old(:,:) = P(:,:) Panom_old(:,:) = Panom(:,:) P(:,:) = 2d0*P(:,:) @@ -166,9 +188,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! Adjust the chemical potential if( abs(trace_1rdm-nO) > thrs_N ) then - call fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB_) - endif ! Extract P and Panom from R @@ -289,15 +309,3 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, & ! if(level_shift > 0d0 .and. Conv > thresh) then ! call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) ! endif - -!block -!integer::iorb1 -!write(*,*) 'Tr[1D] and chem_pot',trace_1rdm,chem_pot -!write(*,*) 'iter',nSCF -!do iorb1=1,nOrb2 -! write(*,*) eHFB_(iorb1) -!enddo -!do iorb1=1,nOrb -! write(*,'(7f10.5)') P(iorb1,:) -!enddo -!end block diff --git a/src/QuAcK/BQuAcK.f90 b/src/QuAcK/BQuAcK.f90 index 01ff9dda..836ee9ac 100644 --- a/src/QuAcK/BQuAcK.f90 +++ b/src/QuAcK/BQuAcK.f90 @@ -1,6 +1,6 @@ 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) + guess_type,mix,temperature) ! Restricted branch of QuAcK @@ -19,6 +19,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, integer,intent(in) :: nV integer,intent(in) :: nR double precision,intent(in) :: ENuc + double precision,intent(in) :: temperature double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart) @@ -93,7 +94,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, 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,Delta) + FHF,Delta,temperature) call wall_time(end_HF) t_HF = end_HF - start_HF diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 0342309a..22d607b9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -73,6 +73,8 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest + double precision :: temperature + character(len=256) :: working_dir ! Check if the right number of arguments is provided @@ -134,7 +136,7 @@ program QuAcK maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, & maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,doppBSE,dBSE,dTDA) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,temperature) !------------------! ! Hardware ! @@ -290,7 +292,8 @@ program QuAcK !--------------------------! if(doBQuAcK) & 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) + S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix, & + temperature) !-----------! ! Stop Test ! diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 8ba40339..11af97ea 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -7,7 +7,7 @@ subroutine read_options(working_dir, maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, & maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,doppBSE,dBSE,dTDA) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA,temperature) ! Read desired methods @@ -72,6 +72,8 @@ subroutine read_options(working_dir, logical,intent(out) :: dBSE logical,intent(out) :: dTDA + double precision,intent(out) :: temperature + ! Local variables character(len=1) :: ans1,ans2,ans3,ans4,ans5 @@ -217,6 +219,13 @@ subroutine read_options(working_dir, if(ans4 == 'T') dBSE = .true. if(ans5 == 'F') dTDA = .false. + ! Options for dynamical Fermi-Dirac occupancies + + temperature = 0d0 + + read(1,*) + read(1,*) temperature + endif ! Close file with options From aa7f5e079a9dcd58a86c271f4ec1a6243be5b706 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 30 Jan 2025 13:45:05 +0100 Subject: [PATCH 12/32] Incl FD distrib occ --- src/HF/HFB.f90 | 12 ++++-- src/HF/fermi_dirac_occ.f90 | 85 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+), 4 deletions(-) create mode 100644 src/HF/fermi_dirac_occ.f90 diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 9c8c1622..c35bda5b 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -126,7 +126,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, allocate(Occ(nOrb)) Occ(:) = 0d0 Occ(1:nO) = 1d0 -! call fermi_dirac_occ(nOrb,chem_pot,Occ,eHF) + call fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) P(:,:) = 0d0 Panom(:,:) = 0d0 do iorb=1,nOrb @@ -148,8 +148,12 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, !------------------------------------------------------------------------ ! Main SCF loop !------------------------------------------------------------------------ - - do while(Conv > thresh .and. nSCF < maxSCF) + + write(*,*) + write(*,*) 'Enterning HFB SCF procedure' + write(*,*) + do while(Conv > thresh .and. nSCF < 2) + !do while(Conv > thresh .and. nSCF < maxSCF) ! Increment @@ -217,7 +221,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Anomalous energy - EL = -0.25d0*trace_matrix(nBas,matmul(Panom,Delta)) + EL = 0.25d0*trace_matrix(nBas,matmul(Panom,Delta)) ! Total energy diff --git a/src/HF/fermi_dirac_occ.f90 b/src/HF/fermi_dirac_occ.f90 new file mode 100644 index 00000000..e34ccd68 --- /dev/null +++ b/src/HF/fermi_dirac_occ.f90 @@ -0,0 +1,85 @@ +subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) + +! Use Fermi Dirac distribution to set up fractional Occs numbers and adjust the chemical potential to integrate to the N for 2N electron systems + + implicit none + +! Input variables + + integer,intent(in) :: nO,nOrb + double precision,intent(in) :: thrs_N + double precision,intent(in) :: temperature + +! Local variables + + integer :: iorb + integer :: isteps + double precision :: delta_chem_pot + double precision :: chem_pot_change + double precision :: grad_electrons + double precision :: trace_1rdm + +! Output variables + + double precision,intent(inout):: chem_pot + double precision,intent(inout):: eHF(nOrb) + double precision,intent(inout):: Occ(nOrb) + + ! Initialize variables + + isteps=0 + delta_chem_pot = 1.0d-3 + trace_1rdm = -1.0d0 + chem_pot_change = 0.0d0 + + write(*,*) + write(*,*)' Fermi-Dirac distribution for the occ numbers' + write(*,*) + write(*,*)'-------------------------------------' + write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') & + '|','Tr[1D]','|','Chem. Pot.','|' + write(*,*)'-------------------------------------' + + + Occ(:) = fermi_dirac(eHF,chem_pot,temperature) + trace_1rdm=sum(Occ(:)) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|' + + + do while( abs(trace_1rdm-nO) > 1.0d-8 .and. isteps <= 100 ) + isteps = isteps + 1 + chem_pot = chem_pot + chem_pot_change + Occ(:) = fermi_dirac(eHF,chem_pot,temperature) + trace_1rdm = sum(Occ(:)) + grad_electrons = ( sum(fermi_dirac(eHF,chem_pot+delta_chem_pot,temperature)) & + - sum(fermi_dirac(eHF,chem_pot-delta_chem_pot,temperature)) )/(2.0d0*delta_chem_pot) + chem_pot_change = -(trace_1rdm-nO)/grad_electrons + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|' + enddo + write(*,*)'-------------------------------------' + write(*,*) + + write(*,*) + write(*,*) ' Initial occ. numbers ' + write(*,*) + do iorb=1,nOrb + write(*,'(3X,F16.10)') Occ(iorb) + enddo + +contains + +function fermi_dirac(eHF,chem_pot,temperature) + implicit none + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: chem_pot + double precision,intent(in) :: temperature + double precision :: fermi_dirac(nOrb) + + fermi_dirac(:) = 1d0 / ( 1d0 + exp((eHF(:) - chem_pot ) / temperature ) ) + +end function fermi_dirac + +end subroutine + From 708964f277681acfe61f4af46038f7c925c1e5d7 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 30 Jan 2025 15:17:29 +0100 Subject: [PATCH 13/32] Works for HFB + Delta --- src/HF/HFB.f90 | 7 +- src/HF/fermi_dirac_occ.f90 | 4 +- src/HF/fix_chem_pot.f90 | 150 +++++++++++++------------------------ 3 files changed, 56 insertions(+), 105 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index c35bda5b..b6b22646 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -111,7 +111,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Initialization - thrs_N = 1d-6 + thrs_N = 1d-8 n_diis = 0 F_diis(:,:) = 0d0 err_diis(:,:) = 0d0 @@ -122,7 +122,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Use Fermi-Dirac occupancies to compute P, Panom, and chem_pot - if(abs(temperature)>1e-4) then ! TODO + if(abs(temperature)>1e-4) then allocate(Occ(nOrb)) Occ(:) = 0d0 Occ(1:nO) = 1d0 @@ -152,8 +152,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*) write(*,*) 'Enterning HFB SCF procedure' write(*,*) - do while(Conv > thresh .and. nSCF < 2) - !do while(Conv > thresh .and. nSCF < maxSCF) + do while(Conv > thresh .and. nSCF < maxSCF) ! Increment diff --git a/src/HF/fermi_dirac_occ.f90 b/src/HF/fermi_dirac_occ.f90 index e34ccd68..e48ce874 100644 --- a/src/HF/fermi_dirac_occ.f90 +++ b/src/HF/fermi_dirac_occ.f90 @@ -47,14 +47,14 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) '|',trace_1rdm,'|',chem_pot,'|' - do while( abs(trace_1rdm-nO) > 1.0d-8 .and. isteps <= 100 ) + do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 ) isteps = isteps + 1 chem_pot = chem_pot + chem_pot_change Occ(:) = fermi_dirac(eHF,chem_pot,temperature) trace_1rdm = sum(Occ(:)) grad_electrons = ( sum(fermi_dirac(eHF,chem_pot+delta_chem_pot,temperature)) & - sum(fermi_dirac(eHF,chem_pot-delta_chem_pot,temperature)) )/(2.0d0*delta_chem_pot) - chem_pot_change = -(trace_1rdm-nO)/grad_electrons + chem_pot_change = -(trace_1rdm-nO)/(grad_electrons+1d-10) write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & '|',trace_1rdm,'|',chem_pot,'|' enddo diff --git a/src/HF/fix_chem_pot.f90 b/src/HF/fix_chem_pot.f90 index 0f4f7559..3847e947 100644 --- a/src/HF/fix_chem_pot.f90 +++ b/src/HF/fix_chem_pot.f90 @@ -11,35 +11,33 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB ! Local variables - logical :: is_up + logical :: backward integer :: iorb integer :: isteps - double precision :: trace_curr,trace_down,trace_up - double precision :: chem_pot_curr,chem_pot_down,chem_pot_up double precision :: delta_chem_pot - double precision :: golden_ratio - double precision :: a - double precision :: c + double precision :: chem_pot_change + double precision :: grad_electrons ! Output variables double precision :: trace_1rdm + double precision :: trace_up + double precision :: trace_down + double precision :: trace_old double precision,intent(inout):: chem_pot double precision,intent(inout):: cp(nOrb2,nOrb2) double precision,intent(inout):: R(nOrb2,nOrb2) double precision,intent(inout):: eHFB_(nOrb2) double precision,intent(inout):: H_hfb(nOrb2,nOrb2) - ! Initialize delta_chem_pot + ! Initialize - is_up=.false. - isteps=0 - golden_ratio = 1.618033988 - delta_chem_pot = 5d-1 - trace_down = 0d0 - trace_up = 0d0 - chem_pot_down = 0d0 - chem_pot_up = 0d0 + backward = .false. + isteps = 0 + delta_chem_pot = 1.0d-1 + chem_pot_change = 0d0 + grad_electrons = 1d0 + trace_1rdm = -1d0 ! Set H_HFB to its non-chemical potential dependent contribution @@ -49,93 +47,47 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB enddo write(*,*) - write(*,*)'-------------------------------------' - write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') & - '|','Tr[1D]','|','Chem. Pot.','|' - write(*,*)'-------------------------------------' - - ! Set interval to search (the good chem_pot is in [chem_pot_down,chem_pot_up]) - call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_curr,H_hfb,cp,R,eHFB_) - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_curr,'|',chem_pot,'|' - chem_pot_curr=chem_pot - if(trace_curr=nO+1 .or. isteps>100) exit ! max 50 au steps for mu (is a lot) - enddo - chem_pot_up=chem_pot - else - ! Decrease chem_pot to occupy less orbs. - do - isteps = isteps+1 - chem_pot = chem_pot - delta_chem_pot - call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_down,H_hfb,cp,R,eHFB_) - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_down,'|',chem_pot,'|' - if(trace_down<=nO-1 .or. isteps>100) exit ! max 50 au steps for mu (is a lot) - enddo - chem_pot_down=chem_pot - endif - if(is_up) then - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_curr,'|',chem_pot_curr,'|' - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_up ,'|',chem_pot_up ,'|' - trace_down=trace_curr - chem_pot_down=chem_pot_curr - else - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_curr,'|',chem_pot_curr,'|' - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_down,'|',chem_pot_down,'|' - trace_up=trace_curr - chem_pot_up=chem_pot_curr - endif - - ! Use Golden-section search algorithm to find chem_pot - isteps=0 - do - isteps = isteps+1 - a=(chem_pot_up-chem_pot_down)/(1d0+golden_ratio) - c=a/golden_ratio - chem_pot_curr=chem_pot_down+a - chem_pot=chem_pot_curr+c - call diag_H_hfb(nOrb,nOrb2,chem_pot_curr,trace_curr,H_hfb,cp,R,eHFB_) - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_curr,'|',chem_pot_curr,'|' - call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_1rdm,'|',chem_pot,'|' - if(abs(trace_1rdm-nO)1000) exit ! 1000 steps for finding mu - if(is_up) then - if(trace_1rdm>=trace_curr .or. abs(trace_1rdm-trace_curr)<1d-8 ) then - chem_pot_down=chem_pot_curr - trace_down=trace_curr - else - chem_pot_up=chem_pot - trace_up=trace_1rdm - endif - else - if(trace_1rdm>=trace_curr .or. abs(trace_1rdm-trace_curr)<1d-8 ) then - chem_pot_up=chem_pot - trace_up=trace_1rdm - else - chem_pot_down=chem_pot_curr - trace_down=trace_curr - endif - endif + write(*,*)'------------------------------------------------------' + write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1A15,2X,A1)') & + '|','Tr[1D]','|','Chem. Pot.','|','Grad N','|' + write(*,*)'------------------------------------------------------' + + ! First approach close the value with an error lower than 1 + + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_old,H_hfb,cp,R,eHFB_) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') & + '|',trace_old,'|',chem_pot,'|',grad_electrons,'|' + do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 ) + isteps = isteps + 1 + chem_pot = chem_pot + delta_chem_pot + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|' + if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then + backward=.true. + chem_pot = chem_pot - 2d0*delta_chem_pot + delta_chem_pot=-delta_chem_pot + endif enddo - write(*,*)'-------------------------------------' + ! Do final search + + isteps = 0 + delta_chem_pot = 1.0d-3 + do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 ) + isteps = isteps + 1 + chem_pot = chem_pot + chem_pot_change + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) + call diag_H_hfb(nOrb,nOrb2,chem_pot+delta_chem_pot,trace_up,H_hfb,cp,R,eHFB_) + call diag_H_hfb(nOrb,nOrb2,chem_pot-delta_chem_pot,trace_down,H_hfb,cp,R,eHFB_) + grad_electrons = (trace_up-trace_down)/(2.0d0*delta_chem_pot) + chem_pot_change = -(trace_1rdm-nO)/(grad_electrons+1d-10) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|' + enddo + write(*,*)'------------------------------------------------------' write(*,*) - + ! Reset H_HFB to its chemical potential version do iorb=1,nOrb From b69451ba382fd59a80ff3bb0832453e3d1468cbf Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 30 Jan 2025 15:51:59 +0100 Subject: [PATCH 14/32] Corrected energy factor Eanom --- src/HF/HFB.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index b6b22646..2c13369a 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -165,7 +165,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) - + ! Diagonalize H_HFB matrix H_hfb(:,:) = 0d0 @@ -220,7 +220,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Anomalous energy - EL = 0.25d0*trace_matrix(nBas,matmul(Panom,Delta)) + EL = trace_matrix(nBas,matmul(Panom,Delta)) ! Total energy From 2744c97f2ee472075f2b528b22360598fe2b62e6 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 30 Jan 2025 16:05:29 +0100 Subject: [PATCH 15/32] TODO print output info add DIIS --- src/HF/fermi_dirac_occ.f90 | 37 +++++++++++++++++++++++++++++++------ src/HF/fix_chem_pot.f90 | 7 ++++--- 2 files changed, 35 insertions(+), 9 deletions(-) diff --git a/src/HF/fermi_dirac_occ.f90 b/src/HF/fermi_dirac_occ.f90 index e48ce874..7cdfc213 100644 --- a/src/HF/fermi_dirac_occ.f90 +++ b/src/HF/fermi_dirac_occ.f90 @@ -12,12 +12,14 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) ! Local variables + logical :: backward integer :: iorb integer :: isteps double precision :: delta_chem_pot double precision :: chem_pot_change double precision :: grad_electrons double precision :: trace_1rdm + double precision :: trace_old ! Output variables @@ -27,10 +29,12 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) ! Initialize variables - isteps=0 - delta_chem_pot = 1.0d-3 - trace_1rdm = -1.0d0 - chem_pot_change = 0.0d0 + backward = .false. + isteps = 0 + delta_chem_pot = 1.0d-1 + chem_pot_change = 0d0 + grad_electrons = 1d0 + trace_1rdm = -1d0 write(*,*) write(*,*)' Fermi-Dirac distribution for the occ numbers' @@ -40,13 +44,34 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) '|','Tr[1D]','|','Chem. Pot.','|' write(*,*)'-------------------------------------' + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|' + + ! First approach close the value with an error lower than 1 Occ(:) = fermi_dirac(eHF,chem_pot,temperature) - trace_1rdm=sum(Occ(:)) + trace_old=sum(Occ(:)) write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_1rdm,'|',chem_pot,'|' + '|',trace_old,'|',chem_pot,'|' + do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 ) + isteps = isteps + 1 + chem_pot = chem_pot + delta_chem_pot + Occ(:) = fermi_dirac(eHF,chem_pot,temperature) + trace_1rdm=sum(Occ(:)) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|' + if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then + backward=.true. + chem_pot = chem_pot - 2d0*delta_chem_pot + delta_chem_pot=-delta_chem_pot + endif + enddo + ! Do final search + write(*,*)'-------------------------------------' + isteps=0 + delta_chem_pot = 1.0d-1 do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 ) isteps = isteps + 1 chem_pot = chem_pot + chem_pot_change diff --git a/src/HF/fix_chem_pot.f90 b/src/HF/fix_chem_pot.f90 index 3847e947..a781d7a1 100644 --- a/src/HF/fix_chem_pot.f90 +++ b/src/HF/fix_chem_pot.f90 @@ -17,13 +17,13 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB double precision :: delta_chem_pot double precision :: chem_pot_change double precision :: grad_electrons + double precision :: trace_up + double precision :: trace_down + double precision :: trace_old ! Output variables double precision :: trace_1rdm - double precision :: trace_up - double precision :: trace_down - double precision :: trace_old double precision,intent(inout):: chem_pot double precision,intent(inout):: cp(nOrb2,nOrb2) double precision,intent(inout):: R(nOrb2,nOrb2) @@ -72,6 +72,7 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB ! Do final search + write(*,*)'------------------------------------------------------' isteps = 0 delta_chem_pot = 1.0d-3 do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 ) From 8dbcd68920c9ef87c5bd1e65f1e6948dd34308c0 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 30 Jan 2025 16:54:12 +0100 Subject: [PATCH 16/32] Remove unused MOM --- src/HF/MOM_overlap.f90 | 52 ------------------------------------------ 1 file changed, 52 deletions(-) delete mode 100644 src/HF/MOM_overlap.f90 diff --git a/src/HF/MOM_overlap.f90 b/src/HF/MOM_overlap.f90 deleted file mode 100644 index 4eb66ecc..00000000 --- a/src/HF/MOM_overlap.f90 +++ /dev/null @@ -1,52 +0,0 @@ -subroutine MOM_overlap(nBas,nO,S,cG,c,ON) - -! Compute overlap between old and new MO coefficients - - implicit none - -! Input variables - - integer,intent(in) :: nBas,nO - double precision,intent(in) :: S(nBas,nBas),cG(nBas,nBas),c(nBas,nBas) - -! Local variables - - integer :: i,j,ploc - double precision,allocatable :: Ov(:,:),pOv(:) - -! Output variables - - double precision,intent(inout):: ON(nBas) - - allocate(Ov(nBas,nBas),pOv(nBas)) - - Ov = matmul(transpose(cG),matmul(S,c)) - - pOv(:) = 0d0 - - do i=1,nBas - do j=1,nBas - pOv(j) = pOv(j) + Ov(i,j)**2 - end do - end do - - pOv(:) = sqrt(pOV(:)) - -! print*,'--- MOM overlap ---' -! call matout(nBas,1,pOv) - - ON(:) = 0d0 - - do i=1,nO - ploc = maxloc(pOv,nBas) - ON(ploc) = 1d0 - pOv(ploc) = 0d0 - end do - - -! print*,'--- Occupation numbers ---' -! call matout(nBas,1,ON) - - - -end subroutine From c1a10548d5149fcf2eb4a7fa97aa9db17ed0cb26 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 30 Jan 2025 17:14:51 +0100 Subject: [PATCH 17/32] Incl. printing function, TODO: DIIS --- src/HF/HFB.f90 | 30 ++++++++++++++++-------------- src/HF/print_HFB.f90 | 29 ++++------------------------- 2 files changed, 20 insertions(+), 39 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 2c13369a..02c0ee11 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -52,14 +52,14 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision :: trace_1rdm double precision :: thrs_N double precision,external :: trace_matrix - double precision,allocatable :: eHFB_(:) + double precision,allocatable :: eigVAL(:) double precision,allocatable :: Occ(:) 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 :: eigVEC(:,:) double precision,allocatable :: H_hfb(:,:) double precision,allocatable :: R(:,:) double precision,allocatable :: P_old(:,:) @@ -95,12 +95,12 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, allocate(err(nBas,nBas)) - allocate(cp(nOrb2,nOrb2)) + allocate(eigVEC(nOrb2,nOrb2)) allocate(H_hfb(nOrb2,nOrb2)) allocate(R(nOrb2,nOrb2)) allocate(P_old(nBas,nBas)) allocate(Panom_old(nBas,nBas)) - allocate(eHFB_(nOrb2)) + allocate(eigVAL(nOrb2)) allocate(err_diis(nBas_Sq,max_diis)) allocate(F_diis(nBas_Sq,max_diis)) @@ -174,15 +174,15 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, H_hfb(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) H_hfb(nOrb+1:nOrb2,1:nOrb ) = transpose(H_hfb(1:nOrb,nOrb+1:nOrb2)) - cp(:,:) = H_hfb(:,:) - call diagonalize_matrix(nOrb2,cp,eHFB_) + eigVEC(:,:) = H_hfb(:,:) + call diagonalize_matrix(nOrb2,eigVEC,eigVAL) ! Build R and extract P and Panom trace_1rdm = 0d0 R(:,:) = 0d0 do iorb=1,nOrb - R(:,:) = R(:,:) + matmul(cp(:,iorb:iorb),transpose(cp(:,iorb:iorb))) + R(:,:) = R(:,:) + matmul(eigVEC(:,iorb:iorb),transpose(eigVEC(:,iorb:iorb))) enddo do iorb=1,nOrb trace_1rdm = trace_1rdm + R(iorb,iorb) @@ -191,7 +191,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Adjust the chemical potential if( abs(trace_1rdm-nO) > thrs_N ) then - call fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB_) + call fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,eigVEC,R,eigVAL) endif ! Extract P and Panom from R @@ -266,18 +266,20 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,err,cp,H_hfb,R,P_old,Panom_old,eHFB_,err_diis,F_diis) + deallocate(J,K,err,eigVEC,H_hfb,R,P_old,Panom_old,eigVAL,err_diis,F_diis) stop end if -! Compute dipole moments - - eHF(1:nOrb)=eHF(1:nOrb)-chem_pot ! TODO fix it +! Compute dipole moments and occupation numbers + eigVEC(:,:) = 0d0 + eigVEC(:,:) = R(1:nOrb,1:nOrb) + call diagonalize_matrix(nOrb2,eigVEC,eigVAL) + eigVAL(:) = 2d0*eigVAL(:) 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,chem_pot,dipole) + call print_HFB(nBas,nOrb,nO,eigVAL,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole) ! Testing zone @@ -292,7 +294,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Memory deallocation - deallocate(J,K,err,cp,H_hfb,R,P_old,Panom_old,eHFB_,err_diis,F_diis) + deallocate(J,K,err,eigVEC,H_hfb,R,P_old,Panom_old,eigVAL,err_diis,F_diis) end subroutine diff --git a/src/HF/print_HFB.f90 b/src/HF/print_HFB.f90 index ef25a8ed..903c7718 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, chem_pot, dipole) +subroutine print_HFB(nBas, nOrb, nO, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_pot, dipole) ! Print one-electron energies and other stuff for G0W0 @@ -12,8 +12,7 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, c 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) :: Occ(2*nOrb) double precision,intent(in) :: ENuc double precision,intent(in) :: ET double precision,intent(in) :: EV @@ -27,18 +26,9 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, c ! Local variables integer :: ixyz - integer :: HOMO - integer :: LUMO - double precision :: Gap logical :: dump_orb = .false. -! HOMO and LUMO - - HOMO = nO - LUMO = HOMO + 1 - Gap = eHF(LUMO)-eHF(HOMO) - ! Dump results write(*,*) @@ -58,10 +48,6 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, c 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.10,A3)') ' Chemical potential = ',chem_pot,' au' write(*,'(A50)') '---------------------------------------' write(*,'(A36)') ' Dipole moment (Debye) ' @@ -72,17 +58,10 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, c ! 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)') ' HFB occupation numbers ' write(*,'(A50)') '---------------------------------------' - call vecout(nOrb, eHF) + call vecout(2*nOrb, Occ) write(*,*) end subroutine From fa0981508e61fa4bc6bb53c87c93bcfa3dc10d5f Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 30 Jan 2025 22:27:59 +0100 Subject: [PATCH 18/32] Solved 1 bug chem pot fix --- src/HF/HFB.f90 | 19 +++++++++++++------ src/HF/fix_chem_pot.f90 | 10 ++++++++-- src/HF/print_HFB.f90 | 16 ++++++++++++++-- 3 files changed, 35 insertions(+), 10 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 02c0ee11..c11d785e 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -51,6 +51,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision :: rcond double precision :: trace_1rdm double precision :: thrs_N + double precision :: norm_anom double precision,external :: trace_matrix double precision,allocatable :: eigVAL(:) double precision,allocatable :: Occ(:) @@ -122,7 +123,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Use Fermi-Dirac occupancies to compute P, Panom, and chem_pot - if(abs(temperature)>1e-4) then + if(abs(temperature)>1d-4) then allocate(Occ(nOrb)) Occ(:) = 0d0 Occ(1:nO) = 1d0 @@ -190,9 +191,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Adjust the chemical potential - if( abs(trace_1rdm-nO) > thrs_N ) then + if( abs(trace_1rdm-nO) > thrs_N ) & call fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,eigVEC,R,eigVAL) - endif ! Extract P and Panom from R @@ -272,14 +272,21 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, end if -! Compute dipole moments and occupation numbers +! Compute dipole moments, occupation numbers and || Anomalous density|| eigVEC(:,:) = 0d0 - eigVEC(:,:) = R(1:nOrb,1:nOrb) + eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb) call diagonalize_matrix(nOrb2,eigVEC,eigVAL) eigVAL(:) = 2d0*eigVAL(:) + Panom_old(:,:) = 0d0 + Panom_old(1:nOrb,1:nOrb) = R(1:nOrb,nOrb+1:nOrb2) + Panom_old = matmul(transpose(Panom_old),Panom_old) + norm_anom = 0d0 + do iorb=1,nBas + norm_anom = norm_anom + Panom_old(iorb,iorb) + enddo call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_HFB(nBas,nOrb,nO,eigVAL,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole) + call print_HFB(nBas,nOrb,nO,norm_anom,eigVAL,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole) ! Testing zone diff --git a/src/HF/fix_chem_pot.f90 b/src/HF/fix_chem_pot.f90 index a781d7a1..baa37783 100644 --- a/src/HF/fix_chem_pot.f90 +++ b/src/HF/fix_chem_pot.f90 @@ -20,6 +20,8 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB double precision :: trace_up double precision :: trace_down double precision :: trace_old + double precision,allocatable :: R_tmp(:,:) + double precision,allocatable :: cp_tmp(:,:) ! Output variables @@ -38,6 +40,8 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB chem_pot_change = 0d0 grad_electrons = 1d0 trace_1rdm = -1d0 + allocate(R_tmp(nOrb2,nOrb2)) + allocate(cp_tmp(nOrb2,nOrb2)) ! Set H_HFB to its non-chemical potential dependent contribution @@ -79,8 +83,8 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB isteps = isteps + 1 chem_pot = chem_pot + chem_pot_change call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) - call diag_H_hfb(nOrb,nOrb2,chem_pot+delta_chem_pot,trace_up,H_hfb,cp,R,eHFB_) - call diag_H_hfb(nOrb,nOrb2,chem_pot-delta_chem_pot,trace_down,H_hfb,cp,R,eHFB_) + call diag_H_hfb(nOrb,nOrb2,chem_pot+delta_chem_pot,trace_up,H_hfb,cp_tmp,R_tmp,eHFB_) + call diag_H_hfb(nOrb,nOrb2,chem_pot-delta_chem_pot,trace_down,H_hfb,cp_tmp,R_tmp,eHFB_) grad_electrons = (trace_up-trace_down)/(2.0d0*delta_chem_pot) chem_pot_change = -(trace_1rdm-nO)/(grad_electrons+1d-10) write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') & @@ -95,6 +99,8 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB H_hfb(iorb,iorb)=H_hfb(iorb,iorb)-chem_pot H_hfb(iorb+nOrb,iorb+nOrb)=H_hfb(iorb+nOrb,iorb+nOrb)+chem_pot enddo + + deallocate(R_tmp,cp_tmp) end subroutine diff --git a/src/HF/print_HFB.f90 b/src/HF/print_HFB.f90 index 903c7718..77258ed8 100644 --- a/src/HF/print_HFB.f90 +++ b/src/HF/print_HFB.f90 @@ -1,7 +1,7 @@ ! --- -subroutine print_HFB(nBas, nOrb, nO, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_pot, dipole) +subroutine print_HFB(nBas, nOrb, nO, N_anom, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_pot, dipole) ! Print one-electron energies and other stuff for G0W0 @@ -21,11 +21,14 @@ subroutine print_HFB(nBas, nOrb, nO, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_p double precision,intent(in) :: EL double precision,intent(in) :: ERHF double precision,intent(in) :: chem_pot + double precision,intent(in) :: N_anom double precision,intent(in) :: dipole(ncart) ! Local variables + integer :: iorb integer :: ixyz + double precision :: trace_occ logical :: dump_orb = .false. @@ -49,6 +52,7 @@ subroutine print_HFB(nBas, nOrb, nO, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_p write(*,'(A33,1X,F16.10,A3)') ' HFB energy = ',ERHF + ENuc,' au' write(*,'(A50)') '---------------------------------------' write(*,'(A33,1X,F16.10,A3)') ' Chemical potential = ',chem_pot,' au' + write(*,'(A33,1X,F16.10,A3)') ' | Anomalous dens | = ',N_anom,' ' write(*,'(A50)') '---------------------------------------' write(*,'(A36)') ' Dipole moment (Debye) ' write(*,'(10X,4A10)') 'X','Y','Z','Tot.' @@ -61,7 +65,15 @@ subroutine print_HFB(nBas, nOrb, nO, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_p write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' HFB occupation numbers ' write(*,'(A50)') '---------------------------------------' - call vecout(2*nOrb, Occ) + trace_occ=0d0 + do iorb=1,2*nOrb + if(abs(Occ(2*nOrb-iorb))>1d-8) then + write(*,'(I7,10F15.8)') iorb,Occ(2*nOrb-iorb) + endif + trace_occ=trace_occ+Occ(iorb) + enddo + write(*,*) + write(*,'(A33,1X,F16.10,A3)') ' Trace [ 1D ] = ',trace_occ,' ' write(*,*) end subroutine From d07bcb4581b9fe2383299ca3be8f4cd8ea1c6b4e Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 3 Feb 2025 13:38:17 +0100 Subject: [PATCH 19/32] 5 point stencil derv N(mu) --- src/HF/HFB.f90 | 97 +++++++++++++++++++++++------------------ src/HF/fix_chem_pot.f90 | 13 ++++-- 2 files changed, 65 insertions(+), 45 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index c11d785e..8cd439ac 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -37,7 +37,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, integer :: iorb integer :: nSCF integer :: nOrb2 - integer :: nBas_Sq + integer :: nOrb_Sq integer :: n_diis double precision :: ET double precision :: EV @@ -57,14 +57,13 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision,allocatable :: Occ(:) double precision,allocatable :: err(:,:) double precision,allocatable :: err_diis(:,:) - double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: H_hfb_diis(:,:) double precision,allocatable :: J(:,:) double precision,allocatable :: K(:,:) double precision,allocatable :: eigVEC(:,:) double precision,allocatable :: H_hfb(:,:) double precision,allocatable :: R(:,:) - double precision,allocatable :: P_old(:,:) - double precision,allocatable :: Panom_old(:,:) + double precision,allocatable :: R_old(:,:) ! Output variables @@ -86,7 +85,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Useful quantities - nBas_Sq = nBas*nBas + nOrb_Sq = nOrb*nOrb nOrb2 = nOrb+nOrb ! Memory allocation @@ -94,17 +93,16 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, allocate(J(nBas,nBas)) allocate(K(nBas,nBas)) - allocate(err(nBas,nBas)) + allocate(err(nOrb2,nOrb2)) allocate(eigVEC(nOrb2,nOrb2)) allocate(H_hfb(nOrb2,nOrb2)) allocate(R(nOrb2,nOrb2)) - allocate(P_old(nBas,nBas)) - allocate(Panom_old(nBas,nBas)) + allocate(R_old(nOrb2,nOrb2)) allocate(eigVAL(nOrb2)) - allocate(err_diis(nBas_Sq,max_diis)) - allocate(F_diis(nBas_Sq,max_diis)) + allocate(err_diis(nOrb_Sq,max_diis)) + allocate(H_hfb_diis(nOrb_Sq,max_diis)) ! Guess chem. pot. @@ -114,7 +112,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, thrs_N = 1d-8 n_diis = 0 - F_diis(:,:) = 0d0 + H_hfb_diis(:,:) = 0d0 err_diis(:,:) = 0d0 rcond = 0d0 @@ -135,13 +133,15 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) + R_old(iorb,iorb) = Occ(iorb) + R_old(iorb+nOrb,iorb+nOrb) = 1d0-Occ(iorb) + R_old(iorb ,iorb+nOrb) = sqrt(Occ(iorb)*(1d0-Occ(iorb))) + R_old(iorb+nOrb,iorb ) = sqrt(Occ(iorb)*(1d0-Occ(iorb))) enddo deallocate(Occ) endif - P_old(:,:) = P(:,:) - Panom_old(:,:) = Panom(:,:) - P(:,:) = 2d0*P(:,:) + P(:,:) = 2d0*P(:,:) Conv = 1d0 nSCF = 0 @@ -178,21 +178,53 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, eigVEC(:,:) = H_hfb(:,:) call diagonalize_matrix(nOrb2,eigVEC,eigVAL) - ! Build R and extract P and Panom + ! Build R - trace_1rdm = 0d0 R(:,:) = 0d0 do iorb=1,nOrb R(:,:) = R(:,:) + matmul(eigVEC(:,iorb:iorb),transpose(eigVEC(:,iorb:iorb))) enddo + trace_1rdm = 0d0 do iorb=1,nOrb - trace_1rdm = trace_1rdm + R(iorb,iorb) + trace_1rdm = trace_1rdm+R(iorb,iorb) enddo ! Adjust the chemical potential if( abs(trace_1rdm-nO) > thrs_N ) & - call fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,eigVEC,R,eigVAL) + call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_hfb,eigVEC,R,eigVAL) + + ! DIIS extrapolation TODO + + !if(max_diis > 1 .and. .false.) then + if(max_diis > 1) then + + write(*,*) ' Doing DIIS' + + err = R_old - R + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nOrb_Sq,nOrb_Sq,n_diis,err_diis,H_hfb_diis,err,H_HFB) + + eigVEC(:,:) = H_hfb(:,:) + call diagonalize_matrix(nOrb2,eigVEC,eigVAL) + + ! Build R and check trace + + trace_1rdm = 0d0 + R(:,:) = 0d0 + do iorb=1,nOrb + R(:,:) = R(:,:) + matmul(eigVEC(:,iorb:iorb),transpose(eigVEC(:,iorb:iorb))) + enddo + do iorb=1,nOrb + trace_1rdm = trace_1rdm + R(iorb,iorb) + enddo + + ! Adjust the chemical potential + + if( abs(trace_1rdm-nO) > thrs_N ) & + call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_hfb,eigVEC,R,eigVAL) + + end if ! Extract P and Panom from R @@ -230,12 +262,9 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, if(nSCF > 1) then - err = P_old - P + err = R_old - R Conv = maxval(abs(err)) - err = Panom_old - Panom - Conv = Conv + maxval(abs(err)) - P_old(:,:) = P(:,:) - Panom_old(:,:) = Panom(:,:) + R_old(:,:) = R(:,:) endif @@ -266,7 +295,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,err,eigVEC,H_hfb,R,P_old,Panom_old,eigVAL,err_diis,F_diis) + deallocate(J,K,err,eigVEC,H_hfb,R,R_old,eigVAL,err_diis,H_hfb_diis) stop @@ -278,13 +307,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb) call diagonalize_matrix(nOrb2,eigVEC,eigVAL) eigVAL(:) = 2d0*eigVAL(:) - Panom_old(:,:) = 0d0 - Panom_old(1:nOrb,1:nOrb) = R(1:nOrb,nOrb+1:nOrb2) - Panom_old = matmul(transpose(Panom_old),Panom_old) - norm_anom = 0d0 - do iorb=1,nBas - norm_anom = norm_anom + Panom_old(iorb,iorb) - enddo + norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2))) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) call print_HFB(nBas,nOrb,nO,norm_anom,eigVAL,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole) @@ -301,21 +324,11 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Memory deallocation - deallocate(J,K,err,eigVEC,H_hfb,R,P_old,Panom_old,eigVAL,err_diis,F_diis) + deallocate(J,K,err,eigVEC,H_hfb,R,R_old,eigVAL,err_diis,H_hfb_diis) end subroutine - - ! ! DIIS extrapolation TODO check and adapt - ! - ! 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 TODO check and adapt ! ! if(level_shift > 0d0 .and. Conv > thresh) then diff --git a/src/HF/fix_chem_pot.f90 b/src/HF/fix_chem_pot.f90 index baa37783..fff18ac2 100644 --- a/src/HF/fix_chem_pot.f90 +++ b/src/HF/fix_chem_pot.f90 @@ -1,4 +1,4 @@ -subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB_) +subroutine fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB_) ! Fix the chemical potential to integrate to the N for 2N electron systems @@ -6,7 +6,7 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB ! Input variables - integer,intent(in) :: nO,nOrb,nOrb2 + integer,intent(in) :: nO,nOrb,nOrb2,nSCF double precision,intent(in) :: thrs_N ! Local variables @@ -17,8 +17,10 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB double precision :: delta_chem_pot double precision :: chem_pot_change double precision :: grad_electrons + double precision :: trace_2up double precision :: trace_up double precision :: trace_down + double precision :: trace_2down double precision :: trace_old double precision,allocatable :: R_tmp(:,:) double precision,allocatable :: cp_tmp(:,:) @@ -50,6 +52,8 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB H_hfb(iorb+nOrb,iorb+nOrb)=H_hfb(iorb+nOrb,iorb+nOrb)-chem_pot enddo + write(*,*) + write(*,'(a,I5)') ' Fixing the Tr[1D] at iteration ',nSCF write(*,*) write(*,*)'------------------------------------------------------' write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1A15,2X,A1)') & @@ -83,9 +87,12 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB isteps = isteps + 1 chem_pot = chem_pot + chem_pot_change call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) + call diag_H_hfb(nOrb,nOrb2,chem_pot+2d0*delta_chem_pot,trace_2up,H_hfb,cp_tmp,R_tmp,eHFB_) call diag_H_hfb(nOrb,nOrb2,chem_pot+delta_chem_pot,trace_up,H_hfb,cp_tmp,R_tmp,eHFB_) call diag_H_hfb(nOrb,nOrb2,chem_pot-delta_chem_pot,trace_down,H_hfb,cp_tmp,R_tmp,eHFB_) - grad_electrons = (trace_up-trace_down)/(2.0d0*delta_chem_pot) + call diag_H_hfb(nOrb,nOrb2,chem_pot-2d0*delta_chem_pot,trace_2down,H_hfb,cp_tmp,R_tmp,eHFB_) +! grad_electrons = (trace_up-trace_down)/(2d0*delta_chem_pot) + grad_electrons = (-trace_2up+8d0*trace_up-8d0*trace_down+trace_2down)/(12d0*delta_chem_pot) chem_pot_change = -(trace_1rdm-nO)/(grad_electrons+1d-10) write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') & '|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|' From a79bcd106b1aa54ffc8f6a67ea2defbe42770520 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 3 Feb 2025 14:37:29 +0100 Subject: [PATCH 20/32] Doings DIIS MO basis... --- src/HF/HFB.f90 | 36 +++++++++--------------------------- 1 file changed, 9 insertions(+), 27 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 8cd439ac..7ba4d35a 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -133,10 +133,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) - R_old(iorb,iorb) = Occ(iorb) - R_old(iorb+nOrb,iorb+nOrb) = 1d0-Occ(iorb) - R_old(iorb ,iorb+nOrb) = sqrt(Occ(iorb)*(1d0-Occ(iorb))) - R_old(iorb+nOrb,iorb ) = sqrt(Occ(iorb)*(1d0-Occ(iorb))) enddo deallocate(Occ) endif @@ -154,7 +150,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*) 'Enterning HFB SCF procedure' write(*,*) do while(Conv > thresh .and. nSCF < maxSCF) - + + ! Increment nSCF = nSCF + 1 @@ -197,11 +194,11 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! DIIS extrapolation TODO !if(max_diis > 1 .and. .false.) then - if(max_diis > 1) then + if(max_diis > 1 .and. nSCF>1) then write(*,*) ' Doing DIIS' - err = R_old - R + err = matmul(H_HFB,R_old) - matmul(R_old,H_HFB) n_diis = min(n_diis+1,max_diis) call DIIS_extrapolation(rcond,nOrb_Sq,nOrb_Sq,n_diis,err_diis,H_hfb_diis,err,H_HFB) @@ -233,41 +230,32 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, P(:,:) = 2d0*matmul(X,matmul(R(1:nOrb,1:nOrb),transpose(X))) Panom(:,:) = matmul(X,matmul(R(1:nOrb,nOrb+1:nOrb2),transpose(X))) - ! 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)) - ! Anomalous energy - EL = trace_matrix(nBas,matmul(Panom,Delta)) - ! Total energy - EHFB = ET + EV + EJ + EK + EL ! Check convergence if(nSCF > 1) then - err = R_old - R + err = matmul(H_HFB,R_old) - matmul(R_old,H_HFB) Conv = maxval(abs(err)) - R_old(:,:) = R(:,:) endif + ! Update R_old + + R_old(:,:) = R(:,:) + ! Dump results write(*,*)'-----------------------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1A16,1X,A1,1X,A10,2X,A1,1X)') & @@ -328,9 +316,3 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, end subroutine - - ! ! Level shift TODO check and adapt - ! - ! if(level_shift > 0d0 .and. Conv > thresh) then - ! call level_shifting(level_shift,nBas,nOrb,nO,S,c,F) - ! endif From 674ceb658cef1fbb3e8b5f9be06452da5bc6ff97 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 3 Feb 2025 16:41:09 +0100 Subject: [PATCH 21/32] TODO DIIS Atomic --- src/HF/HFB.f90 | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 7ba4d35a..5a2f17ce 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -57,11 +57,11 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision,allocatable :: Occ(:) double precision,allocatable :: err(:,:) double precision,allocatable :: err_diis(:,:) - double precision,allocatable :: H_hfb_diis(:,:) + double precision,allocatable :: H_HFB_diis(:,:) double precision,allocatable :: J(:,:) double precision,allocatable :: K(:,:) double precision,allocatable :: eigVEC(:,:) - double precision,allocatable :: H_hfb(:,:) + double precision,allocatable :: H_HFB(:,:) double precision,allocatable :: R(:,:) double precision,allocatable :: R_old(:,:) @@ -96,13 +96,13 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, allocate(err(nOrb2,nOrb2)) allocate(eigVEC(nOrb2,nOrb2)) - allocate(H_hfb(nOrb2,nOrb2)) + allocate(H_HFB(nOrb2,nOrb2)) allocate(R(nOrb2,nOrb2)) allocate(R_old(nOrb2,nOrb2)) allocate(eigVAL(nOrb2)) allocate(err_diis(nOrb_Sq,max_diis)) - allocate(H_hfb_diis(nOrb_Sq,max_diis)) + allocate(H_HFB_diis(nOrb_Sq,max_diis)) ! Guess chem. pot. @@ -112,7 +112,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, thrs_N = 1d-8 n_diis = 0 - H_hfb_diis(:,:) = 0d0 + H_HFB_diis(:,:) = 0d0 err_diis(:,:) = 0d0 rcond = 0d0 @@ -166,13 +166,13 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Diagonalize H_HFB matrix - H_hfb(:,:) = 0d0 - H_hfb(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X)) - H_hfb(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_hfb(1:nOrb,1:nOrb) - H_hfb(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) - H_hfb(nOrb+1:nOrb2,1:nOrb ) = transpose(H_hfb(1:nOrb,nOrb+1:nOrb2)) + H_HFB(:,:) = 0d0 + H_HFB(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X)) + H_HFB(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_HFB(1:nOrb,1:nOrb) + H_HFB(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) + H_HFB(nOrb+1:nOrb2,1:nOrb ) = transpose(H_HFB(1:nOrb,nOrb+1:nOrb2)) - eigVEC(:,:) = H_hfb(:,:) + eigVEC(:,:) = H_HFB(:,:) call diagonalize_matrix(nOrb2,eigVEC,eigVAL) ! Build R @@ -189,7 +189,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Adjust the chemical potential if( abs(trace_1rdm-nO) > thrs_N ) & - call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_hfb,eigVEC,R,eigVAL) + call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_HFB,eigVEC,R,eigVAL) ! DIIS extrapolation TODO @@ -200,9 +200,9 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, err = matmul(H_HFB,R_old) - matmul(R_old,H_HFB) n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nOrb_Sq,nOrb_Sq,n_diis,err_diis,H_hfb_diis,err,H_HFB) + call DIIS_extrapolation(rcond,nOrb_Sq,nOrb_Sq,n_diis,err_diis,H_HFB_diis,err,H_HFB) - eigVEC(:,:) = H_hfb(:,:) + eigVEC(:,:) = H_HFB(:,:) call diagonalize_matrix(nOrb2,eigVEC,eigVAL) ! Build R and check trace @@ -219,7 +219,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Adjust the chemical potential if( abs(trace_1rdm-nO) > thrs_N ) & - call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_hfb,eigVEC,R,eigVAL) + call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_HFB,eigVEC,R,eigVAL) end if @@ -283,7 +283,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,err,eigVEC,H_hfb,R,R_old,eigVAL,err_diis,H_hfb_diis) + deallocate(J,K,err,eigVEC,H_HFB,R,R_old,eigVAL,err_diis,H_HFB_diis) stop @@ -312,7 +312,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Memory deallocation - deallocate(J,K,err,eigVEC,H_hfb,R,R_old,eigVAL,err_diis,H_hfb_diis) + deallocate(J,K,err,eigVEC,H_HFB,R,R_old,eigVAL,err_diis,H_HFB_diis) end subroutine From a75d5b27cfd37a33da5e961d4804159253e5460e Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Mon, 3 Feb 2025 17:30:13 +0100 Subject: [PATCH 22/32] Move 1 space --- src/HF/HFB.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 5a2f17ce..62fd7beb 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -112,7 +112,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, thrs_N = 1d-8 n_diis = 0 - H_HFB_diis(:,:) = 0d0 + H_HFB_diis(:,:) = 0d0 err_diis(:,:) = 0d0 rcond = 0d0 From cad43b2e1b71fc5fa21a1443f7da25dd916e3908 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Tue, 4 Feb 2025 00:43:47 +0100 Subject: [PATCH 23/32] TODO DIIS AOs --- src/HF/HFB.f90 | 51 +++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 46 insertions(+), 5 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 62fd7beb..b8740a9f 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -34,6 +34,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Local variables + integer :: nBas2 integer :: iorb integer :: nSCF integer :: nOrb2 @@ -65,6 +66,11 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision,allocatable :: R(:,:) double precision,allocatable :: R_old(:,:) + double precision,allocatable :: err_ao(:,:) + double precision,allocatable :: S_ao(:,:) + double precision,allocatable :: R_ao_old(:,:) + double precision,allocatable :: H_HFB_ao(:,:) + ! Output variables double precision,intent(out) :: EHFB @@ -87,6 +93,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, nOrb_Sq = nOrb*nOrb nOrb2 = nOrb+nOrb + nBas2 = nBas+nBas ! Memory allocation @@ -94,13 +101,17 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, allocate(K(nBas,nBas)) allocate(err(nOrb2,nOrb2)) - allocate(eigVEC(nOrb2,nOrb2)) allocate(H_HFB(nOrb2,nOrb2)) allocate(R(nOrb2,nOrb2)) allocate(R_old(nOrb2,nOrb2)) allocate(eigVAL(nOrb2)) + allocate(err_ao(nBas2,nBas2)) + allocate(S_ao(nBas2,nBas2)) + allocate(R_ao_old(nBas2,nBas2)) + allocate(H_HFB_ao(nBas2,nBas2)) + allocate(err_diis(nOrb_Sq,max_diis)) allocate(H_HFB_diis(nOrb_Sq,max_diis)) @@ -138,6 +149,9 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, endif P(:,:) = 2d0*P(:,:) + S_ao(:,:) = 0d0 + S_ao(1:nBas ,1:nBas ) = S(1:nBas,1:nBas) + S_ao(nBas+1:nBas2,nBas+1:nBas2) = S(1:nBas,1:nBas) Conv = 1d0 nSCF = 0 @@ -160,7 +174,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, call Hartree_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) - call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) + call anomalous_matrix_AO_basis(nBas,-2d0*Panom,ERI,Delta) ! TODO recover usual Panom F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) @@ -170,7 +184,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, H_HFB(1:nOrb ,1:nOrb ) = matmul(transpose(X),matmul(F,X)) H_HFB(nOrb+1:nOrb2,nOrb+1:nOrb2) = -H_HFB(1:nOrb,1:nOrb) H_HFB(1:nOrb ,nOrb+1:nOrb2) = matmul(transpose(X),matmul(Delta,X)) - H_HFB(nOrb+1:nOrb2,1:nOrb ) = transpose(H_HFB(1:nOrb,nOrb+1:nOrb2)) + H_HFB(nOrb+1:nOrb2,1:nOrb ) = H_HFB(1:nOrb,nOrb+1:nOrb2) eigVEC(:,:) = H_HFB(:,:) call diagonalize_matrix(nOrb2,eigVEC,eigVAL) @@ -193,8 +207,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! DIIS extrapolation TODO - !if(max_diis > 1 .and. .false.) then - if(max_diis > 1 .and. nSCF>1) then + if(max_diis > 1 .and. .false.) then + !if(max_diis > 1 .and. nSCF>1) then write(*,*) ' Doing DIIS' @@ -250,12 +264,37 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, err = matmul(H_HFB,R_old) - matmul(R_old,H_HFB) Conv = maxval(abs(err)) +write(*,*) 'Err MO ',maxval(abs(err)) +do iorb=1,nOrb2 + write(*,'(*(F10.5))') err(iorb,:) +enddo + + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) + H_HFB_ao(:,:) = 0d0 + H_HFB_ao(1:nBas ,1:nBas ) = F(1:nBas,1:nBas) + H_HFB_ao(nBas+1:nBas2,nBas+1:nBas2) = -F(1:nBas,1:nBas) + H_HFB_ao(1:nBas ,nBas+1:nBas2) = Delta(1:nBas,1:nBas) + H_HFB_ao(nBas+1:nBas2,1:nBas ) = Delta(1:nBas,1:nBas) + err_ao = matmul(H_HFB_ao,matmul(R_ao_old,S_ao)) - matmul(matmul(S_ao,R_ao_old),H_HFB_ao) + +write(*,*) 'Err AO ',maxval(abs(err)) +do iorb=1,nBas2 + write(*,'(*(F10.5))') err(iorb,:) +enddo + endif ! Update R_old R_old(:,:) = R(:,:) + R_ao_old(:,:) = 0d0 + R_ao_old(1:nBas ,1:nBas ) = P(1:nBas,1:nBas) + R_ao_old(nBas+1:nBas2,nBas+1:nBas2) = matmul(X(1:nBas,1:nOrb), transpose(X(1:nBas,1:nOrb)))-P(1:nBas,1:nBas) + R_ao_old(1:nBas ,nBas+1:nBas2) = Panom(1:nBas,1:nBas) + R_ao_old(nBas+1:nBas2,1:nBas ) = Panom(1:nBas,1:nBas) + + ! Dump results write(*,*)'-----------------------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1A16,1X,A1,1X,A10,2X,A1,1X)') & @@ -284,6 +323,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*) deallocate(J,K,err,eigVEC,H_HFB,R,R_old,eigVAL,err_diis,H_HFB_diis) + deallocate(err_ao,S_ao,R_ao_old,H_HFB_ao) stop @@ -313,6 +353,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Memory deallocation deallocate(J,K,err,eigVEC,H_HFB,R,R_old,eigVAL,err_diis,H_HFB_diis) + deallocate(err_ao,S_ao,R_ao_old,H_HFB_ao) end subroutine From b260a7a70deb26961c984478923b8caa6eeefe1b Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Tue, 4 Feb 2025 14:02:08 +0100 Subject: [PATCH 24/32] Incl DIIS method --- src/HF/HFB.f90 | 64 +++++++++++++++++++++++--------------------------- 1 file changed, 29 insertions(+), 35 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index b8740a9f..1b2405ee 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -38,7 +38,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, integer :: iorb integer :: nSCF integer :: nOrb2 - integer :: nOrb_Sq + integer :: nBas2_Sq integer :: n_diis double precision :: ET double precision :: EV @@ -56,7 +56,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision,external :: trace_matrix double precision,allocatable :: eigVAL(:) double precision,allocatable :: Occ(:) - double precision,allocatable :: err(:,:) double precision,allocatable :: err_diis(:,:) double precision,allocatable :: H_HFB_diis(:,:) double precision,allocatable :: J(:,:) @@ -64,10 +63,10 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision,allocatable :: eigVEC(:,:) double precision,allocatable :: H_HFB(:,:) double precision,allocatable :: R(:,:) - double precision,allocatable :: R_old(:,:) double precision,allocatable :: err_ao(:,:) double precision,allocatable :: S_ao(:,:) + double precision,allocatable :: X_ao(:,:) double precision,allocatable :: R_ao_old(:,:) double precision,allocatable :: H_HFB_ao(:,:) @@ -91,29 +90,28 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Useful quantities - nOrb_Sq = nOrb*nOrb nOrb2 = nOrb+nOrb nBas2 = nBas+nBas + nBas2_Sq = nBas2*nBas2 ! Memory allocation allocate(J(nBas,nBas)) allocate(K(nBas,nBas)) - allocate(err(nOrb2,nOrb2)) allocate(eigVEC(nOrb2,nOrb2)) allocate(H_HFB(nOrb2,nOrb2)) allocate(R(nOrb2,nOrb2)) - allocate(R_old(nOrb2,nOrb2)) allocate(eigVAL(nOrb2)) allocate(err_ao(nBas2,nBas2)) allocate(S_ao(nBas2,nBas2)) + allocate(X_ao(nBas2,nOrb2)) allocate(R_ao_old(nBas2,nBas2)) allocate(H_HFB_ao(nBas2,nBas2)) - allocate(err_diis(nOrb_Sq,max_diis)) - allocate(H_HFB_diis(nOrb_Sq,max_diis)) + allocate(err_diis(nBas2_Sq,max_diis)) + allocate(H_HFB_diis(nBas2_Sq,max_diis)) ! Guess chem. pot. @@ -152,6 +150,9 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, S_ao(:,:) = 0d0 S_ao(1:nBas ,1:nBas ) = S(1:nBas,1:nBas) S_ao(nBas+1:nBas2,nBas+1:nBas2) = S(1:nBas,1:nBas) + X_ao(:,:) = 0d0 + X_ao(1:nBas ,1:nOrb ) = X(1:nBas,1:nOrb) + X_ao(nBas+1:nBas2,nOrb+1:nOrb2) = X(1:nBas,1:nOrb) Conv = 1d0 nSCF = 0 @@ -174,7 +175,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, call Hartree_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) - call anomalous_matrix_AO_basis(nBas,-2d0*Panom,ERI,Delta) ! TODO recover usual Panom + call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) @@ -205,17 +206,24 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, if( abs(trace_1rdm-nO) > thrs_N ) & call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_HFB,eigVEC,R,eigVAL) - ! DIIS extrapolation TODO + ! DIIS extrapolation - if(max_diis > 1 .and. .false.) then - !if(max_diis > 1 .and. nSCF>1) then + if(max_diis > 1 .and. nSCF>1) then write(*,*) ' Doing DIIS' - err = matmul(H_HFB,R_old) - matmul(R_old,H_HFB) + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) + H_HFB_ao(:,:) = 0d0 + H_HFB_ao(1:nBas ,1:nBas ) = F(1:nBas,1:nBas) + H_HFB_ao(nBas+1:nBas2,nBas+1:nBas2) = -F(1:nBas,1:nBas) + H_HFB_ao(1:nBas ,nBas+1:nBas2) = Delta(1:nBas,1:nBas) + H_HFB_ao(nBas+1:nBas2,1:nBas ) = Delta(1:nBas,1:nBas) + err_ao = matmul(H_HFB_ao,matmul(R_ao_old,S_ao)) - matmul(matmul(S_ao,R_ao_old),H_HFB_ao) + n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nOrb_Sq,nOrb_Sq,n_diis,err_diis,H_HFB_diis,err,H_HFB) + call DIIS_extrapolation(rcond,nBas2_Sq,nBas2_Sq,n_diis,err_diis,H_HFB_diis,err_ao,H_HFB_ao) + H_HFB = matmul(transpose(X_ao),matmul(H_HFB_ao,X_ao)) eigVEC(:,:) = H_HFB(:,:) call diagonalize_matrix(nOrb2,eigVEC,eigVAL) @@ -261,14 +269,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, if(nSCF > 1) then - err = matmul(H_HFB,R_old) - matmul(R_old,H_HFB) - Conv = maxval(abs(err)) - -write(*,*) 'Err MO ',maxval(abs(err)) -do iorb=1,nOrb2 - write(*,'(*(F10.5))') err(iorb,:) -enddo - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) H_HFB_ao(:,:) = 0d0 H_HFB_ao(1:nBas ,1:nBas ) = F(1:nBas,1:nBas) @@ -276,21 +276,15 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, H_HFB_ao(1:nBas ,nBas+1:nBas2) = Delta(1:nBas,1:nBas) H_HFB_ao(nBas+1:nBas2,1:nBas ) = Delta(1:nBas,1:nBas) err_ao = matmul(H_HFB_ao,matmul(R_ao_old,S_ao)) - matmul(matmul(S_ao,R_ao_old),H_HFB_ao) - -write(*,*) 'Err AO ',maxval(abs(err)) -do iorb=1,nBas2 - write(*,'(*(F10.5))') err(iorb,:) -enddo + Conv = maxval(abs(err_ao)) endif ! Update R_old - R_old(:,:) = R(:,:) - R_ao_old(:,:) = 0d0 - R_ao_old(1:nBas ,1:nBas ) = P(1:nBas,1:nBas) - R_ao_old(nBas+1:nBas2,nBas+1:nBas2) = matmul(X(1:nBas,1:nOrb), transpose(X(1:nBas,1:nOrb)))-P(1:nBas,1:nBas) + R_ao_old(1:nBas ,1:nBas ) = 0.5d0*P(1:nBas,1:nBas) + R_ao_old(nBas+1:nBas2,nBas+1:nBas2) = matmul(X(1:nBas,1:nOrb), transpose(X(1:nBas,1:nOrb)))-0.5d0*P(1:nBas,1:nBas) R_ao_old(1:nBas ,nBas+1:nBas2) = Panom(1:nBas,1:nBas) R_ao_old(nBas+1:nBas2,1:nBas ) = Panom(1:nBas,1:nBas) @@ -322,8 +316,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,err,eigVEC,H_HFB,R,R_old,eigVAL,err_diis,H_HFB_diis) - deallocate(err_ao,S_ao,R_ao_old,H_HFB_ao) + deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis) + deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao) stop @@ -352,8 +346,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Memory deallocation - deallocate(J,K,err,eigVEC,H_HFB,R,R_old,eigVAL,err_diis,H_HFB_diis) - deallocate(err_ao,S_ao,R_ao_old,H_HFB_ao) + deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis) + deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao) end subroutine From c05631e1ba8249a5b26d8b1318e776e66a83ec57 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Tue, 4 Feb 2025 14:34:31 +0100 Subject: [PATCH 25/32] Incl sigma to change sign HFB --- src/AOtoMO/anomalous_matrix_AO_basis.f90 | 5 +++-- src/HF/HFB.f90 | 6 +++--- src/QuAcK/BQuAcK.f90 | 6 +++--- src/QuAcK/QuAcK.f90 | 7 ++++--- src/QuAcK/read_options.f90 | 8 +++++--- 5 files changed, 18 insertions(+), 14 deletions(-) diff --git a/src/AOtoMO/anomalous_matrix_AO_basis.f90 b/src/AOtoMO/anomalous_matrix_AO_basis.f90 index 4f2583ef..7a7de116 100644 --- a/src/AOtoMO/anomalous_matrix_AO_basis.f90 +++ b/src/AOtoMO/anomalous_matrix_AO_basis.f90 @@ -1,4 +1,4 @@ -subroutine anomalous_matrix_AO_basis(nBas,Pa,ERI,L) +subroutine anomalous_matrix_AO_basis(nBas,sigma,Pa,ERI,L) ! Compute anomalous L matrix in the AO basis @@ -8,6 +8,7 @@ subroutine anomalous_matrix_AO_basis(nBas,Pa,ERI,L) ! Input variables integer,intent(in) :: nBas + double precision,intent(in) :: sigma double precision,intent(in) :: Pa(nBas,nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -25,7 +26,7 @@ subroutine anomalous_matrix_AO_basis(nBas,Pa,ERI,L) do si=1,nBas do la=1,nBas do mu=1,nBas - L(mu,nu) = L(mu,nu) + Pa(la,si)*ERI(la,si,mu,nu) + L(mu,nu) = L(mu,nu) + sigma*Pa(la,si)*ERI(la,si,mu,nu) end do end do end do diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 1b2405ee..64dc908e 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -1,6 +1,6 @@ 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,Delta, & - temperature) + temperature,sigma) ! Perform Hartree-Fock Bogoliubov calculation @@ -23,7 +23,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - double precision,intent(in) :: temperature + double precision,intent(in) :: temperature,sigma double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) @@ -175,7 +175,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, call Hartree_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) - call anomalous_matrix_AO_basis(nBas,Panom,ERI,Delta) + call anomalous_matrix_AO_basis(nBas,sigma,Panom,ERI,Delta) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:) diff --git a/src/QuAcK/BQuAcK.f90 b/src/QuAcK/BQuAcK.f90 index 836ee9ac..572118b7 100644 --- a/src/QuAcK/BQuAcK.f90 +++ b/src/QuAcK/BQuAcK.f90 @@ -1,6 +1,6 @@ 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,temperature) + guess_type,mix,temperature,sigma) ! Restricted branch of QuAcK @@ -19,7 +19,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, integer,intent(in) :: nV integer,intent(in) :: nR double precision,intent(in) :: ENuc - double precision,intent(in) :: temperature + double precision,intent(in) :: temperature,sigma double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart) @@ -94,7 +94,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, 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,Delta,temperature) + FHF,Delta,temperature,sigma) call wall_time(end_HF) t_HF = end_HF - start_HF diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 22d607b9..af50d219 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -73,7 +73,7 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest - double precision :: temperature + double precision :: temperature,sigma character(len=256) :: working_dir @@ -136,7 +136,8 @@ program QuAcK maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, & maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,doppBSE,dBSE,dTDA,temperature) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & + temperature,sigma) !------------------! ! Hardware ! @@ -293,7 +294,7 @@ program QuAcK if(doBQuAcK) & 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, & - temperature) + temperature,sigma) !-----------! ! Stop Test ! diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 11af97ea..4fb3e4d0 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -7,7 +7,8 @@ subroutine read_options(working_dir, maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, & maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,doppBSE,dBSE,dTDA,temperature) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & + temperature,sigma) ! Read desired methods @@ -73,6 +74,7 @@ subroutine read_options(working_dir, logical,intent(out) :: dTDA double precision,intent(out) :: temperature + double precision,intent(out) :: sigma ! Local variables @@ -219,12 +221,12 @@ subroutine read_options(working_dir, if(ans4 == 'T') dBSE = .true. if(ans5 == 'F') dTDA = .false. - ! Options for dynamical Fermi-Dirac occupancies + ! Options for Hartree-Fock Bogoliubov temperature = 0d0 read(1,*) - read(1,*) temperature + read(1,*) temperature,sigma endif From 502517b1655a8772da264a76c88f8dde77a7b199 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Tue, 4 Feb 2025 15:07:46 +0100 Subject: [PATCH 26/32] Ignore calcs files --- .gitignore | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.gitignore b/.gitignore index 4a6ab19f..6974196e 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,8 @@ __pycache__ .ninja_deps +calcs/ +calcs/input +calcs/int +calcs/mol +calcs/water.out From 2b8c9c8bd2cc3578e00d3cb215853e01891178eb Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Tue, 4 Feb 2025 16:58:22 +0100 Subject: [PATCH 27/32] Allow HF chem_pot init --- src/HF/HFB.f90 | 4 +++- src/QuAcK/BQuAcK.f90 | 5 +++-- src/QuAcK/QuAcK.f90 | 5 +++-- src/QuAcK/read_options.f90 | 9 +++++++-- 4 files changed, 16 insertions(+), 7 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 64dc908e..6b64c1b7 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -1,6 +1,6 @@ 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,Delta, & - temperature,sigma) + temperature,sigma,chem_pot_hf) ! Perform Hartree-Fock Bogoliubov calculation @@ -34,6 +34,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Local variables + logical :: chem_pot_hf integer :: nBas2 integer :: iorb integer :: nSCF @@ -135,6 +136,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, Occ(:) = 0d0 Occ(1:nO) = 1d0 call fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) + if(chem_pot_hf) chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1)) P(:,:) = 0d0 Panom(:,:) = 0d0 do iorb=1,nOrb diff --git a/src/QuAcK/BQuAcK.f90 b/src/QuAcK/BQuAcK.f90 index 572118b7..5bf5637e 100644 --- a/src/QuAcK/BQuAcK.f90 +++ b/src/QuAcK/BQuAcK.f90 @@ -1,6 +1,6 @@ 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,temperature,sigma) + guess_type,mix,temperature,sigma,chem_pot_hf) ! Restricted branch of QuAcK @@ -13,6 +13,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, logical,intent(in) :: doHFB + logical,intent(in) :: chem_pot_hf integer,intent(in) :: nNuc,nBas,nOrb integer,intent(in) :: nC integer,intent(in) :: nO @@ -94,7 +95,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, 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,Delta,temperature,sigma) + FHF,Delta,temperature,sigma,chem_pot_hf) call wall_time(end_HF) t_HF = end_HF - start_HF diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index af50d219..5081fe30 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -73,6 +73,7 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest + logical :: chem_pot_hf double precision :: temperature,sigma character(len=256) :: working_dir @@ -137,7 +138,7 @@ program QuAcK maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & - temperature,sigma) + temperature,sigma,chem_pot_hf) !------------------! ! Hardware ! @@ -294,7 +295,7 @@ program QuAcK if(doBQuAcK) & 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, & - temperature,sigma) + temperature,sigma,chem_pot_hf) !-----------! ! Stop Test ! diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 4fb3e4d0..778e305a 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -8,7 +8,7 @@ subroutine read_options(working_dir, maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & - temperature,sigma) + temperature,sigma,chem_pot_hf) ! Read desired methods @@ -73,6 +73,7 @@ subroutine read_options(working_dir, logical,intent(out) :: dBSE logical,intent(out) :: dTDA + logical,intent(out) :: chem_pot_hf double precision,intent(out) :: temperature double precision,intent(out) :: sigma @@ -224,9 +225,13 @@ subroutine read_options(working_dir, ! Options for Hartree-Fock Bogoliubov temperature = 0d0 + sigma = 1d0 + chem_pot_hf = .false. read(1,*) - read(1,*) temperature,sigma + read(1,*) temperature,sigma,ans1 + + if(ans1 == 'T') chem_pot_hf = .true. endif From de424256839e7e19326ba34675bce7420ddb7a1c Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 6 Feb 2025 14:52:12 +0100 Subject: [PATCH 28/32] Allow restart in HFB --- src/HF/HFB.f90 | 32 +++++++--- src/HF/print_HFB.f90 | 12 ++-- src/HF/read_restart_HFB.f90 | 116 +++++++++++++++++++++++++++++++++++ src/HF/write_restart_HFB.f90 | 37 +++++++++++ src/QuAcK/BQuAcK.f90 | 5 +- src/QuAcK/QuAcK.f90 | 5 +- src/QuAcK/read_options.f90 | 7 ++- 7 files changed, 195 insertions(+), 19 deletions(-) create mode 100644 src/HF/read_restart_HFB.f90 create mode 100644 src/HF/write_restart_HFB.f90 diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 6b64c1b7..0279f67b 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -1,6 +1,6 @@ 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,Delta, & - temperature,sigma,chem_pot_hf) + temperature,sigma,chem_pot_hf,restart_hfb) ! Perform Hartree-Fock Bogoliubov calculation @@ -35,6 +35,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Local variables logical :: chem_pot_hf + logical :: restart_hfb integer :: nBas2 integer :: iorb integer :: nSCF @@ -97,6 +98,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Memory allocation + allocate(Occ(nOrb)) + allocate(J(nBas,nBas)) allocate(K(nBas,nBas)) @@ -132,7 +135,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Use Fermi-Dirac occupancies to compute P, Panom, and chem_pot if(abs(temperature)>1d-4) then - allocate(Occ(nOrb)) Occ(:) = 0d0 Occ(1:nO) = 1d0 call fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) @@ -145,7 +147,20 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) enddo - deallocate(Occ) + endif + + ! Read restart file + + if(restart_hfb) then + call read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot) + P(:,:) = 0d0 + Panom(:,:) = 0d0 + do iorb=1,nOrb + P(:,:) = P(:,:) + Occ(iorb) * & + matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) + Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & + matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) + enddo endif P(:,:) = 2d0*P(:,:) @@ -318,7 +333,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis) + deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis,Occ) deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao) stop @@ -326,14 +341,17 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, end if ! Compute dipole moments, occupation numbers and || Anomalous density|| +! also print the restart file eigVEC(:,:) = 0d0 eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb) call diagonalize_matrix(nOrb2,eigVEC,eigVAL) - eigVAL(:) = 2d0*eigVAL(:) + Occ(1:nOrb) = eigVAL(nOrb+1:nOrb2) + c(1:nBas,1:nOrb) = matmul(X(1:nBas,1:nOrb),eigVEC(1:nOrb,nOrb+1:nOrb2)) norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2))) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_HFB(nBas,nOrb,nO,norm_anom,eigVAL,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole) + call write_restart_HFB(nBas,nOrb,Occ,c,chem_pot) + call print_HFB(nBas,nOrb,nO,norm_anom,Occ,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole) ! Testing zone @@ -348,7 +366,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Memory deallocation - deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis) + deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis,Occ) deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao) end subroutine diff --git a/src/HF/print_HFB.f90 b/src/HF/print_HFB.f90 index 77258ed8..ea61f0c1 100644 --- a/src/HF/print_HFB.f90 +++ b/src/HF/print_HFB.f90 @@ -3,7 +3,7 @@ subroutine print_HFB(nBas, nOrb, nO, N_anom, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_pot, dipole) -! Print one-electron energies and other stuff for G0W0 +! Print one-electron energies and other stuff implicit none include 'parameters.h' @@ -12,7 +12,7 @@ subroutine print_HFB(nBas, nOrb, nO, N_anom, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO - double precision,intent(in) :: Occ(2*nOrb) + double precision,intent(in) :: Occ(nOrb) double precision,intent(in) :: ENuc double precision,intent(in) :: ET double precision,intent(in) :: EV @@ -66,11 +66,11 @@ subroutine print_HFB(nBas, nOrb, nO, N_anom, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF write(*,'(A50)') ' HFB occupation numbers ' write(*,'(A50)') '---------------------------------------' trace_occ=0d0 - do iorb=1,2*nOrb - if(abs(Occ(2*nOrb-iorb))>1d-8) then - write(*,'(I7,10F15.8)') iorb,Occ(2*nOrb-iorb) + do iorb=1,nOrb + if(abs(Occ(nOrb+1-iorb))>1d-8) then + write(*,'(I7,10F15.8)') iorb,2d0*Occ(nOrb+1-iorb) endif - trace_occ=trace_occ+Occ(iorb) + trace_occ=trace_occ+2d0*Occ(iorb) enddo write(*,*) write(*,'(A33,1X,F16.10,A3)') ' Trace [ 1D ] = ',trace_occ,' ' diff --git a/src/HF/read_restart_HFB.f90 b/src/HF/read_restart_HFB.f90 new file mode 100644 index 00000000..8d523c2f --- /dev/null +++ b/src/HF/read_restart_HFB.f90 @@ -0,0 +1,116 @@ + +! --- + +subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot) + +! Write the binary file used to restart calcs + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + double precision,intent(in) :: S(nBas,nBas) + +! Local variables + + integer :: ibas,iorb,iorb1,iunit=667 + + integer :: nBas_ + integer :: nOrb_ + double precision :: chem_pot_ + double precision :: max_diff + double precision,allocatable :: eigVAL(:) + double precision,allocatable :: c_tmp(:,:) + double precision,allocatable :: S_mol(:,:) + double precision,allocatable :: X_mol(:,:) + +! Output variables + + double precision,intent(out) :: chem_pot + double precision,intent(out) :: Occ(nOrb) + double precision,intent(out) :: c(nBas,nOrb) + +! Dump results + + allocate(eigVAL(nOrb),c_tmp(nBas,nOrb),S_mol(nOrb,nOrb),X_mol(nOrb,nOrb)) + + open(unit=iunit,form='unformatted',file='hfb_bin',status='old') + read(iunit) nBas_,nOrb_ + read(iunit) chem_pot_ + do iorb=1,nOrb + read(iunit) c_tmp(:,iorb) + enddo + do iorb=1,nOrb + read(iunit) eigVAL(iorb) + enddo + close(iunit) + + if(nBas==nBas_ .and. nOrb==nOrb_) then + write(*,*) + write(*,*)' Reading restart file' + write(*,*) + + chem_pot = chem_pot_ + Occ(:) = eigVAL(:) + c(:,:) = c_tmp(:,:) + + ! Check the orthonormality + + max_diff=-1d0 + S_mol = matmul(transpose(c),matmul(S,c)) + do iorb=1,nOrb + do iorb1=1,iorb + if(iorb==iorb1) then + if(abs(S_mol(iorb,iorb)-1d0)>1d-8) max_diff=abs(S_mol(iorb,iorb)-1d0) + else + if(abs(S_mol(iorb,iorb1))>1d-8) max_diff=abs(S_mol(iorb,iorb1)) + endif + enddo + enddo + if(max_diff>1d-8) then + write(*,*) ' ' + write(*,'(a)') ' Orthonormality violations. Applying Lowdin orthonormalization.' + write(*,*) ' ' + X_mol = S_mol + S_mol = 0d0 + call diagonalize_matrix(nOrb,X_mol,eigVAL) + do iorb=1,nOrb + S_mol(iorb,iorb) = 1d0/(sqrt(eigVAL(iorb)) + 1d-10) + enddo + X_mol = matmul(X_mol,matmul(S_mol,transpose(X_mol))) + c = matmul(c,X_mol) + max_diff=-1d0 + S_mol = matmul(transpose(c),matmul(S,c)) + do iorb=1,nOrb + do iorb1=1,iorb + if(iorb==iorb1) then + if(abs(S_mol(iorb,iorb)-1d0)>1d-8) max_diff=abs(S_mol(iorb,iorb)-1d0) + else + if(abs(S_mol(iorb,iorb1))>1d-8) max_diff=abs(S_mol(iorb,iorb1)) + endif + enddo + enddo + if(max_diff<=1d-8) then + write(*,*) ' ' + write(*,'(a)') ' No orthonormality violations.' + write(*,*) ' ' + else + write(*,*) ' ' + write(*,'(a)') ' Error in Lowdin orthonormalization.' + write(*,*) ' ' + stop + endif + else + write(*,*) ' ' + write(*,'(a)') ' No orthonormality violations.' + write(*,*) ' ' + endif + + endif + + deallocate(eigVAL,c_tmp,S_mol,X_mol) + +end subroutine diff --git a/src/HF/write_restart_HFB.f90 b/src/HF/write_restart_HFB.f90 new file mode 100644 index 00000000..bcd7ef71 --- /dev/null +++ b/src/HF/write_restart_HFB.f90 @@ -0,0 +1,37 @@ + +! --- + +subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot) + +! Write the binary file used to restart calcs + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + double precision,intent(in) :: chem_pot + double precision,intent(in) :: Occ(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + +! Local variables + + integer :: ibas,iorb,iunit=666 + +! Dump results + + open(unit=iunit,form='unformatted',file='hfb_bin') + write(iunit) nBas,nOrb + write(iunit) chem_pot + do iorb=1,nOrb + write(iunit) c(:,nOrb+1-iorb) + enddo + do iorb=1,nOrb + write(iunit) Occ(nOrb+1-iorb) + enddo + write(iunit) iunit + close(iunit) + +end subroutine diff --git a/src/QuAcK/BQuAcK.f90 b/src/QuAcK/BQuAcK.f90 index 5bf5637e..871af92a 100644 --- a/src/QuAcK/BQuAcK.f90 +++ b/src/QuAcK/BQuAcK.f90 @@ -1,6 +1,6 @@ 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,temperature,sigma,chem_pot_hf) + guess_type,mix,temperature,sigma,chem_pot_hf,restart_hfb) ! Restricted branch of QuAcK @@ -13,6 +13,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, logical,intent(in) :: doHFB + logical,intent(in) :: restart_hfb logical,intent(in) :: chem_pot_hf integer,intent(in) :: nNuc,nBas,nOrb integer,intent(in) :: nC @@ -95,7 +96,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc, 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,Delta,temperature,sigma,chem_pot_hf) + FHF,Delta,temperature,sigma,chem_pot_hf,restart_hfb) call wall_time(end_HF) t_HF = end_HF - start_HF diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 5081fe30..446194d2 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -74,6 +74,7 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest logical :: chem_pot_hf + logical :: restart_hfb double precision :: temperature,sigma character(len=256) :: working_dir @@ -138,7 +139,7 @@ program QuAcK maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & - temperature,sigma,chem_pot_hf) + temperature,sigma,chem_pot_hf,restart_hfb) !------------------! ! Hardware ! @@ -295,7 +296,7 @@ program QuAcK if(doBQuAcK) & 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, & - temperature,sigma,chem_pot_hf) + temperature,sigma,chem_pot_hf,restart_hfb) !-----------! ! Stop Test ! diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 778e305a..65b93540 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -8,7 +8,7 @@ subroutine read_options(working_dir, maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA, & - temperature,sigma,chem_pot_hf) + temperature,sigma,chem_pot_hf,restart_hfb) ! Read desired methods @@ -74,6 +74,7 @@ subroutine read_options(working_dir, logical,intent(out) :: dTDA logical,intent(out) :: chem_pot_hf + logical,intent(out) :: restart_hfb double precision,intent(out) :: temperature double precision,intent(out) :: sigma @@ -227,11 +228,13 @@ subroutine read_options(working_dir, temperature = 0d0 sigma = 1d0 chem_pot_hf = .false. + restart_hfb = .false. read(1,*) - read(1,*) temperature,sigma,ans1 + read(1,*) temperature,sigma,ans1,ans2 if(ans1 == 'T') chem_pot_hf = .true. + if(ans2 == 'T') restart_hfb = .true. endif From 1f8757afb6c1b0567ddf27b2807fa8bdfb04cbe2 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 6 Feb 2025 16:32:08 +0100 Subject: [PATCH 29/32] Fixed silly bug search chem_pot done N2 --- src/HF/fermi_dirac_occ.f90 | 2 +- src/HF/fix_chem_pot.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/HF/fermi_dirac_occ.f90 b/src/HF/fermi_dirac_occ.f90 index 7cdfc213..4af8cbfc 100644 --- a/src/HF/fermi_dirac_occ.f90 +++ b/src/HF/fermi_dirac_occ.f90 @@ -62,7 +62,7 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) '|',trace_1rdm,'|',chem_pot,'|' if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then backward=.true. - chem_pot = chem_pot - 2d0*delta_chem_pot + chem_pot = chem_pot - delta_chem_pot delta_chem_pot=-delta_chem_pot endif enddo diff --git a/src/HF/fix_chem_pot.f90 b/src/HF/fix_chem_pot.f90 index fff18ac2..7b7c9bf0 100644 --- a/src/HF/fix_chem_pot.f90 +++ b/src/HF/fix_chem_pot.f90 @@ -73,7 +73,7 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R '|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|' if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then backward=.true. - chem_pot = chem_pot - 2d0*delta_chem_pot + chem_pot = chem_pot - delta_chem_pot delta_chem_pot=-delta_chem_pot endif enddo From 799d44324fb25b111e9ebb27dcea1b7d6d06e396 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Thu, 6 Feb 2025 22:41:35 +0100 Subject: [PATCH 30/32] Solved RESTART reading bug --- src/HF/HFB.f90 | 14 ++++++++------ src/HF/read_restart_HFB.f90 | 10 +++++++++- src/HF/write_restart_HFB.f90 | 11 +++++++++-- 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 0279f67b..518c88da 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -144,7 +144,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, do iorb=1,nOrb P(:,:) = P(:,:) + Occ(iorb) * & matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) - Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & + Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb)))) * & matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) enddo endif @@ -158,7 +158,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, do iorb=1,nOrb P(:,:) = P(:,:) + Occ(iorb) * & matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) - Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * & + Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb)))) * & matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) enddo endif @@ -343,11 +343,13 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Compute dipole moments, occupation numbers and || Anomalous density|| ! also print the restart file - eigVEC(:,:) = 0d0 + deallocate(eigVEC,eigVAL) + allocate(eigVEC(nOrb,nOrb),eigVAL(nOrb)) + eigVEC(1:nOrb,1:nOrb) = 0d0 eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb) - call diagonalize_matrix(nOrb2,eigVEC,eigVAL) - Occ(1:nOrb) = eigVAL(nOrb+1:nOrb2) - c(1:nBas,1:nOrb) = matmul(X(1:nBas,1:nOrb),eigVEC(1:nOrb,nOrb+1:nOrb2)) + call diagonalize_matrix(nOrb,eigVEC,eigVAL) + Occ(1:nOrb) = eigVAL(1:nOrb) + c = matmul(X,eigVEC) norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2))) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) call write_restart_HFB(nBas,nOrb,Occ,c,chem_pot) diff --git a/src/HF/read_restart_HFB.f90 b/src/HF/read_restart_HFB.f90 index 8d523c2f..5fb09ffb 100644 --- a/src/HF/read_restart_HFB.f90 +++ b/src/HF/read_restart_HFB.f90 @@ -22,6 +22,7 @@ subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot) integer :: nOrb_ double precision :: chem_pot_ double precision :: max_diff + double precision :: val_read double precision,allocatable :: eigVAL(:) double precision,allocatable :: c_tmp(:,:) double precision,allocatable :: S_mol(:,:) @@ -37,11 +38,18 @@ subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot) allocate(eigVAL(nOrb),c_tmp(nBas,nOrb),S_mol(nOrb,nOrb),X_mol(nOrb,nOrb)) + c_tmp=0d0 + S_mol=0d0 + X_mol=0d0 + open(unit=iunit,form='unformatted',file='hfb_bin',status='old') read(iunit) nBas_,nOrb_ read(iunit) chem_pot_ do iorb=1,nOrb - read(iunit) c_tmp(:,iorb) + do ibas=1,nBas + read(iunit) val_read + c_tmp(ibas,iorb) = val_read + enddo enddo do iorb=1,nOrb read(iunit) eigVAL(iorb) diff --git a/src/HF/write_restart_HFB.f90 b/src/HF/write_restart_HFB.f90 index bcd7ef71..d5b4f97d 100644 --- a/src/HF/write_restart_HFB.f90 +++ b/src/HF/write_restart_HFB.f90 @@ -19,6 +19,7 @@ subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot) ! Local variables integer :: ibas,iorb,iunit=666 + double precision :: val_print ! Dump results @@ -26,10 +27,16 @@ subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot) write(iunit) nBas,nOrb write(iunit) chem_pot do iorb=1,nOrb - write(iunit) c(:,nOrb+1-iorb) + do ibas=1,nBas + val_print = c(ibas,nOrb+1-iorb) + if(abs(val_print)<1d-8) val_print=0d0 + write(iunit) val_print + enddo enddo do iorb=1,nOrb - write(iunit) Occ(nOrb+1-iorb) + val_print = Occ(nOrb+1-iorb) + if(abs(val_print)<1d-8) val_print=0d0 + write(iunit) val_print enddo write(iunit) iunit close(iunit) From 539009a05bf1668431ea5adbb25138359826cb0d Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Fri, 7 Feb 2025 10:20:20 +0100 Subject: [PATCH 31/32] Removed unused test --- src/HF/HFB.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 518c88da..35b3c1ba 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -360,8 +360,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, 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 From bc56e7f88ea4bce341c2eb9a16c3384585985c30 Mon Sep 17 00:00:00 2001 From: Mauricio Rodriguez-Mayorga Date: Wed, 12 Feb 2025 22:04:18 +0100 Subject: [PATCH 32/32] eHF is intent_in --- src/HF/HFB.f90 | 2 +- src/HF/fermi_dirac_occ.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index 35b3c1ba..673a1bc7 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -24,6 +24,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc double precision,intent(in) :: temperature,sigma + double precision,intent(in) :: eHF(nOrb) double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: T(nBas,nBas) double precision,intent(in) :: V(nBas,nBas) @@ -75,7 +76,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Output variables double precision,intent(out) :: EHFB - double precision,intent(inout):: eHF(nOrb) double precision,intent(inout):: c(nBas,nOrb) double precision,intent(out) :: P(nBas,nBas) double precision,intent(out) :: Panom(nBas,nBas) diff --git a/src/HF/fermi_dirac_occ.f90 b/src/HF/fermi_dirac_occ.f90 index 4af8cbfc..f6ed6ab7 100644 --- a/src/HF/fermi_dirac_occ.f90 +++ b/src/HF/fermi_dirac_occ.f90 @@ -9,6 +9,7 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) integer,intent(in) :: nO,nOrb double precision,intent(in) :: thrs_N double precision,intent(in) :: temperature + double precision,intent(in) :: eHF(nOrb) ! Local variables @@ -24,7 +25,6 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) ! Output variables double precision,intent(inout):: chem_pot - double precision,intent(inout):: eHF(nOrb) double precision,intent(inout):: Occ(nOrb) ! Initialize variables