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

Incl up to reshufle after MOM

Browse Source
This commit is contained in:
Mauricio Rodriguez-Mayorga 2025-01-28 16:52:17 +01:00
parent fde7d42202
commit f688dd6cec
2 changed files with 168 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

113
src/HF/HFB.f90
View File

@ -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
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
! ! 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
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb,cp,eHF)
c = matmul(X,cp)
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(:,:) = -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

90
src/HF/reshufle_hfb.f90 Normal file
View File

@ -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)<val) then
val=eHFB_tmp(iorb2)
ipos=iorb2
endif
enddo
eHFB_tmp(ipos)=2d10
eHFB(iorb)=val
project(iorb)=project_tmp(ipos)
c_hfb(:,iorb)=c_tmp(:,ipos)
enddo
do iorb=nOrb2/2+1,nOrb2 ! hole 1-RDM
val=1d10
ipos=0
do iorb2=nOrb2/2+1,nOrb2
if(eHFB_tmp(iorb2)<val) then
val=eHFB_tmp(iorb2)
ipos=iorb2
endif
enddo
eHFB_tmp(ipos)=2d10
eHFB(iorb)=val
project(iorb)=project_tmp(ipos)
c_hfb(:,iorb)=c_tmp(:,ipos)
enddo
!write(*,*) 'e_HFB'
!do iorb=1,nOrb2/2
!write(*,*) eHFB(iorb),project(iorb)
!enddo
!write(*,*) '------'
!do iorb=nOrb2/2+1,nOrb2
!write(*,*) eHFB(iorb),project(iorb)
!enddo
deallocate(c_tmp,project_tmp,eHFB_tmp)
end subroutine