LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-03-09 18:22:25 +01:00
Code Releases Activity

Todo golden-section search for mu

Browse Source
This commit is contained in:
Mauricio Rodriguez-Mayorga 2025-01-29 15:52:04 +01:00
parent c85fa22a98
commit 94062e66cd
2 changed files with 164 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
src/HF/HFB.f90
View File

@ -134,6 +134,8 @@ do while(Conv > 1e-7 .and. nSCF < maxSCF) ! TODO
nSCF = nSCF + 1
! TODO remove
chem_pot=-10
! Build Fock and Delta matrices
@ -169,7 +171,7 @@ do while(Conv > 1e-7 .and. nSCF < maxSCF) ! TODO
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 @@ do while(Conv > 1e-7 .and. nSCF < maxSCF) ! TODO
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 @@ end subroutine
! 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

150
src/HF/fix_chem_pot.f90 Normal file
View File

@ -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) then
! Increase chem_pot to occupy more orbs.
do
chem_pot = chem_pot + delta_chem_pot
call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_up,H_hfb,cp,R,eHFB_)
if(trace_up>=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