mirror of
https://github.com/pfloos/quack
synced 2025-04-02 06:51:37 +02:00
Starting Panomalous
This commit is contained in:
Mauricio Rodriguez-Mayorga
parent
cf1ad337e5
commit
79fada93da
@ -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
|
||||
|
||||
37
src/AOtoMO/anomalous_matrix_AO_basis.f90
Normal file
37
src/AOtoMO/anomalous_matrix_AO_basis.f90
Normal file
@ -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
|
||||
|
||||
! ---
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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)') ' <Sz> = ',S
|
||||
write(*,'(A33,1X,F16.6)') ' <S^2> = ',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
|
||||
|
||||
@ -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)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user