mirror of
https://github.com/pfloos/quack
synced 2025-03-09 18:22:25 +01:00
Building R instead of MOM
This commit is contained in:
Mauricio Rodriguez-Mayorga
parent
f7b8a7ce4e
commit
2cf4aa901b
187
src/HF/HFB.f90
187
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)
|
||||
|
||||
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:)
|
||||
write(*,*)
|
||||
write(*,*)'-------------------------------------'
|
||||
write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') &
|
||||
'|','Tr[1D]','|','Chem. Pot.','|'
|
||||
write(*,*)'-------------------------------------'
|
||||
|
||||
! Check convergence
|
||||
! 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
|
||||
|
||||
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 @@ eHF(:)=eHF(:)-chem_pot!TODO remove
|
||||
|
||||
! 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
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user