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

Solved 1 bug chem pot fix

Browse Source
This commit is contained in:
Mauricio Rodriguez-Mayorga 2025-01-30 22:27:59 +01:00
parent c1a10548d5
commit fa0981508e
3 changed files with 35 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
src/HF/HFB.f90
View File

@ -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

10
src/HF/fix_chem_pot.f90
View File

@ -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

16
src/HF/print_HFB.f90
View File

@ -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