PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-01-03 01:55:57 +01:00
Code Releases Activity

restructuration QuAcK for restricted and unrestricted branches

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-10-25 21:38:20 +02:00
parent 4d8b25dfea
commit 0b892dbe31
29 changed files with 1752 additions and 926 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
input/methods
View File

@ -1,5 +1,5 @@
# RHF UHF ROHF RMOM UMOM KS # RHF UHF GHF ROHF
F T F F F F T F F F
# MP2* MP3 # MP2* MP3
F F F F
# CCD pCCD DCD CCSD CCSD(T) # CCD pCCD DCD CCSD CCSD(T)
@ -13,7 +13,7 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3 # G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F F F F F F
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
F T F F F F T F F F F F
# G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh # G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
F F F F F F F F F F F F
# * unrestricted version available # * unrestricted version available

24
input/options
View File

@ -1,17 +1,17 @@
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability # HF: maxSCF thresh DIIS guess mix_guess level_shift stability
5000 0.0000001 T 5 1 1 T 0.0 T 5000 0.0000001 5 1 0.0 0.0 T
# MP: reg # MP: reg
F F
# CC: maxSCF thresh DIIS n_diis # CC: maxSCF thresh DIIS
64 0.0000001 T 5 64 0.0000001 5
# spin: TDA singlet triplet spin_conserved spin_flip # spin: TDA spin_conserved spin_flip
F T T T T F T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm reg # GF: maxSCF thresh DIIS lin eta renorm reg
256 0.00001 T 5 T 0.0 0 F 256 0.00001 5 F 0.0 0 F
# GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg # GW: maxSCF thresh DIIS lin eta TDA_W reg
256 0.00001 T 5 F 0.0 F F 256 0.00001 5 F 0.0 F F
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg # GT: maxSCF thresh DIIS lin eta TDA_T reg
256 0.00001 T 5 F 0.0 F F 256 0.00001 5 F 0.0 F F
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# BSE: phBSE phBSE2 ppBSE dBSE dTDA # BSE: phBSE phBSE2 ppBSE dBSE dTDA

4
src/CC/CC.f90 → src/CC/RCC.f90
View File

@ -1,5 +1,5 @@
subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
! Coupled-cluster module ! Coupled-cluster module

44
src/CI/CID.f90
View File

@ -1,4 +1,4 @@
subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,Fin,E0) subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,eHFin,E0)
! Perform configuration interaction with doubles ! Perform configuration interaction with doubles
@ -14,7 +14,7 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi
integer,intent(in) :: nOin integer,intent(in) :: nOin
integer,intent(in) :: nVin integer,intent(in) :: nVin
integer,intent(in) :: nRin integer,intent(in) :: nRin
double precision,intent(in) :: Fin(nBasin,nBasin) double precision,intent(in) :: eHFin(nBasin)
double precision,intent(in) :: ERIin(nBasin,nBasin,nBasin,nBasin) double precision,intent(in) :: ERIin(nBasin,nBasin,nBasin,nBasin)
double precision,intent(in) :: E0 double precision,intent(in) :: E0
@ -26,7 +26,7 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi
integer :: nV integer :: nV
integer :: nR integer :: nR
double precision,allocatable :: F(:,:) double precision,allocatable :: eHF(:)
double precision,allocatable :: sERI(:,:,:,:) double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: ERI(:,:,:,:) double precision,allocatable :: ERI(:,:,:,:)
@ -61,9 +61,9 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi
nV = 2*nVin nV = 2*nVin
nR = 2*nRin nR = 2*nRin
allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas)) allocate(eHF(nBas),sERI(nBas,nBas,nBas,nBas))
call spatial_to_spin_fock(nBasin,Fin,nBas,F) call spatial_to_spin_MO_energy(nBasin,eHFin,nBas,eHF)
call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI) call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI)
! Antysymmetrize ERIs ! Antysymmetrize ERIs
@ -144,22 +144,22 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi
kcld = kcld + 1 kcld = kcld + 1
tmp = & tmp = &
E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ F(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + eHF(j)*Kronecker_delta(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &
- F(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) & - eHF(j)*Kronecker_delta(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
- F(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) & - eHF(j)*Kronecker_delta(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+ F(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) & + eHF(j)*Kronecker_delta(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
- F(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) & - eHF(i)*Kronecker_delta(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+ F(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) & + eHF(i)*Kronecker_delta(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+ F(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) & + eHF(i)*Kronecker_delta(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
- F(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) & - eHF(i)*Kronecker_delta(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
+ F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + eHF(a)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
- F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & - eHF(a)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
- F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - eHF(a)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+ F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + eHF(a)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & - eHF(b)*Kronecker_delta(b,d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + eHF(b)*Kronecker_delta(b,c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + eHF(b)*Kronecker_delta(b,d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - eHF(b)*Kronecker_delta(b,c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) & - ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) &
+ ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &

56
src/CI/CISD.f90
View File

@ -1,4 +1,4 @@
subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,Fin,E0) subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIin,eHFin,E0)
! Perform configuration interaction with singles and doubles ! Perform configuration interaction with singles and doubles
@ -14,7 +14,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
integer,intent(in) :: nOin integer,intent(in) :: nOin
integer,intent(in) :: nVin integer,intent(in) :: nVin
integer,intent(in) :: nRin integer,intent(in) :: nRin
double precision,intent(in) :: Fin(nBasin,nBasin) double precision,intent(in) :: eHFin(nBasin)
double precision,intent(in) :: ERIin(nBasin,nBasin,nBasin,nBasin) double precision,intent(in) :: ERIin(nBasin,nBasin,nBasin,nBasin)
double precision,intent(in) :: E0 double precision,intent(in) :: E0
@ -26,7 +26,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
integer :: nV integer :: nV
integer :: nR integer :: nR
double precision,allocatable :: F(:,:) double precision,allocatable :: eHF(:)
double precision,allocatable :: sERI(:,:,:,:) double precision,allocatable :: sERI(:,:,:,:)
double precision,allocatable :: ERI(:,:,:,:) double precision,allocatable :: ERI(:,:,:,:)
@ -62,9 +62,9 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
nV = 2*nVin nV = 2*nVin
nR = 2*nRin nR = 2*nRin
allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas)) allocate(eHF(nBas),sERI(nBas,nBas,nBas,nBas))
call spatial_to_spin_fock(nBasin,Fin,nBas,F) call spatial_to_spin_MO_energy(nBasin,eHFin,nBas,eHF)
call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI) call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI)
! Antysymmetrize ERIs ! Antysymmetrize ERIs
@ -111,7 +111,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
do a=1,nV-nR do a=1,nV-nR
ia = ia + 1 ia = ia + 1
tmp = F(i,nO+a) tmp = 0d0
H(ishift+1,jshift+ia) = tmp H(ishift+1,jshift+ia) = tmp
H(jshift+ia,ishift+1) = tmp H(jshift+ia,ishift+1) = tmp
@ -159,9 +159,9 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
do c=1,nV-nR do c=1,nV-nR
kc = kc + 1 kc = kc + 1
tmp = E0*Kronecker_delta(i,k)*Kronecker_delta(a,c) & tmp = E0*Kronecker_delta(i,k)*Kronecker_delta(a,c) &
- F(i,k)*Kronecker_delta(a,c) & - eHF(i)*Kronecker_Delta(i,k)*Kronecker_delta(a,c) &
+ F(nO+a,nO+c)*Kronecker_delta(i,k) & + eHF(a)*Kronecker_delta(a,c)*Kronecker_delta(i,k) &
- ERI(nO+a,k,nO+c,i) - ERI(nO+a,k,nO+c,i)
H(ishift+ia,jshift+kc) = tmp H(ishift+ia,jshift+kc) = tmp
@ -192,11 +192,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
do d=c+1,nV-nR do d=c+1,nV-nR
kcld = kcld + 1 kcld = kcld + 1
tmp = - F(l,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k) & tmp = - ERI(k,l,nO+d,i)*Kronecker_delta(a,c) &
+ F(l,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k) &
- F(k,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l) &
+ F(k,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l) &
- ERI(k,l,nO+d,i)*Kronecker_delta(a,c) &
+ ERI(k,l,nO+c,i)*Kronecker_delta(a,d) & + ERI(k,l,nO+c,i)*Kronecker_delta(a,d) &
- ERI(nO+a,l,nO+c,nO+d)*Kronecker_delta(i,k) & - ERI(nO+a,l,nO+c,nO+d)*Kronecker_delta(i,k) &
+ ERI(nO+a,k,nO+c,nO+d)*Kronecker_delta(i,l) + ERI(nO+a,k,nO+c,nO+d)*Kronecker_delta(i,l)
@ -236,22 +232,22 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI
kcld = kcld + 1 kcld = kcld + 1
tmp = & tmp = &
E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & E0*Kronecker_delta(i,k)*Kronecker_delta(j,l)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ F(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + eHF(j)*Kronecker_delta(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &
- F(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) & - eHF(j)*Kronecker_delta(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) &
- F(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) & - eHF(j)*Kronecker_delta(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) &
+ F(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) & + eHF(j)*Kronecker_delta(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) &
- F(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) & - eHF(i)*Kronecker_delta(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) &
+ F(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) & + eHF(i)*Kronecker_delta(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) &
+ F(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) & + eHF(i)*Kronecker_delta(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) &
- F(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) & - eHF(i)*Kronecker_delta(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) &
+ F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + eHF(a)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
- F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & - eHF(a)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
- F(nO+a,nO+d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - eHF(a)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+ F(nO+a,nO+c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + eHF(a)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & - eHF(b)*Kronecker_delta(b,d)*Kronecker_delta(a,c)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) & + eHF(b)*Kronecker_delta(b,c)*Kronecker_delta(a,d)*Kronecker_delta(i,l)*Kronecker_delta(j,k) &
+ F(nO+b,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + eHF(b)*Kronecker_delta(b,d)*Kronecker_delta(a,c)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - eHF(b)*Kronecker_delta(b,c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
- ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) & - ERI(k,l,i,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c) &
+ ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + ERI(k,l,i,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
+ ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) &

101
src/CI/RCI.f90 Normal file
View File

@ -0,0 +1,101 @@
subroutine RCI(doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, &
epsHF,EHF,cHF,S)
! Configuration interaction module
implicit none
include 'parameters.h'
! Input variables
logical :: doCIS
logical :: doCIS_D
logical :: doCID
logical :: doCISD
logical :: doFCI
logical,intent(in) :: singlet
logical,intent(in) :: triplet
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
double precision :: start_CI ,end_CI ,t_CI
!------------------------------------------------------------------------
! Compute CIS excitations
!------------------------------------------------------------------------
if(doCIS) then
call wall_time(start_CI)
call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CI,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute CID excitations
!------------------------------------------------------------------------
if(doCID) then
call wall_time(start_CI)
call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CI,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute CISD excitations
!------------------------------------------------------------------------
if(doCISD) then
call wall_time(start_CI)
call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI,epsHF,EHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CI,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute FCI
!------------------------------------------------------------------------
if(doFCI) then
call wall_time(start_CI)
write(*,*) ' FCI is not yet implemented! Sorry.'
! call FCI(nBas,nC,nO,nV,nR,ERI,epsHF)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_CI,' seconds'
write(*,*)
end if
end subroutine

21
src/CI/CI.f90 → src/CI/UCI.f90
View File

@ -1,6 +1,5 @@
subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin_conserved,spin_flip, & subroutine UCI(doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS, &
nBas,nC,nO,nV,nR,nS,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int,dipole_int_aa,dipole_int_bb, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,EHF,cHF,S,F)
epsHF,EHF,cHF,S,F)
! Configuration interaction module ! Configuration interaction module
@ -14,10 +13,7 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
logical :: doCID logical :: doCID
logical :: doCISD logical :: doCISD
logical :: doFCI logical :: doFCI
logical :: unrestricted
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: spin_conserved logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip logical,intent(in) :: spin_flip
integer,intent(in) :: nBas integer,intent(in) :: nBas
@ -31,11 +27,9 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: F(nBas,nBas) double precision,intent(in) :: F(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart,nspin) double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart,nspin)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart,nspin) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart,nspin)
@ -50,12 +44,8 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
if(doCIS) then if(doCIS) then
call wall_time(start_CI) call wall_time(start_CI)
if(unrestricted) then call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb, &
call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb, & ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
else
call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF)
end if
call wall_time(end_CI) call wall_time(end_CI)
t_CI = end_CI - start_CI t_CI = end_CI - start_CI
@ -71,7 +61,6 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
if(doCID) then if(doCID) then
call wall_time(start_CI) call wall_time(start_CI)
call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI,F,EHF)
call wall_time(end_CI) call wall_time(end_CI)
t_CI = end_CI - start_CI t_CI = end_CI - start_CI
@ -87,7 +76,6 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
if(doCISD) then if(doCISD) then
call wall_time(start_CI) call wall_time(start_CI)
call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI,F,EHF)
call wall_time(end_CI) call wall_time(end_CI)
t_CI = end_CI - start_CI t_CI = end_CI - start_CI
@ -104,7 +92,6 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin
call wall_time(start_CI) call wall_time(start_CI)
write(*,*) ' FCI is not yet implemented! Sorry.' write(*,*) ' FCI is not yet implemented! Sorry.'
! call FCI(nBas,nC,nO,nV,nR,ERI,epsHF)
call wall_time(end_CI) call wall_time(end_CI)
t_CI = end_CI - start_CI t_CI = end_CI - start_CI

148
src/GF/RGF.f90 Normal file
View File

@ -0,0 +1,148 @@
subroutine RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis, &
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
! Green's function module
implicit none
include 'parameters.h'
! Input variables
logical :: doG0F2
logical :: doevGF2
logical :: doqsGF2
logical :: doG0F3
logical :: doevGF3
integer :: renorm
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
logical,intent(in) :: dophBSE
logical,intent(in) :: doppBSE
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
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,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
double precision :: start_GF ,end_GF ,t_GF
!------------------------------------------------------------------------
! Compute G0F2 electronic binding energies
!------------------------------------------------------------------------
if(doG0F2) then
call wall_time(start_GF)
call G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute evGF2 electronic binding energies
!------------------------------------------------------------------------
if(doevGF2) then
call wall_time(start_GF)
call evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI,dipole_int,epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform qsGF2 calculation
!------------------------------------------------------------------------
if(doqsGF2) then
call wall_time(start_GF)
call qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute G0F3 electronic binding energies
!------------------------------------------------------------------------
if(doG0F3) then
call wall_time(start_GF)
call G0F3(renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute evGF3 electronic binding energies
!------------------------------------------------------------------------
if(doevGF3) then
call wall_time(start_GF)
call evGF3(maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds'
write(*,*)
end if
end subroutine

59
src/GF/GF.f90 → src/GF/UGF.f90
View File

@ -1,7 +1,7 @@
subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,thresh,max_diis, & subroutine UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis, &
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,spin_conserved,spin_flip,linearize,eta,regularize, & dophBSE,doppBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regularize, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_AO,dipole_int,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
! Green's function module ! Green's function module
@ -15,7 +15,6 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
logical :: doqsGF2 logical :: doqsGF2
logical :: doG0F3 logical :: doG0F3
logical :: doevGF3 logical :: doevGF3
logical :: unrestricted
integer :: renorm integer :: renorm
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
@ -26,8 +25,6 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
logical,intent(in) :: TDA logical,intent(in) :: TDA
logical,intent(in) :: dBSE logical,intent(in) :: dBSE
logical,intent(in) :: dTDA logical,intent(in) :: dTDA
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: spin_conserved logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip logical,intent(in) :: spin_flip
logical,intent(in) :: linearize logical,intent(in) :: linearize
@ -56,12 +53,10 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
@ -76,14 +71,9 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doG0F2) then if(doG0F2) then
call wall_time(start_GF) call wall_time(start_GF)
if(unrestricted) then call UG0F2(dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regularize, &
call UG0F2(dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regularize, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,epsHF)
dipole_int_aa,dipole_int_bb,epsHF)
else
call G0F2(dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
@ -99,15 +89,9 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doevGF2) then if(doevGF2) then
call wall_time(start_GF) call wall_time(start_GF)
if(unrestricted) then call evUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
call evUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip, & eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,cHF,epsHF)
dipole_int_aa,dipole_int_bb,cHF,epsHF)
else
call evGF2(dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
ERI,dipole_int,epsHF)
end if
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
@ -123,14 +107,9 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doqsGF2) then if(doqsGF2) then
call wall_time(start_GF) call wall_time(start_GF)
if(unrestricted) then call qsUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,regularize, &
call qsUGF2(maxSCF,thresh,max_diis,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,regularize, & nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
else
call qsGF2(maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc, &
nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
end if
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
@ -146,11 +125,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doG0F3) then if(doG0F3) then
call wall_time(start_GF) call wall_time(start_GF)
if(unrestricted) then print*,'Unrestricted version of G0F3 not yet implemented! Sorry.'
print*,'Unrestricted version of G0F3 not yet implemented! Sorry.'
else
call G0F3(renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
end if
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF
@ -166,11 +141,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t
if(doevGF3) then if(doevGF3) then
call wall_time(start_GF) call wall_time(start_GF)
if(unrestricted) then print*,'Unrestricted version of evGF3 not yet implemented! Sorry.'
print*,'Unrestricted version of evGF3 not yet implemented! Sorry.'
else
call evGF3(maxSCF,thresh,max_diis,renorm,nBas,nC,nO,nV,nR,ERI,epsHF)
end if
call wall_time(end_GF) call wall_time(end_GF)
t_GF = end_GF - start_GF t_GF = end_GF - start_GF

174
src/GT/RGT.f90 Normal file
View File

@ -0,0 +1,174 @@
subroutine RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
! T-matrix module
implicit none
include 'parameters.h'
! Input variables
logical :: doG0T0pp
logical :: doevGTpp
logical :: doqsGTpp
logical :: doG0T0eh
logical :: doevGTeh
logical :: doqsGTeh
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: dophBSE
logical,intent(in) :: dophBSE2
logical,intent(in) :: TDA_T
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: doppBSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
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,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
double precision :: start_GT ,end_GT ,t_GT
!------------------------------------------------------------------------
! Perform G0T0pp calculatiom
!------------------------------------------------------------------------
if(doG0T0pp) then
call wall_time(start_GT)
call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform evGTpp calculatiom
!------------------------------------------------------------------------
if(doevGTpp) then
call wall_time(start_GT)
call evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_GT,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform qsGTpp calculation
!------------------------------------------------------------------------
if(doqsGTpp) then
call wall_time(start_GT)
call qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, &
PHF,cHF,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_GT,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform G0T0eh calculatiom
!------------------------------------------------------------------------
if(doG0T0eh) then
call wall_time(start_GT)
call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform evGTeh calculation
!------------------------------------------------------------------------
if(doevGTeh) then
call wall_time(start_GT)
call evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_GT,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform qsGTeh calculation
!------------------------------------------------------------------------
if(doqsGTeh) then
call wall_time(start_GT)
call qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, &
PHF,cHF,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GT,' seconds'
write(*,*)
end if
end subroutine

72
src/GT/GT.f90 → src/GT/UGT.f90
View File

@ -1,8 +1,7 @@
subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted,maxSCF,thresh,max_diis,doACFDT, & subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,spin_conserved,spin_flip, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int,dipole_int_aa,dipole_int_bb, & ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
PHF,cHF,epsHF)
! T-matrix module ! T-matrix module
@ -17,7 +16,6 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
logical :: doG0T0eh logical :: doG0T0eh
logical :: doevGTeh logical :: doevGTeh
logical :: doqsGTeh logical :: doqsGTeh
logical :: unrestricted
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
@ -32,8 +30,6 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
logical,intent(in) :: dBSE logical,intent(in) :: dBSE
logical,intent(in) :: dTDA logical,intent(in) :: dTDA
logical,intent(in) :: doppBSE logical,intent(in) :: doppBSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: spin_conserved logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip logical,intent(in) :: spin_flip
logical,intent(in) :: linearize logical,intent(in) :: linearize
@ -62,12 +58,10 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(inout):: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(inout):: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(inout):: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
@ -83,14 +77,9 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
if(doG0T0pp) then if(doG0T0pp) then
call wall_time(start_GT) call wall_time(start_GT)
if(unrestricted) then call UG0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
call UG0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,cHF,epsHF)
dipole_int_aa,dipole_int_bb,cHF,epsHF)
else
call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -106,14 +95,9 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
if(doevGTpp) then if(doevGTpp) then
call wall_time(start_GT) call wall_time(start_GT)
if(unrestricted) then call evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
call evUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb, & cHF,epsHF)
cHF,epsHF)
else
call evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -129,15 +113,9 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
if(doqsGTpp) then if(doqsGTpp) then
call wall_time(start_GT) call wall_time(start_GT)
if(unrestricted) then call qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
call qsUGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, & eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
else
call qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, &
PHF,cHF,epsHF)
end if
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -153,14 +131,7 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
if(doG0T0eh) then if(doG0T0eh) then
call wall_time(start_GT) call wall_time(start_GT)
if(unrestricted) then print*,'Unrestricted version of G0T0eh not yet implemented! Sorry.'
print*,'Unrestricted version of G0T0eh not yet implemented! Sorry.'
else
call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
! call soG0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
! linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -176,12 +147,7 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
if(doevGTeh) then if(doevGTeh) then
call wall_time(start_GT) call wall_time(start_GT)
if(unrestricted) then print*,'Unrestricted version of evGTeh not yet implemented! Sorry.'
print*,'Unrestricted version of evGTeh not yet implemented! Sorry.'
else
call evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT
@ -197,13 +163,7 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted
if(doqsGTeh) then if(doqsGTeh) then
call wall_time(start_GT) call wall_time(start_GT)
if(unrestricted) then print*,'Unrestricted version of qsGTeh not yet implemented! Sorry.'
print*,'Unrestricted version of qsGTeh not yet implemented! Sorry.'
else
call qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI,dipole_int_AO,dipole_int, &
PHF,cHF,epsHF)
end if
call wall_time(end_GT) call wall_time(end_GT)
t_GT = end_GT - start_GT t_GT = end_GT - start_GT

172
src/GW/RGW.f90 Normal file
View File

@ -0,0 +1,172 @@
subroutine RGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI,dipole_int_AO,dipole_int,PHF,cHF,epsHF)
! GW module
implicit none
include 'parameters.h'
! Input variables
logical :: doG0W0
logical :: doevGW
logical :: doqsGW
logical :: doufG0W0
logical :: doufGW
logical :: doSRGqsGW
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: doXBS
logical,intent(in) :: dophBSE
logical,intent(in) :: dophBSE2
logical,intent(in) :: TDA_W
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: doppBSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: linearize
double precision,intent(in) :: eta
logical,intent(in) :: regularize
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
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,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
double precision :: start_GW ,end_GW ,t_GW
!------------------------------------------------------------------------
! Perform G0W0 calculatiom
!------------------------------------------------------------------------
if(doG0W0) then
call wall_time(start_GW)
call G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform evGW calculation
!------------------------------------------------------------------------
if(doevGW) then
call wall_time(start_GW)
call evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform qsGW calculation
!------------------------------------------------------------------------
if(doqsGW) then
call wall_time(start_GW)
call qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform SRG-qsGW calculation
!------------------------------------------------------------------------
if(doSRGqsGW) then
call wall_time(start_GW)
call SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, &
singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGW = ',t_GW,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform ufG0W0 calculatiom
!------------------------------------------------------------------------
if(doufG0W0) then
call wall_time(start_GW)
call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF,TDA_W)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Perform ufGW calculatiom
!------------------------------------------------------------------------
if(doufGW) then
call wall_time(start_GW)
call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds'
write(*,*)
end if
end subroutine

59
src/GW/GW.f90 → src/GW/UGW.f90
View File

@ -1,8 +1,7 @@
subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF,thresh,max_diis,doACFDT, & subroutine UGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet,spin_conserved,spin_flip, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & linearize,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int,dipole_int_aa,dipole_int_bb, & ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
PHF,cHF,epsHF)
! GW module ! GW module
@ -17,7 +16,6 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
logical :: doufG0W0 logical :: doufG0W0
logical :: doufGW logical :: doufGW
logical :: doSRGqsGW logical :: doSRGqsGW
logical :: unrestricted
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
@ -32,8 +30,6 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
logical,intent(in) :: dBSE logical,intent(in) :: dBSE
logical,intent(in) :: dTDA logical,intent(in) :: dTDA
logical,intent(in) :: doppBSE logical,intent(in) :: doppBSE
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: spin_conserved logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip logical,intent(in) :: spin_flip
logical,intent(in) :: linearize logical,intent(in) :: linearize
@ -62,12 +58,10 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
double precision,intent(in) :: Hc(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas) double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(inout):: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(inout):: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(inout):: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(inout):: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
@ -82,14 +76,9 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
if(doG0W0) then if(doG0W0) then
call wall_time(start_GW) call wall_time(start_GW)
if(unrestricted) then call UG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
call UG0W0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,cHF,epsHF)
dipole_int_aa,dipole_int_bb,cHF,epsHF)
else
call G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, &
linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_GW) call wall_time(end_GW)
t_GW = end_GW - start_GW t_GW = end_GW - start_GW
@ -105,14 +94,9 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
if(doevGW) then if(doevGW) then
call wall_time(start_GW) call wall_time(start_GW)
if(unrestricted) then call evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA, &
call evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA, & spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S, &
spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,epsHF)
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,epsHF)
else
call evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_GW) call wall_time(end_GW)
t_GW = end_GW - start_GW t_GW = end_GW - start_GW
@ -128,16 +112,9 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
if(doqsGW) then if(doqsGW) then
call wall_time(start_GW) call wall_time(start_GW)
if(unrestricted) then call qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
call qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
else
call qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, &
singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
end if
call wall_time(end_GW) call wall_time(end_GW)
t_GW = end_GW - start_GW t_GW = end_GW - start_GW
@ -153,13 +130,7 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
if(doSRGqsGW) then if(doSRGqsGW) then
call wall_time(start_GW) call wall_time(start_GW)
if(unrestricted) then print*,'Unrestricted version of SRG-qsGW not yet implemented! Sorry.'
print*,'Unrestricted version of SRG-qsGW not yet implemented! Sorry.'
else
call SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, &
singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI, &
dipole_int_AO,dipole_int,PHF,cHF,epsHF)
end if
call wall_time(end_GW) call wall_time(end_GW)
t_GW = end_GW - start_GW t_GW = end_GW - start_GW
@ -175,7 +146,6 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
if(doufG0W0) then if(doufG0W0) then
call wall_time(start_GW) call wall_time(start_GW)
call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF,TDA_W)
call wall_time(end_GW) call wall_time(end_GW)
t_GW = end_GW - start_GW t_GW = end_GW - start_GW
@ -191,7 +161,6 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF
if(doufGW) then if(doufGW) then
call wall_time(start_GW) call wall_time(start_GW)
call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,epsHF)
call wall_time(end_GW) call wall_time(end_GW)
t_GW = end_GW - start_GW t_GW = end_GW - start_GW

150
src/HF/HF.f90
View File

@ -1,150 +0,0 @@
subroutine HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_diis,guess_type,mix,level_shift, &
nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
! Hartree-Fock module
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: doRHF
logical,intent(in) :: doUHF
logical,intent(in) :: doROHF
logical,intent(in) :: doRMOM
logical,intent(in) :: doUMOM
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
logical,intent(in) :: mix
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
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) :: 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,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
double precision :: start_HF ,end_HF ,t_HF
! Output variables
logical,intent(out) :: unrestricted
double precision,intent(out) :: EHF
double precision,intent(out) :: epsHF(nBas,nspin)
double precision,intent(out) :: cHF(nBas,nBas,nspin)
double precision,intent(out) :: PHF(nBas,nBas,nspin)
double precision,intent(out) :: F(nBas,nBas,nspin)
!------------------------------------------------------------------------
! Compute RHF energy
!------------------------------------------------------------------------
if(doRHF) then
! Check that RHF calculation is worth doing...
if(nO(1) /= nO(2)) then
write(*,*) ' !!! The system does not appear to be closed shell !!!'
write(*,*)
stop
end if
call wall_time(start_HF)
call RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute UHF energy
!------------------------------------------------------------------------
if(doUHF) then
! Switch on the unrestricted flag
unrestricted = .true.
call wall_time(start_HF)
call UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UHF = ',t_HF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute ROHF energy
!------------------------------------------------------------------------
if(doROHF) then
! Switch on the unrestricted flag
unrestricted = .true.
call wall_time(start_HF)
call ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ROHF = ',t_HF,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Restricted maximum overlap method
!------------------------------------------------------------------------
if(doRMOM) then
! Check that RMOM calculation is worth doing...
if(nO(1) /= nO(2)) then
write(*,*) ' !!! The system does not appear to be closed shell !!!'
write(*,*)
stop
end if
! call RMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, &
! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
end if
!------------------------------------------------------------------------
! Unrestricted maximum overlap method
!------------------------------------------------------------------------
if(doUMOM) then
! Switch on the unrestricted flag
unrestricted = .true.
! call UMOM(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc, &
! nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,epsHF,cHF,PHF)
end if
end subroutine

8
src/HF/RHF.f90
View File

@ -1,5 +1,5 @@
subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,EHF,eps,c,P) nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,eps,c,P)
! Perform restricted Hartree-Fock calculation ! Perform restricted Hartree-Fock calculation
@ -49,6 +49,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc
double precision,allocatable :: J(:,:) double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:) double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:) double precision,allocatable :: cp(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: Fp(:,:) double precision,allocatable :: Fp(:,:)
! Output variables ! Output variables
@ -57,7 +58,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc
double precision,intent(out) :: eps(nBas) double precision,intent(out) :: eps(nBas)
double precision,intent(out) :: c(nBas,nBas) double precision,intent(out) :: c(nBas,nBas)
double precision,intent(out) :: P(nBas,nBas) double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: F(nBas,nBas)
! Hello world ! Hello world
@ -73,8 +73,8 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc
! Memory allocation ! Memory allocation
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),cp(nBas,nBas),Fp(nBas,nBas), & allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), &
error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) Fp(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and density matrix ! Guess coefficients and density matrix

10
src/HF/ROHF.f90
View File

@ -1,5 +1,5 @@
subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,P) nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,Ptot)
! Perform restricted open-shell Hartree-Fock calculation ! Perform restricted open-shell Hartree-Fock calculation
@ -11,7 +11,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
integer,intent(in) :: guess_type integer,intent(in) :: guess_type
logical,intent(in) :: mix double precision,intent(in) :: mix
double precision,intent(in) :: level_shift double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: nBas integer,intent(in) :: nBas
@ -48,7 +48,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
double precision,allocatable :: F(:,:,:) double precision,allocatable :: F(:,:,:)
double precision,allocatable :: Fp(:,:) double precision,allocatable :: Fp(:,:)
double precision,allocatable :: Ftot(:,:) double precision,allocatable :: Ftot(:,:)
double precision,allocatable :: Ptot(:,:) double precision,allocatable :: P(:,:,:)
double precision,allocatable :: K(:,:,:) double precision,allocatable :: K(:,:,:)
double precision,allocatable :: err(:,:) double precision,allocatable :: err(:,:)
double precision,allocatable :: err_diis(:,:) double precision,allocatable :: err_diis(:,:)
@ -62,7 +62,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
double precision,intent(out) :: EHF double precision,intent(out) :: EHF
double precision,intent(out) :: e(nBas) double precision,intent(out) :: e(nBas)
double precision,intent(out) :: c(nBas,nBas) double precision,intent(out) :: c(nBas,nBas)
double precision,intent(out) :: P(nBas,nBas,nspin) double precision,intent(out) :: Ptot(nBas,nBas)
! Hello world ! Hello world
@ -79,7 +79,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
! Memory allocation ! Memory allocation
allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), & allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), &
Ptot(nBas,nBas),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), & P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), &
err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and demsity matrices ! Guess coefficients and demsity matrices

4
src/HF/UHF.f90
View File

@ -11,7 +11,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
integer,intent(in) :: maxSCF integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis integer,intent(in) :: max_diis
integer,intent(in) :: guess_type integer,intent(in) :: guess_type
logical,intent(in) :: mix double precision,intent(in) :: mix
double precision,intent(in) :: level_shift double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh double precision,intent(in) :: thresh
integer,intent(in) :: nBas integer,intent(in) :: nBas
@ -181,7 +181,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
! Mix guess for UHF solution in singlet states ! Mix guess for UHF solution in singlet states
if(mix .and. nSCF == 1) call mix_guess(nBas,nO,c) if(nSCF == 1) call mix_guess(nBas,nO,mix,c)
! Compute density matrix ! Compute density matrix

9
src/HF/mix_guess.f90
View File

@ -1,4 +1,4 @@
subroutine mix_guess(nBas,nO,c) subroutine mix_guess(nBas,nO,mix,c)
! Guess mixing for UHF calculations on singlet states ! Guess mixing for UHF calculations on singlet states
@ -9,15 +9,16 @@ subroutine mix_guess(nBas,nO,c)
integer,intent(in) :: nBas integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin) integer,intent(in) :: nO(nspin)
double precision,intent(in) :: mix
! Local variables ! Local variables
double precision,parameter :: th = 0.1d0
double precision,allocatable :: HOMO(:), LUMO(:) double precision,allocatable :: HOMO(:), LUMO(:)
double precision,allocatable :: HOMOa(:),LUMOa(:) double precision,allocatable :: HOMOa(:),LUMOa(:)
double precision,allocatable :: HOMOb(:),LUMOb(:) double precision,allocatable :: HOMOb(:),LUMOb(:)
double precision :: th
! Output variables ! Output variables
double precision,intent(inout):: c(nBas,nBas,nspin) double precision,intent(inout):: c(nBas,nBas,nspin)
@ -30,6 +31,8 @@ subroutine mix_guess(nBas,nO,c)
! Perform HOMO and LUMO rotation for guess mixing ! Perform HOMO and LUMO rotation for guess mixing
th = 2d0*pi*mix
HOMO(:) = c(:,nO(1),1) HOMO(:) = c(:,nO(1),1)
LUMO(:) = c(:,nO(1)+1,1) LUMO(:) = c(:,nO(1)+1,1)

62
src/MP/RMP.f90 Normal file
View File

@ -0,0 +1,62 @@
subroutine RMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
! Moller-Plesset module
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: doMP2
logical,intent(in) :: doMP3
logical,intent(in) :: regularize
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision :: start_MP ,end_MP ,t_MP
! Output variables
!------------------------------------------------------------------------
! Compute MP3 energy
!------------------------------------------------------------------------
if(doMP2) then
call wall_time(start_MP)
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute MP3 energy
!------------------------------------------------------------------------
if(doMP3) then
call wall_time(start_MP)
write(*,*) 'MP3 NYI for UHF reference'
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds'
write(*,*)
end if
end subroutine

20
src/MP/MP.f90 → src/MP/UMP.f90
View File

@ -1,4 +1,4 @@
subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF) subroutine UMP(doMP2,doMP3,regularize,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
! Moller-Plesset module ! Moller-Plesset module
@ -9,7 +9,6 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,
logical,intent(in) :: doMP2 logical,intent(in) :: doMP2
logical,intent(in) :: doMP3 logical,intent(in) :: doMP3
logical,intent(in) :: unrestricted
logical,intent(in) :: regularize logical,intent(in) :: regularize
integer,intent(in) :: nBas integer,intent(in) :: nBas
@ -19,8 +18,7 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,
integer,intent(in) :: nR(nspin) integer,intent(in) :: nR(nspin)
double precision,intent(in) :: ENuc double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas) double precision,intent(in) :: epsHF(nBas,nspin)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
@ -38,11 +36,7 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,
if(doMP2) then if(doMP2) then
call wall_time(start_MP) call wall_time(start_MP)
if(unrestricted) then call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
else
call MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
end if
call wall_time(end_MP) call wall_time(end_MP)
t_MP = end_MP - start_MP t_MP = end_MP - start_MP
@ -58,13 +52,7 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa,
if(doMP3) then if(doMP3) then
call wall_time(start_MP) call wall_time(start_MP)
write(*,*) 'MP3 NYI for UHF reference'
if(unrestricted) then
write(*,*) 'MP3 NYI for UHF reference'
stop
else
call MP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF)
end if
call wall_time(end_MP) call wall_time(end_MP)
t_MP = end_MP - start_MP t_MP = end_MP - start_MP

412
src/QuAcK/QuAcK.f90
View File

@ -3,35 +3,27 @@ program QuAcK
implicit none implicit none
include 'parameters.h' include 'parameters.h'
logical :: unrestricted = .false. logical :: doRQuAcK,doUQuAcK,doGQuAcK
logical :: doHF,doRHF,doUHF,doROHF,doRMOM,doUMOM logical :: doRHF,doUHF,doGHF,doROHF
logical :: dostab logical :: dostab
logical :: doKS logical :: doMP2,doMP3
logical :: doMP,doMP2,doMP3 logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical :: doCC,doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical :: dodrCCD,dorCCD,docrCCD,dolCCD logical :: dodrCCD,dorCCD,docrCCD,dolCCD
logical :: doCI,doCIS,doCIS_D,doCID,doCISD,doFCI logical :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical :: doRPA,dophRPA,dophRPAx,docrRPA,doppRPA logical :: dophRPA,dophRPAx,docrRPA,doppRPA
logical :: doGF,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical :: doGW,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW
logical :: doGT,doG0T0pp,doevGTpp,doqsGTpp logical :: doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh
logical :: doG0T0eh,doevGTeh,doqsGTeh
integer :: nNuc,nBas integer :: nNuc,nBas
integer :: nEl(nspin)
integer :: nC(nspin) integer :: nC(nspin)
integer :: nO(nspin) integer :: nO(nspin)
integer :: nV(nspin) integer :: nV(nspin)
integer :: nR(nspin) integer :: nR(nspin)
integer :: nS(nspin) integer :: nS(nspin)
double precision :: ENuc,EHF double precision :: ENuc
double precision,allocatable :: ZNuc(:),rNuc(:,:) double precision,allocatable :: ZNuc(:),rNuc(:,:)
double precision,allocatable :: cHF(:,:,:),epsHF(:,:),PHF(:,:,:)
logical :: doACFDT
logical :: exchange_kernel
logical :: doXBS
double precision,allocatable :: S(:,:) double precision,allocatable :: S(:,:)
double precision,allocatable :: T(:,:) double precision,allocatable :: T(:,:)
@ -39,65 +31,41 @@ program QuAcK
double precision,allocatable :: Hc(:,:) double precision,allocatable :: Hc(:,:)
double precision,allocatable :: X(:,:) double precision,allocatable :: X(:,:)
double precision,allocatable :: dipole_int_AO(:,:,:) double precision,allocatable :: dipole_int_AO(:,:,:)
double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: dipole_int_aa(:,:,:)
double precision,allocatable :: dipole_int_bb(:,:,:)
double precision,allocatable :: F_AO(:,:,:)
double precision,allocatable :: F_MO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:) double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: ixyz
double precision,allocatable :: ERI_MO_aaaa(:,:,:,:)
double precision,allocatable :: ERI_MO_aabb(:,:,:,:)
double precision,allocatable :: ERI_MO_bbbb(:,:,:,:)
double precision :: start_QuAcK ,end_QuAcK ,t_QuAcK double precision :: start_QuAcK,end_QuAcK,t_QuAcK
double precision :: start_int ,end_int ,t_int double precision :: start_int ,end_int ,t_int
double precision :: start_HF ,end_HF ,t_HF
double precision :: start_stab ,end_stab ,t_stab
double precision :: start_KS ,end_KS ,t_KS
double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO
double precision :: start_MP ,end_MP ,t_MP
double precision :: start_CC ,end_CC ,t_CC
double precision :: start_CI ,end_CI ,t_CI
double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_GF ,end_GF ,t_GF
double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT
integer :: maxSCF_HF,max_diis_HF integer :: maxSCF_HF,max_diis_HF
double precision :: thresh_HF,level_shift double precision :: thresh_HF,level_shift,mix
logical :: DIIS_HF,guess_type,ortho_type,mix logical :: guess_type
logical :: reg_MP logical :: reg_MP
integer :: maxSCF_CC,max_diis_CC integer :: maxSCF_CC,max_diis_CC
double precision :: thresh_CC double precision :: thresh_CC
logical :: DIIS_CC
logical :: singlet
logical :: triplet
logical :: spin_conserved logical :: spin_conserved
logical :: spin_flip logical :: spin_flip
logical :: TDA logical :: TDA
integer :: maxSCF_GF,max_diis_GF,renorm_GF integer :: maxSCF_GF,max_diis_GF,renorm_GF
double precision :: thresh_GF double precision :: thresh_GF
logical :: DIIS_GF,lin_GF,reg_GF logical :: lin_GF,reg_GF
double precision :: eta_GF double precision :: eta_GF
integer :: maxSCF_GW,max_diis_GW integer :: maxSCF_GW,max_diis_GW
double precision :: thresh_GW double precision :: thresh_GW
logical :: DIIS_GW,TDA_W,lin_GW,reg_GW logical :: TDA_W,lin_GW,reg_GW
double precision :: eta_GW double precision :: eta_GW
integer :: maxSCF_GT,max_diis_GT integer :: maxSCF_GT,max_diis_GT
double precision :: thresh_GT double precision :: thresh_GT
logical :: DIIS_GT,TDA_T,lin_GT,reg_GT logical :: TDA_T,lin_GT,reg_GT
double precision :: eta_GT double precision :: eta_GT
logical :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA logical :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA
logical :: doACFDT,exchange_kernel,doXBS
!-------------! !-------------!
! Hello World ! ! Hello World !
@ -123,31 +91,31 @@ program QuAcK
! Method selection ! ! Method selection !
!------------------! !------------------!
call read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & call read_methods(doRHF,doUHF,doGHF,doROHF, &
doMP2,doMP3, & doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD, & dodrCCD,dorCCD,docrCCD,dolCCD, &
doCIS,doCIS_D,doCID,doCISD,doFCI, & doCIS,doCIS_D,doCID,doCISD,doFCI, &
dophRPA,dophRPAx,docrRPA,doppRPA, & dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2, & doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, & doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW,doSRGqsGW, & doG0W0,doevGW,doqsGW,doSRGqsGW, &
doufG0W0,doufGW, & doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp, & doG0T0pp,doevGTpp,doqsGTpp, &
doG0T0eh,doevGTeh,doqsGTeh) doG0T0eh,doevGTeh,doqsGTeh)
!--------------------------! !--------------------------!
! Read options for methods ! ! Read options for methods !
!--------------------------! !--------------------------!
call read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, & call read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab, &
reg_MP, & reg_MP, &
maxSCF_CC,thresh_CC,DIIS_CC,max_diis_CC, & maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, & TDA,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, & maxSCF_GF,thresh_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, &
maxSCF_GW,thresh_GW,DIIS_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, & maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,DIIS_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, & doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA) dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
!------------------------------------------------! !------------------------------------------------!
@ -163,7 +131,7 @@ program QuAcK
! = nO*nV ! ! = nO*nV !
!------------------------------------------------! !------------------------------------------------!
call read_molecule(nNuc,nEl,nO,nC,nR) call read_molecule(nNuc,nO,nC,nR)
allocate(ZNuc(nNuc),rNuc(nNuc,ncart)) allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
! Read geometry ! Read geometry
@ -183,9 +151,8 @@ program QuAcK
! Memory allocation for one- and two-electron integrals ! Memory allocation for one- and two-electron integrals
allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin),S(nBas,nBas),T(nBas,nBas), & allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), &
V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart), & ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart))
dipole_int_MO(nBas,nBas,ncart),F_AO(nBas,nBas,nspin))
! Read integrals ! Read integrals
@ -203,281 +170,56 @@ program QuAcK
! Compute orthogonalization matrix ! Compute orthogonalization matrix
call orthogonalization_matrix(ortho_type,nBas,S,X) call orthogonalization_matrix(nBas,S,X)
!---------------------! !---------------------!
! Hartree-Fock module ! ! Choose QuAcK branch !
!---------------------! !---------------------!
doHF = doRHF .or. doUHF .or. doROHF .or. doRMOM .or. doUMOM doRQuAcK = .false.
if(doRHF) doRQuAcK = .true.
if(doHF) then doUQuAcK = .false.
if(doUHF) doUQuAcK = .true.
call wall_time(start_HF) doGQuAcK = .false.
call HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,max_diis_HF, & if(doGHF) doGQuAcK = .true.
guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO,ERI_AO, &
dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for HF = ',t_HF,' seconds'
write(*,*)
end if
!------------------!
! Kohn-Sham module !
!------------------!
if(doKS) then
! Switch on the unrestricted flag
unrestricted = .true.
call wall_time(start_KS)
write(*,*)
write(*,*) 'KS module has been disabled for now! Sorry.'
write(*,*)
! call eDFT(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, &
! nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
! max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EHF,epsHF,cHF,PHF,Vxc)
call wall_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for KS = ',t_KS,' seconds'
write(*,*)
end if
!----------------------------------!
! AO to MO integral transformation !
!----------------------------------!
call wall_time(start_AOtoMO)
write(*,*)
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
if(unrestricted) then
! Read and transform dipole-related integrals
allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart))
dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:)
dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz))
end do
! Memory allocation
allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas))
! 4-index transform for (aa|aa) block
call AOtoMO_integral_transform(1,1,1,1,nBas,cHF,ERI_AO,ERI_MO_aaaa)
! 4-index transform for (aa|bb) block
call AOtoMO_integral_transform(1,1,2,2,nBas,cHF,ERI_AO,ERI_MO_aabb)
! 4-index transform for (bb|bb) block
call AOtoMO_integral_transform(2,2,2,2,nBas,cHF,ERI_AO,ERI_MO_bbbb)
else
! Memory allocation
allocate(ERI_MO(nBas,nBas,nBas,nBas),F_MO(nBas,nBas,nspin))
! Read and transform dipole-related integrals
dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz))
end do
! 4-index transform
call AOtoMO_integral_transform(1,1,1,1,nBas,cHF,ERI_AO,ERI_MO)
F_MO(:,:,1) = F_AO(:,:,1)
call AOtoMO_transform(nBas,cHF,F_MO)
end if
call wall_time(end_AOtoMO)
t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*)
!-----------------------------------!
! Stability analysis of HF solution !
!-----------------------------------!
if(dostab) then
call wall_time(start_stab)
if(unrestricted) then
call UHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb)
else
call RHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO)
end if
call wall_time(end_stab)
t_stab = end_stab - start_stab
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds'
write(*,*)
end if
!-----------------------!
! Moller-Plesset module !
!-----------------------!
doMP = doMP2 .or. doMP3
if(doMP) then
call wall_time(start_MP)
call MP(doMP2,doMP3,unrestricted,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ENuc,EHF,epsHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP = ',t_MP,' seconds'
write(*,*)
end if
!------------------------!
! Coupled-cluster module !
!------------------------!
doCC = doCCD .or. dopCCD .or. doDCD .or. doCCSD .or. doCCSDT .or. &
dodrCCD .or. dorCCD .or. docrCCD .or. dolCCD
if(doCC) then
call wall_time(start_CC)
call CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds'
write(*,*)
end if
!----------------------------------!
! Configuration interaction module !
!----------------------------------!
doCI = doCIS .or. doCID .or. doCISD .or. doFCI
if(doCI) then
call wall_time(start_CI)
call CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin_conserved,spin_flip, &
nBas,nC,nO,nV,nR,nS,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_MO,dipole_int_aa,dipole_int_bb, &
epsHF,EHF,cHF,S,F_MO)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CI = ',t_CI,' seconds'
write(*,*)
end if
!-----------------------------------!
! Random-phase approximation module !
!-----------------------------------!
doRPA = dophRPA .or. dophRPAx .or. docrRPA .or. doppRPA
if(doRPA) then
call wall_time(start_RPA)
call RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, &
TDA,doACFDT,exchange_kernel,singlet,triplet,spin_conserved,spin_flip, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, &
dipole_int_MO,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
write(*,*)
end if
!-------------------------! !-------------------------!
! Green's function module ! ! Restricted QuAcK branch !
!-------------------------! !-------------------------!
doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doG0F3 .or. doevGF3 if(doRQuAcK) &
call RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
if(doGF) then !---------------------------!
! Unrestricted QuAcK branch !
!---------------------------!
call wall_time(start_GF) if(doUQuAcK) &
call GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & call UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,spin_conserved,spin_flip,lin_GF,eta_GF,reg_GF, & dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
dipole_int_AO,dipole_int_MO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc, &
call wall_time(end_GF) S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
t_GF = end_GF - start_GF !--------------------------!
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' ! Generalized QuAcK branch !
write(*,*) !--------------------------!
end if ! if(doGQuAcK) call GQuAcK()
!-----------!
! GW module !
!-----------!
doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW .or. doSRGqsGW
if(doGW) then
call wall_time(start_GW)
call GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet,spin_conserved,spin_flip, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO,dipole_int_MO,dipole_int_aa,dipole_int_bb, &
PHF,cHF,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GW = ',t_GW,' seconds'
write(*,*)
end if
!-----------------!
! T-matrix module !
!-----------------!
doGT = doG0T0pp .or. doevGTpp .or. doqsGTpp .or. doG0T0eh .or. doevGTeh .or. doqsGTeh
if(doGT) then
call wall_time(start_GT)
call GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,spin_conserved,spin_flip, &
lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO,dipole_int_MO,dipole_int_aa,dipole_int_bb, &
PHF,cHF,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GT = ',t_GT,' seconds'
write(*,*)
end if
!--------------! !--------------!
! End of QuAcK ! ! End of QuAcK !

322
src/QuAcK/RQuAcK.f90 Normal file
View File

@ -0,0 +1,322 @@
subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
! Restricted branch of QuAcK
implicit none
include 'parameters.h'
logical,intent(in) :: doRHF,doROHF
logical,intent(in) :: dostab
logical,intent(in) :: doMP2,doMP3
logical,intent(in) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(in) :: dodrCCD,dorCCD,docrCCD,dolCCD
logical,intent(in) :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical,intent(in) :: dophRPA,dophRPAx,docrRPA,doppRPA
logical,intent(in) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical,intent(in) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW
logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh
integer,intent(in) :: nNuc,nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: ENuc
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,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix
logical,intent(in) :: guess_type
logical,intent(in) :: reg_MP
integer,intent(in) :: maxSCF_CC,max_diis_CC
double precision,intent(in) :: thresh_CC
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: TDA
integer,intent(in) :: maxSCF_GF,max_diis_GF,renorm_GF
double precision,intent(in) :: thresh_GF
logical,intent(in) :: lin_GF,reg_GF
double precision,intent(in) :: eta_GF
integer,intent(in) :: maxSCF_GW,max_diis_GW
double precision,intent(in) :: thresh_GW
logical,intent(in) :: TDA_W,lin_GW,reg_GW
double precision,intent(in) :: eta_GW
integer,intent(in) :: maxSCF_GT,max_diis_GT
double precision,intent(in) :: thresh_GT
logical,intent(in) :: TDA_T,lin_GT,reg_GT
double precision,intent(in) :: eta_GT
logical,intent(in) :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA
logical,intent(in) :: doACFDT,exchange_kernel,doXBS
! Local variables
logical :: doMP,doCC,doCI,doRPA,doGF,doGW,doGT
double precision :: start_HF ,end_HF ,t_HF
double precision :: start_stab ,end_stab ,t_stab
double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO
double precision :: start_MP ,end_MP ,t_MP
double precision :: start_CC ,end_CC ,t_CC
double precision :: start_CI ,end_CI ,t_CI
double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_GF ,end_GF ,t_GF
double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT
double precision,allocatable :: cHF(:,:),epsHF(:),PHF(:,:)
double precision :: EHF
double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: ixyz
write(*,*)
write(*,*) '******************************'
write(*,*) '* Restricted Branch of QuAcK *'
write(*,*) '******************************'
write(*,*)
!-------------------!
! Memory allocation !
!-------------------!
allocate(cHF(nBas,nBas),epsHF(nBas),PHF(nBas,nBas), &
dipole_int_MO(nBas,nBas,ncart),ERI_MO(nBas,nBas,nBas,nBas))
!---------------------!
! Hartree-Fock module !
!---------------------!
if(doRHF) then
call wall_time(start_HF)
call RHF(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds'
write(*,*)
end if
if(doROHF) then
call wall_time(start_HF)
call ROHF(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ROHF = ',t_HF,' seconds'
write(*,*)
end if
!----------------------------------!
! AO to MO integral transformation !
!----------------------------------!
call wall_time(start_AOtoMO)
write(*,*)
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
! Read and transform dipole-related integrals
dipole_int_MO(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF,dipole_int_MO(:,:,ixyz))
end do
! 4-index transform
call AOtoMO_integral_transform(1,1,1,1,nBas,cHF,ERI_AO,ERI_MO)
call wall_time(end_AOtoMO)
t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*)
!-----------------------------------!
! Stability analysis of HF solution !
!-----------------------------------!
if(dostab) then
call wall_time(start_stab)
call RHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_MO)
call wall_time(end_stab)
t_stab = end_stab - start_stab
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds'
write(*,*)
end if
!-----------------------!
! Moller-Plesset module !
!-----------------------!
doMP = doMP2 .or. doMP3
if(doMP) then
call wall_time(start_MP)
call RMP(doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP = ',t_MP,' seconds'
write(*,*)
end if
!------------------------!
! Coupled-cluster module !
!------------------------!
doCC = doCCD .or. dopCCD .or. doDCD .or. doCCSD .or. doCCSDT .or. &
dodrCCD .or. dorCCD .or. docrCCD .or. dolCCD
if(doCC) then
call wall_time(start_CC)
call RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, &
maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds'
write(*,*)
end if
!----------------------------------!
! Configuration interaction module !
!----------------------------------!
doCI = doCIS .or. doCID .or. doCISD .or. doFCI
if(doCI) then
call wall_time(start_CI)
call RCI(doCIS,doCIS_D,doCID,doCISD,doFCI,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, &
epsHF,EHF,cHF,S)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CI = ',t_CI,' seconds'
write(*,*)
end if
!-----------------------------------!
! Random-phase approximation module !
!-----------------------------------!
doRPA = dophRPA .or. dophRPAx .or. docrRPA .or. doppRPA
if(doRPA) then
call wall_time(start_RPA)
call RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
write(*,*)
end if
!-------------------------!
! Green's function module !
!-------------------------!
doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doG0F3 .or. doevGF3
if(doGF) then
call wall_time(start_GF)
call RGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, &
dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,lin_GF,eta_GF,reg_GF, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO, &
dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
!-----------!
! GW module !
!-----------!
doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW .or. doSRGqsGW
if(doGW) then
call wall_time(start_GW)
call RGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GW = ',t_GW,' seconds'
write(*,*)
end if
!-----------------!
! T-matrix module !
!-----------------!
doGT = doG0T0pp .or. doevGTpp .or. doqsGTpp .or. doG0T0eh .or. doevGTeh .or. doqsGTeh
if(doGT) then
call wall_time(start_GT)
call RGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, &
lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GT = ',t_GT,' seconds'
write(*,*)
end if
end subroutine

338
src/QuAcK/UQuAcK.f90 Normal file
View File

@ -0,0 +1,338 @@
subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, &
doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, &
maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, &
TDA_W,lin_GW,reg_GW,eta_GW,maxSCF_GT,max_diis_GT,thresh_GT,TDA_T,lin_GT,reg_GT,eta_GT, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,doACFDT,exchange_kernel,doXBS)
implicit none
include 'parameters.h'
logical,intent(in) :: doUHF
logical,intent(in) :: dostab
logical,intent(in) :: doMP2,doMP3
logical,intent(in) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(in) :: dodrCCD,dorCCD,docrCCD,dolCCD
logical,intent(in) :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical,intent(in) :: dophRPA,dophRPAx,docrRPA,doppRPA
logical,intent(in) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical,intent(in) :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW
logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh
integer,intent(in) :: nNuc,nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: ENuc
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,nBas)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas)
integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix
logical,intent(in) :: guess_type
logical,intent(in) :: reg_MP
integer,intent(in) :: maxSCF_CC,max_diis_CC
double precision,intent(in) :: thresh_CC
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
logical,intent(in) :: TDA
integer,intent(in) :: maxSCF_GF,max_diis_GF,renorm_GF
double precision,intent(in) :: thresh_GF
logical,intent(in) :: lin_GF,reg_GF
double precision,intent(in) :: eta_GF
integer,intent(in) :: maxSCF_GW,max_diis_GW
double precision,intent(in) :: thresh_GW
logical,intent(in) :: TDA_W,lin_GW,reg_GW
double precision,intent(in) :: eta_GW
integer,intent(in) :: maxSCF_GT,max_diis_GT
double precision,intent(in) :: thresh_GT
logical,intent(in) :: TDA_T,lin_GT,reg_GT
double precision,intent(in) :: eta_GT
logical,intent(in) :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA
logical,intent(in) :: doACFDT,exchange_kernel,doXBS
! Local variables
logical :: doMP,doCC,doCI,doRPA,doGF,doGW,doGT
double precision :: start_HF ,end_HF ,t_HF
double precision :: start_stab ,end_stab ,t_stab
double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO
double precision :: start_MP ,end_MP ,t_MP
double precision :: start_CC ,end_CC ,t_CC
double precision :: start_CI ,end_CI ,t_CI
double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_GF ,end_GF ,t_GF
double precision :: start_GW ,end_GW ,t_GW
double precision :: start_GT ,end_GT ,t_GT
double precision,allocatable :: cHF(:,:,:),epsHF(:,:),PHF(:,:,:)
double precision :: EHF
double precision,allocatable :: dipole_int_aa(:,:,:),dipole_int_bb(:,:,:)
double precision,allocatable :: ERI_aaaa(:,:,:,:),ERI_aabb(:,:,:,:),ERI_bbbb(:,:,:,:)
integer :: ixyz
write(*,*)
write(*,*) '********************************'
write(*,*) '* Unrestricted Branch of QuAcK *'
write(*,*) '********************************'
write(*,*)
!-------------------!
! Memory allocation !
!-------------------!
allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin), &
dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart), &
ERI_aaaa(nBas,nBas,nBas,nBas),ERI_aabb(nBas,nBas,nBas,nBas),ERI_bbbb(nBas,nBas,nBas,nBas))
!---------------------!
! Hartree-Fock module !
!---------------------!
if(doUHF) then
call wall_time(start_HF)
call UHF(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UHF = ',t_HF,' seconds'
write(*,*)
end if
!------------------!
! Kohn-Sham module !
!------------------!
! if(doKS) then
! call wall_time(start_KS)
! write(*,*)
! write(*,*) 'KS module has been disabled for now! Sorry.'
! write(*,*)
! call eDFT(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nC, &
! nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, &
! max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EHF,epsHF,cHF,PHF,Vxc)
! call wall_time(end_KS)
! t_KS = end_KS - start_KS
! write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for KS = ',t_KS,' seconds'
! write(*,*)
! end if
!----------------------------------!
! AO to MO integral transformation !
!----------------------------------!
call wall_time(start_AOtoMO)
write(*,*)
write(*,*) 'AO to MO transformation... Please be patient'
write(*,*)
! Read and transform dipole-related integrals
dipole_int_aa(:,:,:) = dipole_int_AO(:,:,:)
dipole_int_bb(:,:,:) = dipole_int_AO(:,:,:)
do ixyz=1,ncart
call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz))
call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz))
end do
! 4-index transform for (aa|aa) block
call AOtoMO_integral_transform(1,1,1,1,nBas,cHF,ERI_AO,ERI_aaaa)
! 4-index transform for (aa|bb) block
call AOtoMO_integral_transform(1,1,2,2,nBas,cHF,ERI_AO,ERI_aabb)
! 4-index transform for (bb|bb) block
call AOtoMO_integral_transform(2,2,2,2,nBas,cHF,ERI_AO,ERI_bbbb)
call wall_time(end_AOtoMO)
t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*)
!-----------------------------------!
! Stability analysis of HF solution !
!-----------------------------------!
if(dostab) then
call wall_time(start_stab)
call UHF_stability(nBas,nC,nO,nV,nR,nS,epsHF,ERI_aaaa,ERI_aabb,ERI_bbbb)
call wall_time(end_stab)
t_stab = end_stab - start_stab
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for stability analysis = ',t_stab,' seconds'
write(*,*)
end if
!-----------------------!
! Moller-Plesset module !
!-----------------------!
doMP = doMP2 .or. doMP3
if(doMP) then
call wall_time(start_MP)
call UMP(doMP2,doMP3,reg_MP,nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EHF,epsHF)
call wall_time(end_MP)
t_MP = end_MP - start_MP
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP = ',t_MP,' seconds'
write(*,*)
end if
!------------------------!
! Coupled-cluster module !
!------------------------!
doCC = doCCD .or. dopCCD .or. doDCD .or. doCCSD .or. doCCSDT .or. &
dodrCCD .or. dorCCD .or. docrCCD .or. dolCCD
if(doCC) then
call wall_time(start_CC)
call wall_time(end_CC)
t_CC = end_CC - start_CC
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CC = ',t_CC,' seconds'
write(*,*)
end if
!----------------------------------!
! Configuration interaction module !
!----------------------------------!
doCI = doCIS .or. doCID .or. doCISD .or. doFCI
if(doCI) then
call wall_time(start_CI)
call UCI(doCIS,doCIS_D,doCID,doCISD,doFCI,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS, &
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,EHF,cHF,S)
call wall_time(end_CI)
t_CI = end_CI - start_CI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CI = ',t_CI,' seconds'
write(*,*)
end if
!-----------------------------------!
! Random-phase approximation module !
!-----------------------------------!
doRPA = dophRPA .or. dophRPAx .or. docrRPA .or. doppRPA
if(doRPA) then
call wall_time(start_RPA)
call URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
write(*,*)
end if
!-------------------------!
! Green's function module !
!-------------------------!
doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doG0F3 .or. doevGF3
if(doGF) then
call wall_time(start_GF)
call UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, &
dophBSE,doppBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,lin_GF,eta_GF,reg_GF, &
nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
call wall_time(end_GF)
t_GF = end_GF - start_GF
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds'
write(*,*)
end if
!-----------!
! GW module !
!-----------!
doGW = doG0W0 .or. doevGW .or. doqsGW .or. doufG0W0 .or. doufGW .or. doSRGqsGW
if(doGW) then
call wall_time(start_GW)
call UGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
call wall_time(end_GW)
t_GW = end_GW - start_GW
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GW = ',t_GW,' seconds'
write(*,*)
end if
!-----------------!
! T-matrix module !
!-----------------!
doGT = doG0T0pp .or. doevGTpp .or. doqsGTpp .or. doG0T0eh .or. doevGTeh .or. doqsGTeh
if(doGT) then
call wall_time(start_GT)
call UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF_GT,thresh_GT,max_diis_GT,doACFDT, &
exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip, &
lin_GT,eta_GT,reg_GT,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, &
ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF)
call wall_time(end_GT)
t_GT = end_GT - start_GT
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GT = ',t_GT,' seconds'
write(*,*)
end if
end subroutine

132
src/QuAcK/read_methods.f90
View File

@ -1,14 +1,14 @@
subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & subroutine read_methods(doRHF,doUHF,doGHF,doROHF, &
doMP2,doMP3, & doMP2,doMP3, &
doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & doCCD,dopCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_crCCD,do_lCCD, & do_drCCD,do_rCCD,do_crCCD,do_lCCD, &
doCIS,doCIS_D,doCID,doCISD,doFCI, & doCIS,doCIS_D,doCID,doCISD,doFCI, &
dophRPA,dophRPAx,docrRPA,doppRPA, & dophRPA,dophRPAx,docrRPA,doppRPA, &
doG0F2,doevGF2,doqsGF2, & doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, & doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW,doSRGqsGW, & doG0W0,doevGW,doqsGW,doSRGqsGW, &
doufG0W0,doufGW, & doufG0W0,doufGW, &
doG0T0pp,doevGTpp,doqsGTpp, & doG0T0pp,doevGTpp,doqsGTpp, &
doG0T0eh,doevGTeh,doqsGTeh) doG0T0eh,doevGTeh,doqsGTeh)
! Read desired methods ! Read desired methods
@ -17,7 +17,7 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
! Input variables ! Input variables
logical,intent(out) :: doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS logical,intent(out) :: doRHF,doUHF,doGHF,doROHF
logical,intent(out) :: doMP2,doMP3 logical,intent(out) :: doMP2,doMP3
logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT logical,intent(out) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD logical,intent(out) :: do_drCCD,do_rCCD,do_crCCD,do_lCCD
@ -30,7 +30,7 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
! Local variables ! Local variables
character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6 character(len=1) :: ans1,ans2,ans3,ans4,ans5,ans6
! Open file with method specification ! Open file with method specification
@ -40,19 +40,15 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
doRHF = .false. doRHF = .false.
doUHF = .false. doUHF = .false.
doGHF = .false.
doROHF = .false. doROHF = .false.
doRMOM = .false.
doUMOM = .false.
doKS = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 read(1,*) ans1,ans2,ans3,ans4
if(answer1 == 'T') doRHF = .true. if(ans1 == 'T') doRHF = .true.
if(answer2 == 'T') doUHF = .true. if(ans2 == 'T') doUHF = .true.
if(answer3 == 'T') doROHF = .true. if(ans3 == 'T') doGHF = .true.
if(answer4 == 'T') doRMOM = .true. if(ans4 == 'T') doROHF = .true.
if(answer5 == 'T') doUMOM = .true.
if(answer6 == 'T') doKS = .true.
! Read MPn methods ! Read MPn methods
@ -60,9 +56,9 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
doMP3 = .false. doMP3 = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2 read(1,*) ans1,ans2
if(answer1 == 'T') doMP2 = .true. if(ans1 == 'T') doMP2 = .true.
if(answer2 == 'T') doMP3 = .true. if(ans2 == 'T') doMP3 = .true.
! Read CC methods ! Read CC methods
@ -73,12 +69,12 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
doCCSDT = .false. doCCSDT = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4,answer5 read(1,*) ans1,ans2,ans3,ans4,ans5
if(answer1 == 'T') doCCD = .true. if(ans1 == 'T') doCCD = .true.
if(answer2 == 'T') dopCCD = .true. if(ans2 == 'T') dopCCD = .true.
if(answer3 == 'T') doDCD = .true. if(ans3 == 'T') doDCD = .true.
if(answer4 == 'T') doCCSD = .true. if(ans4 == 'T') doCCSD = .true.
if(answer5 == 'T') doCCSDT = .true. if(ans5 == 'T') doCCSDT = .true.
! Read weird CC methods ! Read weird CC methods
@ -88,11 +84,11 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
do_lCCD = .false. do_lCCD = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4 read(1,*) ans1,ans2,ans3,ans4
if(answer1 == 'T') do_drCCD = .true. if(ans1 == 'T') do_drCCD = .true.
if(answer2 == 'T') do_rCCD = .true. if(ans2 == 'T') do_rCCD = .true.
if(answer3 == 'T') do_crCCD = .true. if(ans3 == 'T') do_crCCD = .true.
if(answer4 == 'T') do_lCCD = .true. if(ans4 == 'T') do_lCCD = .true.
! Read CI methods ! Read CI methods
@ -103,12 +99,12 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
doFCI = .false. doFCI = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4,answer5 read(1,*) ans1,ans2,ans3,ans4,ans5
if(answer1 == 'T') doCIS = .true. if(ans1 == 'T') doCIS = .true.
if(answer2 == 'T') doCIS_D = .true. if(ans2 == 'T') doCIS_D = .true.
if(answer3 == 'T') doCID = .true. if(ans3 == 'T') doCID = .true.
if(answer4 == 'T') doCISD = .true. if(ans4 == 'T') doCISD = .true.
if(answer5 == 'T') doFCI = .true. if(ans5 == 'T') doFCI = .true.
if(doCIS_D) doCIS = .true. if(doCIS_D) doCIS = .true.
! Read RPA methods ! Read RPA methods
@ -119,11 +115,11 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
doppRPA = .false. doppRPA = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4 read(1,*) ans1,ans2,ans3,ans4
if(answer1 == 'T') dophRPA = .true. if(ans1 == 'T') dophRPA = .true.
if(answer2 == 'T') dophRPAx = .true. if(ans2 == 'T') dophRPAx = .true.
if(answer3 == 'T') docrRPA = .true. if(ans3 == 'T') docrRPA = .true.
if(answer4 == 'T') doppRPA = .true. if(ans4 == 'T') doppRPA = .true.
! Read Green's function methods ! Read Green's function methods
@ -134,12 +130,12 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
doevGF3 = .false. doevGF3 = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4,answer5 read(1,*) ans1,ans2,ans3,ans4,ans5
if(answer1 == 'T') doG0F2 = .true. if(ans1 == 'T') doG0F2 = .true.
if(answer2 == 'T') doevGF2 = .true. if(ans2 == 'T') doevGF2 = .true.
if(answer3 == 'T') doqsGF2 = .true. if(ans3 == 'T') doqsGF2 = .true.
if(answer4 == 'T') doG0F3 = .true. if(ans4 == 'T') doG0F3 = .true.
if(answer5 == 'T') doevGF3 = .true. if(ans5 == 'T') doevGF3 = .true.
! Read GW methods ! Read GW methods
@ -151,13 +147,13 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
doufGW = .false. doufGW = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 read(1,*) ans1,ans2,ans3,ans4,ans5,ans6
if(answer1 == 'T') doG0W0 = .true. if(ans1 == 'T') doG0W0 = .true.
if(answer2 == 'T') doevGW = .true. if(ans2 == 'T') doevGW = .true.
if(answer3 == 'T') doqsGW = .true. if(ans3 == 'T') doqsGW = .true.
if(answer4 == 'T') doSRGqsGW = .true. if(ans4 == 'T') doSRGqsGW = .true.
if(answer5 == 'T') doufG0W0 = .true. if(ans5 == 'T') doufG0W0 = .true.
if(answer6 == 'T') doufGW = .true. if(ans6 == 'T') doufGW = .true.
! Read GT methods ! Read GT methods
@ -169,13 +165,13 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, &
doqsGTeh = .false. doqsGTeh = .false.
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 read(1,*) ans1,ans2,ans3,ans4,ans5,ans6
if(answer1 == 'T') doG0T0pp = .true. if(ans1 == 'T') doG0T0pp = .true.
if(answer2 == 'T') doevGTpp = .true. if(ans2 == 'T') doevGTpp = .true.
if(answer3 == 'T') doqsGTpp = .true. if(ans3 == 'T') doqsGTpp = .true.
if(answer4 == 'T') doG0T0eh = .true. if(ans4 == 'T') doG0T0eh = .true.
if(answer5 == 'T') doevGTeh = .true. if(ans5 == 'T') doevGTeh = .true.
if(answer6 == 'T') doqsGTeh = .true. if(ans6 == 'T') doqsGTeh = .true.
! Close file with geometry specification ! Close file with geometry specification

96
src/QuAcK/read_options.f90
View File

@ -1,11 +1,11 @@
subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, & subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab, &
reg_MP, & reg_MP, &
maxSCF_CC,thresh_CC,DIIS_CC,max_diis_CC, & maxSCF_CC,thresh_CC,max_diis_CC, &
TDA,singlet,triplet,spin_conserved,spin_flip, & TDA,spin_conserved,spin_flip, &
maxSCF_GF,thresh_GF,DIIS_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, & maxSCF_GF,thresh_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, &
maxSCF_GW,thresh_GW,DIIS_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, & maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,DIIS_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, & doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA) dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
! Read desired methods ! Read desired methods
@ -16,11 +16,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho
integer,intent(out) :: maxSCF_HF integer,intent(out) :: maxSCF_HF
double precision,intent(out) :: thresh_HF double precision,intent(out) :: thresh_HF
logical,intent(out) :: DIIS_HF
integer,intent(out) :: max_diis_HF integer,intent(out) :: max_diis_HF
integer,intent(out) :: guess_type integer,intent(out) :: guess_type
integer,intent(out) :: ortho_type double precision,intent(out) :: mix
logical,intent(out) :: mix
double precision,intent(out) :: level_shift double precision,intent(out) :: level_shift
logical,intent(out) :: dostab logical,intent(out) :: dostab
@ -28,18 +26,14 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho
integer,intent(out) :: maxSCF_CC integer,intent(out) :: maxSCF_CC
double precision,intent(out) :: thresh_CC double precision,intent(out) :: thresh_CC
logical,intent(out) :: DIIS_CC
integer,intent(out) :: max_diis_CC integer,intent(out) :: max_diis_CC
logical,intent(out) :: TDA logical,intent(out) :: TDA
logical,intent(out) :: singlet
logical,intent(out) :: triplet
logical,intent(out) :: spin_conserved logical,intent(out) :: spin_conserved
logical,intent(out) :: spin_flip logical,intent(out) :: spin_flip
integer,intent(out) :: maxSCF_GF integer,intent(out) :: maxSCF_GF
double precision,intent(out) :: thresh_GF double precision,intent(out) :: thresh_GF
logical,intent(out) :: DIIS_GF
integer,intent(out) :: max_diis_GF integer,intent(out) :: max_diis_GF
logical,intent(out) :: lin_GF logical,intent(out) :: lin_GF
integer,intent(out) :: renorm_GF integer,intent(out) :: renorm_GF
@ -48,7 +42,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho
integer,intent(out) :: maxSCF_GW integer,intent(out) :: maxSCF_GW
double precision,intent(out) :: thresh_GW double precision,intent(out) :: thresh_GW
logical,intent(out) :: DIIS_GW
integer,intent(out) :: max_diis_GW integer,intent(out) :: max_diis_GW
logical,intent(out) :: TDA_W logical,intent(out) :: TDA_W
logical,intent(out) :: lin_GW logical,intent(out) :: lin_GW
@ -57,7 +50,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho
integer,intent(out) :: maxSCF_GT integer,intent(out) :: maxSCF_GT
double precision,intent(out) :: thresh_GT double precision,intent(out) :: thresh_GT
logical,intent(out) :: DIIS_GT
integer,intent(out) :: max_diis_GT integer,intent(out) :: max_diis_GT
logical,intent(out) :: TDA_T logical,intent(out) :: TDA_T
logical,intent(out) :: lin_GT logical,intent(out) :: lin_GT
@ -86,22 +78,16 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho
maxSCF_HF = 64 maxSCF_HF = 64
thresh_HF = 1d-6 thresh_HF = 1d-6
DIIS_HF = .false. max_diis_HF = 1
max_diis_HF = 5
guess_type = 1 guess_type = 1
ortho_type = 1 mix = 0d0
mix = .false.
level_shift = 0d0 level_shift = 0d0
dostab = .false. dostab = .false.
read(1,*) read(1,*)
read(1,*) maxSCF_HF,thresh_HF,ans1,max_diis_HF,guess_type,ortho_type,ans2,level_shift,ans3 read(1,*) maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,ans1
if(ans1 == 'T') DIIS_HF = .true. if(ans1 == 'T') dostab = .true.
if(ans2 == 'T') mix = .true.
if(ans3 == 'T') dostab = .true.
if(.not.DIIS_HF) max_diis_HF = 1
! Read MPn options ! Read MPn options
@ -115,91 +101,73 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho
maxSCF_CC = 64 maxSCF_CC = 64
thresh_CC = 1d-5 thresh_CC = 1d-5
DIIS_CC = .false. max_diis_CC = 1
max_diis_CC = 5
read(1,*) read(1,*)
read(1,*) maxSCF_CC,thresh_CC,ans1,max_diis_CC read(1,*) maxSCF_CC,thresh_CC,max_diis_CC
if(ans1 == 'T') DIIS_CC = .true.
if(.not.DIIS_CC) max_diis_CC = 1
! Read excited state options ! Read excited state options
TDA = .false. TDA = .false.
singlet = .false.
triplet = .false.
spin_conserved = .false. spin_conserved = .false.
spin_flip = .false. spin_flip = .false.
read(1,*) read(1,*)
read(1,*) ans1,ans2,ans3,ans4,ans5 read(1,*) ans1,ans2,ans3
if(ans1 == 'T') TDA = .true. if(ans1 == 'T') TDA = .true.
if(ans2 == 'T') singlet = .true. if(ans2 == 'T') spin_conserved = .true.
if(ans3 == 'T') triplet = .true. if(ans3 == 'T') spin_flip = .true.
if(ans4 == 'T') spin_conserved = .true.
if(ans5 == 'T') spin_flip = .true.
! Read GF options ! Read GF options
maxSCF_GF = 64 maxSCF_GF = 64
thresh_GF = 1d-5 thresh_GF = 1d-5
DIIS_GF = .false. max_diis_GF = 1
max_diis_GF = 5
lin_GF = .false. lin_GF = .false.
eta_GF = 0d0 eta_GF = 0d0
renorm_GF = 0 renorm_GF = 0
reg_GF = .false. reg_GF = .false.
read(1,*) read(1,*)
read(1,*) maxSCF_GF,thresh_GF,ans1,max_diis_GF,ans2,eta_GF,renorm_GF,ans3 read(1,*) maxSCF_GF,thresh_GF,max_diis_GF,ans1,eta_GF,renorm_GF,ans2
if(ans1 == 'T') DIIS_GF = .true. if(ans1 == 'T') lin_GF = .true.
if(ans2 == 'T') lin_GF = .true. if(ans2 == 'T') reg_GF = .true.
if(ans3 == 'T') reg_GF = .true.
if(.not.DIIS_GF) max_diis_GF = 1
! Read GW options ! Read GW options
maxSCF_GW = 64 maxSCF_GW = 64
thresh_GW = 1d-5 thresh_GW = 1d-5
DIIS_GW = .false. max_diis_GW = 1
max_diis_GW = 5
lin_GW = .false. lin_GW = .false.
eta_GW = 0d0 eta_GW = 0d0
reg_GW = .false. reg_GW = .false.
TDA_W = .false. TDA_W = .false.
read(1,*) read(1,*)
read(1,*) maxSCF_GW,thresh_GW,ans1,max_diis_GW,ans2,eta_GW,ans3,ans4 read(1,*) maxSCF_GW,thresh_GW,max_diis_GW,ans1,eta_GW,ans2,ans3
if(ans1 == 'T') DIIS_GW = .true. if(ans1 == 'T') lin_GW = .true.
if(ans2 == 'T') lin_GW = .true. if(ans2 == 'T') TDA_W = .true.
if(ans3 == 'T') TDA_W = .true. if(ans3 == 'T') reg_GW = .true.
if(ans4 == 'T') reg_GW = .true.
if(.not.DIIS_GW) max_diis_GW = 1
! Read GT options ! Read GT options
maxSCF_GT = 64 maxSCF_GT = 64
thresh_GT = 1d-5 thresh_GT = 1d-5
DIIS_GT = .false. max_diis_GT = 1
max_diis_GT = 5
lin_GT = .false. lin_GT = .false.
eta_GT = 0d0 eta_GT = 0d0
reg_GT = .false. reg_GT = .false.
TDA_T = .false. TDA_T = .false.
read(1,*) read(1,*)
read(1,*) maxSCF_GT,thresh_GT,ans1,max_diis_GT,ans2,eta_GT,ans3,ans4 read(1,*) maxSCF_GT,thresh_GT,max_diis_GT,ans1,eta_GT,ans2,ans3
if(ans1 == 'T') DIIS_GT = .true. if(ans1 == 'T') lin_GT = .true.
if(ans2 == 'T') lin_GT = .true. if(ans2 == 'T') TDA_T = .true.
if(ans3 == 'T') TDA_T = .true. if(ans3 == 'T') reg_GT = .true.
if(ans4 == 'T') reg_GT = .true.
if(.not.DIIS_GT) max_diis_GT = 1
! Options for adiabatic connection ! Options for adiabatic connection

103
src/RPA/RRPA.f90 Normal file
View File

@ -0,0 +1,103 @@
subroutine RRPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,singlet,triplet, &
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF,cHF,S)
! Random-phase approximation module
implicit none
include 'parameters.h'
! Input variables
logical :: dophRPA
logical :: dophRPAx
logical :: docrRPA
logical :: doppRPA
logical,intent(in) :: TDA
logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet
logical,intent(in) :: triplet
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: epsHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
double precision :: start_RPA ,end_RPA ,t_RPA
!------------------------------------------------------------------------
! Compute (direct) RPA excitations
!------------------------------------------------------------------------
if(dophRPA) then
call wall_time(start_RPA)
call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute RPAx (RPA with exchange) excitations
!------------------------------------------------------------------------
if(dophRPAx) then
call wall_time(start_RPA)
call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPA,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute crRPA excitations
!------------------------------------------------------------------------
if(docrRPA) then
call wall_time(start_RPA)
call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute ppRPA excitations
!------------------------------------------------------------------------
if(doppRPA) then
call wall_time(start_RPA)
call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_RPA,' seconds'
write(*,*)
end if
end subroutine

39
src/RPA/RPA.f90 → src/RPA/URPA.f90
View File

@ -1,7 +1,5 @@
subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, & subroutine URPA(dophRPA,dophRPAx,docrRPA,doppRPA,TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip, &
TDA,doACFDT,exchange_kernel,singlet,triplet,spin_conserved,spin_flip, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, &
dipole_int,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
! Random-phase approximation module ! Random-phase approximation module
@ -14,13 +12,10 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
logical :: dophRPAx logical :: dophRPAx
logical :: docrRPA logical :: docrRPA
logical :: doppRPA logical :: doppRPA
logical :: unrestricted
logical,intent(in) :: TDA logical,intent(in) :: TDA
logical,intent(in) :: doACFDT logical,intent(in) :: doACFDT
logical,intent(in) :: exchange_kernel logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet
logical,intent(in) :: triplet
logical,intent(in) :: spin_conserved logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip logical,intent(in) :: spin_flip
integer,intent(in) :: nBas integer,intent(in) :: nBas
@ -34,11 +29,9 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
double precision,intent(in) :: epsHF(nBas,nspin) double precision,intent(in) :: epsHF(nBas,nspin)
double precision,intent(in) :: cHF(nBas,nBas,nspin) double precision,intent(in) :: cHF(nBas,nBas,nspin)
double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
@ -53,12 +46,8 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
if(dophRPA) then if(dophRPA) then
call wall_time(start_RPA) call wall_time(start_RPA)
if(unrestricted) then call phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
call phURPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
else
call phRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
@ -74,12 +63,8 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
if(dophRPAx) then if(dophRPAx) then
call wall_time(start_RPA) call wall_time(start_RPA)
if(unrestricted) then call phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, &
call phURPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,epsHF,cHF,S)
else
call phRPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
@ -95,11 +80,7 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
if(docrRPA) then if(docrRPA) then
call wall_time(start_RPA) call wall_time(start_RPA)
if(unrestricted) then write(*,*) 'Unrestricted version of crRPA not yet implemented! Sorry.'
write(*,*) 'Unrestricted version of crRPA not yet implemented! Sorry.'
else
call crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA
@ -115,11 +96,7 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted,
if(doppRPA) then if(doppRPA) then
call wall_time(start_RPA) call wall_time(start_RPA)
if(unrestricted) then call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF)
call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF)
else
call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF)
end if
call wall_time(end_RPA) call wall_time(end_RPA)
t_RPA = end_RPA - start_RPA t_RPA = end_RPA - start_RPA

5
src/utils/orthogonalization_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine orthogonalization_matrix(ortho_type,nBas,S,X) subroutine orthogonalization_matrix(nBas,S,X)
! Compute the orthogonalization matrix X ! Compute the orthogonalization matrix X
@ -6,7 +6,7 @@ subroutine orthogonalization_matrix(ortho_type,nBas,S,X)
! Input variables ! Input variables
integer,intent(in) :: nBas,ortho_type integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: S(nBas,nBas)
! Local variables ! Local variables
@ -14,6 +14,7 @@ subroutine orthogonalization_matrix(ortho_type,nBas,S,X)
logical :: debug logical :: debug
double precision,allocatable :: UVec(:,:),Uval(:) double precision,allocatable :: UVec(:,:),Uval(:)
double precision,parameter :: thresh = 1d-6 double precision,parameter :: thresh = 1d-6
integer,parameter :: ortho_type = 1
integer :: i integer :: i

28
src/utils/read_molecule.f90
View File

@ -1,4 +1,4 @@
subroutine read_molecule(nNuc,nEl,nO,nC,nR) subroutine read_molecule(nNuc,nO,nC,nR)
! Read number of atoms and number of electrons ! Read number of atoms and number of electrons
@ -6,19 +6,18 @@ subroutine read_molecule(nNuc,nEl,nO,nC,nR)
include 'parameters.h' include 'parameters.h'
! Input variables
integer,intent(out) :: nNuc
integer,intent(out) :: nEl(nspin)
integer,intent(out) :: nO(nspin)
integer,intent(out) :: nC(nspin)
integer,intent(out) :: nR(nspin)
! Local variables ! Local variables
integer :: nCore integer :: nCore
integer :: nRyd integer :: nRyd
! Output variables
integer,intent(out) :: nNuc
integer,intent(out) :: nO(nspin)
integer,intent(out) :: nC(nspin)
integer,intent(out) :: nR(nspin)
! Open file with geometry specification ! Open file with geometry specification
open(unit=1,file='input/molecule') open(unit=1,file='input/molecule')
@ -26,7 +25,7 @@ subroutine read_molecule(nNuc,nEl,nO,nC,nR)
! Read number of atoms and number of electrons ! Read number of atoms and number of electrons
read(1,*) read(1,*)
read(1,*) nNuc,nEl(1),nEl(2),nCore,nRyd read(1,*) nNuc,nO(1),nO(2),nCore,nRyd
if(mod(nCore,2) /= 0 .or. mod(nRyd,2) /= 0) then if(mod(nCore,2) /= 0 .or. mod(nRyd,2) /= 0) then
@ -35,9 +34,8 @@ subroutine read_molecule(nNuc,nEl,nO,nC,nR)
end if end if
nO(:) = nEl(:)
nC(:) = nCore/2 nC(:) = nCore/2
nR(:) = nRyd nR(:) = nRyd/2
! Print results ! Print results
@ -46,9 +44,9 @@ subroutine read_molecule(nNuc,nEl,nO,nC,nR)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,*) write(*,*)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of spin-up electrons',nEl(1) write(*,'(A28,1X,I16)') 'Number of spin-up electrons',nO(1)
write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nEl(2) write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nO(2)
write(*,'(A28,1X,I16)') ' Total number of electrons',sum(nEl(:)) write(*,'(A28,1X,I16)') ' Total number of electrons',sum(nO(:))
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'
write(*,*) write(*,*)
write(*,'(A28)') '----------------------' write(*,'(A28)') '----------------------'