PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-01 06:21:37 +02:00
Code Releases Activity

Merge pull request #12 from marm314/master

Browse Source 
HFB SCF
This commit is contained in:
Pierre-Francois Loos 2025-02-12 22:12:31 +01:00 committed by GitHub
commit 810ab896d2
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
17 changed files with 1142 additions and 69 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored
View File

@ -6,3 +6,8 @@
__pycache__
.ninja_deps
calcs/
calcs/input
calcs/int
calcs/mol
calcs/water.out

2
PyDuck.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python
#!/usr/bin/env python3
import os, sys
import argparse

6
src/AOtoMO/Hartree_matrix_AO_basis.f90
View File

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

38
src/AOtoMO/anomalous_matrix_AO_basis.f90 Normal file
View File

@ -0,0 +1,38 @@
subroutine anomalous_matrix_AO_basis(nBas,sigma,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) :: sigma
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) + sigma*Pa(la,si)*ERI(la,si,mu,nu)
end do
end do
end do
end do
end subroutine
! ---

3
src/AOtoMO/exchange_matrix_AO_basis.f90
View File

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

373
src/HF/HFB.f90 Normal file
View File

@ -0,0 +1,373 @@
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,Panom,F,Delta, &
temperature,sigma,chem_pot_hf,restart_hfb)
! Perform Hartree-Fock Bogoliubov calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
double precision,intent(in) :: temperature,sigma
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
logical :: chem_pot_hf
logical :: restart_hfb
integer :: nBas2
integer :: iorb
integer :: nSCF
integer :: nOrb2
integer :: nBas2_Sq
integer :: n_diis
double precision :: ET
double precision :: EV
double precision :: EJ
double precision :: EK
double precision :: EL
double precision :: chem_pot
double precision :: dipole(ncart)
double precision :: Conv
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(:)
double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: H_HFB_diis(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: eigVEC(:,:)
double precision,allocatable :: H_HFB(:,:)
double precision,allocatable :: R(:,:)
double precision,allocatable :: err_ao(:,:)
double precision,allocatable :: S_ao(:,:)
double precision,allocatable :: X_ao(:,:)
double precision,allocatable :: R_ao_old(:,:)
double precision,allocatable :: H_HFB_ao(:,:)
! Output variables
double precision,intent(out) :: EHFB
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)
double precision,intent(out) :: Delta(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'*****************************'
write(*,*)'* HF Bogoliubov Calculation *'
write(*,*)'*****************************'
write(*,*)
! Useful quantities
nOrb2 = nOrb+nOrb
nBas2 = nBas+nBas
nBas2_Sq = nBas2*nBas2
! Memory allocation
allocate(Occ(nOrb))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(eigVEC(nOrb2,nOrb2))
allocate(H_HFB(nOrb2,nOrb2))
allocate(R(nOrb2,nOrb2))
allocate(eigVAL(nOrb2))
allocate(err_ao(nBas2,nBas2))
allocate(S_ao(nBas2,nBas2))
allocate(X_ao(nBas2,nOrb2))
allocate(R_ao_old(nBas2,nBas2))
allocate(H_HFB_ao(nBas2,nBas2))
allocate(err_diis(nBas2_Sq,max_diis))
allocate(H_HFB_diis(nBas2_Sq,max_diis))
! Guess chem. pot.
chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1))
! Initialization
thrs_N = 1d-8
n_diis = 0
H_HFB_diis(:,:) = 0d0
err_diis(:,:) = 0d0
rcond = 0d0
P(:,:) = matmul(c(:,1:nO), transpose(c(:,1:nO)))
Panom(:,:) = 0d0
! Use Fermi-Dirac occupancies to compute P, Panom, and chem_pot
if(abs(temperature)>1d-4) then
Occ(:) = 0d0
Occ(1:nO) = 1d0
call fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
if(chem_pot_hf) chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1))
P(:,:) = 0d0
Panom(:,:) = 0d0
do iorb=1,nOrb
P(:,:) = P(:,:) + Occ(iorb) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb)))) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
enddo
endif
! Read restart file
if(restart_hfb) then
call read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
P(:,:) = 0d0
Panom(:,:) = 0d0
do iorb=1,nOrb
P(:,:) = P(:,:) + Occ(iorb) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb)))) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
enddo
endif
P(:,:) = 2d0*P(:,:)
S_ao(:,:) = 0d0
S_ao(1:nBas ,1:nBas ) = S(1:nBas,1:nBas)
S_ao(nBas+1:nBas2,nBas+1:nBas2) = S(1:nBas,1:nBas)
X_ao(:,:) = 0d0
X_ao(1:nBas ,1:nOrb ) = X(1:nBas,1:nOrb)
X_ao(nBas+1:nBas2,nOrb+1:nOrb2) = X(1:nBas,1:nOrb)
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*) 'Enterning HFB SCF procedure'
write(*,*)
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Build Fock and Delta matrices
call Hartree_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
call anomalous_matrix_AO_basis(nBas,sigma,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 ) = H_HFB(1:nOrb,nOrb+1:nOrb2)
eigVEC(:,:) = H_HFB(:,:)
call diagonalize_matrix(nOrb2,eigVEC,eigVAL)
! Build R
R(:,:) = 0d0
do iorb=1,nOrb
R(:,:) = R(:,:) + matmul(eigVEC(:,iorb:iorb),transpose(eigVEC(:,iorb:iorb)))
enddo
trace_1rdm = 0d0
do iorb=1,nOrb
trace_1rdm = trace_1rdm+R(iorb,iorb)
enddo
! Adjust the chemical potential
if( abs(trace_1rdm-nO) > thrs_N ) &
call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_HFB,eigVEC,R,eigVAL)
! DIIS extrapolation
if(max_diis > 1 .and. nSCF>1) then
write(*,*) ' Doing DIIS'
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:)
H_HFB_ao(:,:) = 0d0
H_HFB_ao(1:nBas ,1:nBas ) = F(1:nBas,1:nBas)
H_HFB_ao(nBas+1:nBas2,nBas+1:nBas2) = -F(1:nBas,1:nBas)
H_HFB_ao(1:nBas ,nBas+1:nBas2) = Delta(1:nBas,1:nBas)
H_HFB_ao(nBas+1:nBas2,1:nBas ) = Delta(1:nBas,1:nBas)
err_ao = matmul(H_HFB_ao,matmul(R_ao_old,S_ao)) - matmul(matmul(S_ao,R_ao_old),H_HFB_ao)
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas2_Sq,nBas2_Sq,n_diis,err_diis,H_HFB_diis,err_ao,H_HFB_ao)
H_HFB = matmul(transpose(X_ao),matmul(H_HFB_ao,X_ao))
eigVEC(:,:) = H_HFB(:,:)
call diagonalize_matrix(nOrb2,eigVEC,eigVAL)
! Build R and check trace
trace_1rdm = 0d0
R(:,:) = 0d0
do iorb=1,nOrb
R(:,:) = R(:,:) + matmul(eigVEC(:,iorb:iorb),transpose(eigVEC(:,iorb:iorb)))
enddo
do iorb=1,nOrb
trace_1rdm = trace_1rdm + R(iorb,iorb)
enddo
! Adjust the chemical potential
if( abs(trace_1rdm-nO) > thrs_N ) &
call fix_chem_pot(nO,nOrb,nOrb2,nSCF,thrs_N,trace_1rdm,chem_pot,H_HFB,eigVEC,R,eigVAL)
end if
! 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)))
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
! Anomalous energy
EL = trace_matrix(nBas,matmul(Panom,Delta))
! Total energy
EHFB = ET + EV + EJ + EK + EL
! Check convergence
if(nSCF > 1) then
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) - chem_pot*S(:,:)
H_HFB_ao(:,:) = 0d0
H_HFB_ao(1:nBas ,1:nBas ) = F(1:nBas,1:nBas)
H_HFB_ao(nBas+1:nBas2,nBas+1:nBas2) = -F(1:nBas,1:nBas)
H_HFB_ao(1:nBas ,nBas+1:nBas2) = Delta(1:nBas,1:nBas)
H_HFB_ao(nBas+1:nBas2,1:nBas ) = Delta(1:nBas,1:nBas)
err_ao = matmul(H_HFB_ao,matmul(R_ao_old,S_ao)) - matmul(matmul(S_ao,R_ao_old),H_HFB_ao)
Conv = maxval(abs(err_ao))
endif
! Update R_old
R_ao_old(:,:) = 0d0
R_ao_old(1:nBas ,1:nBas ) = 0.5d0*P(1:nBas,1:nBas)
R_ao_old(nBas+1:nBas2,nBas+1:nBas2) = matmul(X(1:nBas,1:nOrb), transpose(X(1:nBas,1:nOrb)))-0.5d0*P(1:nBas,1:nBas)
R_ao_old(1:nBas ,nBas+1:nBas2) = Panom(1:nBas,1:nBas)
R_ao_old(nBas+1:nBas2,1:nBas ) = Panom(1:nBas,1:nBas)
! 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(*,*)
end do
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis,Occ)
deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao)
stop
end if
! Compute dipole moments, occupation numbers and || Anomalous density||
! also print the restart file
deallocate(eigVEC,eigVAL)
allocate(eigVEC(nOrb,nOrb),eigVAL(nOrb))
eigVEC(1:nOrb,1:nOrb) = 0d0
eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb)
call diagonalize_matrix(nOrb,eigVEC,eigVAL)
Occ(1:nOrb) = eigVAL(1:nOrb)
c = matmul(X,eigVEC)
norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2)))
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call write_restart_HFB(nBas,nOrb,Occ,c,chem_pot)
call print_HFB(nBas,nOrb,nO,norm_anom,Occ,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole)
! Testing zone
if(dotest) then
call dump_test_value('R','HFB energy',EHFB)
call dump_test_value('R','HFB dipole moment',norm2(dipole))
end if
! Memory deallocation
deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis,Occ)
deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao)
end subroutine

52
src/HF/MOM_overlap.f90
View File

@ -1,52 +0,0 @@
subroutine MOM_overlap(nBas,nO,S,cG,c,ON)
! Compute overlap between old and new MO coefficients
implicit none
! Input variables
integer,intent(in) :: nBas,nO
double precision,intent(in) :: S(nBas,nBas),cG(nBas,nBas),c(nBas,nBas)
! Local variables
integer :: i,j,ploc
double precision,allocatable :: Ov(:,:),pOv(:)
! Output variables
double precision,intent(inout):: ON(nBas)
allocate(Ov(nBas,nBas),pOv(nBas))
Ov = matmul(transpose(cG),matmul(S,c))
pOv(:) = 0d0
do i=1,nBas
do j=1,nBas
pOv(j) = pOv(j) + Ov(i,j)**2
end do
end do
pOv(:) = sqrt(pOV(:))
! print*,'--- MOM overlap ---'
! call matout(nBas,1,pOv)
ON(:) = 0d0
do i=1,nO
ploc = maxloc(pOv,nBas)
ON(ploc) = 1d0
pOv(ploc) = 0d0
end do
! print*,'--- Occupation numbers ---'
! call matout(nBas,1,ON)
end subroutine

110
src/HF/fermi_dirac_occ.f90 Normal file
View File

@ -0,0 +1,110 @@
subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
! Use Fermi Dirac distribution to set up fractional Occs numbers and adjust the chemical potential to integrate to the N for 2N electron systems
implicit none
! Input variables
integer,intent(in) :: nO,nOrb
double precision,intent(in) :: thrs_N
double precision,intent(in) :: temperature
double precision,intent(in) :: eHF(nOrb)
! Local variables
logical :: backward
integer :: iorb
integer :: isteps
double precision :: delta_chem_pot
double precision :: chem_pot_change
double precision :: grad_electrons
double precision :: trace_1rdm
double precision :: trace_old
! Output variables
double precision,intent(inout):: chem_pot
double precision,intent(inout):: Occ(nOrb)
! Initialize variables
backward = .false.
isteps = 0
delta_chem_pot = 1.0d-1
chem_pot_change = 0d0
grad_electrons = 1d0
trace_1rdm = -1d0
write(*,*)
write(*,*)' Fermi-Dirac distribution for the occ numbers'
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,'|'
! First approach close the value with an error lower than 1
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_old=sum(Occ(:))
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_old,'|',chem_pot,'|'
do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + delta_chem_pot
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_1rdm=sum(Occ(:))
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|'
if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then
backward=.true.
chem_pot = chem_pot - delta_chem_pot
delta_chem_pot=-delta_chem_pot
endif
enddo
! Do final search
write(*,*)'-------------------------------------'
isteps=0
delta_chem_pot = 1.0d-1
do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + chem_pot_change
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_1rdm = sum(Occ(:))
grad_electrons = ( sum(fermi_dirac(eHF,chem_pot+delta_chem_pot,temperature)) &
- sum(fermi_dirac(eHF,chem_pot-delta_chem_pot,temperature)) )/(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,A1)') &
'|',trace_1rdm,'|',chem_pot,'|'
enddo
write(*,*)'-------------------------------------'
write(*,*)
write(*,*)
write(*,*) ' Initial occ. numbers '
write(*,*)
do iorb=1,nOrb
write(*,'(3X,F16.10)') Occ(iorb)
enddo
contains
function fermi_dirac(eHF,chem_pot,temperature)
implicit none
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: chem_pot
double precision,intent(in) :: temperature
double precision :: fermi_dirac(nOrb)
fermi_dirac(:) = 1d0 / ( 1d0 + exp((eHF(:) - chem_pot ) / temperature ) )
end function fermi_dirac
end subroutine

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

@ -0,0 +1,159 @@
subroutine fix_chem_pot(nO,nOrb,nOrb2,nSCF,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,nSCF
double precision,intent(in) :: thrs_N
! Local variables
logical :: backward
integer :: iorb
integer :: isteps
double precision :: delta_chem_pot
double precision :: chem_pot_change
double precision :: grad_electrons
double precision :: trace_2up
double precision :: trace_up
double precision :: trace_down
double precision :: trace_2down
double precision :: trace_old
double precision,allocatable :: R_tmp(:,:)
double precision,allocatable :: cp_tmp(:,:)
! 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
backward = .false.
isteps = 0
delta_chem_pot = 1.0d-1
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
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(*,'(a,I5)') ' Fixing the Tr[1D] at iteration ',nSCF
write(*,*)
write(*,*)'------------------------------------------------------'
write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1A15,2X,A1)') &
'|','Tr[1D]','|','Chem. Pot.','|','Grad N','|'
write(*,*)'------------------------------------------------------'
! First approach close the value with an error lower than 1
call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_old,H_hfb,cp,R,eHFB_)
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') &
'|',trace_old,'|',chem_pot,'|',grad_electrons,'|'
do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + delta_chem_pot
call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_)
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|'
if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then
backward=.true.
chem_pot = chem_pot - delta_chem_pot
delta_chem_pot=-delta_chem_pot
endif
enddo
! Do final search
write(*,*)'------------------------------------------------------'
isteps = 0
delta_chem_pot = 1.0d-3
do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 )
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+2d0*delta_chem_pot,trace_2up,H_hfb,cp_tmp,R_tmp,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_)
call diag_H_hfb(nOrb,nOrb2,chem_pot-2d0*delta_chem_pot,trace_2down,H_hfb,cp_tmp,R_tmp,eHFB_)
! grad_electrons = (trace_up-trace_down)/(2d0*delta_chem_pot)
grad_electrons = (-trace_2up+8d0*trace_up-8d0*trace_down+trace_2down)/(12d0*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)') &
'|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|'
enddo
write(*,*)'------------------------------------------------------'
write(*,*)
! 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
deallocate(R_tmp,cp_tmp)
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

79
src/HF/print_HFB.f90 Normal file
View File

@ -0,0 +1,79 @@
! ---
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
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
double precision,intent(in) :: Occ(nOrb)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET
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) :: 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.
! Dump results
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A33)') ' Summary '
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',ET + EV,' au'
write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',ET,' au'
write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',EV,' au'
write(*,'(A50)') '---------------------------------------'
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'
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.'
write(*,'(10X,4F10.4)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
write(*,'(A50)') '---------------------------------------'
write(*,*)
! Print results
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' HFB occupation numbers '
write(*,'(A50)') '---------------------------------------'
trace_occ=0d0
do iorb=1,nOrb
if(abs(Occ(nOrb+1-iorb))>1d-8) then
write(*,'(I7,10F15.8)') iorb,2d0*Occ(nOrb+1-iorb)
endif
trace_occ=trace_occ+2d0*Occ(iorb)
enddo
write(*,*)
write(*,'(A33,1X,F16.10,A3)') ' Trace [ 1D ] = ',trace_occ,' '
write(*,*)
end subroutine

124
src/HF/read_restart_HFB.f90 Normal file
View File

@ -0,0 +1,124 @@
! ---
subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
! Write the binary file used to restart calcs
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
double precision,intent(in) :: S(nBas,nBas)
! Local variables
integer :: ibas,iorb,iorb1,iunit=667
integer :: nBas_
integer :: nOrb_
double precision :: chem_pot_
double precision :: max_diff
double precision :: val_read
double precision,allocatable :: eigVAL(:)
double precision,allocatable :: c_tmp(:,:)
double precision,allocatable :: S_mol(:,:)
double precision,allocatable :: X_mol(:,:)
! Output variables
double precision,intent(out) :: chem_pot
double precision,intent(out) :: Occ(nOrb)
double precision,intent(out) :: c(nBas,nOrb)
! Dump results
allocate(eigVAL(nOrb),c_tmp(nBas,nOrb),S_mol(nOrb,nOrb),X_mol(nOrb,nOrb))
c_tmp=0d0
S_mol=0d0
X_mol=0d0
open(unit=iunit,form='unformatted',file='hfb_bin',status='old')
read(iunit) nBas_,nOrb_
read(iunit) chem_pot_
do iorb=1,nOrb
do ibas=1,nBas
read(iunit) val_read
c_tmp(ibas,iorb) = val_read
enddo
enddo
do iorb=1,nOrb
read(iunit) eigVAL(iorb)
enddo
close(iunit)
if(nBas==nBas_ .and. nOrb==nOrb_) then
write(*,*)
write(*,*)' Reading restart file'
write(*,*)
chem_pot = chem_pot_
Occ(:) = eigVAL(:)
c(:,:) = c_tmp(:,:)
! Check the orthonormality
max_diff=-1d0
S_mol = matmul(transpose(c),matmul(S,c))
do iorb=1,nOrb
do iorb1=1,iorb
if(iorb==iorb1) then
if(abs(S_mol(iorb,iorb)-1d0)>1d-8) max_diff=abs(S_mol(iorb,iorb)-1d0)
else
if(abs(S_mol(iorb,iorb1))>1d-8) max_diff=abs(S_mol(iorb,iorb1))
endif
enddo
enddo
if(max_diff>1d-8) then
write(*,*) ' '
write(*,'(a)') ' Orthonormality violations. Applying Lowdin orthonormalization.'
write(*,*) ' '
X_mol = S_mol
S_mol = 0d0
call diagonalize_matrix(nOrb,X_mol,eigVAL)
do iorb=1,nOrb
S_mol(iorb,iorb) = 1d0/(sqrt(eigVAL(iorb)) + 1d-10)
enddo
X_mol = matmul(X_mol,matmul(S_mol,transpose(X_mol)))
c = matmul(c,X_mol)
max_diff=-1d0
S_mol = matmul(transpose(c),matmul(S,c))
do iorb=1,nOrb
do iorb1=1,iorb
if(iorb==iorb1) then
if(abs(S_mol(iorb,iorb)-1d0)>1d-8) max_diff=abs(S_mol(iorb,iorb)-1d0)
else
if(abs(S_mol(iorb,iorb1))>1d-8) max_diff=abs(S_mol(iorb,iorb1))
endif
enddo
enddo
if(max_diff<=1d-8) then
write(*,*) ' '
write(*,'(a)') ' No orthonormality violations.'
write(*,*) ' '
else
write(*,*) ' '
write(*,'(a)') ' Error in Lowdin orthonormalization.'
write(*,*) ' '
stop
endif
else
write(*,*) ' '
write(*,'(a)') ' No orthonormality violations.'
write(*,*) ' '
endif
endif
deallocate(eigVAL,c_tmp,S_mol,X_mol)
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

44
src/HF/write_restart_HFB.f90 Normal file
View File

@ -0,0 +1,44 @@
! ---
subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot)
! Write the binary file used to restart calcs
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
double precision,intent(in) :: chem_pot
double precision,intent(in) :: Occ(nOrb)
double precision,intent(in) :: c(nBas,nOrb)
! Local variables
integer :: ibas,iorb,iunit=666
double precision :: val_print
! Dump results
open(unit=iunit,form='unformatted',file='hfb_bin')
write(iunit) nBas,nOrb
write(iunit) chem_pot
do iorb=1,nOrb
do ibas=1,nBas
val_print = c(ibas,nOrb+1-iorb)
if(abs(val_print)<1d-8) val_print=0d0
write(iunit) val_print
enddo
enddo
do iorb=1,nOrb
val_print = Occ(nOrb+1-iorb)
if(abs(val_print)<1d-8) val_print=0d0
write(iunit) val_print
enddo
write(iunit) iunit
close(iunit)
end subroutine

84
src/QuAcK/BQuAcK.f90
View File

@ -1,4 +1,6 @@
subroutine BQuAcK(working_dir,dotest,doHFB)
subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,temperature,sigma,chem_pot_hf,restart_hfb)
! Restricted branch of QuAcK
@ -11,10 +13,43 @@ subroutine BQuAcK(working_dir,dotest,doHFB)
logical,intent(in) :: doHFB
logical,intent(in) :: restart_hfb
logical,intent(in) :: chem_pot_hf
integer,intent(in) :: nNuc,nBas,nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: temperature,sigma
double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix
integer,intent(in) :: guess_type
! Local variables
double precision :: start_HF ,end_HF ,t_HF
double precision :: start_int, end_int, t_int
double precision,allocatable :: eHF(:)
double precision,allocatable :: cHF(:,:)
double precision,allocatable :: PHF(:,:)
double precision,allocatable :: PanomHF(:,:)
double precision,allocatable :: FHF(:,:)
double precision,allocatable :: Delta(:,:)
double precision :: ERHF,EHFB
double precision,allocatable :: ERI_AO(:,:,:,:)
write(*,*)
write(*,*) '******************************'
write(*,*) '* Bogoliubov Branch of QuAcK *'
@ -25,14 +60,43 @@ subroutine BQuAcK(working_dir,dotest,doHFB)
! Memory allocation !
!-------------------!
!---------------------!
! Hartree-Fock module !
!---------------------!
allocate(eHF(nOrb))
allocate(cHF(nBas,nOrb))
allocate(PHF(nBas,nBas))
allocate(PanomHF(nBas,nBas))
allocate(FHF(nBas,nBas))
allocate(Delta(nBas,nBas))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
call wall_time(end_int)
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds'
write(*,*)
!--------------------------------!
! Hartree-Fock Bogoliubov module !
!--------------------------------!
if(doHFB) then
! Run first a RHF calculation
call wall_time(start_HF)
! call HFB(dotest)
call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
write(*,*)
! 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,PanomHF, &
FHF,Delta,temperature,sigma,chem_pot_hf,restart_hfb)
call wall_time(end_HF)
t_HF = end_HF - start_HF
@ -41,4 +105,14 @@ subroutine BQuAcK(working_dir,dotest,doHFB)
end if
! Memory deallocation
deallocate(eHF)
deallocate(cHF)
deallocate(PHF)
deallocate(PanomHF)
deallocate(FHF)
deallocate(Delta)
deallocate(ERI_AO)
end subroutine

13
src/QuAcK/QuAcK.f90
View File

@ -73,6 +73,10 @@ program QuAcK
logical :: dotest,doRtest,doUtest,doGtest
logical :: chem_pot_hf
logical :: restart_hfb
double precision :: temperature,sigma
character(len=256) :: working_dir
! Check if the right number of arguments is provided
@ -108,7 +112,7 @@ program QuAcK
!------------------!
call read_methods(working_dir, &
doRHF,doUHF,doGHF,doROHF, &
doRHF,doUHF,doGHF,doROHF,doHFB, &
doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD, &
@ -134,7 +138,8 @@ program QuAcK
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
dophBSE,dophBSE2,doppBSE,dBSE,dTDA, &
temperature,sigma,chem_pot_hf,restart_hfb)
!------------------!
! Hardware !
@ -289,7 +294,9 @@ program QuAcK
! Bogoliubov QuAcK branch !
!--------------------------!
if(doBQuAcK) &
call BQuAcK(working_dir,doGtest,doHFB)
call BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix, &
temperature,sigma,chem_pot_hf,restart_hfb)
!-----------!
! Stop Test !

8
src/QuAcK/read_methods.f90
View File

@ -1,5 +1,5 @@
subroutine read_methods(working_dir, &
doRHF,doUHF,doGHF,doROHF, &
doRHF,doUHF,doGHF,doROHF,doHFB, &
doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
@ -22,7 +22,7 @@ subroutine read_methods(working_dir, &
! Output variables
logical,intent(out) :: doRHF,doUHF,doGHF,doROHF
logical,intent(out) :: doRHF,doUHF,doGHF,doROHF,doHFB
logical,intent(out) :: doMP2,doMP3
logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
@ -58,13 +58,15 @@ subroutine read_methods(working_dir, &
doUHF = .false.
doGHF = .false.
doROHF = .false.
doHFB = .false.
read(1,*)
read(1,*) ans1,ans2,ans3,ans4
read(1,*) ans1,ans2,ans3,ans4,ans5
if(ans1 == 'T') doRHF = .true.
if(ans2 == 'T') doUHF = .true.
if(ans3 == 'T') doGHF = .true.
if(ans4 == 'T') doROHF = .true.
if(ans5 == 'T') doHFB = .true.
! Read MPn methods

21
src/QuAcK/read_options.f90
View File

@ -7,7 +7,8 @@ subroutine read_options(working_dir,
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
dophBSE,dophBSE2,doppBSE,dBSE,dTDA, &
temperature,sigma,chem_pot_hf,restart_hfb)
! Read desired methods
@ -72,6 +73,11 @@ subroutine read_options(working_dir,
logical,intent(out) :: dBSE
logical,intent(out) :: dTDA
logical,intent(out) :: chem_pot_hf
logical,intent(out) :: restart_hfb
double precision,intent(out) :: temperature
double precision,intent(out) :: sigma
! Local variables
character(len=1) :: ans1,ans2,ans3,ans4,ans5
@ -217,6 +223,19 @@ subroutine read_options(working_dir,
if(ans4 == 'T') dBSE = .true.
if(ans5 == 'F') dTDA = .false.
! Options for Hartree-Fock Bogoliubov
temperature = 0d0
sigma = 1d0
chem_pot_hf = .false.
restart_hfb = .false.
read(1,*)
read(1,*) temperature,sigma,ans1,ans2
if(ans1 == 'T') chem_pot_hf = .true.
if(ans2 == 'T') restart_hfb = .true.
endif
! Close file with options